20241015 Aanvullende info Rapport Bodemdalingsprognose_Flevoland.PDF

Figuur 7 Toepassing op bodemkolommen binnen een 3D voxel model van de ondergrond.

De berekeningen zijn onderverdeeld in ‘stress periods’ en tijdstappen. Dit is rechtsonder in Figuur 6 weergegeven. Deze berekeningen vinden plaats binnen het witte kader in Figuur 6. Binnen een stress period heeft de GLG een vast peil. Door bodemdaling treedt dus daling op van voxels ten opzichte van de GLG (‘verdrinking’). Dit remt de bodemdalingsprocessen. Op de overgang naar een volgende stress period kan het GLG peil worden gewijzigd. Over het algemeen betreft dit GLG verlaging. Dit is in Figuur 6 weergegeven met het blok ‘aandrijving’. Hierin kunnen o.a. scenario’s voor de invloed van waterbeheer en klimaatverandering op de GLG worden gemodelleerd. Zo kan GLG peilverlaging bijvoorbeeld worden gekoppeld aan een scenario voor oppervlaktewaterpeilverlaging13. Ook kan de GLG geheel of gedeeltelijk worden “geïndexeerd” voor de bodemdaling in de voorgaande stress periode. Peilverlaging of GLG-indexactie voor bodemdaling kan worden gekoppeld aan “peilgebieden” (links in Figuur 6). Het GLG-peil kan ook worden aangepast volgens een gespecificeerde andere trend, bijvoorbeeld om de invloed van een klimaatontwikkeling te simuleren. Al deze invloeden kunnen worden gecombineerd. In Atlantis wordt de GLG om diverse redenen niet direct berekend via een grondwatermodel.

—————————————— 13

32 van 114

Oppervlaktewaterpeil zelf is echter niet expliciet aanwezig in Atlantis.

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

7

AHN analyse Dit hoofdstuk beschrijft de werkzaamheden en resultaten binnen het onderdeel ‘bodemdalingsmetingen’ (Figuur 1). ‘Bodemdalingsmetingen’ staat hier tussen aanhalingstekens omdat het feitelijk gaat om verwerking van hoogtedata uit het Actueel hoogtebestand Nederland (AHN). Deze hoogtedata zijn en worden niet verkregen voor bodemdalingstoepassingen. Daadwerkelijke (geodetische) bodemdalingsmetingen richten zich formeel op metingen van hoogteverandering, wat meer omvat dan individuele hoogtemetingen. Bij gebrek aan bodemdalingsmetingen worden de hoogtemetingen van het AHN hier wel benut om informatie uit af te leiden over hoogteverandering. De analyse bestaat uit twee onderdelen. In het eerste onderdeel zijn hoogteschattingen voor AHN3 en AHN4 toegevoegd aan bestaande datasets van ZZL voor de rijpingsterreinen in Zuidelijk en Oostelijk Flevoland. In het tweede onderdeel is een provincie dekkende kaart gemaakt van de ‘recente bodemdalingssnelheid’ (2009-2022) op basis van AHN2, AHN3 en AHN4.

7.1

Aangevulde tijdreeksen maaivelddaling rijpingsterreinen De datasets en de grafieken die met de Bodemdalingskaart Flevoland (De Lange e.a., 2012) zijn opgeleverd zijn uitgebreid met de maaivelddalingsschattingen voor edities 3 & 4 van het AHN. AHN3 en AHN4 zijn zowel in Oostelijk als in Zuidelijke Flevoland respectievelijk ingewonnen in de jaren 2016 en 2020.

7.1.1

Methode Om een consistente dataset te behouden is de oorspronkelijke werkwijze uit 2012 toegepast: - Gebruik van de 5m resolutie rasterkaarten AHN. Per rijpingsterrein (oppervlak ca. 50 m x 60 m): - Selectie van alle cellen waarvan het middelpunt zich in de polygoon van het rijpingsterrein bevindt. - Bepaling van de gemiddelde, minimale en maximale NAP hoogte alsmede variantie en standaarddeviatie van de NAP hoogte. - Bepaling van de daling voor de gemiddelde, minimale en maximale NAP hoogte ten opzichte van de start NAP hoogte zoals vermeld in de oorspronkelijke dataset (2012). Voor een aantal rijpingsterreinen is de werkwijze ook toegepast op AHN2 en zijn resultaten vergeleken met de reeds bestaande waarden in de dataset. Daaruit bleek dat de verschillen verwaarloosbaar zijn (kleiner dan 1 mm).

7.1.2

Resultaten De resultaten zijn opgenomen in digitale bijlage 19A.2. Figuur 8 toont als voorbeeld de grafiek voor rijpingsterrein OF 1. De datapunten van de RIJP metingen en AHN1 en AHN2 zijn overgenomen uit de oorspronkelijke dataset. De datapunten voor AHN3 en AHN4 zijn de nieuwe toevoegingen. Het symbool voor de AHN datapunten correspondeert met de gemiddelde hoogte van het rijpingsterrein. Belangrijk: De bandbreedte die is aangegeven voor elk van de AHN dalingsschattingen is géén maat voor de onzekerheid in de daling. De bandbreedte correspondeert met het bereik van de hoogte binnen het rijpingsterrein (minimale en maximale waarde). De onzekerheid van de daling wordt bepaald door de onzekerheid/fout in de AHN hoogte. Voor agrarisch land

33 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

is die onzekerheid relatief groot ten opzichte van nauwkeurig definieerbare oppervlak zoals wegdek (maaiveldhoogte van bouwland heeft een lage idealisatie precisie). Voor agrarisch land en met name bouwland moet gedacht worden aan een standaardafwijking van de hoogte van 1 a 2 decimeter (vakgroep Geoscience and Remote Sensing TUD, pers. comm.).

Figuur 8 Grafische weergave van het resultaat voor rijpingsterrein OF 1. “RIJP metingen” zijn gebaseerd op Flevobericht 388. De bandbreedte aangeven bij de AHN dalingsschattingen zijn o.b.v. de minimale en maximale hoogtewaarden (5 m x 5 m) binnen het rijpingsterrein.

Er is geen uitgebreide studie gedaan naar wat de trends in de data aan informatie verschaffen. De tijdreeksen voor diverse terreinen wekken de indruk dat de daling tussen 2016 en 2020 wat sterker was dan in de jaren daarvoor (OF6, OF17, OF14, ZF18, ZF19, ZF29, ZF34), wat zou kunnen wijzen op een invloed van de recente droge zomers, maar dat kan niet met zekerheid worden vastgesteld gezien de mogelijke fout in de individuele hoogteschattingen. Op diverse plaatsen lijken ook recente ophogingen het dalingsbeeld te verstoren zoals nabij Zeewolde (ZF10) en Almere (ZF26). Een opvallende grote dalingsversnelling lijkt te hebben plaatsgevonden in het noordoostelijk deel van de Oostvaardersplassen (ZF33). Mogelijk hangt dit samen met (tijdelijke) verlaging van het waterpeil om rietgroei te beperken.

7.2

Maaivelddalingssnelheidskaart op basis van AHN2,3,4

7.2.1

Methode Voor de productie van de snelheidskaart is gebruik gemaakt van de 5 m rasters van ANH2 (2009), AHN3 (2016) en AHN4 (ZF en OF: 2020; NOP: 2022). De totale periode is dus 11 tot 13 jaar. Voor elk rasterpunt levert dat hoogteschattingen voor drie verschillende jaren (Figuur 9). Per rasterpunt is vervolgens een lineaire regressie uitgevoerd om de lineaire snelheid van maaiveldhoogteverandering én de fout daarin te schatten (Figuur 10). De foutschatting is essentieel om de kaartinformatie te kunnen gebruiken (significantie). Met advisering door de vakgroep Geosciences and Remote Sensing van de TU Delft is de onderstaande rekenwijze opgesteld en gebruikt. Box 2: Rekenwijze voor de productie van de snelheidskaart De notatie van de parameters voor de ‘meetdata’ voor elk gridpunt in kaartbeeld is weergegeven in Figuur 9. Aan elke AHN hoogte is een standaarddeviatie  toegekend.

34 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 9 Gebruikte notaties voor de AHN hoogte (h), standaarddeviatie van de hoogte (s) en tijd (t). Tijd is de tijd sinds de AHN2 opname (2009). Voor de schatting van de lineaire snelheid wordt een ‘simpele’ kleinste-kwadraten methode gebruikt (regressie). De standaarddeviaties worden meegenomen in een covariantiematrix. Het lineaire model, weergegeven me de rode lijn in Figuur 10, is:

hˆ0 = hˆ0 hˆ = hˆ + vt 1

0

1

hˆ2 = hˆ0 + vt2

Figuur 10 Notatie van de hoogten voor het lineaire model (rode lijn) dat wordt gefit aan de meetdata. In matrixvorm wordt het model geschreven als

hˆ = Axˆ

 hˆ0    hˆ =  hˆ1     hˆ2   

 0 1   A =  t1 1  t 1  2 

v xˆ =    hˆ   0

In deze notatie worden de meetdata (de drie ANH hoogtes) geschreven als

 h0    h =  h1  h   2

35 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

De covariantiematrix bij de meetdata:

  02 0 0    Qh =  0  12 0   0 0  22    De kleinste-kwadraten methode levert als oplossing: xˆ = Qx AT Qhinv h

Qx =  AT Qhinv A

inv

De gezochte dalingssnelheid v is de eerste component in de oplossingsvector x̂ . De gezochte standaardafwijking van v is de tweede-machtswortel van de eerste component van de variantiematrix van de oplossing Qx . De matrixvergelijkingen zijn opgenomen en toegepast in een Python script. De methode is daarmee in de toekomst makkelijk is uit te breiden met extra metingen, bijvoorbeeld AHN5.

7.2.2

Resultaat Figuur 11 en Figuur 12 tonen het snelheidsbeeld voor respectievelijk Oostelijk en Zuidelijk Flevoland en voor de Noordoostpolder in m/jaar. Digitale versies van de kaarten zijn beschikbaar als digitale bijlage 19A.4. De 1x fout in de snelheid bedraagt ten minste 0,01 m/jaar (werkwaarde 0,015 m/jaar). Dit wordt toegelicht in Box 3. Box 3: Foutschatting voor de snelheidskaart De 1x (standaarddeviatie) van de snelheid volgt uit de standaarddeviaties van de drie hoogtewaarden van AHN2, AHN3 en AHN4. AHN vermeldt een maximale systematische fout van 5 cm en een maximale stochastische fout van 5 cm voor deze hoogtedata. Als werkgetal voor de standaarddeviatie zou dan een waarde tussen 5 en 10 cm voor de ANH hoogten passend zijn voor de berekening van standaarddeviatie van de snelheid. Deze hoogtefouten zijn echter geschat op basis van validatietests voor vrij ideale oppervlakken zoals parkeerplaatsen of wegdek. Er moet vanuit worden gegaan dat de fout voor veel minder ideale oppervlakken zoals akkerland aanzienlijk groter is (geringere idealisatieprecisie). Tabel 3 toont de standaarddeviatie van de snelheid voor vier waarden voor de standaarddeviatie van de hoogte. De standaarddeviatie van de snelheid is dus behoorlijk onzeker, maar groter dan 0,01 m/jaar.

Als werkwaarde voor duiding van de kaarten wordt 0,015 m/jaar voorgesteld Tabel 3 Fout in de snelheid als functie van de fout in de hoogte

36 van 114

Standaarddeviatie hoogte

Standaarddeviatie bodemdalingssnelheid

7,5 cm

1,0 cm/jaar

10 cm

1,4 cm/jaar

15 cm

2,1 cm/jaar

20 cm

2,8 cm/jaar

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 11 Lineaire maaiveldsnelheid (m/jaar) op basis van AHN2, AHN3 en AHN4 voor Zuidelijk en Oostelijk Flevoland (beperkte resolutie weergave). Negatieve waarden betreffen maaivelddaling. Positieve waarden (stijging) zijn weergegeven als geen bodemdaling (geel) om dalingsschattingen duidelijk in beeld te krijgen. De 1x fout in de kaartwaarden wordt geschat op 0,015 m/jaar.

Figuur 12 Lineaire maaiveldsnelheid (m/jaar) op basis van AHN2, AHN3 en AHN4 voor de Noordoostpolder (beperkte resolutie weergave). Negatieve waarden betreffen maaivelddaling. Positieve waarden (stijging) zijn weergegeven als geen bodemdaling (geel) om dalingsschattingen duidelijk in beeld te krijgen. De 1x fout in de kaartwaarden wordt geschat op 0,015 m/jaar.

7.2.3

Duiding o

37 van 114

Alleen waar v > 0,015 m/jaar (donkerblauw) is kan met vrij grote waarschijnlijkheid (kans ca. 83%) worden vastgesteld dat er sprake is van bodemdaling. Dit geldt dus slechts voor

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

een beperkt deel van het totale areaal. Tegelijkertijd is er overal 50% kans dat de bodemdalingssnelheid groter is (of de bodemstijging minder is) dan de kaartwaarden. o

Donkerblauwe tinten signaleren gebiedsdelen met relatief sterke bodemdaling van meer dan 0,01 m/jaar voor de meetperiode. Deze gebieden komen voor in het centrale deel van ZF en aan de zuidzijde van de Oostvaardersplassen. In de NOP komen deze gebieden voor ten NW van Schokland en, iets minder overtuigend, ten NW van Emmeloord.

o

De snelheidswaarden in agrarisch gebied worden sterk bepaald door landbewerking (ploegen) en/of veranderingen in landbedekking (gewas). Dat kan worden opgemaakt uit het feit dat snelheden sterk kunnen verschillen tussen naast elkaar gelegen landbouwpercelen. Of er wel/niet geploegd is kort voor de AHN opname of er bij één of twee opnamen een gewas heeft gestaan kan veel uitmaken in de verkregen AHN hoogte. Deze ‘toevallige’ hoogtevariaties maken geen onderdeel uit van de ‘gezochte’ bodemdalingstrend door rijping, oxidatie en compactie en beperken daarmee de informatiewaarde voor de bodemdalingsprognose. Daarnaast wordt verwacht dat krimpzwel bewegingen zorgen voor variatie in maaiveldhoogte gedurende het jaar. Zwel en bodemstijging ontstaat wanneer uitgedroogde klei weer vocht opneemt. Deze variaties zijn naar verwachting echter klein in vergelijking met invloeden van landbewerking zoals door ploegen.

In stedelijk gebied zijn invloeden van ophoging (drukt maaivelddalingssnelheid of levert maaiveldstijging, ook bij sterke bodemdaling), of verwijdering van voorbelastingen (snelle maaivelddaling niet representatief voor bodemdaling) aanwezig.

38 van 114

o

In de NOP zijn de tochten veelal duidelijk herkenbaar door de (donker)blauwe tinten. Dit wordt geïnterpreteerd als een artefact in de hoogte informatie en snelheidsschatting zonder betekenis voor bodemdaling.

o

De kaarten bevatten vanzelfsprekend geen informatie over de diepte-intervallen in de ondergrond die bijdragen aan de bodemdaling.

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

8

Ondergrondmodel

8.1

Inleiding Het ondergrondmodel is een voxelmodel, d.w.z., een ruimtelijk model dat is opgebouwd uit horizontale lagen die elk bestaan uit aaneengesloten identieke blokken. Elke laag heeft een dikte van 0,2 m, elk blok meet 100 m bij 100 m in het horizontale vlak en is 0,2 m dik. De basis voor het voxelmodel is het ondergrondmodel van TNO (Gunnink, 2021). De bovenste 1,6 m van dit model is geactualiseerd op basis van recente bodemgegevens van WENR, beschikbaar via de BRO. In dit hoofdstuk wordt eerst het diepere deel van het voxelmodel beschreven (paragraaf 8.2), vervolgens wordt toegelicht hoe het bovenste deel van het voxelmodel, d.w.z. het bodemdeel, tot stand is gekomen (paragraaf 8.3). Vervolgens wordt beschreven hoe de NAP hoogten zijn toegekend aan het voxelmodel (paragraaf 8.4).

8.2

Productiewijze ‘diepe deel’ Voor de diepere ondergrond (vanaf 1,60 m diepte) is gebruik gemaakt van een bestaand voxelmodel van de ondergrond. Dit model is in 2021 gemaakt t.b.v. een studie van de hydraulische weerstand van de deklaag in de Flevopolders (Gunnink, 2021). Voor een gedetailleerde beschrijving van de werkwijze wordt verwezen naar bovengenoemd rapport. Hieronder volgt een beknopte beschrijving. Het voxelmodel van de ondergrond is gebaseerd op boringen die in de DINO database en BRO database zijn opgenomen, stand van zaken maart 2021. Dit zijn voor een deel zgn. RIJP boringen, die zijn gezet t.b.v. onderzoek naar de bodemgesteldheid na inpoldering. Boringen van Alterra zijn tevens gebruikt en ruim 4200 sonderingen. De sonderingen zijn bewerkt om tot een inschatting van de lithoklasse te komen. Ook zijn alle andere gegevens bewerkt om een eenduidige lithoklasse uit de boorbeschrijving te extraheren. Omdat het model is gemaakt voor de berekening van de deklaagweerstand zijn de verschillende zandklassen niet onderscheiden, omdat is aangenomen dat deze een dusdanig hoge verticale doorlatendheid hebben dat deze niet van belang zijn in het bepalen van de deklaagweerstand. De lithoklassen die zijn onderscheiden zijn veen, klei, zandige klei en zand (ongedifferentieerd). Het voxelmodel is mede gebaseerd op een lagenmodel van de lithostratigrafische eenheden in het gebied.

De basis van het model is de top van het Pleistocene zand. Vervolgens komt de Formatie van Nieuwkoop, Basisveenlaag (soms verslagen) voor en daarop de Formatie van Naaldwijk, Laagpakket van Wormer. Daarop is het Hollandveen Laagpakket van de Formatie van Nieuwkoop gevormd en daarna het laagpakket van Walcheren (Formatie van Naaldwijk). Tevens is de Almere Laag (Formatie van Naaldwijk, Laagpakket van Walcheren) gemodelleerd. De basis van iedere eenheid is in de boringen aangegeven en vervolgens is een vlak dekkend model van de basis van de eenheid berekend. De verschillende lagen zijn vervolgens “op elkaar gestapeld”, waardoor een consistent lagenmodel ontstaat. Binnen iedere laag zijn de voxels geselecteerd die in die laag vallen en is een schatting berekend van de kans op iedere lithoklasse d.m.v. 3D kriging interpolatie van de boorgegevens. De kansen op iedere lithoklasse in elke voxel is bewerkt om tot een meest waarschijnlijke lithoklasse te komen volgens het algoritme dat is beschreven in Soares (1992). De meest waarschijnlijke lithoklasse is gebruikt voor het Atlantis model. Verder zijn in het voxelmodel voor Atlantis de voxeldikten verlaagd naar 0,2 m.

39 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

De boringen die gebruikt zijn in de modellering zijn deels van vlak na de inpoldering. Lagen in de boringen zijn zeer waarschijnlijk gecompacteerd en veen is geoxideerd. Daardoor is de toenmalige beschrijving niet een accurate beschrijving van de huidige situatie van de bodemopbouw. In het voxelmodel dat in 2021 gemaakt is voor de deklaagweerstand is hiermee geen rekening gehouden. Het ondiepe deel van het model tot 1,6m diepte is daarom vernieuwd, gebaseerd op meer recente boringen. De verwachting is dat de lagen dieper dan 1,60 m niet in sterke mate gecompacteerd zijn en dat het model dieper dan 1,60 m een voldoende accurate weergave van de huidige opbouw geeft.

8.3

Productiewijze bodemdeel De gebruikte bodemgegevens zijn verzameld in de periode 2012 tot en met 2018. Ze zijn dus van recente datum. Voor processen als rijping en veenoxidatie is dat essentieel anders representeren ze niet de huidige situatie. De bodemgegevens bestaan uit 1461 profielbeschrijvingen die verricht zijn door ervaren veldbodemkundigen. Voor meer informatie over het bodemkundig onderzoek dat tot deze dataset heeft geleid wordt verwezen naar Brouwer et al. (2018). De bodemprofielen maken deel uit van de Basisregistratie Ondergrond (broloket.nl). De bodemprofielen worden opgedeeld in acht laagjes van 0,20 m dikte. Het bovenste laagje begint aan maaiveld, het onderste laagje eindigt op een diepte van 1,6 m. Van elk laagje is de lithoklasse bepaald, d.w.z., of het gaat om klei, veen, zand of zavel. Deze classificatie is uitgevoerd conform de richtlijnen in ten Cate et al. (1995). Van elk laagje is tevens de rijpingsklasse bekend, die varieert van 1 (ongerijpt) tot 5 (volledig gerijpt). Veen en zand hebben rijpingsklasse “0”, een code die aangeeft dat rijping niet van toepassing is. Andere bodemprofielkenmerken zijn: het gehalte aan organische stof (uitgedrukt als massapercentage), het lutumgehalte (uitgedrukt als massapercentage), en indien veen aanwezig is, de veensoort. Vervolgens worden deze bodemprofielen gebruikt om de bovenste 1,6 m van het voxelmodel te vullen.

Hiervoor hebben we gebruik gemaakt van een geostatistische methode die bekend staat als Bayesian Maximum Entropy (BME). BME berekent voor elke voxel de kans op een bepaalde bodemeigenschap gegeven de waargenomen bodemeigenschappen in de omliggende bodemprofielen. BME kan zowel gebruikt worden voor nominale variabelen, zoals lithoklassen, als continue variabelen, zoals het massapercentage lutum en het massapercentage organische stof. Deze methode biedt voordelen t.o.v. alternatieve ruimtelijke interpolatiemethoden. Er is bijvoorbeeld geen nabewerking nodig om valide kansen te krijgen. Voor alternatieve methoden, zoals indicator kriging (Journel, 1983), zijn dergelijke nabewerkingen vaak wel nodig. Daarnaast stelt BME minder restricties aan de modellen die bij het interpoleren mogen worden gebruikt en biedt het mogelijkheden om aanvullende informatie (zoals fysische principes) in te brengen. Nadeel van BME is dat de rekentijden aanzienlijk kunnen zijn. Voor meer informatie over BME zie Christakos (2000) en Bogaert (2002). BME resulteert in een kansverdeling voor elk attribuut. Voor nominale variabelen zoals de lithoklassen zijn dit discrete verdelingen, d.w.z., de kans op het voorkomen van elke lithoklasse. We kiezen dan de lithoklasse met de grootste kans op voorkomen.

Omdat we beschikken over de kansverdeling, kunnen we ook uitspraken doen over de nauwkeurigheid van de interpolatie. Voor continue variabelen zoals het massapercentage organische stof en lutum kiezen we voor de gemiddelde waarde volgens de kansverdeling.

40 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 13 geeft de geïnterpoleerde lithoklasse voor laag 1 en laag 6. Voor de overige lagen wordt verwezen naar de bijlage 19E. Figuur 14 geeft de bijbehorende onzekerheid van deze kaarten. We zien dat we het meest onzeker zijn bij overgangen in lithologie.

Figuur 13 De geïnterpoleerde lithoklasse voor laag 1 (links) en laag 6 (rechts) van het voxelmodel.

Figuur 14 Onze onzekerheid over de geïnterpoleerde lithoklasse voor laag 1 (links) en laag 6 (rechts) van het voxelmodel. De onzekerheid is hier uitgedrukt als een genormaliseerde entropie: 0 betekent dat we heel zeker zijn over de lithoklasse, 1 is dat we heel onzeker zijn.

Om de consistentie tussen lithoklasse, GLG, en rijpingsklasse te waarborgen is na interpolatie een nabewerking uitgevoerd: - Veen en zand krijgen altijd rijpingsklasse “0” (= niet van toepassing); - Voxels beneden GLG-niveau krijgen rijpingsklasse ”1” (ongerijpt); Hoewel in zeer droge jaren de grondwaterstand onder de GLG kan zakken wordt verondersteld dat dit geen extra rijping van de lagen onder de GLG veroorzaakt. - Klei- en zavelvoxels met rijpingsklasse “0” krijgen rijpingsklasse “5” (volledig gerijpt) als het de bovenste laag betreft (0 - 0.2 m-mv). Anders wordt de rijpingsklasse

41 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

afhankelijk gesteld van de GLG: de laag net boven GLG-niveau krijgt rijpingsklasse “1-3”, die daarboven “4”, en de lagen daarboven “5”. Figuur 15 geeft de geïnterpoleerde rijpingsklasse voor laag 1 en laag 6. Voor de overige lagen wordt verwezen naar bijlage 19E..

Figuur 15 De geïnterpoleerde rijpingsklasse voor laag 1 (links) en laag 6 (rechts) van het voxelmodel.

Figuur 16 geeft de geïnterpoleerde massapercentages lutum voor laag 1 en laag 6. Voor de overige lagen wordt verwezen naar bijlage 19E.

Figuur 16 De geïnterpoleerde massapercentages lutum voor laag 1 (links) en laag 6 (rechts) van het voxelmodel.

Figuur 17 geeft de geïnterpoleerde massapercentages organische stof voor laag 1 en laag 6. Voor de overige lagen wordt verwezen naar bijlage 19E.

42 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 17 De geïnterpoleerde massapercentages organische stof voor laag 1 (links) en laag 6 (rechts) van het voxelmodel.

In de bodemkunde en de geologie is het gebruikelijk om het lutumgehalte uit te drukken als massapercentage op de minerale delen. Voor Atlantis is echter het massapercentage lutum nodig uitgedrukt op de totale massa grond (dus inclusief organische stof). We rekenen het lutumgehalte uitgedrukt als massapercentage op de minerale delen (L) om naar het lutumgehalte uitgedrukt als massapercentage op de totale massa grond (L’) met de volgende expressie: L’ = L(1 - H/100) waarin H het massapercentage organische stof is. Het boven- en ondergrondmodel zijn gecombineerd door de bovenste 1,60 m van het ondergrondmodel te vervangen door die van het bodemmodel. Door de TNO lithoklassen en de WENR lithoklassen te aggregeren naar een kleinere set lithoklassen bestaande uit de hoofdgroepen klei, veen, zand en zavel proberen we abrupte overgangen tussen het bodemmodel en het diepere ondergrondmodel te reduceren. Voor het ondergrondmodel gebruiken we daarvoor de volgende aggregaties: veen (vnl Hollandveen) en veen op pleistoceen worden samengevoegd tot veen, klei blijft klei, zandige klei wordt zavel, en zand blijft zand.

De aggregatie van het bodemmodel is als volgt: veen, zandig veen, kleiig veen, venig zand, en venige klei worden samengevoegd tot veen, lichte klei en matig zware klei worden samengevoegd tot klei, kleiig silt, matig lichte zavel, zandige leem, zeer lichte zavel, en zware zavel worden samengevoegd tot zavel, en grind, kleiarm zand, kleiig zand, leemarm zand, sterk lemig zand, zeer sterk lemig zand, en zwak lemig zand worden samengevoegd tot zand.

43 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

8.4

NAP hoogte Het vervaardigde voxelmodel van de ondergrond is in eerste instantie geproduceerd t.o.v. maaiveld. In dit geval betekent dit dat het voxelmodel overal een vaste basis heeft van 10 meter onder maaiveld, waarbij iedere verticale kolom is opgebouwd uit een stapeling van 50 voxels met een dikte van 20 cm. Echter, voor het gebruik in Atlantis moet de hoogte in m NAP worden toegekend aan de voxels. Deze informatie is toegevoegd door de basis van het model (voorheen 10 m-mv) een absolute hoogte van in totaal 10 meter onder de maaiveldhoogte in m NAP te geven. Hiervoor is het 5x5 meter rasterbestand van het Actueel Hoogtebestand Nederland 4 (AHN4) gebruikt. Dit is verwerkt tot een 100x100 meter grid, door de mediane maaiveldhoogte van alle 5x5 meter cellen binnen een 100x100 meter cel, toe te wijzen aan die cel.

44 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

9

Kartering rijpingsgraad

9.1

Werkwijze

In Flevoland komen op diverse locaties ongerijpte of onvolledig gerijpte bodemhorizonten voor op relatief ondiepe bodemniveaus (boven de GLG). De rijpingsgraad is op twee manieren beschikbaar: 1. Eenvoudige veldtest; 2. Berekend op basis van een set van gemeten/geschatte bodemparameters die resulteren in de n-waarde (zie bijlage 19C). Uit vooronderzoek is gebleken dat de dataset (2) op basis van metingen erg beperkt is en niet bruikbaar voor een gebiedsdekkend beeld. Er is daarop alleen gebruik gemaakt van de veldmetingen. In 2018 is de bodemkaart vernieuwd (Brouwer et al, 2018). Voor het actualiseren van de bodemkaart zijn nieuwe boringen uitgevoerd en is gekeken of oude databestanden nog geschikt waren: • Voor de Noordoostpolder vormen de 477 nieuwe boormonsterbeschrijvingen, samen met 76 (alleen vanaf 2005) boorgegevens uit BIS Nederland en 3351 boorgegevens uit het DINO Loket (alleen gebruikt voor de bepaling pleistocene zanddiepte en maaiveldhoogte) de kalibratieset voor de modellen waarmee de ruimtelijke voorspellingen is uitgevoerd voor het deelgebied van de Noordoostpolder.

Validatie heeft plaatsvonden met 35 aparte boringen • Voor Oostelijk en Zuidelijk Flevoland vormen de 730 nieuwe boormonsterbeschrijvingen, samen met 45 boorgegevens uit BIS Nederland en 67 boorgegevens uit het DINO Loket de kalibratieset voor de modellen waarmee de ruimtelijke voorspellingen is uitgevoerd voor het deelgebied Oostelijk en Zuidelijk Flevoland. Validatie heeft plaatsgevonden met 77 aparte boringen • De nieuwe boringen hebben zich gericht op het verzamelen van actuele informatie over de laagopbouw van de ondergrond tussen 1,2 en 3,5 à 4 m –mv. met speciale aandacht voor de begindiepte van de pleistocene ondergrond en slappe kleilagen, en de dikte van het minerale dek en de veenlagen. De diverse puntdata leveren o.a. data over de diepteligging, de bodemkundige en geologische indeling, het type veen en de mate van rijping. Verder is van elke locatie bekend wat de maaiveldhoogte is en is een schatting gemaakt van de GHG en GLG. Deze data zijn te vinden op www.bodemdata.nl. Van alle puntdata is bekend of er klei voorkomt en in welke mate deze klei is gerijpt (Tabel 4). Om een verklaring te vinden voor het al dan niet aanwezig zijn van ongerijpte klei is o.a. gekeken naar landgebruik, kwel/infiltratie; GLG en ontstaansgeschiedenis. Tenslotte is een vergelijking gemaakt tussen de diepte waarop ongerijpte klei wordt aangetroffen en de op locatie gemeten GLG.

Tabel 4 Voorbeeld van puntdata met informatie over de aanwezigheid en diepte van ongerijpt klei.

45 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Tabel 5 rijpingsklasses, n-waarde en beschrijving rijping kleilagen

n-waarde 0,7 0.85 1,2 1,7 2,2 0

9.2

rijpingsklasse

5 4 3 2 1 0

beschrijving gerijpt vrijwel gerijpt half gerijpt vrijwel ongerijpt ongerijpt geen klei

Resultaten Bij circa 100 boorlocaties is in het profiel een bepaalde mate van ongerijpte klei aangetroffen. De rijpingsgraad is verdeeld in vijf klasses (Tabel 5). Per laag van 20 cm is per boorlocatie de mate van rijping aangegeven (zie Figuur 18 als voorbeeld, volledig overzicht in bijlage 19E). In het voxelmodel (hoofdstuk 8) zijn per voxellaag de puntdata geïnterpoleerd naar een vlakdekkend geheel. Om hier indicatie van te geven is de rijpingsklasse van het diepste ondiepe voxel (1,5 m-mv) weergeven (Figuur 19). De kleien van Flevoland zijn voor het grootste gedeelte gerijpt. Belangrijkste gebieden met nog ongerijpte gronden zijn: een gebied ten westen van Emmeloord; een gebied aan de zuidoostelijke rand van de Noordoostpolder, een gebied in het centrum van de Zuidelijk Flevoland en een gebied ten noorden van Almere. Het landgebruik, de kwel/infiltratie en de GLG geven niet een duidelijke verklaring voor het voorkomen van ongerijpte klei. De belangrijkste verklarende parameter is het afzettingsmilieu. Zo was er ten westen van Emmeloord een binnenmeer; deze is opgevuld met “Union klei” en met “Cardiumklei” (Figuur 20 Verspreiding Union klei (Wiggers, 1955)). Dit zijn (vrij zware) zeekleien en zouden daarom veel langzamer rijpen. In Zuidelijk Flevoland lijken de niet-gerijpte kleien met name voor te komen ter hoogte van kleiaanvulling van de Eembedding (Figuur 21).

Omdat de ongerijpte kleien zowel voorkomen in de oudste gedeeltes van Flevoland als in de nieuwere delen is het tijdstip van drooglegging niet een goede maat voor de mate van rijping of de rijpingssnelheid. Zo komen er in de Noordoostpolder 80 jaar na de drooglegging nog ongerijpte kleien voor. De rijpingsdata bieden daarmee ook onvoldoende basis om een goede inschatting te maken van de snelheid waarmee de ongerijpte kleien naar een volledige gerijpte toestand zullen zijn overgegaan. Mogelijk zullen bepaalde delen nooit volledig rijpen. Vooralsnog wordt uitgegaan van 90 jaar +/- 30 jaar; sterk afhankelijk van uitdrogingssituatie. Om te bepalen of de volcapillaire zone van invloed is op de mate van rijping is een vergelijking gemaakt tussen de diepte waarop ongerijpt (klasses 1-5, 0 = geen klei)) klei wordt aangetroffen en de GLG (Figuur 22). Het blijkt dat voor de niet-gerijpte lagen (1-4) gemiddeld gezien op of onder de GLG liggen. Dit betekent dat de volcapillaire zone geen bescherming lijkt te bieden tegen het rijpingsproces.

46 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 18 De verspreiding ongerijpte klei op 1,5 m -mv . 0 = geen klei, 1 is ongerijpt- 5 is volledig gerijpt (data boorlocaties)

47 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 19 Verspreiding ongerijpte klei op 1,5 m -mv . 0 = geen klei, 1 is ongerijpt- 5 is volledig gerijpt (geïnterpoleerd).

Figuur 20 Verspreiding Union klei (Wiggers, 1955)

48 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 21 Verspreiding oude zeeklei-afzetting (Menk et al, 1998)

49 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 22 box-and-whisker plots voor de Noordoostpolder (NOP) en Zuidelijk en Oostelijk Flevoland (ZOF) geven de relatie weer tussen de rijpingsklasse en de afstand tot de GLG. De horizontale lijn in het midden van elke box is de mediane waarde. De box zelf wordt begrensd door de kwartielen. De verticale lijnen geven de ‘spreiding’ van de gegevens, en de punten potentiële uitbijters (outliers) (0= geen klei, 1 = volledig ongerijpt, 5 = volledig gerijpt)

50 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

10

Kartering veentypen

10.1

Werkwijze Door diverse afzettingsmilieus is er in Flevoland een groot aantal verschillende veentypes ontstaan. Uit het BRO-archief zijn de boringen geselecteerd waar tot een diepte van 1,6 mmv een of meerdere horizonten veen- of moerige lagen voorkomen (n= 2337). Van deze horizonten is bepaald, welk type veen het is, wat de trofiegraad is (eutroof, mesotroof of oligotroof) en hoeveel klei wordt aangetroffen (veldschatting) (Tabel 6). Tabel 6 Karakterisering veenlagen tot 1.6 m-mv Flevoland (afkorting zie Tabel 7) Veensoort

51 van 114

Trofie

Aantal

Lutum (%) Gemiddeld

Min

Max

BA

eutroof

8

31

25

40

B

eutroof

27

28

12

40

BE

eutroof

75

21

5

30

DK

eutroof

66

23

4

40

VV

eutroof

1365

17

1

40

R

eutroof

19

22

10

40

GY

eutroof

97

19

15

22

GL

eu-/mesotroof

3

16

3

30

OV

eu-/mesotroof

26

17

1

30

CR

eu-/mesotroof

25

15

10

30

RC

meso-/eutroof

260

12

1

32

BM

mesotroof

127

14

1

30

C

mesotroof

161

10

1

30

DV

mesotroof

25

11

1

25

DZ

oligotroof

27

6

2

20

S

oligotroof

24

8

1

15

SP

oligotroof

1

12

12

12

WO

oligotroof

1

5

5

5

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Tabel 7 Veensoorten

Vervolgens zijn de diverse veensoorten ingedeeld in vijf klasses op basis van de trofiegraad (Tabel 8). De trofiegraad wordt als belangrijkste bepalende factor gezien voor de oxidatiesnelheid (zie ook 12.2.2)

Tabel 8 Indeling veensoorten naar veenklasses

10.2

Akorting veensoort

Veenklasse

BA", "B", "BE", "DK", "VV", "R", "GY":

eutroof

"GL", "OV", "CR", "RC":

eu-/mesotroof

BM", "C", "DV":

mesotroof

"DZ", "S", "SP", "WO":

oligotroof

"VV":

verslagen

Resultaten Diverse boorlocaties laten een groot aantal verschillende veensoorten zien. De veensoorten zijn verdeeld in vijf klasses (Tabel 8). Per laag van 20 cm is per boorlocatie de veenklasse bepaald. Figuur 23 en Figuur 24 geven een beeld van de veenklasse per datapunt op een diepte van 0,9 en 1,5 m-mv. In het voxelmodel (hoofdstuk 8) zijn per voxellaag de puntdata geïnterpoleerd naar een vlakdekkend geheel. Om het veentype in Flevoland in beeld te vatten is in Figuur 25 een geïnterpoleerde kaart weergeven van de veensoorten op 1,1 m-mv. Ondanks dat de veentypes in de diepte variëren kan deze figuur worden gebruikt om een algemeen beeld te krijgen van de veensoorten in de bovenste 1,6 m. De volledige kaarten per voxellaag van het veentype en geïnterpoleerde veensoorten zijn te vinden in bijlage 19E.

52 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Het grootste gedeelte van het veen in Flevoland bestaat uit verslagen veen. Uitzondering is Schokland, op dit voormalig eiland wordt een min of meer onverstoorde veenondergrond (in situ) bestaande uit broekvenen (eutroof) en zeggevenen (mesotroof) gevonden. Verder worden er in de volgende gebieden andere veentypes gevonden. Hierbij is niet altijd mogelijk een relatie te leggen met de (veelal complexe) ontstaansgeschiedenis: • Rand van Almere en Oostvaardersplassen. Hier wordt Glyde en Gyttja (eutroof) gevonden; • Rand zuidelijk Flevoland. Hier worden mesotrofe venen gevonden; • Midden-west Noordoostpolder. Hier worden zeggevenen (Mesotroof) gevonden; • Verspreid Noordoostpolder. De oligotrofe venen liggen in de buurt van ondiep gelegen dekzand (dekzandopduiking). Dieper dan 1,6 m-mv zijn beperkte boordata met een veentypering beschikbaar. Uit de beperkte diepe boringen blijkt dat meer naar de diepte er meer zegge- en broekveen wordt gevonden en dat het aandeel verslagen veen sterk afneemt. Zegge- en broekveen bevat meer organische stof en minder lutum dan de verslagen venen.

Figuur 23 Veentypering op 0,9 m-mv (puntdata)

53 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 24 Veentypering op 1,5 m-mv (puntdata)

54 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 25 Veentype op 1,1 m-mv (geïnterpoleerd)

55 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

11

Rekenwijze rijping

11.1

Inleiding Rijping vindt plaats bij het droogvallen (drainage/ontwatering) van klei-rijke sedimenten die onder water zijn afgezet. Luchtindringing in de bodem zet fysische en chemische veranderingen in gang waarbij het slappe sediment geleidelijk verandert in een stevige, doorlatende bodem met structuur. Fysische rijping gaat gepaard met een onomkeerbaar verlies aan water en een aanzienlijk volumeverlies, variërend van 10% bij lichte zavel tot wel 30% bij zware klei (Jongmans et al. 2013). Dit wordt ook wel (onomkeerbare) krimp genoemd. Wateronttrekking door begroeiing versnelt het rijpingsproces aanzienlijk. Door het volumeverlies ontstaan rijpingsscheuren in de bodem en treedt maaivelddaling op. In de Flevopolders was rijping (krimp) een hoofdoorzaak van de bodemdaling door de drooglegging van de waterbodems van het IJsselmeer. Grote delen van het minerale dek zijn inmiddels sterk gerijpt waardoor een verdere bijdrage aan bodemdaling afwezig of zeer beperkt is. Echter, op diverse locaties in Flevoland komen desondanks nog ongerijpte of onvolledig gerijpte bodemhorizonten voor op relatief ondiepe bodemniveaus (boven de GLG). Dat betekent dat op deze locaties ook nu nog rekening moet worden gehouden met een rijpingsbijdrage (restrijping).

Bovendien is klei en zavel beneden de diepste grondwaterstand over het algemeen niet aan uitdroging blootgesteld, wat betekent dat in scenario’s waarbij de grondwaterstand (GLG) zou worden verlaagd, rijping weer zal gaan bijdragen. Er is daarom besloten om rijping als apart bodemdalingsproces mee te nemen in de bodemdalingsprognose voor Flevoland (zie uitgangspunt 2 in hoofdstuk 2). Omdat bij aanvang van dit project rijping in Atlantis ontbrak, is hiervoor een nieuwe rekenwijze opgesteld.

11.2

Uitgangspunten De volgende uitgangspunten zijn gebruikt bij het opstellen van de rekenwijze. •

De rekenwijze maakt gebruik van de GLG als de algemene stuurvariabele voor de karakterisatie van grondwatercondities in Atlantis. • Als toestandsvariabele die de mate van rijping karakteriseert word gebruik gemaakt van het bestaande concept n-factor (rijpingsgraad of rijpingsgetal) (Pons en Zonneveld, 1965). Zie bijlage 19B.. • Rijping is compleet als een laagste limietwaarde voor de n-factor is bereikt. Er kan dan geen verdere rijping meer plaatsvinden. • Rijping is onomkeerbaar. • Rijping is een proces (volumeafname van voxels) dat ook bij gelijkblijvende grondwatercondities (hier GLG) plaatsvindt14. Het alternatief zou zijn dat een forcering (GLG verandering ten opzichte van het bodemmateriaal) nodig is om volumeverlies door rijping te bewerkstelligen. Hoewel er geen zekerheid bestaat over de juiste zienswijze, pleiten de volgende overwegingen voor het gekozen uitgangspunt: - Het feit dat rijpingsgraad globaal toeneemt met de ouderdom van diverse andere, oudere polders, zoals in de Dollard. - Bij de drooglegging van Flevoland was natuurlijk sprake van forcering door de drooglegging en ontwatering. Vervolgens was echter ook sprake van een —————————————— 14

56 van 114

Dit wordt ook wel een seculair proces genoemd.

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

-

11.3

seculaire ontwikkeling; rijping als naijling die in ieder geval decennia lang door heeft gewerkt, en mogelijk nog steeds gaande is. Bij bodemkartering rond 2000 is geconstateerd dat in Noord- Zuid Holland (De Vries et al, 2017) de rijpingsgraad van bodemhorizonten was toegenomen in vergelijking met de karting in de 60’ en 70 er jaren. Dit is in ieder geval niet strijdig met een seculair karakter van de rijping, maar forcering door veranderd landgebruik en ingrepen in de waterhuishouding zullen hier ook aan hebben bijgedragen en mogelijk ook dominant zijn geweest.

Rekenwijze De afname rijpingsgraad n wordt beschreven als exponentiële relatie (met een e-folding tijdschaal  die bepaalt hoe snel de eindtoestand wordt ‘bereikt’)

dn ngerijpt − n = dt 

(1)

Bij ongewijzigde omstandigheden volgt dan bij een startwaarde nstart op t = tstart

n(t ) = ngerijpt + (nstart − ngerijpt )  e − (t −tstart )/

(2)

Voor een tijdstap van eindige grootte  t van t naar t + t :

n(t + t ) = n(t ) + ( n(t ) − ngerijpt ) ( e−t / − 1)

(3)

Een afname van rijpingsgraad voor een tijdstap of stress periode die volgt uit (1) wordt vertaald naar een relatieve volumeafname (= volumerek) van een voxel (of deel van een voxel dat binnen het rijpingsbereik ligt). V / Vref = factor  n

(4)

Afname levert hier een negatief getal. Voor Vref wordt uit praktisch oogpunt gekozen voor het grondvolume bij de rijpingsgraad bij de start van de berekening, n0 . Daardoor kunnen de

volumeveranderingen V voor een tijdstap worden berekend uit het start volume van een voxel. In de implementatie is ervoor gekozen om zowel het referentievolume als n0 elke tijdstap te updaten15. De factor voor het model van Pons en Zonneveld (1965) (afleiding in bijlage 19C):

factor =

( L + bH ) / w n0 ( L + bH ) / w + 0,2R / w + C

(5)

R = 100 − L − H C = R / R + H / H + L / L

—————————————— Een belangrijke reden is dat voxels gedurende de berekening kunnen worden ‘opgeknipt’ wanneer een deel ervan geen rijping moet ondergaan; daardoor ontstaat er een nieuw referentievolume. Ook kan het volume afnemen door andere processen als oxidatie van organisch materiaal in organisch-rijke klei. 15

57 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Benodigde parameters zijn:

(massapercentage t.o.v. vaste stof) L, H  w ,  L ,  H ,  R (dichtheid water en vaste stof componenten) (startwaarde rijpingsgetal) n0 Voor de relatieve hoogteverandering van een voxel (of deel van een voxel dat binnen het rijpingsbereik ligt) is verder een geometrie factor r nodig 1

H  V  r = 1 +  −1 H  V 

(6)

Hierin geeft r = 3 uniforme krimp in alle richtingen (isotroop) en r = 1 puur verticale krimp (Rijniersce (1983)) . Er wordt een ondergrens voor het rijpingsbereik in de bodem gedefinieerd dat een hoogte Hv0,rijping boven de GLG ligt (vergelijkbaar met de rekenwijze bij oxidatie). Om te kunnen rekenen moeten waarden voor n, L en H bij aanvang zijn toegekend aan alle voxels. Deze parametrisatie wordt toegelicht in Hoofdstuk 15. Omdat voxels die bij aanvang beneden het rijpingsbereik liggen op een later tijdstip in principe na verloop van tijd binnen het rijpingsbereik kunnen komen en vervolgens aan indroging en rijping worden blootgesteld, moeten alle voxels een n-waarde worden toegekend.

58 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

12

Rekenwijze oxidatie

12.1

Rekenwijze De bestaande rekenwijze in Atlantis voor oxidatie is gebruikt16 (Erkens e.a., 2021). De methode is gebaseerd op de langjarige waarnemingen aan maaivelddaling door veenoxidatie op de locatie van het veenweide innovatiecentrum (VIC) te Zegveld en een aantal andere locaties in Nederland (van den Akker e.a. 2007). Op basis van deze metingen zijn door van den Akker lineaire relaties opgesteld tussen de maaivelddalingssnelheid en GLG. Dat is gedaan voor veen zonder kleidek, veen met een kleidek, en voor veen met kleidek met een toepassing van een correctie voor kleidekdikte. In de methode is de relatie tussen de GLG en maaivelddalingssnelheid door oxidatie omgezet naar een massa-gebaseerde oxidatiesnelheid (in massa organisch materiaal dat wordt geoxideerd per volume bodemmateriaal per tijdseenheid). De hoeveelheid daling die een voxel kan leveren door oxidatie van organisch materiaal wordt bepaald door het organisch stofgehalte van de voxel (massa organische stof per massa droge stof). Dit maakt dat oxidatie niet alleen wordt toegepast op veenlagen, maar op alle lagen waarvoor een organisch stofgehalte groter dan een gedefinieerde grenswaarde aanwezig is. Verder is daardoor de volumeverandering voor een organisch-rijke klei geringer dan voor een schoon veen bij eenzelfde hoeveelheid aan geoxideerd organisch materiaal.

De berekeningen betreffen bodemdaling door aerobe oxidatie, dus door invloed van zuurstofindringing in de bodem. Voor de kwantitatieve beschrijving van de methode wordt verwezen naar Erkens e.a. (2021; hoofdstuk 6). De belangrijkste parameters die worden gebruikt in de methode zijn:



oxidatiesnelheid/afbraaksnelheid

kg/m3/dag

Voor de nationale bodemdalingsvoorspellingskaarten is hiervoor één waarde gehanteerd van 4,6x10-3 kg/m3/dag die hoofdzakelijk is bepaald door bovengenoemde waarnemingen uit het veenweidegebied bij Zegveld.

H v0,oxidatie

hoogte boven de GLG zonder oxidatie

m

Het dieptebereik in de bodemkolom waarbinnen oxidatie wordt berekend strekt zich niet uit van maaiveld tot de GLG, maar tot een ondieper niveau dat daar H v0,oxidatie boven ligt. Voor de nationale bodemdalingsvoorspellingskaarten is hiervoor één waarde van 0,3 m gehanteerd. Deze waarde is hoofdzakelijk bepaald door bovengenoemde waarnemingen uit het veenweidegebied bij Zegveld. De parameter vertegenwoordigt naar verwachting o.a. invloeden van capillaire opstijging die zorgen voor een hoge verzadigingsgraad boven de GLG en een sterke beperking van zuurstofindringing. max org

maximale bulk dichtheid organische stof (schoon veen)

kg/m3

Deze parameter bepaalt in belangrijke mate het volumeverlies van voxels waarin oxidatie optreedt. Een hogere dichtheid, bijvoorbeeld voor meer samengedrukt veen, betekent minder volumeverlies. Voor de nationale bodemdalingsvoorspellingskaarten is één waarde —————————————— Deze methode is in eerdere rapporten aangeduid met ‘voorraad-voxel methode’, omdat de methode expliciet de boekhouding van het organisch materiaal in een voxel bijhoudt. 16

59 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

gehanteerd van 100 kg/m3 die ongeveer overeenkomt met de gemiddelde waarde van een grote nationale dataset aan dichtheidsmetingen (Figuur 26).

Figuur 26 Organische stof dichtheidsmetingen als functie van organisch stof gehalte voor een grote nationale dataset (Erkens e.a., 2016). In Atlantis wordt voor de nationale bodemdalingsprognosekaarten een waarde van 100 kg/m3 gebruikt voor schoon veen (organisch stofgehalte van 100%). Dit komt ongeveer overeen met de gemiddelde waarde aangegeven met de continue horizontale lijn (a) in de grafiek.

12.2

Onderzoek naar gebiedsspecifieke parametrisatie max ) zou Er is onderzocht of de parametrisatie voor de drie parameters (  , H v0,oxidatie ,  org

kunnen worden aangepast voor de specifieke Flevolandse situatie. Daarbij is o.a. gekeken naar mogelijkheden voor ruimtelijke differentiatie van  op basis van veentype en de dikte van het minerale dek (uitgangspunt 2 in hoofdstuk 2). 12.2.1

De invloed van veentypen Veentypen De standaardwaarde voor  in Atlantis wordt is hoofdzakelijk bepaald door waarnemingen in Zegveld. Zegveld kent overwegend veenprofielen met een dun (15-40 cm) kleidek. Het kleidek is gemiddeld matig zware klei en het onderliggende veen bestaat veelal uit bosveen (eutroof) (soms broekveen). Bosveen heeft een organisch stofgehalte van 50-70% o.s. en ongeveer 10-30% lutumbijmenging. Bosveen is meestal slibrijk. Denk aan zacht ooibos dat regelmatig kon worden overstroomd door rivierwater. Hieronder worden resultaten besproken van een kort onderzoek rond de vraag of de venen die in Flevoland worden aangetroffen anders zijn dan die in Zegveld en in hoeverre daarvoor andere oxidatiesnelheden zouden gelden. Oxidatiesnelheid Op basis van de analyse van de diverse veentypes in Flevoland is een vergelijking gemaakt met het veen in Zegveld (Zie ook Hoofdstuk 10). De venen in Flevoland in de ondiepe ondergrond (tot 1,6 m-mv) zijn sterk afwijkend van het veen in Zegveld. Slechts op enkele locatie worden bos- en/of broekvenen aangetroffen. 60 % van de veenhorizonten is verslagen veen (gesedimenteerd veen); waarschijnlijk heeft dit verslagen veen een veel hogere bulkdichtheid dan andere veensoorten. 70 % van het veen is eutroof, oligotroof veen komt nauwelijks voor. Dit betekent dat de trofiegraad vergelijkbaar is

60 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

met Zegveld. De trofiegraad is duidelijk gecorreleerd aan het lutumgehalte. Hoe hoger de trofie, des te hoger is meestal ook het lutumgehalte. Hier moet wel worden opgemerkt dat aan de exacte getallen weinig waarde kan worden gehecht omdat dit schatting uit het veld zijn en geen labmetingen. In veel profielen is het veen aan de bovenkant veraard (DK, DZ en DV) (mondelinge mededeling Brouwer, 2022). In het kader van het Nationaal Onderzoeksprogramma Broeikasgassen Veenweide (NOBV) is onderzoek gedaan naar de potentiële oxidatiesnelheid (basal respiration rate) van diverse veensoorten (Figuur 27). Hierbij is gekeken naar veen dat voorkomt in het volledige geoxideerd gedeelte, het gedeelte rond de grondwaterspiegel en het gedeelte onder de grondwaterspiegel (anaeroob). Opvallend is dat de oxidatiesnelheid tussen de diverse veensoorten weinig lijkt te verschillen. De enige uitzondering zijn de meer veraarde venen in Vegelinsoord waar een hogere oxidatiesnelheid wordt gemeten. De in Flevoland veel voorkomende verslagen venen zijn in het kader van NOBV niet onderzocht. Het vermoeden is dat de verslagen venen een lagere oxidatiesnelheid en een hogere dichtheid hebben. Doordat metingen ontbreken kan dit niet worden gestaafd. Op basis van het bovenstaande kan geen onderbouwd verschil in oxidatiesnelheden tussen de diverse venen worden vastgesteld.

Voor prognose 2023 voor Flevoland is daarom geen ruimtelijke differentiatie voor type veen toegepast. Wel is het duidelijk dat de venen in Flevoland sterk afwijkend zijn van de venen die voorkomen in de veenweidegebieden. Om te weten hoeveel invloed een verschil in oxidatiesnelheid kan hebben op de bodemdaling is dit meegenomen in de gevoeligheidsanalyse. Samen met de vervaardigde veentype kaart (paragraaf 10.2)(eutroof veen, mesotroof veen en verslagen veen (ongedifferentieerd)) geeft dit een beeld van de mogelijke implicaties voor de prognose.

Figuur 27 Basal respiration voor diverse veensoorten in het oxische bereik. KR is kragge, RIB is ribbe in Kalenberg is een “natuurlijk” veensysteem met groeiend Riet-Sphagnum vegetatie ( deels met dopheide). MP is een maatregel plot dus onderwaterdrainage of andere watermaatregel. RF is referentie om te vergelijken met de MP plot. Moerige gronden op Vegelinsoord is een spalterveen laag van 20-30 cm dik. Demmerik is een toemaakdek op bosveen Andere bosveen sites zijn: Zegveld, Langeweide. Rouveen is een puur ( heel mooi) zeggeveen in een kwellocatie. Aldeboarn in een oligotroof veen met heide, veenmos en wollegras dus van oorsprong regenwater afhankelijk maar onder een klei-dek van 40 cm Assendelft is een zegge-rietveen met zoutwater invloeden (heel sulfide rijk mogelijk veel anaerobe afbraak en dan stimulans van afbraakproducten door toediening zuurstof). (bron: Weidner et al, 2023).

61 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

12.2.2

Invloed minerale dek Een mineraal dek heeft invloed op de aeratie van het er onderliggende veen. Zuurstofdiffusie is in (natte) klei lager dan in zand. Echter door uitdroging in de zomer gaat klei scheuren en kan er makkelijk zuurstof naar beneden worden getransporteerd. De dikte van het minerale dek speelt ook een rol. Er wordt dus verwacht dat deze factoren van invloed zijn op de oxidatiesnelheid van veen onder een mineraal dek. Literatuur Diverse onderzoeken wijzen erop dat de totale hoeveelheid bodemdaling door veenoxidatie in aanwezigheid van een kleidek vooral wordt beperkt door het feit dat het kleidek zelf niet (of nauwelijks) oxideert (o.a. Van den Akker e.a. 2007). Er zijn uit de literatuur echter geen aanwijzingen dat het onderliggende veen minder bijdraagt (in vergelijking met wanneer het kleidek niet aanwezig is). Goede info over de invloed van het zuurstoftransport van een mineraaldek op de oxidatiesnelheid ontbreekt. Voorprognose 2023 voor Flevoland is daarom geen ruimtelijke differentiatie in de oxidatiesnelheid  toegepast in verband met het minerale dek. Kartering veraard veen in relatie tot het minerale dek Op basis van beschikbare boorbeschrijvingen van de bodemkaart is de begindiepte waarop veraard veen wordt beschreven gekarteerd en vergeleken met de dikte van het minerale dek. (Figuur 28). Daaruit wordt het volgende opgemaakt: • Tot ca.

80 cm biedt het mineraal dek weinig tot geen bescherming tegen oxidatie. Dit houdt waarschijnlijk verband met het feit dat de veelal kleiige deklaag, mede door de sterke ontwatering, sterk is gescheurd waardoor zuurstof makkelijk tot grote diepte kan doordringen. • Een uitzondering lijkt het gebied ten noordwesten van Emmeloord waar het minerale dek vrij dun is (de GLG over het algemeen dieper ligt dan 110 cm -mv) en desondanks vrijwel geen veraard veen wordt beschreven.

62 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 28 Begindiepte veraard veen (op boorlocaties) geprojecteerd op de dikte van het minerale dek in Flevoland.

12.2.3

Hoogte boven de GLG zonder oxidatie De parameter H v0,oxidatie bepaalt tot welke hoogte boven de GLG geen oxidatie wordt berekend. De standaardwaarde van 0.3 m die wordt gebruikt in Atlantis is geschat op basis van correlaties van GLG- en langjarige bodemdalingsmetingen (door oxidatie) in Zegveld. Vergelijkbare metingen zijn niet beschikbaar voor Flevoland. Er is daarom gekeken welk inzicht kan worden verkregen uit (1) schattingen van water- en luchtgehalte van de bodem boven de grondwaterstand (GLG) en (2) de diepte tot waar veraard veen wordt aangetroffen in Flevoland in vergelijking met de GLG. Water- en luchtgehalte van de bodem ten opzichte van de GLG In het kader van het Spaarwater-project is door Acacia (2019) voor een klei op veenprofiel gemeten aan bodemvocht in relatie tot de afstand tot de grondwaterspiegel. Onder de aanname dat het profiel hydrostatisch is, kan de hoogte boven de grondwaterspiegel gelijkgesteld worden aan de zuigspanning (Figuur 29) en volgt dat de bodem op ca. 30 cm boven de grondwaterspiegel een waterverzadiging heeft van ca. 85 % (dit wordt dan de volcapillaire zone genoemd). Dit is deels te verklaren doordat de planten met name vocht onttrekken aan de topkleilaag en niet uit de onderliggende veenlaag.

63 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 29 Gemeten bodemvochtprofielen Spaarwater (Acacia, 2019)

Met behulp van de BOFEK-systematiek3 (https://www.wur.nl/nl/show/bodemfysischeeenhedenkaart-bofek2020.htm) kan worden bepaald welk gestandaardiseerde pF-curve past bij welk bodemprofiel. In Zegveld is dit de BOFEK-bouwsteen O17. Op basis van dat standaard profiel kan dan worden vastgesteld dat bij Zegveld er bij een drukhoogte van -30 cm 4 % lucht in het profiel voorkomt; bij – 70 cm is dat 10 % (Figuur 30). Om een indicatie te krijgen of het lucht-/zuurstofprofiel sterk afhankelijk is van de bodemopbouw, is het verloop ook berekend voor een ander profiel uit Flevoland. Dit profiel (BOFEK-bouwsteen O2) bestaat uit een zandige toplaag met een moerige tussenlaag en daaronder een zwaklemig zand. Dit profiel geeft aan dat bij een luchtgehalte van 4 %, de drukhoogte ca. – 35 cm bedraagt (Figuur 31). Dit betekent dat dit profiel tot iets hoger in het profiel verzadigd is.

Figuur 30 Verloop lucht- en watergehalte Zegveld-profiel (O17)

64 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 31 Verloop lucht- en watergehalte Zandprofiel (O2)

In het bovenstaande rekenvoorbeeld is er op basis van Zegvelddata vanuit gegaan dat er geen oxidatie optreedt bij een luchtpercentage lager dan 4 %. Er is echter niet gemeten aan het daadwerkelijk luchtgehalte en of dit al dan niet oxidatie veroorzaakt. Stel dat de grens niet op 4 % maar op 2 % ligt dan daalt de drukhoogte voor het Zegveld-profiel tot – 15 cm. Dit betekent dat H v0,oxidatie halveert en daarmee dus veel meer veen zou kunnen oxideren. Hierbij moet ook worden meegenomen dat in werkelijkheid de grondwaterstand in de zomerperiode tot onder de GLG (de drie laagste jaarlijkse grondwaterstanden gemiddeld over acht jaar) kan uitzakken waardoor het luchtgehalte op het betreffende niveau boven de GLG gedurende bepaalde perioden ook nog wat hoger kan zijn. In het kader van het NOBV is op de locatie Vlist (klei op veen) gemeten aan de zuurstofconcentratie in het profiel. Hieruit blijkt dat de zuurstofconcentratie tot op ca. 25 cm boven de gemeten grondwaterstand nog vrijwel nul is (Figuur 32). Wel is te zien dat de grondwaterstand bij droge zomers tot onder de GLG kan zakken en dat daarmee ook de zuurstofgrens op bepaalde momenten veel dichter bij de GLG kan komen te liggen. Op basis van indicatieve berekeningen lijkt er niet veel onderscheid te zijn tussen de dikte van de zone waarin de diverse profielen nog vrijwel volledig verzadigd zijn.

Op basis van berekening en metingen lijkt een volcapillaire zone van 25/30 cm een redelijke inschatting. Gezien de onzekerheid m.b.t. een passende waarde voor H v0,oxidatie is deze parameter meegenomen in de gevoeligheidsanalyse.

65 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 32 Gemeten zuurstofconcentratie en grondwaterstanden en de GLG en GLZ (gemiddelde diepte waarop er geen zuurstof wordt aangetroffen) locatie Vlist (klei (40 cm) op veen). (bron: Wageningen Environmental Research, niet gepubliceerd)

Diepte van veraard veen in Flevoland t.o.v. de GLG Op basis van beschikbare boorbeschrijvingen van de bodemkaart is een kaart gemaakt van de diepte tot waarop veraard veen wordt beschreven in de boringen en de GLG (Figuur 33). Daaruit wordt het volgende opgemaakt: • Met name in Zuidelijke en Oostelijk Flevoland en rond Schokland wordt tot dicht (040 cm) boven de GLG veraard veen aangetroffen. • Dit is minder zichtbaar voor de Noordoostpolder (NOP). In het noordelijke gedeelte van NOP houdt dit verband met het feit dat het veen alleen heel hoog in het profiel wordt aangetroffen en dus niet dieper veraard kan zijn. In het middendeel van de NOP zit de GLG relatief diep (> 150 cm), wat kan bijdragen aan een grotere afstand tussen GLG en de veraardingsdiepte (in Atlantis wordt een maximale diepte van oxidatie gehanteerd van 120 cm onder maaiveld, ongeacht GLG en H v0,oxidatie ). •

De kaart levert indicaties dat H v0,oxidatie kleiner dan 30 cm kan zijn.

Gezien de onzekerheid m.b.t. een passende waarde voor H v0,oxidatie is deze parameter meegenomen in de gevoeligheidsanalyse.

66 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

Figuur 33 Verschil tussen GLG en onderkant veraard veen in boorbeschrijvingen

12.2.4

Organische stofdichtheid veen Analyses van veen gerapporteerd in het kader vanprognose 2012 (Grontmij rapport 2012; bijlagen 4 t/m 6) bevatten indirecte informatie over de organische stofdichtheid van veen in Flevoland. De tabellen bevatten data voor 28 locaties. Eerst een selectie gemaakt van analyses met organisch stofgehalte > 20% waarvoor ook een bulk dichtheid bepaling aanwezig is zodat organische stofdichtheid berekend kan worden. Dit leverde resultaten op voor 9 locaties (4 in NOP, 4 in OF en 1 in ZF). De organische stofdichtheidsbepalingen zijn weergegeven Figuur 34 ten opzichte van het organisch stofgehalte. Uit deze beperkte data wordt het volgende opgemaakt: • De spreiding van de dichtheid is zeer vergelijkbaar met die van de bredere nationale dataset van Erkens et al. (2016).

67 van 114

Bodemdalingsprognose Flevoland 11208645-002-BGS-0003, 22 november 2023

•

•

Er is geen basis om ruimtelijk te differentiëren binnen Flevoland. Ook zijn er geen aanwijzingen voor een relatie met het bovengewicht (o.a. dikte minerale dek) of diepte onder maaiveld (analyse niet getoond). De standaardwaarde van 100 kg/m3 die is gehanteerd in Atlantis is ook voor Flevoland een passende waarde.

Figuur 34 Organische stofdichtheid van veen in Flevoland als functie van het organisch stofgehalte.

12.3

Verwerking van de temperatuurontwikkeling in het klimaatscenario Aangezien de afbraaksnelheid van organische materiaal (oxidatiesnelheid) gecorreleerd is met (bodem)temperatuur, zal een structurele stijging van de gemiddelde luchttemperatuur als gevolg van klimaatverandering ook de oxidatiesnelheid beïnvloeden. De oxidatiesnelheid (𝛼) wordt in de prognoses daarom stapsgewijs gecorrigeerd voor de verwachte temperatuurstijging uit KNMI’14 scenario GL. Hiervoor is gebruik gemaakt van de door Hendriks en Vermeulen (1997) ontwikkelde empirische relatie tussen de oxidatiesnelheid van organisch materiaal in Nederlandse veenbodems en de (lucht)temperatuur op basis van laboratoriumexperimenten. De empirische vergelijking toont een exponentiële ontwikkeling van de CO2-uitstoot met de (lucht)temperatuur: De relatieve oxidatiesnelheid : 𝑘𝑟𝑒𝑙 = 0.264 𝑒 0.127𝑇 ; hierin is T de bodetemperatuur. Een aantal relevante beperkingen van de studie van Hendriks en Vermeulen zijn: • De monsters zijn enkel blootgesteld aan 6 verschillende temperaturen (19,5; 10,5; 5,5; 3,2; 1,7 en 0°C) in een gecontroleerde opstelling. Dit betekent dat het toepasbaarheidsbereik van de relatie tussen 0°C en 19,5°C ligt. • De empirische vergelijking is gebaseerd op 8 veenmonsters van verschillende dieptes afkomstig van drie locaties in Nederland.