Lisa 3

JÄÄKSOODE VEEREŽIIMI TAASTAMISE KOMPLEKSUURINGU LÕPPARUANNE 1. PROJEKTI KESTUS Algus

(kuu/aasta): 24.04.2017 Lõpp:

(kuu/aasta) 01.09.2023

2. PROJEKTI TAOTLEJA (teadusasutus): Tartu Ülikool

Telefon: +372 7 375826

Aadress: Ülikooli 18, 50090 Tartu

Registrikood: 74001073

Panga rekvisiidid: SEB Pank AS, Tornimäe 2, 15010 TALLINN, arvelduskonto (IBAN): EE281010102000234007 , SWIFT/BIC: EEUHEE2X , käibemaksukohustuslase nr (VAT number): EE100030417 , tehingupartneri kood (TP kood): 605201 3. PROJEKTI JUHT: Ain Kull

(Ees- ja perekonnanimi) kaasprofessor, PhD (Amet, teaduskraad)

4. PROJEKTI PÕHITÄITJAD ARUANDEPERIOODI VÄLTEL Projekti põhitäitjad:

Ees- ja perekonnanimi Teaduskraad Ametikoht 1. Ain Kull PhD loodusgeograafia kaasprofessor 2. Valentina Sagris PhD geoinformaatika teadur 3. Edgar Karofeld PhD rakendusökoloogia kaasprofessor 4. Kai Vellak PhD taimeökoloogia kaasprofessor 5. Alar Läänelaid PhD maastikuökoloogia emeriitdotsent 6. Gert Veber PhD loodusgeograafia teadur 7. Marko Kohv PhD rakendusgeoloogia teadur 8. Mae Uri Dipl./BSc spetsialist (keemik) 9. Edgar Sepp MSc, doktorant geoinformaatika spetsialist 10. Martin Maddison PhD keskkonnatehnoloogia kaasprofessor 11. Ivika Ostonen-Märtin PhD juureökoloogia professor 12. Kristina Sohar PhD loodusgeograafia teadur 13. Iuliia Burdun PhD doktorant, kaitsnud PhD 2020 14. Tauri Tampuu PhD doktorant, kaitsnud PhD 2022 Projektiga seotud abitööjõud: 1. Kärt Erikson BSc/MSc magistrant (2023) / doktorant (2023 sept.) 2. Birgit Viru MSc/PhD doktorant, kaitsnud PhD 2020 5. PROJEKTI KULUD ARUANDEPERIOODIL 2023.a. 69267,18 eurot

Kokku

Töötasud (põhitäitjad +abitööjõud) 40702.70 Sotsiaalmaks 13408.58 Töötuskindlustusmaks 325.08 Ostetud teenused 4879.82 Lähetuskulud 4750.46 Materjalid, tarvikud, masinad, seadmed 5087.24 Muud kulud 113.30 Kokku 69267.18

Ostetud teenuste selgitus 4879.82 Mulla- ja veekeemia analüüsid biogeokeemia

laboris Lähetuskulude selgitus 4750.46 Kõik lähetused on seotud uurimisaladel

gaasi- ning veeproovide regulaarse kogumisega, drooniseire ja taimkatteseirega

Materjalide, tarvikute, masinate ja seadmete selgitus

5087.24 Fotosünteetiliselt aktiivse kiirguse mõõtmise andurid ja temperatuuriandurid. Soetati mõõteseadmetele patareisid ja akusid, seadmete hooldusmaterjale, mõõdulinte, teipe jmt. tarvikuid

Muude kulude selgitus 113.3 Kummikud, töökindad välitöödeks. Kulurida ei kajasta Tartu Ülikooli üldkulueraldist (20%) RMK-lt 2023.a. esitatud aruannete eest (arvestuslik summa 16328.22), mis kajastub eelarves pärast aruande heakskiitmist ja lepingutasu laekumist tartu Ülikoolile.

6. PROJEKTI TÄITMISE LÕPPARUANNE Rakendusuuringu „Ammendatud turbamaardlate vee-režiimi taastamise kompleksuuringu metoodika väljatöötamine ja uuringu läbiviimine“ eesmärgiks oli perioodil 2017 – 2023 luua jääksoode seisundi ja korrastamisjärgsete muutuste seiramise metoodika, rajada valimisse kuuluvas viies jääksoos seirealad ning viia läbi kogu perioodi hõlmav kompleksseire. Lõpparuandes antakse ülevaade projekti raames 2017.a. aprillist kuni 2023.a. septembrini Laiuse, Kõima, Maima, Kildemaa ja Ess-soo jääksoodes läbi viidud tegevustest ja esmastest tulemustest ning tuuakse välja peamised seire käigus tehtud tähelepanekud korrastamistööde edukust mõjutavatest teguritest. Koondaruandes käsitletud teemade detailsem analüüs (eeskätt metoodika osas) on esitatud aruande lõpus viidatud lisana esitatud artiklites ja lisamaterjalides. Seirealade rajamine, seire ja korrastamistööde ajajoon Eelneva ruumianalüüsi ning välitööde tulemuste põhjal valiti RMK poolt korrastatavate jääksoode hulgast seiratavateks aladeks Laiuse, Ess-soo, Maima, Kõima ja Kildemaa jääksood. Neist Maima ja Ess-soo moodustavad väga sarnase paari, kus on esindatud mahajäetud freesturbaväljad ja iseseisvalt taimestunud turbavõtuaugud ning nende vahel kuivemad metsastuvad tervikud ning Ess-soos ka freesturbavälja laiendamiseks eelkuivendusega kuid algse rabataimestiku eemaldamiseta ala. Selle paari puhul oli eesmärgiks korrastamistööde käigus võrrelda Lääne-Eesti ja Kagu-Eesti erinevusi (ilmastik, pealiskorraga seotud mõjutused turba- ning veeomadustele) järgnevate töötlustega aladel:

a) võrdlusalad (korrastamistööde käigus mõjutamata veerežiim ja taimestik); b) madalaveeline veekogu loodusliku taimestumisega; c) alad pinnaspaisudega stabiliseeritava veerežiimiga, kus taimestik areneb iseseisvalt; d) alad lausalise kraavide täitmisega stabiliseeritava veerežiimiga, kus taimestik areneb iseseisvalt; e) alad pinnaspaisudega stabiliseeritava veerežiimiga, kus turbasambla fragmentide siirdamisega

kiirendatakse taimestumist; f) alad lausalise kraavide täitmisega stabiliseeritava veerežiimiga, kus turbasambla fragmentide

siirdamisega kiirendatakse taimestumist. Laiuse jääksoo puhul oli kobraste tegevuse tulemusena lõunapoolses jääksoo osas kujunenud madalaveeline veekogu, põhjapoolses osas aga mahajäetud freesturbaväljal suhteliselt suure ida- läänesuunalise kõrgusgradiendiga vähe kuni mõõdukalt taimestunud ala. Korrastamistööde käigus säilitati loodusliku taimestumise teel soostuv veekogu, põhjapoolne freesturbaala jagati aga neljaks erineva veetasemega osaks, kus idapoolses osas on maapinna suhtes kõige sügavam veetase ning läänepoolses osas veetase maapinnale lähedane ning keskel võrdlusala. Kõikjal peale võrdlusala eemaldati puurinne ning veetase stabiliseeriti kraavidele rajatud pinnaspaisudega. Läänepoolsel ala rajati veetaseme tõstmise tõttu laienenud veepeegliga kraavidega alal ka raiejääkidest lainerahusti. Taimestiku puhul eeldati looduslikku arengut. Kildemaa jääksoo hõlmas nii mahajäetud freesturbavälja kui rabapoolses osas ka freesturbavälja laiendamiseks eelkuivendusega kuid algse rabataimestiku eemaldamiseta ala, mis on lausaliselt

puurindega kaetud (sarnane Ess-soo vastavale alale). Korrastamismeetmena oli kavandatud kraavide sulgemine pinnaspaisudega, tihedama puistu raadamine ning taimestiku iseseisev areng stabiliseeritud veerežiimi tingimustes. Kõima jääksoos oli korrastamisalal nii turbavõtuaukudega ala (sarnane Maimaja Ess-soo vastava tüübiga) kui ka eelkuivendusega ala (sarnane Kildemaa ja Ess-soo vastavale tüübile kuid oluliselt kõrgema veetasemega ning vähem metsastunud). Kõima uurimisala oli lausaliselt rabaliikidega taimestunud ja vajas korrastamistööde käigus vaid pinnaspaisudega kraavide sulgemist veetaseme taastamiseks ning tervikutel ja kraavide servades puurinde eemaldamist/harvendamist. 2017.a. suvel viidi seiratavates jääksoodes läbi turbalasundi sondeerimine, mille käigus hinnati turbalasundi tüseduse ruumilist varieeruvust ning määrati organoleptiliselt turba tüüp ning lagunemisaste (joonis 1). Jääklasundi omaduste ning taimestumise iseloomu järgi valiti igas jääksoos (v.a. Laiuse) kaks algseisu kõige paremini esindavat piirkonda võrdlusaladeks, kus kogu uuringu raames muudatusi ei tehta ja mille suhtes võrreldakse korrastatavate alade muutusi.

Joonis 1. Üldvaade võrdlusaladele Kõima (vasakul) ja Maima (paremal) jääksoos ning vastavalt kumbagi jääksoo võrdlusalade A ja B turbaprofiilidele. 2017.a. augustis rajati kõigil uuringualadel võrdlusalad (alad mis jäävad muutumatuks ka korrastamistööde käigus ehk referentsalad). Neile aladele installeeriti turbaveevaatluskaevud (Ess-soo 3 tk, Laiuse 2 tk, Maima 2 tk, Kõima 2 tk, Kildemaa 2 tk), veetaseme mõõtekaevud ning gaasivoogude mõõtmise püsivahendid. Samal kuul alustati igakuiste vee- ja gaasiproovide kogumist. Igakuiselt (Laiuse ja Ess-soo puhul ka kaks korda kuus) koguti võrdlusaladelt gaasiproovid (CO2, N2O, CH4), mõõdeti vaatluskaevudes ning kraavides veetase, portatiivsete seadmetega O2 sisaldus (mg/l) ning küllastatustase (O2%), pH, konduktiivsus (µS/cm), ORP (mV) ja koguti veeproovid laboratoorseteks analüüsideks. Laboratoorselt määrati igakuiselt vaatluskaevudest ning võrdlusaladega piirnevatest kraavidest kogutud veeproovidest üldsüsiniku ja üldlämmastiku, lahustunud üldsüsiniku, lahustunud orgaanilise süsiniku, lahustunud anorgaanilise süsiniku ning lahustunud üldlämmastiku sisaldus. 2018.a. suvel alustati taimkatte maapealse seire välitöid kõigis viies jääksoos (Kõima, Maima, Laiuse, Kildema ja Ess-soo), kus võrdlusaladel ja erineva planeeritava korrastamismeetodiga aladel märgistati 1x1 m püsiruudud (kokku 156), võeti nende nurgapostide koordinaadid (RTK, kasutati ka drooniseirel ankurpunktidena). Taimkatteseire püsiruudud fotografeeriti ja neil teostati taimkatte analüüs (üldkatvus, eri rinnete ja taimeliikide esinemine ja katvus). 2017-2020. a. osaleti jääksoode korrastamisprojektide koostamisel ja anti sisend projekteerijatele. 2019.a. suvel-sügisel korrastati Laiuse ja Kõima jääksood. 2020.a. rajati korrastatud Laiuse ja Kõima jääksoodes täiendavad püsiproovialad, installeeriti korrastatud aladele täiendavad piesomeetrid ning vaatluskaevud ja rajati täiendavad gaasivoogude mõõtmise alad (sh. täiendavad ujuvkambrid kraavidele ning veekogule). Neil aladel alustati igakuist seiret ning proovialadele rajatud ülevooludel veeseiret. 2020.a suvel-sügisel korrastati Maima jääksoo, oktoobris installeeriti korrastatud aladele täiendavad piesomeetrid ning vaatluskaevud ja rajati täiendavad gaasivoogude mõõtmise alad. 2021.a. suvel parandati Laiuse jääksoos eraldusvalli kõrge veetasemega ala ja kontrollala vahel ning korrigeeriti kahel alal (reguleerimatute) ülevoolude kõrgust. 2021.a. sügisel-talvel korrastati Ess-soo jääksoo ning alates 2021.a. novembrist alustati

korrastamisjärgset seiret värskelt rajatud ülevooludest. 2022.a. kevadel installeeriti Ess-soos korrastatud aladele täiendavad piesomeetrid, vaatluskaevud ja rajati täiendavad gaasivoogude mõõtmise alad ning alustati korralist igakuist seiret. 2023.a. augustis alustati Kildemaa jääksoo korrastamist. Uuringuperioodi ilmastiku ülevaade Ilmastik mõjutab oluliselt soode veerežiimi aastatevahelist muutlikkust ning seeläbi nii veega toitainete ärakannet, ökosüsteemi gaasivahetust kui ka taimestiku arengut. Looduslikus seisundis sood on suhteliselt suure puhverdamisvõimega, jääksood aga vähese puhverdamisvõimega ning eriti tundlikud on ilmastiku suhtes värskelt korrastatud alad. Kuna aastate lõikes on ilmastik olnud väga erinev, tuleb seda silmas pidada ka aastatevahelisi veetaseme, gaasivoo ning ärakande väärtusi võrreldes ning korrastamistööde üldise edukuse hindamisel. Kui korrastamiseelsel perioodil oli väga põuane vaid 2018 aasta (mis järgnes keskmisest vihmasemale 2017 lõpule), siis korrastamisjärgne periood oli oluliselt kuivem nii 2021, 2022 kui ka 2023 juulini. Uurimisperioodi iseloomustab pikaajalisest keskmisest kõrgem õhutemperatuur (joonis 2). Eriti soojad olid 2019 ja 2020 aastad, mil talvekuudel oli kuu keskmine temperatuur normist isegi kuni 6 kraadi soojem. Soojemad talvekuud tähendasid sagedasi sulaperioode ja kevadel väiksemat lumeveevaru, mistõttu soo veetase sõltus juba varakevadest alates peamiselt sademete hulgast ning päikesekiirguse intensiivsusest.

Joonis 2. Kuu keskmise õhutemperatuuri (joonisel heleroheliste tulpadena) erinevus ( ̊C) uurimisperioodil võrreldes pikaajalise keskmise kliimanormiga (joongraafik). Samblafragmentidega jääksoode korrastamine toimus 2020.a. (Maima) ja 2021.a. (Ess-soo), seetõttu on eriti oluline tähelepanu pöörata perioodi 2021-2023 ilmastikule. Talv algas suhteliselt varakult 2021.a. keskmisest külmema detsembriga kuid maapind külmus varase lumikatte (novembris) tõttu vaid osaliselt ja keskmisest soojem talve jätk (joonis 2) soodustas nii lume sulamist kui külmumata pinnasest gaasivoo eritumist. Korrastamisalade seisukohast oli aga kõige olulisem sulailmadega lumeveevaru kahanemine ja (pool)külmunud pinnasel tekkinud lombid, mis tuule tekitatud lainetusega uhtusid samblafargmentide katteks laotatud põhu vaaludesse. Talvistele keskmiselähedastele sajuhulkadele järgnes aga 2022.a. erakordselt kuiv märts (joonis 3) ning kogu järgneva aasta jooksul oli iga kuu sademete hulk ligi 40% väiksem pikaajalisest kuu sademete normist. Vaid kahel järgneval talvekuul oli sademeid keskmisest enam, kuid 2023.a. kevad algas taas väga tugeva põuaga kui sademeid langes kuude lõikes vaid 30-40 pikaajalisest normist.

Joonis 3. Kuu keskmise sademete summa (joonisel tumesiniste tulpadena) erinevus (%)uuringuperioodil võrreldes pikaajalise kuu keskmise sademete summaga (joongraafik; kuu sademete summa millimeetrites). Keskmisest kõrgem õhutemperatuur, erakordselt väike sademete hulk (132 mm normist vähem 2022.a.) ja keskmisest päikeselisem ilm (eriti 2023 aprillist juunini; joonis 4) tingis intensiivse evapotranspiratsiooni tõttu maist alates kiire veetaseme alanemise (auramine ületas sademete hulka juba märtsist) ja kuivastressi 2021.a. samblafragmentide laotamisega korrastatud uurimisaladel Ess-soo jääksoos, aga ka 2020.a. sarnaselt korrastatud Maima uurimisaladel.

Joonis 4. Kuu keskmine päikesepaiste kestuse summa (joonisel kollaste tulpadena) erinevus (%) uuringuperioodil võrreldes pikaajalise kuu keskmise päikesepaistega tundide summaga (joongraafik; kuu päikesepaistega tundide summa). Kuivale 2022. aastale järgnenud väga kuiv, päikeseline ja suhteliselt tuuline 2023.a. kevad tõi kaasa väga kiire veetaseme alanemise ning samblafragmentide kasvu pidurdumise või isegi kohati hukkumise. Maima jääksoos Ala 9 (kood P) kannatas 2020.a. suhteliselt hästi kasvama läinud turbasammal tugeva talvise külmakohrutuse all. Külmakohrutusest on Maimas iga talv olnud tugevalt mõjutatud ka võrdlusala 1 (Ala 6-2) ja Ala 8 (kood O). Ess-soos on külmakohrutuse mõju olnud väiksem, kõige enam on seal mõjutatud olnud ala 4 (kood F).

TULEMUSED Pinna- ja turbavesi Nii turbavees kui jääksoo kraavides on süsinik ja lämmastik valdavalt lahustunud vormis (vastavalt 92% ja 92% üldsüsinikust ja üldlämmastikust), lahustunud ja lahustumata vormid aga omavahel tugevalt korreleeritud. Lämmastiku ja süsinikusisaldus oli mõõtmisperioodil kõrgem turbavees, kraavides oli pindmise äravoolu ja sademetevee tõttu toimunud vähesel määral lahjendumine. Küll aga on nii kraavi- kui turbavees lahustunud süsiniku (DC) ja lahustunud lämmastiku (DN) suhe väga heaks turba mineraliseerumise indikaatoriks: DC/DN suhe on seda kõrgem ja regressioonvõrrandi seos tugevam, mida enam on ala häiritud (joonis 5).

Maima

y = 11.953Ln(x) + 27.244 R2 = 0.6754

10

20

30

40

50

60

70

0.0 1.0 2.0 3.0 4.0 5.0

Dissolved N (mg/l)

D is

so lv

ed C

(m g/

l)

Kildemaa

y = 19.46Ln(x) + 39.445 R2 = 0.7029

0

10

20

30

40

50

60

70

0.0 1.0 2.0 3.0 4.0 5.0 Dissolved N (mg/l)

D is

so lv

ed C

(m g/

l)

Joonis 5. Lahustunud süsiniku ja lahustunud üldlämmastiku vaheline seos Kildemaa ja Maima jääksoo näitel. Veekvaliteedi ruumilise autokorrelatsiooni hindamiseks koguti veeproovid nii tootmisväljakutevahelistest kraavidest, piirdekraavidest kui väljavooludest 2018 ja 2019 mais ning 2019 septembris. Veeproovide tulemused näitavad tugevat autokorrelatsiooni lahustunud süsiniku kontsentratsiooni (joonis 6) ja mõõdukat korrelatsiooni lahustunud üldlämmastiku (joonis 7) osas. Lahustunud süsinik on kõigi seirealade puhul orgaanilise süsinikuna. Vaid Maima jääksoos esines kevadel ning sügisel pikema vihmaperioodi järel anorgaanilist süsinikku (karbonaadina) Ala 5 (kood B) ja 6-1 (võrdlusala 2) seirekaevudes. Korrastamisjärgselt on 2021-2023 aastatel Ala 5 seirekaevus põhjaveeline toide suurenenud ning anorgaanilise süsiniku (DIC) ja lahustunud lämmastiku (DN) osakaal suurenenud. Eeldatavasti on see seotud pinnaspaisude jaoks liiga sügavalt ekskavaatoriga materjali võtmisel veepidemeks olnud hästilagunenud turbakihi rikkumisest.

Joonis 6. Lahustunud süsiniku (DC) kontsentratsioon (mg/l C) jääksoo kraavi- ja turbavees Ess-soo ning Kildemaa uurimisalade näitel.

Joonis 7. Lahustunud lämmastiku (DN) kontsentratsioon jääksoo kraavivees. Lahustunud lämmastiku sisaldus on väga madal (joonis 7) nii kraavides kui turbavees. Kui DC kontsentratsioon sõltus taimestumisest mõõdukalt ja olulisem oli kraavituse seisukord, siis toitainevaeses keskkonnas on lahustunud lämmastiku kontsentratsioon selgelt madalam taimestunud kraavides ning eriti madal piirkondades kus maapind on lausaliselt kaetud taimestikuga. Korrastamistööde mõju toitainete ja lahustunud süsiniku ärakandele on lühiajaline kui tööde käigus ei mõjutata põhjaveetoitelisust. Maima jääksoos suurenes lahustunud lämmastiku kontsentratsioon (joonis 8) kõige enam alal 5 (Maima B) ja 10 (Maima D). P ja E alal võib osaliselt kõrgemat DN fooni korrastamisjärgselt selgitada ka samblafragmentide laotamise järel laguneva orgaanilise ainega (kattepõhk ning surnud fragmendid).

Joonis 8. Lahustunud lämmastiku kontsentratsiooni muutus uuringuperioodi jooksul Maima jääksoos. Kuigi ka Laiuse jääksoo on toitainerikkam ja osaliselt põhjaveelise toitega, on seal korrastamisjärgsel ajal (pärast 2019 sügist) lahustunud lämmastiku sisaldus kerges langustrendis (joonis 9) ja selgitatav kiiresti areneva taimkattega ja vähese äravooluga. Kui vesi on pikema viibeajaga, siis tarbitakse lämmastikku nii taimestiku poolt kui sõltuvalt redokspotentsiaalist denitrifikatsiooni/nitrifikatsiooni protsessides. Võrdlusalal kus muutused on olnud väikesed, on ka DN püsinud stabiilsena.

Joonis 9. Lahustunud lämmastiku kontsentratsiooni muutus uuringuperioodi jooksul Maima jääksoos. Ess-soo ja Kildemaa uurimisaladel on DN sisaldus vees madal ja aastatevahelised erinevused statistiliselt ebaolulised. Küll aga on selgelt tuvastatav aastaajaline käik kõrgema kontsentratsiooniga suvekuudel, ning kuivematel aastatel (2018, 2022, 2023).

Joonis 10. Lahustunud lämmastiku kontsentratsiooni muutus uuringuperioodi jooksul Ess-soo ja Kildemaa jääksoos. Lahustunud orgaanilise süsiniku sisaldus on Maima jääksoos võrdlusalal 1 jäänud muutumatuks, võrdlusalal Maima 2 langes korrastamistööde käigus põhjaveelise toite lisandumisel väga madalale (joonis 11), aga teistel aladel on sarnane võrdlusalaga 1.

Joonis 10. Lahustunud orgaanilise süsiniku kontsentratsiooni muutus uuringuperioodi jooksul Maima jääksoos. Sarnase selge aastase käiguga kuid olulise trendita on DOC sisalduse käik ka teiste uurimisalade puhul (Laiuse, Ess-soo, Kildemaa, Kõima). Ess-soost võetud puursüdamiku 8-kuune laboratoorne inkubatsioonikatse näitab, et temperatuuri seos DOC-ga pole erinevalt CO2 voost lineaarne, aga 25 C ületav temperatuur suurendab oluliselt DOC teket turbas, olles samas sõltuvuses veetasemest/aereeritusest (Palviainen et al., 2023). Laiuse jääksoos tuleb esile korrastamistööde järgselt madalaveelise veekogu DOC sisalduse langus ja stabiliseerumine 60 mg/li tasandi lähedal.

Joonis 11. Lahustunud orgaanilise süsiniku kontsentratsiooni muutus uuringuperioodi jooksul Laiuse jääksoos.

Vaatamata suhteliselt kõrgemale kontsentratsioonile nii DN kui DOC osas, ei ole ärakanne korrastatud jääksoost suur kuna vee äravool korrastatud aladel on viimastel aastatel (ülevoolu rajamisest saati) olnud vaid lühikesel perioodil talviste sulade ajal ning kevadel lume sulamise järel, mil kontsentratsioonid on keskmisest madalamad. Laiuse jääksoo ülevoolude puhul on äravool vaid märtsis-aprillis, läänepoolses ülevoolus (madalaveelise veekogu ja Lehtmetsa raba vesi) kuni 4 kuud (märtsist juunini). Sarnane on äravoolu periood ka Maima ning Ess-soo puhul. Kõima edelapoolse kraavi äravoolu pole võimalik hinnata kuna vesi valgub ühtlaselt metsa alla. Kirdepoolses äravoolus liigub vesi novembrist maini. Täpse äravoolu koguse hindamine on takistatud kuna Ess-soos viis 30. augusti sadu ülevoolu kõrvalt pinnase ja mitmel sügiskuul puudus äravoolu mõõtmine, Laiuse läänepoolsel ülevoolul muutis kobras V-ülevoolu kuju ja suurust ning Maimal on suure veetaseme kõikumise tõttu olnud vaja vähemalt kaks korda aastas ülevoolu kõrgust reguleerida. Vooluhulga ja kontsentratsiooni järgi hinnates on süsiniku ärakanne DOC kujul jääksoodest vahemikus 62-87 kg/ha*aastas. Mullastik Korrastamistöödega seotud muutused mulla keemilistes omadustes on väga väikesed ja üldjuhul statistiliselt ebaolulised (joonis 12). Ainus oluline muutus on seotud Maima jääksoo mulla happesusega, kus ilmselt on põhjuseks vettpidava turbakihi häirimine ja selle tulemusena suurem põhjavee sissevool alale (eriti Ala 5 (B), aga ka 2 (L), 10 (D) ning 11 (E). Teiste parameetrite osas olulisi muutusi ei toimunud, aga pinnasetööde tõttu suurenes ruumiline varieeruvus. Samblafragmentide laotamisega alal tõusis pindmises kihis süsinikusisaldus keskmiselt ligi 1% võrra, kuid pole selge kas seda tingis täiendav orgaanilise aine lisandumine (sammal, põhk) või eelnevalt osaliselt mineraliseerunud pinnase koorimine.

Joonis 12. Mulla pH, üldfosfori ja üldlämmastiku sisalduse muutus korrastamistööde käigus. Maima jääksoos mulla pH muutuse ja põhjaveelise toitumuse suurenemise vahelist seost kinnitab ka lahustunud anorgaanilise lämmastiku (DN) sisalduse suurenemine poorivees ning kraavides (joonis 8). DN sisaldus on suurenenud samadel aladel (B, D, E) kus tõusis mulla pH sisaldus ning poorivee karbonaatiooni sisaldus, aga muutus ei avaldu võrdlusalal ega selle juures kraavi vees. Jääksoode mullaanalüüsi andmeid kasutati üleriigilise suuremõõtkavalise mulla fosforisisalduse kaardi koostamisel. Valminud kaart on GIS andmestikuna vabavaraliseks kasutamiseks ja metoodika osas detailsemalt kirjeldatud artiklis: Kull, Anne; Kikas, Tambet; Penu, Priit; Kull, Ain (2023). Modeling Topsoil Phosphorus—From Observation-Based Statistical Approach to Land-Use and Soil-Based High- Resolution Mapping. Agronomy, 13 (5), 1183. DOI: 10.3390/agronomy13051183. Biomassi lagunemiskatsed Jääksoodes viidi läbi standardiseeritud teekotikatse rohelise ning punase (rooibos) teega ning Laiuse ja Ess-soo jääksoodes korrastamisjärgselt maa-aluse ja maapealse biomassi lagunemiskatsed. Standardiseeritud teekottide (punane e. rooibos ja rohelise tee) katse esmased tulemused Laiuse jääksoos alustatud eksperimendist lubavad oodata selget seost nii veetasemega kui taimestikuga (joonis 13). Esimese aasta massikadu on Maima ja Kildemaa jääksoos punase tee puhul sarnane Laiuse jääksoos teekottide massikaoga, kuid erinevus rohelise ja punase tee vahel on väiksem. Kildemaa jääksoos on lagunemine mõnevõrra kiirem kui Maima uurimisaladel, eriti rohelise tee osas.

Joonis 13. Vasakul teekottide paigutuse skeem katsealadel, parempoolsel joonisel punase ja rohelise tee jääkmass 3 kuu, 6 kuu, 1 aasta, 1,5 aasta ja 2 aasta pärast Laiuse võrdlusalal (control), rabametsas (Raba), kuivenduse mõjuga rabametsa servas (Kuivendatud mets) ja pinnaspaisudega korrastatud keskmise veetasemega uurimisalal (Keskmine veetase; korrastamisprojektis Ala 2, uurimisala kood Laiuse E) ning alumisel joonisel jääkmass esimese aasta lõpuks Maima ning Kildemaa võrdlusaladel. Rohelise tee lämmastikusisaldus on kõrgem (3-5%) ja imiteerib peenjuurte lagunemist ning on happelises pinnases suhteliselt suure hajuvusega. Punane tee imiteerib rohkem okaste varist ning selle lagunemine on erineva taimestiku ning veerežiimiga aladel ühtlasem. See viitab ka voortevahelises Laiuse jääksoos (turba pH 2.5-3.5, mediaan 3.1) lagundavate mikroorganismide ühtlast aktiivsust erinevates kooslustes ja rohelise teega võrreldes suhteliselt madalamat leostumiskadu, eriti esimese 6 kuu jooksul. Tulemuste põhjalikum analüüs koos kõigi keskkonnategurite (temperatuur, veetase, mullakeemia, sademed jmt) toimub koostöös Iiri ning Rootsi teadlastega ja võrdluses nende sarnaste katsete andmetega. Sarnaselt 2021.a. varakevadel Laiuse jääksoos alustatud maa-aluse ja maapealse biomassi lagunemiskatsetega laiendati 2022.a. katset värskelt korrastatud Ess-soo alale. Lagunemiskatsesse lisati standardiseeritud teekotikatsele ka eraldi proovid männi ja sookase ning villpea, jõhvika ja mustika/sinika peenjuurte ning varisega. Lagunemiskatsed (vahetult maapinnal ning 5-10 cm sügavusel turbas) rajati kuivemal ja märjemal võrdlusalal, turbasambla fragmentide laotamisega alal, pinnaspaisudega tõstetud veetasemega alal ning suletud kuivenduskraavidega rabametsa alal (joonis 14). Katse on korduste arvu järgi planeeritud kolmeaastasena.

Joonis 14. Lagunemiskatse rajamine Ess-soos uurimisalal 2022. aastal. Vasakpoolsel joonisel proovide paigaldamine alale nr. 11 (kood C) pinnaspaisudega suletud kraavidega alal ning parempoolsel joonisel proovide paigaldamine suletud kuivenduskraavidega rabametsa alal. Kaugseire Arvestades jääksoode suurt pindala, raskesti ligipääsetavust, alasisest suurt heterogeensust ning korrastamistööde puhul ka võimalikku kiiret taimkatte arengu dünaamikat, on kaugseire potentsiaalselt hea vahend seisundi hindamiseks. Käesoleva uuringu raames hinnati nii optilise seire (droon ja satelliit) kui radarkaugseire (satelliit) rakendamise võimalusi. Drooniseire peamiseks eeliseks on väga hea lahutusvõime ja võimalus lennata vastavalt vajadusele ning ilmastikuoludele. RGB kaameraga droonid on praeguseks kujunenud laiatarbekaubaks ja pildi kvaliteet on väga hea. Peamised RGB kaameraga droonide kasutamisega seotud metoodilised küsimused puudutavad erinevate aastate lõikes homogeensete aegridade saavutamist, sest vaatamata päikesekiirgusandurite ja kalibreeritud peegeldusplaatide kasutamisele on drooniseireks liiga suurte (eriti Ess-soo ja Maima) alade puhul probleemiks suur kiirgusspektri ajaline varieeruvus. Lennuaja jooksul muutuvad valgusolud ja kiirgusspekter kahandab piltide põhjal automatiseeritud taimkatteklassifitseerimise edukust erinevate ülelendude vahel, aga ka isegi sama päeva lendude osas kui kiirgusintensiivsus jõuab pika lennuaja jooksul oluliselt muutuda. Paremate sensoritega (kiirgusspektri andurid nii üles kui allasuunatuna) droonid, kalibreeritud peegeldusplaadid, georefereeritud ankurpunktid jmt. muudab aga lendamise kalliks ja töömahukaks (sh. kameraalne järeltöötlus). Maima ning Kõima jääksoo korrastamise eelse drooniandmestiku põhjal hinnati erinevate masinõppe algoritmide rakendatavust ja nende maakatte klassifitseerimise täpsust. MarjanSadat Barekaty leidis oma magistritöös Maima jääksoo põhjal, et nii Random Forest (RF), Support Vector Machine (SVM) ja K-Nearest Neighbours (KNN) meetod annavad suhteliselt sarnase tulemuse RGB kaameraga drooniandmestiku puhul. Kõrgeim kaalutud keskmine F1-skoor saadi RF vaikemudeliga kombineerituna vegetatsiooniindeksitega (0,59), sellele järgnesid KNN (0,58) ja SVM (0,57) kombineerituna vegetatsiooniindeksite ja MinMaxScaleriga. Pildi suurem pikslitihedus ei parandanud klassifitseerimise tulemust. Klassifitseerimist raskendas oluliselt UAV ortofoto kõrgest ruumilisest lahutusest tingitud müra ja maakatteklasside mitte tasakaalus olev koosseis (erinevate liikide/koosluste ruumiline esinemine ebavõrdne, mis on aga looduses tavapärane olukord). Teistele uuringutele tuginedes saaks ilmselt klassifitseerimistulemusi parandada kasutades objektipõhist pildianalüüsi (OBIA), mis töötaks paremini puurinde ning mättaid moodustavate taimede puhul ning lisades kalibreeritud multispektraalsed andmed ning lisatunnused (nt. LIDAR andmed). Sarnaselt droonipiltide töötlemise ja sellelt taimkatte tuvastamise metoodikale on võimalik automatiseeritult tuvastada ka taimestumise osakaalu taimkatteruutude fotode alusel. Uuringu käigus arendati QGIS tarkvara baasil fototuvastussüsteemi, et kõrge lahutusega fotodelt (joonis 15 a ja b) RGB kanalites automatiseeritult eristada kasvama läinud turbasambla fragmentide pindalalist katvust. Selleks fotografeeriti standardselt kõrguselt 1x2 m raami jäävad ruudud (100 tk), neist 33 kasutati õpetusalana ja 67 ala automaattuvastuse alana ning neist omakorda 33 lisaks käsitsi klassifitseeritavate kontrollaladena (joonis 15 c).

Joonis 15. Kasvama läinud turbasambla fragmentide tuvastamiseks kasutatud fotod (a ja b), mis georefereeriti ja transformeeriti ortofotodeks. Käsitsi klassifitseeritud alad kasvavate turbasamblafragmentidega (15c) on kujutatud roostepruunide areaalidega. Näited hästi tuvastatavatest fragmentidest (15 d) ja raskesti tuvastatavatest fragmentidest (15 e). Automaatne klassifitseerimine osutus tõhusaks punaka, lillaka, roheka ja rohekaskollaka tooniga turbasammalde puhul (summaarne tuvastamistõhusus 78%; joonis 15d) kuid tõhusus jäi madalamaks kollakaspruuni tooniga sammalde puhul, kus tuvastamist segasid õlgedele ning lagunevatele taimejäänustele sarnased spektraalsed omadused (joonis 15e). Samuti oli raskusi üksikute väga väikeste hajusalt paiknevate või osaliselt õlgedega kaetud väikeste fragmentide tuvastamisega. Sarnaselt drooniandmestiku töötlemisele on ka tavafotode töötlemise puhul eelduseks pildistamine sarnastes valgusoludes, suur õpetusandmestik ja suhteliselt väike eristatavate klasside arv. Suurem klasside arv või Random Forest/Bagging algoritmide kasutamine tekitab rohkem segaklasse, mille sisu on raskesti tõlgendatav. Lisaks RGB kaamerale katsetati Laiuse testalal ka infrapunakaameraga (IR) drooniseiret, et ühest küljest parandada RGB kaameraga kombineeritult taimkatteklasside eristamise võimet ja teiseks hinnata taimestumise edukust maapinna temperatuuri alusel (suvine kõrge pinnatemperatuur on hüpoteesikohaselt taimestumisele oluline takistus) ning maapinna erineva soojenemise kaudu (kaks ülelendu IR kaameraga hommikul jahtunud maapinnaga ning pärastlõunal maksimaalselt soojenenud maapinnaga ajal) välja töötada maapinna niiskuse arvutamise metoodika. Paralleelselt IR droonilennule viidi läbi ka maapinnal kontaktmeetodil pinnatemperatuuri ja mullaniiskuse (m3/m3) mõõtmine (joonis 16). IR kaameraga testiti ka erineva lennukõrguse mõju 5 m kõrguse muuduga vahemikus 35-150 meetrit, sobivaimaks lennukõrguseks on taimkattestruktuuri määramiseks 70-80 m, maapinna temperatuuri kontrasti järgi niiskuse hindamiseks piisab ka 150 m lennukõrgusest.

a) b)

c) d) e)

Joonis 16. IR kaameraga mõõdetud maapinna temperatuur (23.aug.2018, kl. 15) ja samal ajal maapinnal kontaktmeetodil mõõdetud mullaniiskus (iga lilla ja kollane punkt tähistab mõõtepunkti). Jahutuseta laiatarbe infrapunakaamera droonidele osutus kogu uurimisala katva komposiitpildi koostamiseks ebatäpseks (vt. joonis 15 vasakpoolse kujutise lennusuunast sõltuvat triibulisust) ja mõjutab seeläbi lõpptulemust. Samas temperatuurikontrasti väärtused (pärastlõunasest temperatuuri komposiitpildist lahutatud hommikune temperatuuri komposiitpilt) korreleerusid mõõdetud mullaniiskuse väärtustega. Termokaameraga droon sobib suurepäraselt ka allikaliste kohtade või hilissügisest kevadeni lekkivate pinnaspaisude tuvastamiseks. Arvestades seda, et droonipildi alusel on väga keeruline (ja/või kulukas) koostada aastateülest homogeenset aegrida, on taimestikuseire puhul kõige tõhusam drooniseire kasutamine üldise taimestumise hindamiseks dominantliikide/koosluste alusel ning nende piiride pikemaajalise muutumise jälgimiseks. Lausalise kaardistamise aluseks võiks olla k-means meetodil loodud aluskaart (selle loomine ei eelda eelnevat ala seiret), mille klassidele antakse sisu georefereeritud väliuuringute abil. Eristatavate klasside arv sõltub kasutatud lähteandmestikust, jäädes Sentinel satelliidi optiliste kanalite ja indeksite kasutamisel enamasti 5-7 klassi vahemikku, drooniseire RGB andmete puhul 7-9 klassi ning multispektraalsete kanalite kasutamisel 10-12 klassi piiresse. K-means meetodil loodud dominantklasside arvu määramine on empiiriline, eeskätt ekspertteadmistel põhinev ning vajab reeglina 3-4 erineva versiooni loomist, mille puhul statistiliselt eristunud klassid sisustatakse georefereeritud välitööandmestiku alusel uurimisalal. Neist versioonidest valitakse lõpuks välitingimustes reaalselt tuvastatavate ja looduses eristuvate klasside alusel sobivaima klasside arvuga aluskaart. Seega on drooniseire kõige paremini kasutatav a) ala (visuaalse) eelhinnangu ja seirealade esindusliku paigutuse koostamiseks, b) väiksemate alade detailseks sagedaseks võrdlemiseks (nt. veepiiri või mingi taimestikuareaali aastaajaline dünaamika), c) termokaameraga vee liikumise ja allikaliste kohtade ning paisude lekete tuvastamine, d) sisend ajas dünaamilise ruumilise kasvuhoonegaaside mudeli jaoks taimestiku katvuse muutuse alusel (eeldab vähemalt 3-4 perioodi katmist igal aastal: varakevadine lumesulamine, kevaduvine tärkamine, suvine rohtse biomassi maksimum, sügisene samblarinde seisundi hindamise aeg). Küll aga eeldab selline detailsusaste drooniseire puhul suurt arvutusjõudlust, ajakulu ning arvestatavat rahalist ressurssi. Optiline satelliitseire tagab samaaegselt suure ala katvuse, kuid on väikse ruumilise lahutusega (piksel u. 5-30 m vahemikus) ja ei saa valida ilmastiku järgi sobivat ülelennu aega. Arvestades ülelendude sagedust ja meie laiuskraadil tavapärast pilvisust, on Sentinel-2 missiooni näitel kuu kohta keskmiselt kasutada 1-2 päeva kvaliteetset kujutist (valdavalt pilvevaba) huvipakkuvast alast. Sügisel ja talvel võib esineda kuid, mil kvaliteetset kujutist ei saadagi. Korrastamata jääksoode puhul on see piisav kuna muutused on üldjuhul väikesed (erandiks kevadeti üleujutatavad alad), aga korrastamisjärgseks seireks on see aastaajalise arengu dünaamika hindamiseks ebapiisav. Küll aga sobib selline sagedus pikaajaliseks (paljude aastate üleseks, st. enam kui 10-aastase perioodi muutuste) kindla fenofaasi või aastaaja alusel (madalsoo ja rohundirikka ala puhul kesksuvine, turbasammaldega aladel sügisene periood) hindamiseks. Ülelendude sagedus aga omakorda on seotud kaetava ala suurusega (piksli suurusega) – nii näiteks saab MODIS missiooni Terra ja AQUA satelliitide abil arvutada maapinna ööpäevase temperatuuri amplituudi, aga piksli suurus ulatub kilomeetrini ja huvipakkuva ala sisu kipub hägustuma kuna hõlmab nii freesturbaväljakuid, kraave kui servas ümbritsevat ala (joonis 17).

Jääksoo korrastamine 08-10.20219

Joonis 17. Päevane maapinna temperatuur (°C) Laiuse korrastamisalal (Laiuse 1) ja looduslikus seisundis rabametsas (Laiuse_natural) Terra satelliidi andmestiku alusel aastatel 2017-2021.

Samas on sel viisil aastane pidev temperatuurikäik uuritavalt alalt tagatud ja seda saab kasutada näiteks sisendina mullahingamise (Rsoil) või ökosüsteemi hingamise (Reco) modelleerimiseks nagu näidatud jääksoode näitel Burdun et al., 2021 poolt. Vaatamata madalale ruumilisele lahutusele on selline maapinna temperatuur sisendandmestikuna parem kui lähimas ilmajaamas mõõdetud õhutemperatuuri või maapinna temperatuuri vahetu kasutamine, kuna ilmajaam asub mineraalpinnasel, kus termiline režiim on soomuldadest erinev. Ökosüsteemihingamise modelleerimiseks nii vahetult mõõdetud kui kaugseire andmete alusel on sobilik järgmine valem (Riutta et al., 2007; Järveoja et al., 2016):

Metaanivoo hindamine satelliidi andmetel põhineva maapinna temperatuuri alusel ei anna häid tulemusi kuna metanogenees on seotud sügavama anaeroobse turbakihiga ning aereeritud tsooni temperatuur pigem soosib metaani oksüdeerimist/metanotroofide poolt tarbimist ja kahandab metaanivoogu ning selle seost sügavama kihi termiliste omadustega. Kui looduslikus soos metaanivoog ligikaudu järgib aastast temperatuurikäiku (mõningase ajalise nihkega), siis jääksoodes on seos nõrk ja olulisem on sademete hulk ning poorides vee küllastatu hapnikuga. Neid näitajaid paraku praeguste teadmiste kohaselt pinnakihist sügavamal kaugseire vahenditega piisava ruumilise lahutuse ning ajasammuga ei ole võimalik tuletada. Teataval määral võimaldab seda satelliitradarandmestik (SAR), kuid ka seal on avalikult kasutatava andmestiku lainepikkus sobiv vaid väga õhukese pinnakihi kirjeldamiseks.

Joonis 18. Mõõdetud ja satelliidi maapinnatemperatuuri andmete alusel modelleeritud Reco looduslikes soodes (Männikjärve, Linnusaare), kuivendusega jääksoode osas (Kõima 1, Kildemaa 2) ja jääksoo freesturbaväljadel (allikas: Burdun et al., 2021). Optilise kaugseire abil jääksoode korrastamistööde järgse arengu kirjeldamiseks on tavapärase nähtava valguse spektriosa (RGB) kõrval otstarbekas kasutada erinevate spektriosade alusel koostatud indekseid. Kuna jääksoode seisund, korrastamismeetodid (veekogu, metsastamine, rohttaimedega madalsoo-suunaline korrastamine, samblafragmentide laotamine, pinnaspaisude kasutamine isetaimestumisega jne.) on alade lõikes varieeruvad, on vajalik erinevaid indekseid kasutada. Madalaveeliste taimestuvate veekogude puhul annab parima tulemuse NRG indeks, taimestumata veekogu piiritlemiseks aga NDPI. Avavett ja väga niisket pinnast kajastavad paremini NRG ja NGR indeksid, kuid NGR puuduseks on see, et ei suuda edasi anda infot kuivema taimestumata turbaga piirkondade kohta (mis jääksoo korrastamise seisukohast on oluline määratleda). Rohundirikka taimestikuga jääksoo, metsastunud/metsastatud jääksoo kirjeldamiseks sobib hästi laialt kasutatav taimkatteindeks NDVI. Joonis 19 illustreerib 2020.a. korrastatud Maima jääksoo erineva taimestumismääraga (ja korrastamisviisiga) alade ning seda ümbritseva looduslähedase rabataimestikuga ala näitel erinevate indeksite võimekust seisundit kirjeldada.

Joonis 19. Sentinel-2 satelliidi andmete alusel arvutatud indeksid Maima korrastatud jääksoo näitel (21.09.2023). RGB (Red/Green/Blue) iseloomustab tavapilti nähtavas spektriosas, NRG (nIR/R/G) indeksit kus sinine spektriosa on asendatud lähisinfrapunaga, NDVI (normalized difference vegetation index) taimkatet kajastav indeks, NGR (nIR/G/R) sarnane NRG indeksiga niiskuse kirjeldamiseks, NDPI (Normalized Difference Pond Index; (mIR1- Green)/(mIR1+Green) ja NNR (nIR/nIR/Red).

Satelliidiseire andmete alusel kiire hinnangu andmiseks korrastamise edukuse kohta lühiajalise perioodi alusel (mõned aastad) on takistuseks väheste pilvevabade kaadrite esinemine. Atmosfääri läbipaistvus (eriti pilvisus, veeaur) mõjutab oluliselt kõigi optilise seire kanalite alusel arvutatud indeksite väärtust ja võib mõjutada arvutatud ajalisi trende. On üsna sage, et kogu kuu lõikes pole ühtegi hea lahutusega (piksel 10m või väiksem) pilti kogu uurimisala kohta ning erineva pilvisusega tehtud piltide alusel komposiitpilt lahendab probleemi vaid osaliselt. Joonis 20 iseloomustab Maima jääksoo näitel 2022 sügisest (kuiva pika põuase suvega aasta) ja 2023 sügiseni (kuiva kevadsuvega aasta) näitel ühe aasta jooksul RGB, NDVI ja NRG indeksite aastaajalist dünaamikat. Tähelepanu tuleks pöörata Ala 1 (kood M), 5 (B), 7 (N) kiirele taimestumisele valdavalt pilliroo, villpea ja tarnadega ning samblafragmentide laotamisega kuid kõrge veetaseme all kannatavate alade 3 (K) ja 4 (C) kokkuveoteeäärse tsooni muutustele ning normaaltingimustes sobiliku ala 9 (P) seisundi muutusele.

Joonis 20. Korrastatud Maima jääksoo seisundi muutus iga kuu parima kvaliteediga (pilvevabama) pildi alusel RGB (vasakpoolne veerg), NDVI (keskmine veerg) ja NRG (parempoolne veerg) näitel 2022 sügisest alates kuni 2023.a. sügiseni. NDVI mustja ja punakad toonid iseloomustavad rohelise taimestikuta (ja/või veega ning tehispinnasega alasid, tumeroheline lausaliselt taimestunud alasid).

Joonis 20. järg

Joonis 20. järg Joonis 20 illustreerib hästi kuidas 2023.a. väga kuiva suve järel kahanes juulini veega kaetud ala, aga septembris oli taastunud liiga kõrge veetase peaaegu kevadise seisuni (eriti ilmekas Ala 6-1, 4, 1, 7 ja 11 näitel, eriti NRG indeksiga väljendatuna). Seejuures pilvevabade piltide puudumise tõttu ei tule kaugseire andmetest välja, et muutus toimus lühikese aja jooksul just vahetult pärast augustikuise pildi tegemist ning järgmiste ülelendude ajal oli taevas lausalise pilvkattega.

Satelliidi radarandmestiku (SAR) puhul on eeliseks selle vähene sõltuvus ilmastikust või pilvisusest, kuid ülelendude sagedus on väike ja aluspinna koherentsuse muudu alusel pindalaline lahutusvõime tagasihoidlik (enamasti vajalik hektarile lähenev pindala, et sisulisi muutusi ajas eristada). Kõrgusmuudu kaudu niiskusrežiimi muutumise hindamine DInSAR (järjestikuste kujutiste faaside alusel arvutamise meetod) on soos võimalik (nn. soo hingamise mõõtmine) ja enamasti üsna täpne (mõõdetav millimeetrites), kuid probleemiks on ülelendude sagedus, sest erandlikel juhtudel võib kahe pildi vahelisel perioodil maapinna kõrguse muut sadude tõttu ületada faasi ulatust (Sentinel 1 C-band puhul u. 2.5 cm) ja sel juhul tegelik kõrgusmuut jääb teadmata arvu faaside võrra ekslikuks (Tampuu et al., 2023). SAR andmestikku on võimalik kasutada muutuste tuvastamiseks ka koherentsuse kaudu. Sel juhul on soodes vertikaalne-vertikaalne polarisatsioon muutuste kirjeldamiseks tõhusam kui vertikaalne- horisontaalne polarisatsiooni kasutamine, kuna viimasel on just jääksoodes suurem hajuvus (joonis 21).

Joonis 21. Sünteetilise apertuurradari (SAR) erinevate polarisatsioonide hajuvus (nii tõusva kui laskuva suhtelise orbiidi RON alusel 6-päevase sammuga andmestiku põhjal lumevabal perioodil) loodusliku lageraba, jääksoo ning kasutuses oleva freesturbavälja võrdluses (Tampuu et al., 2020). Maima jääksoo uurimisperioodi hõlmav koherentsuse muutuses endisel freesturbaväljal ja turbavõtuaukudega alal võrreldes loodusliku taustaalaga tuleb väga selgelt esile järsk muutus freesturbaväljal alates 2020 a. lõpust (joonis 22), mil veetase järsult freesturbaväljakutel tõusis ning seejärel kajastuvad 2022.a. kuiv suvine-sügisene ning 2023 kuiv suvine periood kasvava koherentsusena (veega kaetud ala kahaneb). Looduslik ning turbavõtuaukudega ala reageerivad 2022 põuale aga vastandsuunalisena (kuiva turbasambla niiskus ja vastavalt elektrijuhtivus kahaneb).

Joonis 22. SAR kahe suhtelise orbiidi (RON 58 ja 80) alusel vertikaalne-vertikaalne polarisatsiooniga jääksoo korrastamisega seotud muudatuste tuvastamine Maima jääksoos endisel freesturbaväljal (Maima_frees), turbaaukude piirkonnas (serv) ning raba looduslikul taustaalal (looduslik).

Ka Ess-soo uurimisalal on täheldatavad sarnased muutused SAR andmestiku alusel (joonised 23, 24).

Joonis 23. Ess-soo SAR pilt suhteliselt orbiidilt RON 160 kevadel kõrgema veetasemega perioodil 1. märtsil 2022. Sinakad toonid iseloomustavad madalat koherentsust (puurinne, vaba veepind) ning kollakad ja punakad toonid suuremat koherentsust. Mustad piirjooned tähistavad Ess-soo erinevaid korrastamisalasid, millest on välja jäetud eraldavad pinnaspaisud, kokkuveotee ning kraavid ja üleminekulised tsoonid.

Joonis 22. SAR andmete alusel vertikaalne-vertikaalne polarisatsiooniga jääksoo korrastamisega seotud muudatuste tuvastamine Maima ja Ess-soo jääksoos. Ülemine joonis iseloomustab 2020.a lõpus järsu veetaseme tõusu tõttu suurt muutust Maima jääksoos, kuid 2021.a. sügisel Ess-soos sarnast õleujutust ei esinenud ning muutus koherentsuses on tagasihoidlikum. Alumine joonis iseloomustab korrastatud alase väga sarnast sünkroonsust maapinna niiskuse muutuses põua tõttu 2022 ja 2023.a., kuid toob ka välja erineva suhtelise orbiidi (RON) valiku olulisuse niiskuse kirjeldamise seisukohast.

Veetaseme dünaamika Veetaset, kasvuhoonegaaside voogu ning Maimas ja eriti Ess-soos värskelt korrastatud aladel samblafragmentide kasvama minekut (ka laiemalt alade taimestumist) mõjutas väga tugevalt 2022.a. ja 2023.a. ilmastik. Kui Maima jääksoos tõusis pärast korrastamist 2020.a. lõpus ja 2021.a. veetase sammalde kasvuks ebasoodsalt kõrgeks, siis Ess-soo korrastamisele järgnes kaks väga kuiva suve. Talv algas suhteliselt varakult 2021.a. keskmisest külmema detsembriga kuid maapind külmus varase lumikatte (novembris) tõttu vaid osaliselt ja keskmisest soojem talve jätk (joonis 3) soodustas nii lume sulamist kui külmumata pinnasest gaasivoo eritumist. Korrastamisalade seisukohast oli aga kõige olulisem sulailmadega lumeveevaru kahanemine ja (pool)külmunud pinnasel tekkinud lombid, mis tuule tekitatud lainetusega uhtusid samblafargmentide katteks laotatud põhu ning osaliselt ka samblafragmendid vaaludesse. Talvistele keskmiselähedastele sajuhulkadele järgnes aga 2022.a. erakordselt kuiv märts (joonis 3) ning kogu järgneva aasta jooksul oli iga kuu sademete hulk ligi 40% väiksem pikaajalisest kuu sademete normist. 2023.a. kevadsuvi osutus aga veelgi kuivemaks ja veetase alanes taas väga kiiresti, langedes Maima jääksoo võrdlusalal ning samblafragmentide laotamisega pinnaspaisudega alal 9 (kood P) ligi 60 cm sügavusele maapinna suhtes. Juulis alanud sademed küll tõstsid veetaset, aga optimaalse tasemeni (-20 cm) jõudis see alles septembris (joonis 23).

Joonis 23. Kuu keskmise veetaseme dünaamika Maima jääksoo korrastatud aladel ning võrdlusalal. Alade tähises „sph“ näitab turbasamblafragmentide laotamist, „Pais“ ala korrastamist ainult pinnaspaisude rajamisega, „Täis“ lausaliselt pinnasega täidetud kraave. Halli varjutusega ala indikeerib eelistatud veetaseme vahemikku korrastatud alal (veetase maapinna suhtes vahemikus 0...-20 cm). Laiuse jääksool on Lehtmetsa raba näol suur tagamaa madalaveelisel veekogul ning mõningane põhjavee toide, mis koostoimes Lehtmetsa peakraavil toimetavate kobrastega tagasid suhteliselt hea veetaseme stabiilsuse kogu korrastatud ala ulatuses (v.a. kõige kõrgema maapinnaga väike eraldatud idapoolne nurk) ja veetase oli kogu aasta ulatuses vahemikus 0...-40 cm (joonis 24). Sellest tulenevalt algas 2022 aastal ja jätkus 2023.a. jõudsalt ülepinnaline taimestumine Laiuse kesksel korrastamisalal (kood Laiuse E) ning läänepoolsel alal (Laiuse W), kuid jäi puudulikuks kõige kuivemal väikesel idapoolsel alal. Samuti laienes keskmiselt 4.4 meetri võrra veekogu suunas taimestunud vöönd madalaveelise veekogu põhja-, edela- ja lõunaservas, mis on madalamad ja laugema kaldaga. Kõima jääksoos on küll veetase tänu suurele looduslikule puhverdavale tagamaale ning juba algselt lausalisele samblakattele optimaalse lähedal, aga nii 2021. kui 2022.a. on veetase ilmastikust tingituna augustiks langenud madalamale kui eelnevatel aastatel. Seevastu Kõima turbavõtuaukude veetase on oluliselt tõusnud (eriti gradiendiga korrastamisala edelaosa suunas) ja turbavõtuaukude vahelised tervikud on muutunud niiskemaks, veetase kõrgem (Kõima S tervik; joonis 24) kui võrdlusalal ja kvalitatiivselt on märgatav ala lääne- ning edelaosas tervikute servades turbasambla laienemist aukudest tervikule, kanarbiku ja samblike hääbumist ning nokkheina ja villpea lisandumist.

Joonis 24. Kuu keskmise veetaseme dünaamika Kõima ja Lause jääksoo korrastatud aladel ning võrdlusalal. Kasvuhoonegaaside voog Korrastamise käigus saavutatud kõrge veetase on kahandanud turba lagunemise kiirust ja süsihappegaasi lendumist korrastatud aladelt nii Kõima, Laiuse, Ess-soo kui Maima jääksoos. Peamine mõju on Maima ja Ess-soo alal saavutatud turba lagunemise aeglustumise kaudu, Laiuse jääksoos aga ka kiiresti arenema hakanud taimestiku tõttu (peamiselt karusammal, jõhvikas, pilliroog, lääneoas ka turbasammal). Juba algselt lausalise taimkattega Kõima jääksoos gaasivoo osas statistiliselt olulisi muutusi ei ole, pigem on muutused selgitatavad aastate vahelisest ilmastiku erinevusest. Kõima jääksoo puhul on turbavõtuaukudes kõrgema veetaseme tõttu edelapoolses osas lopsakalt arenemas älvestele iseloomulikud turbasambla liigid ning kohati laieneb turbasammal ka madalamatele terviku osadele. Enamasti on tervikud siiski aeglase taimestumisega ja gaasivoogu mõjutab enam turba niiskusrežiimi muutus. Maima jääksoo kontrollala nr. 2 on aastaringselt lausaliselt 30-50 cm paksuse veekihiga kaetud ja jäetud antud analüüsist välja kuna ei vasta enam kontrollala kriteeriumitele. Kontrollala nr. 1 on samuti korrastamistööde järel märjemaks muutunud (eriti kevadel ja sügisel), mistõttu ka põuasel 2022 ja 2023.a. suvel oli seal veetase sarnane uuringuperioodi algusega, aga 2023.a. ei avaldunud see mõju veel taimestiku arengus väljaspool kraavi servasid ning ala 8 (O) piirdevalliga külgnevat peenart, kus on intensiivne jõhvika areaali laienemine. Kogu endise freesturbavälja ulatuses on domineeriv mullahingamine, autotroofne hingamine ja taimede fotosüntees on aastase voo mõttes enamasti tagasihoidlik. Erandi moodustavad pillirooga kattuvad alad (Ala 1 (M), 5 (B), 7 (N) ja turbasamblaga endised turbavõtuaugud (ala 12 (G)), kus keskpäevane ökosüsteemi CO2 sidumine (NEE, Net Ecosystem Exchange) võib ulatuda pilliroo puhul -192 mg CO2-C m2 h-1 ja turbasamblal -77 mg CO2-C m2 h-1. Enamasti jääb siiski aeglase taimestumise, laotatud põhu ja surnud samblafragmentide tõttu NEE isegi suvekuudel Maimas emissiooni poolele. Kui 2021.a. oli samblafragmentidega korrastatud aladel süsihappegaasi emissioon ligi poole väiksem kui kontrollalal ning lausalise kraavide täitmisega alal omakorda väiksem kui pinnaspaisudega suletud kraavidega alal, siis 2022.a. sellist erinevust ei esinenud ja vaid suuremalt jaolt veega üleujutatuks jäänud alad (C ja K) olid teistest väiksema emissiooniga, kuid 2023.a. oli ka neil aladel voog ülejäänuga sarnasemaks muutunud. Sellest tulenevalt on ökosüsteemi hingamine (Reco) jätkuvalt hea indikaator süsihappegaasi emissiooni väljendamiseks (joonis 25), mis toob kombineeritult välja nii mullahingamise kui taimestiku arengu mõju. Ökosüsteemi hingamine jäi vaatamata kahele järjestikusele soojale kuivale suvele valdavalt samale tasemele kui eelnevatel aastatel. 2021.a. veega kaetud aladel aga 2022 ja 2023.a. põuastel suvedel vesi soojenes kiiresti ja veetase alanes, jättes maapinna kohati mudaga kaetuks ja suurendades süsihappegaasi voogu. Erandlik on joonisel ala B (paisudega suletud kraavid, veega osaliselt üleujutatud), kus 2022.a. suvine Reco CO2-C piik on seotud intensiivse pilliroo kasvuga ning taime hingamine kombineerub sooja mudaja pinnase emissiooniga. Lisaks mõjutas üleujutatud alade voogu ka surnud kanarbiku jmt. lagunemine. 2023.a. taimestumise tõttu päevane NEE suurenes sel alal ligi 20 mg CO2-C m2 h-1 võrra ja surnud taimede lagunemine on aeglustunud.

Joonis 25. Ökosüsteemi hingamine (Reco) Maima jääksoos. Märge „veega“ iseloomustab korrastamise järgselt üleujutatud ala, „norm“ tähistab normaalse veerežiimiga ala, kus veetase jäi valdavalt maapinnast sügavamale. Aasta CO2 bilanss oli 2022.a. sambla fragmentidega korrastatud üleujutatud aladel emiteeriv (0.46 t/ha C), kraavidel pinnaspaisudega korrastatud aladel 0.86 t/ha C ning koos fragmentide laotamisega aladel 0.75 t/ha C. Turbaaukudes aga toimus tänu ohtrale päikesekiirgusele ning optimaalse lähedasele veetasemele (kohev sammal liigub sünkroonselt veetaseme muutusega, veetase 0...-15 cm) sidumine NEE -1.04 t/ga C. 2023.a. kuni augustini samblafragmentidega korrastatud aladel päevane NEE suurenes sidumise suunas ligi 10 mg CO2-C m2 h-1 võrra, aga enamikul suvekuudel jäi endiselt emiteerivaks isegi päevasel fotosünteesi toimumise ajal väga hõreda taimestiku tõttu. Kõima jääksoos on võrdlusala emiteeriv (1.9 t/ha C), turbaaukude vaheline tervik emiteerib 2.6 t/ha C, samas kui turbaaugu emissioon on 1.1 t/ha C ning kuiva suve tõttu oli ka looduslähedases seisus rabaosa emiteeriv (0.6 t/ha C) ning 2023.a. suvekuudel emissioon kasvas 2022.a. võrreldes veelgi. Laiuse jääksoos algas 2022.a. suve teises pooles ja jätkus kogu 2023.a. kiire taimkatte levik varasemalt palja turbaga alal. Kasvuala laiendasid kõige jõudsamalt karusammal ja jõhvikas, kraavides pilliroog, tarnad ning valge vesiroos. Läänepoolses osas kus turbaaukudele rajati lainetõkked, laienes kiiresti pilliroo ning hundinuiaga kaetud ala, tervikutel ja madalamates niiskemates lohkudes turbasammal. Kiire taimkatte muutuse tõttu allus mõõtmisandmestik modelleerimisele gaasimõõtmisrõngaste lõikes erinevalt (R2 0.43-0.95). Laiuse 1 (võrdlusala) on läbi kõigi aastate olnud CO2 emiteerija, Laiuse idapoolne (Laiuse E) korrastamisala oli 2021.a. emiteeriv, kuid 2022.a. saavutas sidumise jõhvikaga kaetud alal ning pillirooga taimestunud alal. Kõige märjemal alal (Laiuse W) on 2022.a. süsinikuneutraalsed või siduvad kõik taimestunud alad (joonis 26). Kraavide ning Laiuse madalaveelise veekogu süsinikubilanss on positiivne, keskmine emissioon 0.42 t/ha C. Seisva veega kraavides võib küll suve alguses vetika vohamise tõttu mõnel kuul süsihappegaasi sidumine olla intensiivne, aga suve teises pooles algab tekkinud biomassi lagunemine ja eritub nii süsihappegaasi kui metaani. 2023.a. suurenes päevane NEE sidumine kõigis korrastatud jääksoo osades 10-30 mg CO2-C m2 h-1 võrra võrreldes 2022.a., kõige enam pillirooga taimestunud pinnaspaisuga suletud kraavil.

Joonis 26. CO2 bilanss korrastatud Laiuse jääksoos. Laiuse 1 on võrdlusala, Laiuse E keskne korrastamisala võrdlusalast idas ning Laiuse W kontrollalast läänes paiknev maapinnalähedase veetasemega korrastamisala. Kõrge maapinna temperatuur ning suhteliselt kõrge veetase soodustavad metaani teket. 2022.a. olid pika põua tingimustes Maima jääksoos metaani tekkeks äärmiselt soodsad tingimused. Kuigi veega kaetud korrastatud aladel oli suvel keskmiselt kõrgem metaani emissioon, oli ka nii pinnaspaisude kui täidetud kraavidega korrastatud alasid, kus metaani voog oli suur. Samas pinnaspaisudega ala 10 (D) ja võrdlusala olid endiselt väga madala metaani emissiooniga, aga eelneval 2021.a. suvel oli just ala 10 kõrge vooga kui seal veetase lühiajaliselt väga kiiresti muutus. 2023.a. jäi eelnevate aastatega võrreldes metaanivoog oluliselt väiksemaks oli korrastamismeetodist sõltumatult sarnane.

Joonis 27. Kuu keskmine süsiniku kadu metaanina lendumise kaudu Maima kontrollalal (2017-2023) ja korrastamisjärgselt nii kontrollalal kui korrastatud aladel.

Naerugaasi voog oli 2022.a. sarnaselt eelnevatele aastatele toitainevaestes tingimustes kõigis uuritavates jääksoodes ebaoluliselt väike (joonis 28). Suhteliselt pika kuiva perioodi ja hoovihmadest tingitud veetaseme kiirete kõikumiste tulemusel suurenes N2O voog korrastamise järgselt Maimal juba 2021.a. ning veelgi selgemalt 2022.a., aga ka need vood on väga väikesed. Ainsaks erandiks oli september kui pärast pikka põuaperioodi ja sügavale langenud veetaseme juures algasid intensiivsed sajuhood, mis kiirelt täitsid pinnaspoore ning soodustasid lühiajalist naerugaasi heidet. Sarnane põuajärgne järsk naerugaasi voo lendumine septembris leidis aset ka teistel uurimisaladel. 2023.a. kiire kevadine veetaseme alanemine ja püsimine stabiilsena kuni augustis ohtrate sademete tõttu veetase taas kiirelt tõusis ja stabiliseerus, jäi naerugaasi voog väga madalaks. Teistest aladest eristub juba teist põuast aastat järjest Maima võrdlusala 1, kus korrastamise järgselt on veetase muutunud kõikuvamaks (kevadel ja sügisel naaberalade tõttu veetase kraavides tõuseb) ning see on kaasa toonud võrdlusalal naerugaasi suurema emissiooni, mis absoluutväärtuselt on siiski ebaoluline.

Joonis 28. Naerugaasi emissioon Maima jääksoost perioodil 2017-2023. Märgalade gaasivood on ajaliselt ja ruumiliselt suure varieeruvusega, seetõttu on ennatlik paari korrastamisjärgse aasta ning ühe või kahe ala tulemuste põhjal teha järeldusi korrastamismeetmete tõhususe osas. Maima jääksoos on samblafragmentide abil taimestumise kiirendamine valdavalt ebaõnnestunud liiga kõrge ning kõikuva veetaseme tõttu, aga samas on kõikidel samblafragmentide laotamisega aladel vähemalt mingil määral hajusalt kasvama läinud samblaid ning lisandunud on teisi raba liike. Aladel kuhu samblafragmente ei laotatud ei ole ka sõltumata veetasemest või paiknemisest looduslikuma taimestikuga rabaosa suhtes turbasamblaid iseseisvalt alale ilmunud. Kaks põuast suve on Maimal taimestumist oluliselt kiirendanud, eriti aladel 1 (M), 2 (L), 5 (B) ning 7 (N). Turbavõtuaukude juures on korrastamise mõju vähemärgatav, ilmselt põuaste suvede tõttu, sest kevadel on turbaaukudes veetase tervikute tasapinnani, kuid suveks taandub oluliselt. Sama tähelepanek kehtib ka Ess-soo kohta, kus taimestiku taastumise aeg on olnud oluliselt lühem, aga Ess-soos on just turbavõtuaukude ja metsas eelkuivenduskraavide sulgemisega alal veetase püsinud hästi ka suvedel ning turbasammalde laienemine olnud kiire (sh. ekskavaatori tekitatud aukudes ja pinnaspaisude külgedel niisemates lohkudes). Seevastu Laiuse jääksoos on kõigil korrastatud väljakutel ilmunud vähemalt mõnes piirkonnas ka iseseisvalt turbasamblaid, kohati on turbasammalde areaali laienemine alates 2022.a. suve lõpust muutunud kiireks. Detailne ülevaade taimkatte muutustest seiratavate jääksoode püsiseireruutudes on esitatud aruande II osas „RMK taimestiku seire KOONDARUANNE.pdf“.

Tähelepanekuid ja soovitusi korrastamisalade põhjal Alade jagamine väiksemateks hüdroloogilisteks üksusteks on ennast õigustanud, vähendades nii veelgi ulatuslikumaid üleujutusi või ulatuslikumaid liiga kuivi alasid. Iga eraldusvalli sees tuleb väljaku madalama osa juures tekitada ülevool, ülevoolude puhul tuleks kasutada reguleeritava kõrgusega ülevoolu lahendusi (joonis 29). Need on lihtsad, kuid võimaldavad vähemalt esimestel taimestumise seisukohast kriitilisel aastatel ilmastiku, projekteerimis- või ehitusvigade tõttu tekkinud veetaseme probleeme leevendada.

Joonis 29. Pinnasvall jääksoode eraldajana (vasakul) ning lihtne kuid tõhus reguleeritav ülevool veetaseme reguleerimiseks. Selliseid ülevoole tuleks kasutada iga hüdroloogiliselt eraldatud jääksoo eraldusvalli juures. Pinnaspaisudega suletavate kraavide puhul tuleks rajada igale kraavile pinnaspais iga 30 cm kõrgusmuudu kohta, aga vähemalt 3 pinnaspaisu. See võimaldab hüdroloogiliselt eraldatud üksustes juhtida väljaku madalamas piirkonnas lumesulamise vee serpentiinina läbi ala nii, et pinnaspaisudega kraavid oleks üle ühe ühendatud erineval pool kesksest pinnaspaisude reast (joonis 30). Antud lahendus on väga hästi toiminud Ess-soo põhjapoolsel alal (Alad 5, 7, 11, kood J ja H, C), soodustades vee pikemat säilitamist korrastatud alal, kuid kahandades suuremat üleujutust.

Joonis 30. Serpentiinina ühendatud kraavid Ess-soos.

Laiuse jääksoo Laiuse jääksoo oli esimene mis sai uuringualadest korrastatud 2019.a oktoobriks ning seega on korrastamisjärgseid muutusi saanud jälgida peaaegu 4 aastat. Esimesel kahel aastal oli taimestiku kujunemine aeglane ja kohati mõjutas liiga kõrge veetase, aga samas soodustas see veelindude saabumist uurimisalale, kes levitavad ka taimede seemneid ja eoseid, aga ka väetavad ala. Väetamise efekt on tugev kevadel, mil laudteed on lausaliselt väljaheidetega kaetud ning paigal seismist nõudvate välitööde korral on kasuks vihmavarju või kapuutsiga kummimantli kasutamine sõltumata ilmast. Mõju avaldub selgelt ka madalaveelise veekogu kõrgendatud DN sisalduses, mis tõuseb ka sügisese rändeperioodi ajal, kuid sügisvihmade lahjendava toime tõttu pole sama tuntav kui kevadel. Alates 2022.a. algas kiire taimestiku areng, mis jätkus jõudsalt 2023.a. 2023.a. algas ka madalaveelise veekogu kallastel taimestiku laienemine veekogu suunas. Jääksoo lääneosas laiematele kraavidele rajatud lainerahusti (joonis 31) on oma eesmärki täitnud suurepäraselt ja soodustanud kiiret taimestiku laienemist kraavides. Laienenud on peamiselt pilliroog ja hundinui, aga ka tarnad ja kohati turbasamblad.

Joonis 31. Laiuse jääksoo korrastamisalad. Taimkatte arengut Laiuse jääksoos enne korrastamist, vahetult pärast korrastamist ja uuringu lõpuaastal 2023 septembrikuiste satelliidipiltide alusel iseloomustab joonis 32. Võrreldes algseisuga on oluliselt paremini taimestunud ala loodepoolne osa, madalaveelise veekogu kallastele on kujunenud kuni mõnekümne meetri laiune sootaimedega taimestunud kaldavöönd-õõtsik, vaatamata puurinde eemaldamisele on kirdepoolne osa on NDVI indeksi väärtuse järgi taimestunud juba paremini kui kontrollala, aga kui kontrollala indeksit mõjutab eeskätt puurinne ja villpea, siis loodepoolses osas on domineerivad sootaimed (tarnad, pilliroog, rabakarusammal, jõhvikas jmt). Jõhvika areaal laieneb aastas keskmiselt 40-50 cm võrra peenarde keskosast ääreala suunas ja moodustab kohati lausalise katte.

Joonis 32. Taimkatte muutused Laiuse jääksoos enne korrastamist (2018), korrastanmise ajal (2019) ja uuringuperioodi lõpus (2023) Sentinel-2 satellidipildi alusel NRG ja NDVI indeksitena väljendatuna. Kõima jääksoo Kõima jääksoo korrastati 2019.a. lõpuks. Ala oli juba eelnevalt peaaegu lausaliselt taimestunud ja vaid üksikutes kohtades turbavõtu aukude vahelistel tervikutel oli taimestumata laike. Korrastamistööde käigus eemaldati suuremad puud evapotranspiratsiooni kahandamiseks, pinnaspaisudega suleti kraavid ja väljavool turbavõtu aukudest. Tänu eelnevalt olemasolevale rabataimestikule taastus kogu alal korrastamise käigus paljandunud pinnas kiiresti. Madalamad alad kattusid nii nokkheina kui turbasammaldega (joonis 33), kõrgemad pinnaspaisud peamiselt kanarbiku, karusambla ja villpeaga (joonis 34).

Joonis 33. Pinnaspaisude rajamiseks turba võtmise auk (vasakul) ja endine kirdepoolne kogujakraav (Kõima-N väljavool) on turbasammalde,villpea ning nokkheinaga kattumas.

Joonis 34. Kõrgemad pinnaspaisud kattuvad villpea, kanarbiku, karusambla ja murakaga, madalamad servad ja turba võtmisel tekkinud lohud nokkheinaväljaga.

Joonis 35. Pinnaspaisude tõttu seisva veega kraavid ning turbavõtuaugud täituvad turbasammaldega, tervikutel märjemates piirkondades kanarbik hääbub.

Joonis 36. Edelapoolses osas kus veetase on kõrgem ja püsivalt maapinnale lähedal ka põuastel suvedel (maapinna kalle tagab vee pealevoolu) on ka suuremad pinnaspaisud peaaegu täielikult taimestunud.

Joonis 37. Edelaoas lausaliselt täidetud kogujakraav ning vee liikumist tõkestavad massiivsed pinnaspaisud hoiavad turbavõtuaukudes veetaset kõrgena ka kesksuvel ning tagavad soodsad tingimused kiireks taimestumiseks. Turbasammaldega kaetud areaal on nelja aastaga jõudsalt laienenud.

Joonis 38. Edelapoolne väljavool on aastaringselt kuiv, soost valguv vesi on leidnud endale tee metsa alla, kuhu valgub ühtlaselt laial alal. Maima jääksoo Maima jääksoo korrastamine toimus 2020.a. sügisel ja oli uurimisaladest esimene kus kasutati kõiki erinevaid korrastamisvõtteid (madalaveeline veekogu, pinnaspaisud kraavidel, kraavide lausaline täitmine, pinnaspaisud kraavidel ja turbasambla fragmentide laotamine, kraavide lausaline täitmine ja turbasambla fragmentide laotamine, turbavõtuaukude väljavoolude sulgemine). Kavandatud tegevused osaliselt ebaõnnestusid ebaõige veetaseme tõttu, kuid soovitust kõrgem veetase ei takista soostumist. Turbasammalde areng ja levik on alal liiga kõrge või muutliku veetaseme tõttu piiratud, pilliroo, tarnade, villpea ja nokkheina, üksikutes piirkondades ka jõhvika laienemine on viimase aastaga kiirenenud.

Taimestumist kiirendas kõige enam 2022 ja 2023.a. põuased suved, mis tagas taimede arenguks soodsama veetaseme. 2021.a. lõpus väga edukalt laienenud nokkheina areaali alal 2 (L) ja 10 (D) hävitas peaaegu täielikult sügisrände eel 100-200-pealine sookurgede parv, mis rebis taimed lausaliselt juurtega välja. Uurimisalal on kohatud merikotkast (sageli Ala 4 rabapoolsel küljel kõrgema männi ladvas), kuni 20 luigest koosnevat parve, koovitajaid, põtra, hunti ja pruunkaru. Uurimisala projekteerimisel/korrastamisel tehtud suurim eksimus oli liigse vee äravoolu planeerimine läbi olemasoleva osaliselt täidetud kogujakraavi. Kõrge veetaseme korral täidetud kraavis mudajas mass kerkib koos veetasemega ja takistab vee äravoolu, suvel alaneva veetaseme korral aga alaneb ka mudajas mass äravoolukraavis ja pigem soodustab kiiremat veetaseme alanemist turbaväljal. Maima eksimust võeti arvesse Ess-soos, kus kogujakraav täideti või sulgeti pinnaspaisudega ja liigvee äravooluks kujundati eraldi voolunõva serpentiinina läbi turbaväljade.

Joonis 39. Madalaveelise märgala (Ala 1, kood M) veetase jäi planeeritust madalamaks kuna soovitud veetaseme korral oleks kõik rabapoolsed väljakud veelgi sügavamalt üleujutatud olnud. Taimestumise seisukohast on veetase alal soodne ja pilliroo, hundinuia ning tarnade jõudne levik algas 2022.a. ja 2023.a. sügiseks on taimestumine peaaegu lausaline.

Joonis 40. Kraavide lausalise täitmisega ja samblafragmentide laotamisega alal (ala 3, K) on veetase liiga kõrge, valdava osa aastast on ala veega kaetud ja taimestunud on peamiselt täidetud kraavide kohad, kus maapind kerkib koos veetasemega. Siiski leidub hajusalt ka üksikuid elusaid turbasamblaid.

Joonis 41. Ebasoodsalt kõrge veetaseme puhul toimub taimestumine kiiremini just täidetud kraavide kohal kuna seal maapind kerkib koos veetasemega. Taimestikus domineerivad villpead, hundinui, tarnad, nokkhein ja mätaste vahel üksikuid turbasamblaid, mis on fragmentide laotamisest säilinud. Veega kaetud ala on luikede kasutuses.

Joonis 42. Pinnaspaisudega suletud kraavidega turbasambla fragmentide laotamise alal on taimestumine võimalik vaid soodsa niiskusrežiimiga vööndis. Liiga sügava veega alal toimub aeglane taimestumine pilliroo ning hundinuiaga. Sobivates tingimustes on turbasambla katvus hea ja sammal elujõuline.

Joonis 43. Pinnaspaisudega suletud kraavidega alal 5 (B) toimub looduslik taimestumine kõrge veetaseme tingimustes ja levivad madalsoole iseloomulikud liigid. Taimestumise kiirust toetab sel alal lahustunud lämmastikuga rikastunud põhjavee väljakiildumine.

Joonis 44. Kraavide lausalise täitmisega ja samblafragmentide laotamisega alal (ala 11, E) on veetase liiga kõrge, valdava osa aastast on ala veega kaetud. See ala vajas korrastamisel pinnase suuremamahulist tasandamist ja seetõttu pole pindmine turbakiht tihedalt alumiste kihtidega seotud ning liigub koos veetasemega kaasa. Taimestunud on peamiselt täidetud kraavide kohad, kus maapind kerkib kergemini koos veetasemega, aga taimestumine on ulatuslikum kui sarnaselt töödeldud alal 3 (K).

Joonis 44. Lausaliselt täidetud kraavidega alal (10, D) toimub iseeneslik taimestumine ebaühtlaselt. Kanarbik ja sinikas hääbuvad, villpea, nokkhein, tarnad ja jõhvikas laiendavad areaali. Vaatamata sobivatele niisketele laikudele ala sees ja külgnemisele samblafragmentide laotamise alaga, pole iseseisvalt turbasamblaid ilmunud.

Joonis 45. Pinnaspaisudega suletud kraavidega turbasambla laotamisega alal (9, P) oli esimesel aastal sammalde elulevus väga hea, kuid järgneval talvel kannatas kõrge lumesulavee uhtumise ja tugeva külmakohrutuse all. 2022/2023 kohrutuse kahju kordus. Siiski on kogu ala samblafragmentidega hajusalt

kaetud, püsivad koloniseerimistuumakesed tekkinud ning alal on esindatud paljud tüüpilised rabaliigid. Taimestumine on küll oodatust aeglasem, aga püsiv. Sel alal on niiskemal perioodil sambalaga paremini kattunud gaasirüngastes mõõdetud päevasel ajal ökosüsteemi hingamist ületavaid CO2 sidumise väärtusi.

Joonis 46. Võrdlusala on kõige kehvemini taimestunud. Alustaimestikus domineerivad üksikud hajusalt paiknevad villpeamättad, kraavi kallastel ka jõhvikas, samblikud. Kased ja männid kannatavad mineraliseerumise ning tuuleerosiooni tõttu paljanduvate juurte käes. Kuigi korrastamise käigus võrdlusala veetase tõusis, ei ole see veel oluliselt mõjutanud taimestumist.

Joonis 47. Pinnaspaisudega suletud kraavidega ja turbasambla fragmentide laotamisega ala (7, N), mis oli korrastamise eelselt tugevalt pilliroo ja noorte mändidega kaetud, on esimesest aastast saati olnud kõikuva veetasemega, aga juba esimesel sügisel risoomidest võrsunud varred takistasid lainetusel laotatud samblafragmente ja kattepõhku ära uhtuda ning pilliroo vahel esineb ohtralt turbasammalt, huulheina, kanarbikku. Esimeste aastate tulemus on paljulubav ja samblad elujõulised, kuid ebaselge on kas pikemas perspektiivis hakkab pilliroog turbasammalt varjutama või suudab sammal moodustada tugeva ühtlase katte.

Ess-soo Ess-soo ala korrastati 2021. a. sügisel ja selle käigus tehti võrreldes Maima alaga projektis mitmeid muudatusi. Eeldatavad veetasemed modelleriti iga ala lõikes, äravooluteed planeeriti serpentiinina väljakute keskosa kaudu, kohati säilitati üleujutuste vältimiseks avatud kogujakraavi lõike ning looduslikuma rabaosa ja freesturbavälja vahele rajati kogujakraavile veekogu. Korrastamistööde käigus tehti jooksvalt täiendusi vastavalt nivelleerimise tulemustele, lisati pinnaspaise kraavidele ning rajati madala pinnasvalliga eraldatud terrasseeritud väljak. Kuigi korrastamisest on seireperioodi lõpuks möödunud alla 2 aasta ja mõlemad korrastamisjärgsed suved on olnud äärmiselt põuased ning ebasoodsad turbasambla fragmentide siirdamisega korrastamiseks, on üldtulemused siiski lootustandvad ja metoodilises mõttes võib Ess-soo korrastamist edukaks näiteks pidada.

Joonis 48. Kuigi lausaliselt täidetud kraavidega turbasambla fragmentide laotamise ala nr. 2 (N) on arvestatava pikisuunalise nõlvakaldega ja külgnev avatud äravoolukraaviga, on see põuased suved kõige edukamalt üle elanud ja elujõulisi samblafragmente esineb lausaliselt. Sarnaselt samasugusele korrastamismeetodile Maima jääksoos, on ka siin kiiremini taimestunud just täidetud kraaviga osa. Samas domineerib täidetud kraavi osas villpea, turbasammalt leidub hajusalt kõikjal ning kuivemal osal on enam kanarbikku. Kohati esineb ka nokkheina laike.

Joonis 49. Pinnaspaisudega kraavid hoiavad ka arvestatava nõlvakalde korral edukalt veetaset üleval. Samblafragmentide laotamine sel alal (1, L) pole enamasti sama edukas kui kraavide lausalise täitmisega naaberalal (2, N), kuid sobiva niiskusega piirkonnas on kujunenud ulatuslik lausalise turbasambla katvusega ala. Taimestumine toimub edukamalt ka kraavide kallastel, aga traavidevahelistel väljakutel on taimestik hõre ning elusaid turbasamblaid vähem kui naaberalal. Kas samblafragmendid on säilitanud kahe põuase suve järel elujõu, selgub järgnevatel aastatel.

Joonis 50. Kõikidel aladel kuhu samblafragmente on laotatud, on hajusalt elujõulisi turbasablaga laigukesi ning seisva veega kraavilõikudes sageli ka vohavat turbasammalt (ala 9, E).

Joonis 51. Kuigi pealtnäha mõõdab eddy covariance mast nukralt tühja välja (ala 10, D) CO2 ja CH4 voogu, on siiski kogu alal hajusalt elujõulisi turbasambla laigukesi ja soodsamate aastate saabumisel võib sambla katvus kiiresti laieneda. Sarnaselt Maima jääksoos üleujutatud mudasele väljale (ala 11, E) on ka siin esimese 2 aasta taimestumine väga tagasihoidlik, aga tüüpiliselt toimub taimestumine alguses kiiremini just täidetud kraavide kohal. Pioneerliigiks villpea, kuivematel aladel kanarbik, hajusalt elus turbasamblaid, nokkheina, huulheina.

Joonis 52. Looduslikult kujunenud sootaimestiku säilitamine korrastamise ajal kiirendab veetaseme tõstmisel maapinna kattumist taimedega. Alal 11 (B) on pinnaspaisude abil veetaseme tõstmisel ja stabiliseerimisel jõhvikas laiendanud kaetavat areaali 60-70 cm võrra aastas, kraavides hõljuval mudal laiutavad villpead ning servas laiendavad kasvuala nokkheinad.

Joonis 53. Avatuks jäävatel kraavilõikudel haost või põhupallidest tõkete tekitamine/säilitamine on kasulik nii heljumi kahandamiseks kui kraavi kinnikasvaise kiirendamise seisukohast. Kogujakraav K-17 on tõketevahelisel lõigul täitumas turbasammalde, ubalehtede, soovõhkade, villpeade ja tarnadega.

Joonis 54. Kuigi pinnaspaisud kraavidel ja terrasseerimine madala eraldusvalliga hoiavad veetaset võrdväärselt teiste korrastatud väljakutega, on alal 4 (F) turbasammalde kasvama minek oluliselt kehvem. Üheks põhjuseks oli külmunud kängardes fragmentide laotamine pärast tugevat öökülma külmunud maapinnale, teiseks põhjuseks oli sel alal erandlikult esinev külmakohrutus 2022/2023 talvel ning kolmandaks põhjuseks 2022.a. 30. augustil esinenud erakordselt intensiivne sadu (Korelas mõõdeti 24 h jooksul 84 mm sademeid), mis tulvaveega uhtus ära Ess-soo uurimisala peamise ülevoolu mulde (P3) ning uhtus peenarde kõrgematesse osadesse nii samblafragmendid kui kattepõhu.

Joonis 55. Näide meandreeruvast paisudega suletud kraave ühendavast vooluteest (vasakul) ning liiga kõrget veetaset vältivast voolunõvast (paremal).

Joonis 56. Kogujakraavi võib sulgeda laiade turbaga täidetud lõikudega, kus sulgev lävend on kaetud taimede juurtega tihedalt läbikasvanud mätastega. Selline veekogu aitab hoida freesturbavälja otstele iseloomulikku kõrgemat serva niiskemana ja tagab kiirema taimestumise ning väiksema erosiooni, mis muidu kannaks turvast madalamal paiknevatele laotatud samblafragmentidele. KOKKUVÕTE Korrastamata jääksood olid olulised CO2 allikad. Enne korrastamist oli CO2 emissioon sõltuvalt aasta ilmastikust ja alast 4.7 (3.2 – 8.3) CO2-C t/ha*a. Metaani emissioon oli tagasihoidlik 0.09 t CH4-C t/ha*a. Toitainevaese rabaturbaga jääksoode naerugaasi emissioon oli samuti väike (0.0003 N2O-N t/ha*a) ja korrastamisejärgsel oluliselt ei muutunud. CO2 voog korrastamisjärgselt kahanes ja Laiuse jääksoos neli aastat pärast korrastamist jõudis aastabilansina süsinikuneutraalsuseni. Teistel korrastatud aladel oli aasta bilanss CO2 osas jätkuvalt emiteeriv 0.4-1.9 CO2-C t/ha*a. Kuigi gaasivood on suuremad suvekuudel (v.a. naerugaas, mil puudub selge aastaajaline käik), võivad külmumata pinnasega talvekuud oluliselt mõjutada gaasivoo aastast bilanssi. Süsihappegaasi sidumist mõjutab kõige enam fotosünteetiliselt aktiivne kiirgus (PAR), temperatuur (õhu ja pindmise 10 cm mullatemperatuur). Viimastel aastatel Eestis enam Keskkonnaagentuuri hallatavates ilmajaamades PAR ei mõõdeta ja ainsad teadaolevad pidevad PAR mõõtmised toimuvad hetkel RMK jääksoodes paiknevates mõõtekohtades. Ilma PAR pideva aegreata ei ole ökosüsteemi puhasgaasivahetuse (Net Ecosystem Exchange, NEE) usaldusväärne modelleerimine võimalik. Korrastamisjärgse seire periood vastavalt 4, 3, 2 ja 0 aastat on ebapiisav, et teha järeldusi meetodite tõhususe, taimestumise kiiruse või kasvuhoonegaaside voo kahanemise kohta. Esimestel aastatel mõjutab kasvuhoonegaase samblafragmentidega korrastataval alal põhu ja surnud fragmentide lagunemine. Äärmiselt suur määramatus on seotud ilmastikuga. Taimestumine kiirenes alates kolmandast korrastamisjärgsest aastast, kuid selgusetu on kui suurt rolli selle juures mängisid viimased kaks põuase suvega aastat. Kõikidel aladel kus rakendati turbasamblafragmentide laotamist, on vähemalt hajusalt elusaid turbasambla kogumeid ja vähestel aladel moodustavad ka väiksemaid lausalise katvusega alasid. Meetodi edukust kahandas projekteerimisviga veetaseme osas Maima jääksoos ning vahetult korrastamisele järgnenud 2 väga põuast suve Ess-soos. Aladel kus turbasambla fragmente ei laotatud, iseseisvalt turbasamblaid kasvama hakanud ei ole. Samuti on sambla fragmentide laotamisega aladel

rohkem rabale iseloomulikke liike. Esimeste aastate tulemused näitavad, et samblafragmentide laotamise teel korrastatavate jääksoode puhul taimestuvad nii üleujutatavate kui põuast mõjutatud aladel kiiremine lausaliselt täidetud kraavidega alad, aga pikemas ajaskaalas ei pruugi see kehtida. Ka pinnaspaisudega kraavide kallastel laieneb taimestik. Kriitiline on siiski sobilik veetaseme vahemik ja suvine niiskuse olemasolu, see sõltub aga nii võimalikust külgnevast tagamaast kui konkreetsete aastate ilmastikust. Drooniseire on väga tõhus abivahend korrastatava alaseisundi eelnevaks kaardistamiseks, seirealade optimaalseks valikuks, allikaliste alade tuvastamiseks, pinnaspaisude lekete avastamiseks ning ligikaudseks pinnase niiskuse määramiseks. Taimkatte kaardistamiseks on võimalik kasutada k-means klasterdamisel põhinevat lähenemist koos välitööde käigus klassidele sisu andmisega või suure õpetusandmestiku olemasolul masinõppe meetodil (random forest, bagging jmt). Pikaajalise homogeense aegrea saavutamine taimkatte dünaamika kaardistamiseks on väga kallis (tehniliselt ning tööjõukulult), aeganõudev ja keeruline, mida omakorda mõjutab tehnoloogia kiire areng ning sensorite muutus. Satelliidiseire on jääksoode korrastamise tulemuslikkuse jälgimiseks asjakohane, kuid kasutegur on suurem pika seireperioodi puhul. Lühikese perioodi puhul jääb muutuste suhtes väga tundliku jääksoo dünaamika oluliselt kiiremaks (nt. veetasemete muutus ja üleujutatavate alade ulatus) kui pilvevabade piltide saamine satelliitidelt. Samuti eeldab selline seire suuremate seireruutude rakendamist maapealses seires, et andmestik oleks võrreldav piksli suurusega. Paljude klassikaliste indeksite kasutamise muudab keeruliseks ka jääksoodele sagedane olukord, kus taimede vahelt paistab vesi, mitte maapind. See raskendab ka muidu pilvedest vähem mõjutaud radari andmestiku kasutamist. Dendrokronoloogia abil on võimalik näha turbaväljade rajamisega kaasnevaid mõjusid, raskustega luua kronoloogiaid jääksoos kasvavate puude osas (puud erivanuselised ja seega muutuva nooruskasvuga ning samaaegselt kiire keskkonnatingimuste muutusega), aga veetaseme tõstmise avaldumise tuvastamiseks ei ole 3-4 aastat piisav. Männid jätavad ebasoodsates tingimustes aastarõngaid vahele ja seega nii lühikesed perioodid ei allu kronoloogia loomisele. Lahustunud orgaanilise süsiniku ja lämmastiku kontsentratsioonid korrastamisjärgselt küll kuni kaheks aastaks tõusid, kuid selgusetu on seos ilmastiku (kuumad põuased suved) ja korrastamistööde osakaalu osas. Ärakanne on aga tagasihoidlik kuna vee äravoolu esineb uuritud aladel 2-4 kuud aastas ja needki madalama kontsentratsiooniga hilissügisel ja varakevadel. Vooluhulga ja kontsentratsiooni järgi hinnates on süsiniku ärakanne DOC kujul jääksoodest vahemikus 62-87 kg/ha*aastas. Lagunemiskatse esmased tulemused näitavad, et peamine massikadu toimub väga kiiresti esimese aasta jooksul ning selles mängib omakorda suurimat rolli esimeste kuude jooksul leostumiskadu. Veetase ja taimestik mõjutavad lagunemist oluliselt. Erinevate taimsete materjalide (varis, peenjuured, erinevad liigid) lagunemiskatsete tulemused selguvad kolme aasta pärast.

7. PROJEKTIGA HAAKUVAD TEADUSTEEMAD, GRANDID, DOKTORI- JA MAGISTRITÖÖD, JÄRELDOKTORITE UURIMISTEEMAD, LEPINGUD: Teavitustegevus: Lühiartikkel projekti eesmärkidest Eesti Loodus 8/2017, lk. 5. http://www.eestiloodus.ee/arhiiv/Eesti_Loodus08_2017.pdf ja Rahvusvahelise Märgalade Kaitse Grupi kuukirjas IMCG Bulletin, June, 2017 pp. 13-14 http://www.imcg.net/media/2017/imcg_bulletin_1706.pdf. 04.10.2017 TÜ geograafia osakonna seminar projekti eesmärgi ja hetkeseisu tutvustamiseks ning võimalike täiendavate huviliste (omafinantseeringu korras) kaasamine ülikooli teistest uurimisgruppidest. 18.-19.oktoober 2017 Toilas Keskkonnaministeeriumi turbaümarlaual projekti eesmärgi ja hetkeseisu tutvustamine. Suuline ettekanne: A. Kull & G. Veber, Abandoned peat extraction sites – will future be wetter and better? 10.-12.10.2018 Tartu, 18th Baltic Peat Producers Forum. Jääksoode korrastamisega seonduvat on laiema üldsuse teavitamiseks käsitletud populaarteaduslikus väljaandes "Samblasõber" nr 23, 2020, lk 10-15: https://sisu.ut.ee/sites/default/files/samblasober/files/samblasober_23_0.pdf Esinemine ERR Aktuaalne Kaamera, Osoon ja Vikerraadios intervjuudega.

Magistritööd ja doktoritööd Ott Toomsalu, 2019. Jääksoodes toimuvate muutuste analüüsimine LiDAR andmetel. Magistritöö. Kaitsutud Tartu Ülikooli geograafia osakonnas. https://dspace.ut.ee/handle/10062/65031 MarjanSadat Barekaty, 2021. Compare the performance of applying Machine Learning concepts to landcover classification models using very high-resolution UAV data. Magistritöö. Kaitsutud Tartu Ülikooli geograafia osakonnas. https://dspace.ut.ee/handle/10062/72820 Kärt Erikson, 2022. Veerežiimi häiringute ja ilmastiku mõju hariliku männi (Pinus sylvestris L.) radiaalsele juurdekasvule Lehtmetsa soo näitel. Magistritöö. Kaitsutud Tartu Ülikooli geograafia osakonnas. https://dspace.ut.ee/handle/10062/82873 Tauri Tampuu doktoritöö: Application of spaceborne SAR polarometry and interferometry for landscape ecological studies in bogs (Tartu Ülikool, kaitstud 2022.a. augustis). Artiklid Birgit Viru, Gert Veber, Jaak Jaagus, Ain Kull, Martin Maddison, Mart Muhel, Alar Teemusk, and Ülo Mander, 2017. Winter nitrous oxide and methane emissions from drained peatlands. Geophysical Research Abstracts, Vol. 21, EGU2019-15964. The abstract identification number EGU2019-15964. https://meetingorganizer.copernicus.org/EGU2019/EGU2019-15964.pdf?pdf Tampuu, Tauri; Praks, Jaan; Uiboupin, Rivo; Kull, Ain (2020). Long Term Interferometric Temporal Coherence and DInSAR Phase in Northern Peatlands. Remote Sensing, 12 (10), ARTN 1566. DOI: 10.3390/rs12101566 Tampuu, T.; Praks, J.; Kull, A. (2020). Insar Coherence for Monitoring Water Table Fluctuations in Northern Peatlands. International Geoscience and Remote Sensing Symposium (IGARSS). IGARSS, 4738−4741. DOI: 10.1109/IGARSS39084.2020.9323709 Burdun, Iuliia; Kull, Ain; Maddison, Martin; Veber, Gert; Karasov, Oleksandr; Sagris, Valentina; Mander, Ülo (2021). Remotely Sensed Land Surface Temperature Can Be Used to Estimate Ecosystem Respiration in Intact and Disturbed Northern Peatlands. Journal of Geophysical Research Biogeosciences, 126 (11), e2021JG006411. DOI: 10.1029/2021JG006411 T. Tampuu, J. Praks, A. Kull, R. Uiboupin, T. Tamm, K. Voormansik (2021).Detecting peat extraction related activity with multi-temporal Sentinel-1 InSAR coherence time series. International Journal of Applied Earth Observation and Geoinformation,Vol. 98,102309, https://doi.org/10.1016/j.jag.2021.102309 Tampuu, Tauri; De Zan, Francesco; Shau, Robert; Praks, Jaan; Kohv, Marko; Kull, Ain (2022). Can Bog Breathing be Measured by Synthetic Aperture Radar Interferometry. 2022-July, 16−19. DOI: 10.1109/IGARSS46834.2022.9883421. Kull, Anne; Kikas, Tambet; Penu, Priit; Kull, Ain (2023). Modeling Topsoil Phosphorus—From Observation-Based Statistical Approach to Land-Use and Soil-Based High-Resolution Mapping. Agronomy, 13 (5), 1183. DOI: 10.3390/agronomy13051183 Palviainen, M., Könönen, M., Peltomaa, E., Pumpanen, J., Ojala, A., Hasselquist, E., Laudon, H., Ostonen, I., Renou-Wilson, F., Kull, A., Veber, G., Mosquera, V., and Laurén, A.: Processes affecting lateral carbon fluxes from drained forested peatlands, EGU General Assembly 2023, Vienna, Austria, 24–28 Apr 2023, EGU23-6367, https://doi.org/10.5194/egusphere-egu23-6367, 2023. Tampuu, T.; Praks, J.; De Zan, F.; Kohv, M.; Kull, A. (2023). Relationship between ground levelling measurements and radar satellite interferometric estimates of bog breathing in ombrotrophic northern bogs. Mires and Peat, 29, 1−28. DOI: 10.19189/MaP.2022.OMB.Sc.1999815

8. Projekti juht (nimi): Ain Kull

Allkiri: allkirjastatud digitaalselt

Kuupäev: allkirjastatud digitaalselt

9. Taotleja allkirjaõigusliku esindaja kinnitus aruande õigsuse kohta (nimi, amet): Ain Kull, kaasprofessor

Allkiri: allkirjastatud digitaalselt

Kuupäev: allkirjastatud digitaalselt

NB! Aruanne esitada elektrooniliselt aadressil [email protected]

Agronomy 2023, 13, 1183. https://doi.org/10.3390/agronomy13051183 www.mdpi.com/journal/agronomy

Article

Modeling Topsoil Phosphorus—From Observation-Based

Statistical Approach to Land-Use and Soil-Based

High-Resolution Mapping

Anne Kull 1, Tambet Kikas 2, Priit Penu 2 and Ain Kull 3,*

1 Institute of Agricultural and Environmental Sciences, Estonian University of Life Sciences, Kreutzwaldi 1,

51006 Tartu, Estonia 2 Centre of Estonian Rural Research and Knowledge, Teaduse 4, 75501 Saku, Estonia 3 Institute of Ecology and Earth Sciences, University of Tartu, Vanemuise 46, 51003 Tartu, Estonia

* Correspondence: [email protected]

Abstract: Phosphorus (P) is a macronutrient that often limits the productivity and growth of terres-

trial ecosystems, but it is also one of the main causes of eutrophication in aquatic systems at both

local and global levels. P content in soils can vary largely, but usually, only a small fraction is plant-

available or in an organic form for biological utilization because it is bound in incompletely weath-

ered mineral particles or adsorbed on mineral surfaces. Furthermore, in agricultural ecosystems,

plant-available P content in topsoil is mainly controlled by fertilization and land management. To

understand, model, and predict P dynamics at the landscape level, the availability of detailed ob-

servation-based P data is extremely valuable. We used more than 388,000 topsoil plant-available P

samples from the period 2005 to 2021 to study spatial and temporal variability and land-use effect

on soil P. We developed a mapping approach based on existing databases of soil, land-use, and

fragmentary soil P measurements by land-use classes to provide spatially explicit high-resolution

estimates of topsoil P at the national level. The modeled spatially detailed (1:10,000 scale) GIS da-

taset of topsoil P is useful for precision farming to optimize nutrient application and to increase

productivity; it can also be used as input for biogeochemical models and to assess P load in inland

waters and sea.

Keywords: agricultural land; geographical information system; interpolation; land use; machine

learning bagging model; soil phosphorus mapping

1. Introduction

Phosphorus (P) is an important macronutrient for plant growth, and the primary role

of P in plants is to store energy needed for plant growth and reproduction produced by

photosynthesis. Most phosphorus is relatively immobile in soil and even P from phosphate

fertilizers will readily react with soil minerals, making it less available for plants.

Soil phosphorus consists of two forms: organic (non-plant available) and inorganic

(plant available) P. Approximately 29 to 65 percent of total soil phosphorus is in organic

forms, which is not plant available, while the remaining 35 to 71 percent is in inorganic

forms [1] available to plants mainly in the form of phosphate as labile or occluded forms

of P [2]. Organic forms of phosphorus include dead plant/animal residues and soil micro-

organisms. Soil microorganisms play a key role in processing and transforming these or-

ganic forms of phosphorus into plant-available forms, especially in natural ecosystems.

The inorganic form (plant-available form of P) is highly reactive and can be tied up to

chemical compounds (e.g., iron, aluminum, and calcium) in soils as absorbed phosphorus

[3]. The inorganic P has to be dissolved into a solution (P-solution) for plants to be able to

Citation: Kull, A.; Kikas, T.; Penu, P.;

Kull, A. Modeling Topsoil

Phosphorus—From

Observation-Based Statistical

Approach to Land-Use and

Soil-Based High-Resolution

Mapping. Agronomy 2023, 13, 1183.

https://doi.org/10.3390/agron-

omy13051183

Academic Editor: Alberto San

Bautista

Received: 13 February 2023

Revised: 26 March 2023

Accepted: 12 April 2023

Published: 22 April 2023

Copyright: © 2023 by the authors.

Licensee MDPI, Basel, Switzerland.

This article is an open access article

distributed under the terms and

conditions of the Creative Commons

Attribution (CC BY) license

(https://creativecommons.org/license

s/by/4.0/).

Agronomy 2023, 13, 1183 2 of 24

uptake it. P absorption (incorporation of plant-available P in soil solution into soil clay min-

erals, such as Fe, Al, and Ca oxides) is a fast process, and absorbed P can be released into

soil solution (plant-available P) via desorption. P absorption (retention) is higher in soils that

contain more clay and iron and aluminum oxides [4]. A lower soil pH favors P absorption

by aluminum and iron (pH below 5.5), and a higher pH (over 7.5) enhances plant-available

P fixation by calcium [4,5].

The solution P pool (plant-available P) is very small and ranges from 0.001 mg/L to 1

mg/L [6]. In general, roots absorb phosphorus in the form of orthophosphate but can also

absorb certain forms of organic phosphorus (e.g., glycerophosphate, lecithin, and phytin)

[7]. The active (labile) pool is P in the solid phase that is more readily released into soil

solution, which is the water surrounding soil particles. As plants take up phosphate, the

concentration of phosphate in the solution decreases and some phosphate from the active

P pool is released. An active P pool consists of mainly P absorbed into soil clay minerals

(Fe, Al, and Ca oxides) and easily mineralizing organic P. Labile inorganic P in soils is

predominantly present as specifically adsorbed orthophosphate. This P is in equilibrium

with the soil solution and acts to buffer the soil solution against concentration change [8,9].

Thus, P is desorbed into solution in response to plants’ P uptake or adsorbed from solution

when the P concentration is raised by mineralization or by the addition of a fertilizer ma-

terial.

Soil age is also an important factor influencing P availability, with P becoming in-

creasingly limiting in ancient soils because it gradually disappears through leaching and

erosion [10,11]. In arable lands, P is also partly removed with yield, and additional P fer-

tilizer application is needed to maintain soil fertility. Contrarily to agricultural lands, in

natural ecosystems, P uptaken by plants will mainly remain in place after plants’ death,

and after decomposition and mineralization, it will be available in topsoil for other plants.

Climatic conditions affect P availability as a higher temperature increases the activity of

soil microorganisms and contributes to faster organic matter decomposition and P minerali-

zation. On the other hand, a higher temperature also increases P sorption. In well-aerated soils,

P is released faster than in saturated wet soils. A neutral soil pH (6–7.5) is the best for P avail-

ability, while pH values below 5.5 and between 7.5 and 8.5 limit P availability to plants due to

fixation by aluminum, iron (in case of lower pH), or calcium (in case of higher pH) [5]. Wet

soil conditions decrease soil pH, which increases P sorption to Fe and Al oxides, but flooding

the soil reduces P sorption by increasing the solubility of phosphates that are bound to alumi-

num and iron oxides and amorphous minerals [12]. This aspect has to be considered when

mapping soil P in high latitudes rich in histosols across all land uses and soil types.

Soil phosphorus maps are used to estimate P availability to calculate the need for

fertilization and to model P loss. As only a small amount of P fertilizers used may be read-

ily available for plant uptake due to P sorption and P loss via erosion and runoff, the use

of fertilizers should be based on deeper knowledge about soil properties. Improper P fer-

tilizer use contributes to higher P loss and eutrophication rather than better plant growth.

Detailed and spatially high-resolution soil P measurements covering equally all land

uses and soil types are not usually available for larger areas to produce accurate data-

driven P maps, which are based only on soil sampling data. If high-spatial-resolution P

sampling data are available (e.g., for small areas), geostatistical methods, such as interpo-

lation techniques (e.g., ordinary kriging, co-kriging, or regression kriging) are usually

used to produce soil P maps [13,14]. In the case that soil P sampling grid is not dense

enough, other available environmental variables are also used to predict P content by us-

ing machine learning algorithms, co-kriging, or some other multi-criteria algorithms, in-

cluding multiple regression models. Hybrid geostatistical methods, which incorporate

spatially distributed soil observations and readily available ancillary environmental data

(e.g., topographic variables and satellite data) are recommended instead of univariate

methods, such as ordinary kriging, in cases where natural soil-forming processes are com-

plex or landscapes have high anthropogenic influence [15,16].

Agronomy 2023, 13, 1183 3 of 24

Gaussian process regression (GPR) works principally like co-kriging with a number

of covariates. Ballabio et al. [17] applied GPR to create maps of different LUCAS topsoil

chemical properties (pH, P, N, CaCO3, K, etc.), including indexes calculated from remote

sensing data, meteorological parameters, XY coordinates, and topographic, land-use, and

geological variables.

Recently, machine learning algorithms are increasingly used in predictive mapping

of soil parameters over larger areas as the number of available variables and computing

capacity are increased. Matos-Moreira et al. [18] used a machine learning tool (Cubist) to

develop rule-based predictive models from a calibration dataset. Covariates included pe-

dological, geological, agricultural, terrain and geophysical-related attributes. Rossel and

Bui [19] used a machine learning algorithm (Cubist) to predict total phosphorus in six

different soil layers and ordinary kriging to predict the residuals at each of the standard

depths in Australian soils. To derive the final estimates of total P, the predictions from

Cubist and the kriging estimates of the residuals were summed. The spatial modeling was

performed on 50 bootstraps, and 90 m grid size was used for modeling. Hosseini et al. [20]

used different statistical and machine learning algorithms, such as genetic algorithm (GA),

artificial neural network (ANN), fuzzy inference system (FIS), adaptive neuro-fuzzy in-

ference system (ANFIS), partial least squares (PLS), principal component regression

(PCR), ordinary least squares (OLS), and multiple regression (MR), to determine the best

model to predict soil P in Iran. Four predictive variables (clay, sand, soil organic matter

(SOM), and pH) were used to predict soil P, and the best results were obtained by PLS

(among the regression models) and GA and ANN (among the intelligent models).

Esfandiarpour-Boroujeni et al. [21] used different methods (decision tree (DT), ran-

dom forest (RF), artificial neural network (ANN), and support vector machine (SVM)) for

predicting soil classes (different WRB classification levels) based on environmental varia-

bles (planar and profile curvature, aspect, elevation, slope, catchment area, topographic

wetness index (TWI), LS factor, NDVI, etc.) and expert knowledge from soil scientists

(presence or absence of soil horizons and other soil properties). The SVM algorithm had

the highest overall accuracy for prediction of all qualitative soil properties. The ANN al-

gorithm showed good performance in predicting some quantitative variables (e.g., pure

clay percentage), and the DT algorithm had the lowest uncertainty value.

However, most of these studies that have been conducted to predict soil P with good

results have used detailed measured data for smaller areas [16,18,20–22] or small-scale

maps over larger areas [17,19,23].

The aim of this study was to create a high-resolution (1:10,000) topsoil plant-available

P map for the entire Estonia (area 45,339 km2) covering all land uses and soil types based

on datasets with highly unbalanced data availability and spatial resolution across land-

use categories. The available soil sampling data are usually collected for different pur-

poses and, thus, have different spatial and temporal resolutions, usually covering agricul-

tural areas mainly on mineral soils well while being extremely scarce over natural ecosys-

tems (e.g., forests and wetlands) on histosols and other less fertile soils. Based on the typ-

ical large bias of sampling data along soil types and land-use categories, we hypothesized

that topsoil P content can be sufficiently accurately mapped with geostatistical methods

(e.g., kriging interpolation) in arable lands, but the use of machine learning algorithms

increases prediction accuracy, especially in less intensively managed land-use categories

(e.g., short-term and permanent grasslands), while these established relationships are not

applicable for natural ecosystems (e.g., forests and wetlands), where more robust models

such as the two parameter ordination method should be applied.

Our study objectives were to establish the relationships between observed soil plant-

available P and environmental variables by soil and land-use classes, to determine the

most effective predictive factors related to soil P, establish a cost-effective modeling ap-

proach, assess the accuracy of different modeling methods, and create a high-resolution

topsoil P map at the national level.

Agronomy 2023, 13, 1183 4 of 24

2. Materials and Methods

2.1. Soil Phosphorus Data

Soil phosphorus data available for agricultural soils were obtained from the PANDA

database, which contains regular soil monitoring and voluntary soil sampling data by

farmers. The database contains soil sampling data collected from agricultural lands for

the period 2005–2021 and is managed by the Centre of Estonian Rural Research and

Knowledge. For this period in the PANDA database, there are data about 387,904 compo-

site soil samples all over Estonia, including topsoil phosphorus content (mg/kg). Compo-

site soil samples in the PANDA database were collected by licensed personnel following

the same prescribed sampling strategy (prescribed sampling route and sampling density,

depth, and volume) and analyzed in the same laboratory using the Mehlich-3 method. All

soil samples from arable lands were collected in late summer or autumn at minimum 2–3

months after last fertilization and, in case of winter crops in crop rotation, before the sow-

ing of crops. In the case of application of organic fertilizers, the soil samples were collected

not earlier than 6 months after fertilization. The same field was resampled in each 4–5

years. Up to 2012, the Mehlich-3 colorimetric and, since 2013, the Mehlich-3 inductively

coupled plasma (ICP) analysis method were used to determine soil available phosphorus.

The Mehlich-3 P extraction method is the main method used for estimating plant-available

P in Estonia since 2004. This method is robust, provides the advantage of multielement

analysis [24,25], and, thus, is well-suited for long-term monitoring of Estonian agricultural

soils, where pH is in the range between 5 and 7 for 75.9% of the samples. The PANDA

dataset from 2021 (29,261 soil samples) was not included in model building but used as

an additional independent test dataset to evaluate model performance.

There are no similar comprehensive topsoil phosphorus content databases for other

land-use categories available in Estonia. Comparable values of topsoil phosphorus content

for other land-use types (forest, wetland, peat extraction areas, and quarries) on different

soil types were searched through a literature review of scientific papers [26–44], reports

[45,46], and the Estonian Environmental Monitoring System and supplemented by origi-

nal unpublished datasets of the authors. Therefore, the datasets for these land-use catego-

ries vary by sample size, sampling, and analysis methods. For all samples from agricul-

tural soils, the Mehlich-3 method was applied, while the topsoil P concentration in the

forest and peatland soil samples collated from multiple studies and literature sources

were determined with various methods (Aqua Regia, Olsen, and Kjeldahl), and, thus, con-

version coefficients based on Kulhánek et al. [47] and Wolf and Baker [48] were used prior

to statistical analysis to convert phosphorus content to the same level with the Mehlich-3

method.

2.2. Land Use and Land Cover

Land-use data were combined from the ARIB (Agricultural Registers and Information

Board) 2020 database, where main crop types for agricultural lands are registered. Data on

natural grasslands were taken from the EELIS (Estonian Nature Information System) database

for semi-natural grassland layer, and all other land-use types (wetlands, forests, mining areas,

settlement, waterbodies, etc.) were taken from the ETD (1:10,000, Estonian Topographic Data-

base) (Estonian Land Board, 2020). For topographic information (ground elevation), a 10 m

resolution digital elevation model generated from LiDAR data (Estonian Land Board, 2021)

was used.

Crop data by year were obtained from the ARIB database. In the ARIB database, all

crops that are grown in agricultural massive (complex of neighboring fields) are listed in

alphabetical order. Agricultural massive is an agricultural unit that is in the agricultural

registry and can be applied for EU agriculture support (area-related aid). The ARIB crops

were classified into 8 categories: natural grassland, fallow, cultural grassland, permanent

cultures (e.g., orchards and berries), crop (e.g., rye, wheat, and barley), legumes (e.g., pea,

Agronomy 2023, 13, 1183 5 of 24

bean, and lentil), technical cultures (rapeseed, flax, and hemp), and vegetables (e.g., po-

tato, carrot, and cabbage). The main crop type was selected for each agricultural massive

from the ARIB crop list by prioritizing the crops by their probability of having the largest

share of this agricultural massive (e.g., the area for crop growing is probably larger than

the area for vegetables). From the EELIS database, data on semi-natural grasslands were

obtained and added to the natural grassland category. These areas not covered by the

ARIB data, and EELIS semi-natural grasslands were obtained from the ETD database and

additional categories of wetland, peat mining areas, forest, settlements, and waterbodies

were added. Settlements and waterbodies were omitted from the further analysis as we

did not estimate phosphorus for these areas.

We distinguished 11 land-use classes, with 8 of which being agricultural land catego-

ries (crop types) for which topsoil P values had been obtained from the PANDA database

(permanent grassland, fallow, cultivated grassland, permanent cultures (e.g., orchards and

berries), cereals (e.g., rye, wheat, and barley), legumes (e.g., pea, bean, and lentil), technical

cultures (rapeseed, flax, and hemp), vegetables (e.g., potato, carrot, and cabbage)). Three

non-agricultural land-use types (forest, wetland, and peat extraction area) were included in

the models to allow P value prediction for the entire Estonia (excluding settlement areas,

waterbodies, and mining areas).

2.3. Soil Data

The Estonian Soil Map (1:10,000) in digital form based on extensive field surveys in-

cludes information about soil type, fertility, texture, coarse fragments, and rock content

by layers. Kmoch et al. [49] used the Estonian Soil Map to develop EstSoil-EH: A high-

resolution eco-hydrological modelling parameters dataset for Estonia where soil map in-

formation was made machine-readable and used to predict soil properties (such as clay,

silt and sand content, organic carbon content, and bulk density) by pedotransfer functions

(PTF). In our study, we used the Estonian Soil map data obtained from Kmoch et al. [49].

Kmoch et al. [49] used the SAGA GIS functions (Conrad et al., 2015) to calculate the mean,

median, and standard deviation of several topographic factors (slope, USLE slope length

and steepness factor LS, topographic wetness index (TWI), and topographic ruggedness

index (TRI)) and environmental variables (share of drained area and different land-use

types within the soil polygon) as predictor variables in a random forest model to predict

soil organic carbon content (SOC). The 5 m resolution LiDAR-based DEM provided by the

Estonian Land Board was used for deriving topographic factors. In our study, we modeled

only topsoil phosphorus content; thus, we used only top-layer soil variables (bulk density

(bd1); clay (clay1), silt (silt1), and sand (sand1) contents; soil organic carbon content (soc1);

hydrologic conductivity (k1); and available water content (awc1)) from Kmoch et al. [49].

Quaternary sediment deposit map in a scale of 1:400,000 (Estonian Land Board, 2020)

was used as an indicator of soil parent material. This information was classified into 9

categories (alluvial, glaciolacustrine, glaciofluvial, till, aeolian, wetland sediment, shallow

quaternary sediment, lacustrine, and marine deposits).

2.4. Data Processing

2.4.1. Soil-Type Aggregation

In the Estonian Soil map (1:10,000), there are 119 soil types distinguished according

to the Estonian national soil classification system. The Estonian soil classification and soil

texture types based on the Kachinsky texture system [50] cannot be directly converted to

WRB or FAO classification as the soil types are distinguished based on different principles

[49]. To avoid potential conversion-related issues, we made all our calculations based on

the original texture and soil-type classes.

The original Estonian soil types were in some cases grouped into larger categories

based on genesis and expected similarities in soil properties related to phosphorus content

(soil wetness, gleyic processes, erosion, anthropogenic influences, etc.) to increase soil P

Agronomy 2023, 13, 1183 6 of 24

sample number per soil group. Median phosphorus content by soil types or groups was

calculated (Soil_MedP variable in Supplementary Table S1), and differences between soil

groups were tested by using the Kruskal–Wallis and Dunn tests with Bonferroni multiple

comparison correction in R. There were no statistically significant differences between dif-

ferent Gleysols (Gh, GI, Gk, Go, and Gor). Due to the limited number of samples and the

mixed origin of soils in the anthropogenic soil category (mainly quarries and recultivated

quarries), these soils also had no statistically significant difference from many other soil

categories. Altogether, 26 soil categories were distinguished.

Soil texture was classified into 12 categories based on Estonian soil texture types [50]

as some of these categories could not be directly converted to the USDA texture system. An

approximate reference was used as listed (Estonian texture codes given in brackets): gravel

(kr), gravel limestone (r), sand (s), sandy loam (ls1), loamy sand (ls2, ls3, and sl), fine sand

(pl), clay (s), peat (t), fine sandy loam (tls), and fine loamy sand (tsl). Differences in median

phosphorus content between soil texture categories were calculated (Texture_MedP variable

in Table S1), and differences between categories were tested by using the Kruskal–Wallis

and Dunn tests in R.

2.4.2. Data Analysis and Modeling

Three different approaches were used to model phosphorus content: (1) spatial inter-

polation (kriging), (2) statistical-empirical soil and land-use ordination scalar, and (3) ma-

chine learning bagging (bootstrap aggregation) model.

Spatial interpolation techniques (kriging in our case) create continuous smooth raster

surfaces based on the predictive variable values. Soil and land-use ordination and bagging

methods need some data units (polygons) for which all predictive variables are available.

Therefore, for model building, only variables that were available for all study area (entire

Estonia) could be used. For topsoil P content prediction, we used the unique units created

by intersecting soil and land-use polygons.

The ArcGIS Pro (2.9) tools were used to prepare units (so-called soil–land use poly-

gons) for P prediction. Soil map layer enriched with the predicted soil parameters (EstSoil-

EH) was intersected with land-use/land-cover layer (composed from the ARIB, EELIS, and

ETD data as described above) to create a layer that contains information about both soil

and land-use parameters. Settlements and waterbodies were excluded from the prediction

dataset. The output layer contained over 3 million polygons. Sediment (parent material)

information was added to the soil–land use polygons by spatial join function with the

largest share option. The soil–land use polygons central point X and Y coordinates were

calculated, and the central point elevation was extracted from a 10 m resolution LiDAR

DEM to represent longitude, latitude, and elevation variables in the bagging model.

Kriging

The most widespread approach to create a continuous surface from spatially well-

distributed point measurements is interpolation [51–53]. Kriging is among the most pop-

ular spatial interpolation techniques because it can give the best linear unbiased predic-

tion if suitable parameters are selected. This method is widely used in regional hydrolog-

ical and spatial nutrient runoff models, but also for mapping spatial distribution of nutri-

ents in detailed field-level test areas.

In our study, topsoil phosphorus content for agricultural areas all over the Estonia

was interpolated based on the measured phosphorus values of the PANDA database by

using simple kriging and stable variogram method in the ArcGIS Pro Geostatistical Ana-

lyst extension. Normal score transformation was used to transform the measured phos-

phorus data to follow a univariate standard normal distribution. Kriging is using Gauss-

ian process to estimate the mean value of a dependent variable and thus normal distribu-

tion of data are very important. Kriging assumes normal distribution and stationarity of

data (close points should have quite similar values, with low variance nearby). Spatial

Agronomy 2023, 13, 1183 7 of 24

trend was investigated based on the scatter plots of P values against the X and Y coordi-

nates, but, as no trend was detected, trend removal was not applied in the kriging models.

Two kriging models were built in order to compare the kriging results with other models:

the Kriging1 method, which used all available sample points to build a model to get the

best results, and the Kriging2 method, where the original phosphorus dataset was divided

into training (75%) and test (25%) datasets by using the ArcGIS Geostatistical Analyst ex-

tension tool Subset Features, which created a random subset from the dataset. The

Kriging2 model served only a model validation purpose.

Soil–Land Use Ordination

As the PANDA database covers only agricultural lands, kriging interpolation predic-

tions outside agricultural areas are not valid as the topsoil phosphorus content is largely

depending on anthropogenic factors, such as fertilizer application, pollution, and crop

harvesting, which are different in agricultural lands and cannot be interpolated or extrap-

olated to other land-use categories without taking into consideration other factors.

Ordination can be considered a synonym for multivariate gradient analysis. Ordina-

tion methods are mainly used in community ecology, where multiple spatial axes are used

to arrange or order multiple variables (of species and/or sample units) along gradients.

The most famous ordination methods are principal component analysis (PCA), factor

analysis (FA), detrended correspondence analysis (DCA), and other similar multidimen-

sional statistical analysis methods [54,55].

According to the exploratory data analysis, soil and land-use types were the main

independent factors influencing phosphorus content in topsoil. The median values of P

content of different soil types, land-use types, and their combinations were calculated to

create a land-use/soil-type matrix. Soil and land-use types were ordered ascendingly

based on the median topsoil phosphorus values. Some soil-type/land-use combinations

had no or limited number of measured soil P values. To interpolate topsoil phosphorus

values for these soil–land use combination without sufficient sampling data, the so-called

soil–land use ordination surface was created. Soil type was used as the X axis and land-

use type was used as the Y axis to generate Cartesian coordinates. The X and Y coordinates

representing soil and land-use types and the corresponding class’s median topsoil phos-

phorus content for groups with sufficient measured topsoil P data were used to interpo-

late trend surface in Surfer (ver. 24) by using radial basis multiquadratic (R2 parameter

0.2) gridding. The created ordination model surface represents the P value for different

soil and land-use combinations. This model was used to assign median phosphorus values

for all soil–land use polygons in the entire Estonia according to particular combination of

soil type and land-use class. The soil–land use ordination model results were validated by

calculating the linear relationship between the model results predicted for the soil–land

use polygons all over Estonia and all phosphorus sampling points in the PANDA data-

base. Additionally, the soil–land use polygons’ median P values calculated from the meas-

ured P data from the PANDA database were used for validation. As the P values in the

PANDA database are not normally distributed, a nonparametric correlation analysis was

implemented, and Kendall tau correlation metrics between the measured and predicted

data were calculated.

Bagging

To predict the topsoil phosphorus content for soil–land use/land cover polygons all

over Estonia based on multiple variables, the prediction model was built by using bag-

ging-based machine learning method (bootstrap aggregation and classification and re-

gression trees (CART) method) [56]. Bagging works best with high-variance unstable

models, and it is usually applied to decision tree methods, where a forest of many trees is

created to average the model and predict the best outcome, to add stability and accuracy,

and to reduce overfitting and variance. Bagging decreases the variance in prediction in

high-variance models without increasing the bias. It creates new bootstrap training sets

Agronomy 2023, 13, 1183 8 of 24

by sampling with replacement and then fits a model to each new training set. These mod-

els are combined by averaging the predictions for the regression case. Other machine

learning methods, such as random forest and extreme gradient boosting methods, were

tried, but they gave very similar results; thus, bagging was chosen because the bagging

prediction was faster in case of large dataset (over 3 million soil–land use polygons). Mi-

crosoft R and R Studio (v4.0.2) were used to build machine learning models and for ex-

ploratory data analysis. For bagging, the tidymodels and baguette libraries were used. As

the dependent variable (topsoil phosphorus content (mg/kg)) was a continuous variable,

a regression mode and rpart (recursive partitioning and regression tree) model [57] was

used for bagging. The dataset was randomly divided into training (75%, 268,882 samples)

and test data (25%, 89,625 samples), and high and low phosphorus values were selected

proportionally. A total of 25 bootstrap samples were taken from the training dataset to fit

the regression model. Additionally, 50 and 100 bags (bootstrap samples) were also tried,

but the model improvement was very limited compared to the increase in computing time.

For bagging, individual trees are grown deep and not pruned. These trees have high var-

iance and low bias. Overfitting and high variance are solved by averaging the model over

the number of trees. In this study, bagging fitted 25 regression models and averaged the

results to reduce the variance. A total of 32 variables (Supplementary Table S1) were used

to build the bagging model to predict topsoil P content. The categorical variables (soil

type, land-use/land-cover type, sediments (parent material), soil texture, and soil hydro-

logical group) were replaced with the median phosphorus content value of each category

(Supplementary Table S1).

2.4.3. Model Evaluation

To compare the results from different models, similar evaluation methods and metrics

were used when applicable. The main metrics we used to evaluate the model performance

were coefficient of determination (R2), root mean squared error (RMSE), and mean absolute

error (MAE). R2 is the metric used to describe the performance of any models (including

cross-validation and test set validation). R2 explains to what extent the variance of the pre-

dicted P explains the variance of the measured P (or polygon median P). In our case, R2 also

evaluates the linear regression between the measured and predicted values. RMSE is the

standard deviation of the residuals (prediction errors) in original units. As the effect of each

error (deviation from regression line) is proportional to the size of the squared error, RMSE

is sensitive to outliers. MAE is the average of the absolute values of the errors in original

units.

In the kriging method, prediction errors are evaluated using cross-validation by leav-

ing one sample point out and using all other data points to predict the value for this omit-

ted point, but cross-validation method uses all data points to estimate the trend and auto-

correlation. As the cross-validation differs from validation applied for evaluation of other

models, our original phosphorus dataset was divided into training (75%) and test (25%)

datasets by using the ArcGIS Geostatistical Analyst extension tool Subset Features, which

created a random subset from the dataset. Both training data cross-validation and test data

validation against the measured P values were used to evaluate the kriging model's per-

formance. As the ArcGIS Geostatistical extension kriging evaluation statistics does not in-

volve any other metrics, the cross-validation dataset was exported to Excel, and the eval-

uation metrics (R2, RMSE, and MAE) were calculated separately.

The bagging model was built based on 75% of data (268,882 samples) being used as

the training set. The other 25% of data (89,625 samples) were used for model validation.

To test the bagging model performance, R2, RMSE, and MAE were calculated.

As bagging uses random subsampling with replacement, each bagged predictor was

trained on approximately 63% of training data, and the remaining 37% of data were called

out-of-bag (OOB) observations. OOB error estimates require less computation than cross-

validation and can be used to evaluate model performance in the training process. An-

other benefit of OOB is that it does not require a separate test dataset. Cross-validation

Agronomy 2023, 13, 1183 9 of 24

(CV) and 10-fold random sampling were used to evaluate the model performance on the

training data, and the evaluation metrics R2, RMSE, and MAE, were calculated.

Besides using traditional (validation on test data) and model-specific (OOB error and

cross-validation) validation methods, all model results were also compared with the meas-

ured composite sample P values and the soil–land use polygon median P values. The

model prediction results were extracted at the measured P sample locations (358,947 sam-

ples), and evaluation metrics (MAE, RMSE, and R2) to compare the measured vs. predicted

values were calculated. The topsoil P samples collected in 2021 (29,261 samples) served as

an additional independent validation dataset not used in any other analysis.

The topsoil P samples were also compared to the soil–land use polygon median P

values as the ordination and bagging method predicted values to soil–land use polygons.

Likewise, the topsoil P samples collected in 2021 were also compared to the soil–land use

polygon median P values, which were calculated based on the topsoil P samples taken in

the period 2005–2020.

3. Results

3.1. Long-Term Changes of Topsoil Phosphorus in Agricultural Soils

Plant-available P is relatively stable in natural ecosystems but strongly influenced by

management in agricultural areas. To assess trends in topsoil P content, we used long-

term monitoring data from the PANDA database for agricultural land in the period 2005–

2021.

The lowest number of P samples was collected in 2012 (5243 samples), while 24,857

P samples were collected in 2015. The median P value was lowest in 2005 (47 mg/kg) and

highest in 2016 (82 mg/kg). The P sample locations, sample size, and median P content by

years are presented in the supplementary materials (Figure S1). As the exploratory analy-

sis detected no specific temporal trend for P values collected in different years (Figure 1),

the temporal component was left out from further topsoil P modeling. The P values vary

interannually mostly because different locations were sampled in different years and

resampling usually occurred in four years. Moreover, the topsoil P value is affected by

fertilizer application and crop rotation. Extremely high P values might be measured in

case of recent fertilizer application nearby and were treated as outliers.

Figure 1. Measured topsoil P values (mg/kg) for agricultural lands in 2005–2021 and corresponding

median P values for main land-use classes.

Agronomy 2023, 13, 1183 10 of 24

3.2. Dependence of Topsoil P on Soil Parameters in Agricultural Lands

To understand the dynamics and potential loss of P at the landscape level, it is im-

portant to understand how soil parameters are related to P content. We used exploratory

statistical analysis to check the measured topsoil P content relations with other variables

related to soil, land use, and topography available for the entire Estonia (Figure 2). The

correlation analysis shows that in agricultural ecosystems, plant-available P content in

topsoil is mainly controlled by fertilization and land management. The best explanatory

variables in the single-parameter linear regression to predict P values are soil type (R2 =

0.15) and other soil-related categorical variables, including texture (R2 = 0.08) and soil hy-

dro group (R2 = 0.09). Land-use category contributes less to P prediction in the single-

parameter model (R2 = 0.02). According to the linear regression model built by using soil

and land-use categories, both soil and land-use categories are statistically significant fac-

tors describing the measured P values (R2 = 0.16). The exploratory data analysis shows

that higher median P values are associated with automorphic soils (original soil classes

shown, approximate WRB reference is provided in [49], mainly with K (Kh 116 mg/kg, KI

93 mg/kg, Ko 81 mg/kg, and Kr 75 mg/kg) and L (Lk 126 mg/kg, LP 101 mg/kg, and LkG

77 mg/kg) soils and also D soils (69 mg/kg). Lower P values are common for hydromorphic

soils (AM 16, M 18, T 18, AG 24, G1 27, and Gk 29 mg/kg) (see also Figure 3).

Figure 2. Kendall correlation between topsoil P values (mg/kg) for agricultural lands in 2005–2020

and variables used in the bagging model.

Agronomy 2023, 13, 1183 11 of 24

Quaternary sediments as an indicator of parent material are directly related to soil

properties, and different sediments contribute differently to topsoil P values. The highest

median P values are associated with glaciofluvial deposits (MedP = 90 mg/kg). Shallow

quaternary deposits in the case of Rendzic soils (MedP = 76), aeolian deposits (MedP = 74),

and till (MedP = 74) have medium median P values. The lowest median P values belong

to wetland sediments (peat) where MedP = 35, lacustrine deposits (MedP = 44), glaciola-

custrine deposits (MedP = 51), and marine and alluvial deposits (MedP = 59).

3.3. Topsoil Phosphorous Modeling

3.3.1. Soil–Land Use Ordination Model

Soil–land use/land cover ordination is principally a categorical model that assigns the

same P values for a particular soil–land use category combination. It can distinguish values

for types in a matrix predefined by soil categories (29) multiplied by land-use categories

(11), which leads to a total of 298 soil–land use combinations. Each unique soil–land use

ordination combination corresponds to single predicted P value. The advantage of this

method is that if the soil and land use ordination is known, then further interpolation can

also provide topsoil P values for categories where observed values are missing or limited in

number, or estimate continuous values for transition from one soil/land use category to an-

other. An ordination model was created to predict the median topsoil phosphorus values

for pre-defined soil–land use combinations. The model predicts topsoil P values for pre-

defined soil–land use combinations based on the trend surface interpolated by using all

available P values (including data from the literature), and the soil and land-use categories

are ordinated in an ascending order (Figure 3).

Figure 3. Data gap filling by radial basis multiquadratic spatial interpolation along soil and crop

ordination axes to model topsoil P values for any soil-type and land-use combinations.

The results show that soil and land-use types have influence on topsoil phosphorus

content, but these factors do not describe the majority of variance, which depends on

many other factors. This kind of simple two-parameter ordination and classification

model does not cover enough variance in topsoil P content to use the predicted P values

for, e.g., P loss modeling. To improve the ordination model, the ordination values were

replaced with the soil–land use polygon median phosphorus values (Figure 4) calculated

Agronomy 2023, 13, 1183 12 of 24

based on the PANDA database (2005–2020) in those polygons where phosphorus sam-

pling points existed (there were 108,128 polygons out of 3,063,558 where at least one sam-

pling point was located). In this way, the model quality was improved in agricultural ar-

eas. Model validation against the measured P sample values (2005–2020) gave the follow-

ing results: R2 = 0.74, MAE = 18.8, and RMSE = 36.4, showing the relationship between

soil–land use polygon median topsoil phosphorus content and measured phosphorus

sampling points. The results suggest that there is quite high variance in topsoil phospho-

rus content even within the same agricultural parcel with an identical soil type.

Figure 4. Modeled topsoil median plant-available P content (mg/kg) based on soil–land use ordina-

tion + polygon median P model.

3.3.2. Kriging Model

Simple kriging was used to interpolate topsoil P content in Estonia (Figure 5) based

on the P sampling points (PANDA database) collected in 2005–2020. Two kriging models

were built in order to validate the kriging results by different methods. The first kriging

model (Kriging1) was built traditionally by using all sample points and validated by cross-

validation (by leaving-one-out method). The Kriging1 model cross-validation evaluation

gave the following results: R2 = 0.63, MAE = 27.1, and RMSE = 43.1.

Agronomy 2023, 13, 1183 13 of 24

Figure 5. Modeled topsoil median plant-available P content (mg/kg) based on kriging model. Non-

agricultural land uses are excluded from kriging to avoid extrapolation.

To test the kriging model performance on the test data, Kriging2 model was built by

using only 75% of randomly selected P sample data. The remaining 25% of P samples

(89,661 samples) were used for model validation, which gave the following results: R2 =

0.62, MAE = 27.7, and RMSE = 44.4.

The kriging models represent the topsoil P values well in fields which are uniform in

soil cover and are excessively sampled, but they fail in areas where soil cover is mosaic,

crop rotation has been different over the years, or the distance between existing sampling

points is greater than the average soil mapping unit diameter.

3.3.3. Bagging Model

The bagging model was introduced in order to overcome the disadvantages of the

kriging and ordination models, namely a low share of P variation explained by the two-

parameter soil–land use ordination model and the unreliability of the kriging model re-

sults outside of agricultural areas and far from the P sampling points.

Ensemble machine learning models can give much more accurate predictions than

simple models, but the result is like a Black Box Model which underlying functional form

is so complex that it cannot be written down as a simple formula. We used 32 variables

(Supplementary Table S1) to build the bagging model. Ensemble models, such as bagging

models, that average the results of multiple models are not easy to interpret as each vari-

able may appear in different trees in different positions, if at all. However, bagging gives

the aggregated variable importance scores (Figure 6), which can be used to understand

which variables are more important in predicting the topsoil P values. According to the

bagging model that we used to predict topsoil P values (Figure 7), the two most important

variables were the Y coordinate (latitude) and X coordinate (longitude). The next im-

portant variables were soil-type median phosphorus content (Soil_MedP) and the soil-

type median phosphorus content if the land use was natural grassland (NatSoil_MedP).

The absolute elevation was in the fifth rank of importance. The most important topo-

graphic variables calculated from the digital elevation model were the standard deviation

Agronomy 2023, 13, 1183 14 of 24

of slope in soil polygon (slp_stdev), mean terrain ruggedness index (tri_mean), and stand-

ard deviation of LS factor (ls_stdev). Soil organic carbon (soc1) was also important as ex-

pected, and so was soil polygon size (unit_area).

Figure 6. Standardized variable importance score of factors in the bagging model.

In many cases, the importance of X and Y coordinates is complicated to interpret, but in

our model, they represent a known set of natural phenomena following these directions. The

north–south direction represents the gradient in bedrock (from limestone to sandstone) and

the presence of naturally occurring phosphorite brought up by glacial activity in the northern-

most part of Estonia [58]. The west–east direction coincides with the gradients of maritime to

continental climate, lower to higher topography, and soil formation age after the retreat of the

last continental ice sheet.

Figure 7. Modeled topsoil median plant-available P content (mg/kg) based on the bagging model.

Agronomy 2023, 13, 1183 15 of 24

3.3.4. Hybrid Map of Topsoil Phosphorus

All models created based on different methods have limitations inherited from data

availability or bias by land-use or soil-type classes. Agricultural land is spatially well cov-

ered with sampling data, but only limited soil types and few sampling points are available

in natural ecosystems (e.g., forests and wetlands) or human-altered environment (e.g., quar-

ries), thus severely hindering use of data-driven methods such as machine learning or

kriging. Moreover, the relationship between the explanatory variables and predicted varia-

ble established in a managed environment, such as agricultural land, could not be directly

applied in natural ecosystems. On the other hand, simple land-use/soil-type relationship

with topsoil P in a relatively uniform natural environment does not describe the spatial and

temporal variability of managed ecosystems where the same soil polygon might be divided

by different land uses or experience various fertilization practices or crop rotation affecting

the soil phosphorus status.

To overcome the problems related to different methods, including (a) interpolation

performs well only in densely sampled agricultural lands, (b) the bagging model covers

with sufficient accuracy most land-use categories with meaningful sampling data by land-

use and soil classes, but topsoil P predictions overshoot in natural ecosystems with limited

data availability, and (c) the ordination model’s substantive problem that single variables,

such as soil and land use/land cover, explain only a minor share of variance related to

topsoil P content in managed and fertilized land-use classes, we created a hybrid map

(Figure 8).

Figure 8. Hybrid model of topsoil median plant-available P content (mg/kg) based on the bagging

model and soil–land use polygon median P values if available for agricultural land, and the ordina-

tion model for natural vegetation-covered areas.

The hybrid map was created by combining arable land polygons with the median

value of existing sampling data, the bagging-based predicted values for other agricultural

lands, and the two-parameter soil–land use ordination matrix-based model values for for-

ests, wetlands, and peat extraction areas. The topsoil P values for each soil–land use unit

of eight agricultural land-use categories (permanent grassland, fallow, cultivated grass-

land, permanent cultures, cereals, legumes, technical cultures, and vegetables) were based

Agronomy 2023, 13, 1183 16 of 24

on the observed polygon median topsoil P content value if a particular soil–land use unit

had sufficient number of sampling data to calculate the median value. The bagging model

predicted topsoil P values were used in agricultural land where sampling data were lim-

ited or missing. The remaining three non-agricultural land-use types (forest, wetland, and

peat extraction area) were based on the ordination model results.

The created high-resolution topsoil phosphorus map reveals clear a latitudinal–longitu-

dinal and elevation-related spatial pattern. Higher predicted P values are concentrated in the

northern and western coastal areas. In the northern coast region (Maardu, Toolse, and Rak-

vere), there are phosphorite deposits due to glacial deposits that naturally contribute to the

higher P values in the topsoil of this region. In West Estonian islands and in uplands, the P

content in the topsoil is higher than average, partly due to long-term management and fertili-

zation. In Estonian northern and western coastal areas, thin and calcareous juvenile soils con-

tribute to the higher P content as soil age is one of the factors contributing to the P content.

Low P values are associated with large wetland areas and older soils with lighter texture,

which contain less P due to leaching and erosion.

3.3.5. Validation of Topsoil P Prediction Models

All models were validated against the measured topsoil P values, the median of the

measured P values of the soil–land use unit area, as well against each other.

• Ordination validation

Due to the data-driven approach of the ordination model, classical model validation

methods (test dataset and cross-validation) were not available. All topsoil P sample points

were used to evaluate the ordination model performance (R2 = 0.15, MAE = 42.9, and RMSE

= 66.8) (Table 1, Ordination vs. P).

Table 1. Ordination model validation results.

Model MAE RMSE R2 Explanation

Ordination vs. P 42.9 66.8 0.15 Soil–land use ordination model results validated against

measured P samples (all P samples used for validation).

Ordination vs.

poly MedP 36.7 56.5 0.2

Soil–land use ordination model results validated against

polygon median P values (all P samples used for validation).

Ordination vs. P

2021 46.2 73.4 0.17

Soil–land use ordination model results validated against

P data collected in 2021 and not used for model building.

Ordination +

polygon MedP

vs. P 2021

34.5 57.4 0.48

Soil–land use ordination model values replaced with

soil–land use polygon median P values in polygons

where P sample points measured in 2005–2020 exist, and

evaluated against P samples collected in 2021.

The model results were validated with an independent dataset (P samples collected in

2021), and the results were as follows: R2 = 0.17, MAE = 46.2, and RMSE = 73.4 (in Table 1

Ordination vs. P 2021). The ordination model results were also evaluated against the soil–land

use polygon median P values as the ordination model predicts values based on the soil–land

use types. The ordination model evaluation against the polygon median P values gave slightly

better results (R2 = 0.20, MAE = 36.7, and RMSE = 56.5) (Table 1, Ordination vs. poly MedP)

compared to the evaluation against the sampled P values. The ordination + polygon median

P model results were validated with an independent dataset (P samples collected in 2021), and

the results were as follows: R2 = 0.48, MAE = 34.5, and RMSE = 57.4 (Table 1: Ordination +

polygon MedP vs. P 2021).

• Kriging validation

The Kriging1 model cross-validation evaluation gave the following results: R2 = 0.63,

MAE = 27.1, and RMSE = 43.1 (Table 2, Kriging1 CV vs. P). The R2 between the Kriging1 pre-

dicted values and all measured phosphorus values is 0.82, MAE = 19.8, and RMSE = 31.3 (Table

2, Kriging1 vs. P), but this comparison is not valid as the validation considering the kriging1

Agronomy 2023, 13, 1183 17 of 24

model was built by using the same P sample points, and the kriging interpolator is the exact

interpolator; thus, the result rather suggests that variation in the topsoil P values is quite high

even in a very short distance.

Table 2. Validation of kriging model results.

Model MAE RMSE R2 Explanation

Kriging1 CV vs.

P 27.1 43.1 0.63

Kriging1 (all sample points used for interpolation) cross-

validation (by leaving-one-out method).

Kriging1 vs. P 19.8 31.3 0.82

Kriging1 (all sample points used for interpolation) valida-

tion against all measured P samples (same points that were

used for interpolation).

Kriging1 vs.

poly MedP 15.6 25.1 0.84

Kriging1 (all sample points used for interpolation) validation

against polygon median P values (if there was only one P sam-

ple in the polygon, then the same values were used for inter-

polation).

Kriging1 vs. P

2021 30.7 49.6 0.59

Kriging1 (all sample points used for interpolation) validation

against P samples (collected in 2021) not used for model build-

ing.

Kriging2_test

vs. P 27.7 44.4 0.62

Kriging2 model (surface was interpolated by using only

75% of P sample data) validated against test data (25% of

data). Only 25% of data were used to calculate these met-

rics.

Kriging2_test

vs. poly MedP 19.7 31.6 0.75

Kriging2 model (surface was interpolated by using only

75% of P sample data) validated against test data (25% of

data) polygon median P values. Only 25% of data were

used to calculate these metrics.

The Kriging1 model results were also compared with the soil–land use polygon median

P values (kriging model value in the measured P sample location was compared with the

median P value of soil–land use polygon; Table 2, Kriging1 vs. poly MedP). These results

show very good model performance; however, as the model was built by using all P samples

and it was evaluated using the same data, the good agreement was expected. The Kriging1

model performance based on the data that were not included in model building (P samples

taken in 2021) gave the following results: R2 = 0.59, MAE = 30.7, and RMSE = 49.6 (Table 2,

Kriging1 vs. P 2021), indicating significant interannual and spatial topsoil P variation due to

fertilization and crop rotation.

The Kriging2 model results were evaluated by comparing the test dataset that was

not used in the kriging surface prediction to the soil–land use polygon median P values

(R2 = 0.75, MAE = 19.7, and RMSE = 31.6; Table 2, Kriging2_test vs. poly MedP). With rela-

tively good accordance of this comparison in proportion to the abovementioned compar-

ison between the Kriging1 model and 2021 sampling data, it reflects that interannual var-

iation in arable-land topsoil P content might be a more important factor than spatial vari-

ation.

• Bagging validation

Bagging model performance was validated on the test dataset (bagging model was

built on 75% of randomly selected bootstrap sample and 25% of data were left for model

validation), and the results were as follows: R2 = 0.54, MAE = 30.8, and RMSE = 48.6 (Table

3, Bagging_test vs. P). The evaluation of the bagging model results against all measured P

samples (both training and test data) show moderate prediction capability (R2 = 0.56, MAE

= 30.4, and RMSE = 47.3 (Table 3, Bagging vs. P)).

Agronomy 2023, 13, 1183 18 of 24

Table 3. Validation results of the bagging model.

Model MAE RMSE R2 Explanation

Bagging vs. P 30.4 47.3 0.56 Bagging prediction results validated against measured P val-

ues (all P samples used for validation).

Bagging vs.

poly MedP 20.9 33.6 0.72

Bagging prediction results validated against polygon median

P values (all P samples used for validation).

Bagging vs. P

2021 36.9 58.3 0.43

Bagging prediction results validated against P samples col-

lected in 2021 and not used for model building.

Bagging + pol-

ygon MedP

vs. P 2021

34.3 56.5 0.49

Bagging model values replaced with soil–land use polygon

median P values in polygons where P sample points measured

in 2005–2020 exist, and evaluated against P samples collected

in 2021.

Bagging_test

vs. P 30.8 48.6 0.54

Bagging model results validated on test dataset (25% of data

not used for model building).

Bagging 10-

fold CV 31.3 49.1 0.52

Bagging model performance evaluation by using a random 10-

fold cross-validation dataset (75% of data used for model build-

ing).

Bagging OOB

error 50.3

Bagging model performance evaluation by using out-of-bag

(OOB) samples. From the training dataset (75% of data) se-

lected randomly for each bagged predictor, approximately

63% of data were used for training, and the remaining 37%

were used as OOB samples.

Poly MedP vs.

P 18.8 36.4 0.74

Polygon median P values compared to P values measured in

2005–2020 (polygon median P value was calculated based on

the same sample).

Poly MedP vs.

P2021 32.6 54.6 0.53

Polygon median P values (calculated based on P sample val-

ues in 2005–2020) compared to P values measured in 2021.

The bagging predicted P values against the polygon median P values show good

agreement: R2 = 0.72, MAE = 20.9, and RMSE = 33.6 (Table 3, Bagging vs. poly MedP). The

PANDA database P samples collected in 2021 were used as an additional validation da-

taset to evaluate the bagging model performance. The result shows a weaker prediction

capability like in the case of kriging1 with similar validation: R2 = 0.43, MAE = 36.9, and

RMSE = 58.3 (Table 3, Bagging vs. P 2021). When the bagging prediction results were re-

placed with the soil–land use polygon median P values (in polygons where P samples

were taken in 2005–2020) and the results were compared to the P samples taken in 2021,

the results were only slightly better: R2 = 0.49, MAE = 34.3, and RMSE = 56.5 (Table 3,

Bagging + polygon MedP vs. P 2021).

The bagging model performance was also evaluated by calculating the out-of-bag

estimate, which used the training data that were left out from the random bootstrap sam-

pling with replacement (OOB sample) as a test dataset (RMSE = 50.3; Table 3, Bagging

OOB error). The cross-validation (CV) and 10-fold random sampling gave the following

results: RMSE = 49.1, R2 = 0.52, and MAE = 31.3 (Table 3, Bagging 10-fold CV).

4. Discussion

High-resolution modeling of topsoil plant-available P over large areas is a compli-

cated task as P content is influenced by both natural and human-driven processes; some

of these phenomena are continuous in time and space (e.g., soil development, parent ma-

terial, erosion, climate), while management of agricultural lands can affect P content in

topsoil sharply at the border of land parcels, even in cases of the same soil and crop types,

or over the years, bearing in mind legacy P as well.

Previous studies have shown that soil P content is heterogeneous even at the field

level [59], although land use is also an important factor determining P content. Our results

show that the topsoil P values inside the same soil–land use polygon can have quite high

variance (mean P range inside soil–land use polygon 58 mg/kg, median range 38 mg/kg,

Agronomy 2023, 13, 1183 19 of 24

maximum range of 198 mg/kg). Very high P values are usually related to croplands but

are sometimes observed in pasture areas as well, which can be explained by the presence

of livestock and continuous addition of manure. This phenomenon has been observed in

extensively grazed natural grasslands and wooded meadows in Estonia, and it is also rec-

orded by Roger et al. [13]. Roger et al. [13] showed that due to continuous addition of

manure (on average, 737 mg/kg organic P), permanent grasslands and mountain pasture

areas, compared to croplands, have the lowest quantity of total P, but the highest quantity

of available P.

There is an inherent problem of unbalanced monitoring data by land-use categories

but also by soil types, which affects spatial analysis of soil properties [60,61]. The relative

lack of information from non-agricultural areas results in bias of data by spatial extent and

land-use types. The relationship established between the explanatory variables and pre-

dicted variable in a managed environment could not be directly applied in natural ecosys-

tems. While topsoil P is relatively uniform by soil type in natural environments, the spatial

and temporal variability of managed ecosystems is very high, being affected by land uses,

fertilization practices, crop rotation, as well as the legacy P. Most authors, who have pre-

dicted P values for larger areas, have found that prediction uncertainty is high [23], and

higher uncertainty has been associated with sparse measurements and higher measured P

values [19].

Most commonly used ancillary environmental variables for P prediction models in-

clude soil properties, e.g., clay, sand, soil organic matter, pH, and soil type [13,17,20,62],

different topographic variables, such as slope, elevation, curvature, LS factor and topo-

graphic wetness index [13,17,20], vegetation and remote sensing data mostly in the form of

vegetation indices, e.g., NDVI [16,22], land use [13,17], and climatic variables (e.g., temper-

ature and precipitation indices) [17]. Similar to our results, higher P values are associated

with anthropogenic land use and fertilizer application [19]. According to Hosseini et al. [20]

attributes related to agricultural practices contributed more to plant-available P models,

whereas soil and geological attributes contributed more to TP (total P) models. Roger et al.

[13] found that environmental variables (elevation, slope, wetness index, and plan curvature

derived from a digital elevation model) explained about 20–25% of spatial variance of dif-

ferent P forms. In our study, soil type explained about 15% of P variance, land-use category

explained 13%, and other parameters contributed individually less, although most of them

were statistically significant in the bagging model.

The kriging model shows the best results in model evaluation, but its use is limited

to land-use categories with extensive sampling data available as interpolation methods

are primarily based on spatial autocorrelation. The bagging and ordination models are

not related to certain locations, although the bagging model includes also spatial variables

(latitude and longitude), and, therefore, the values that these models predict are not di-

rectly connected to the sample point values in particular locations.

5. Conclusions

Three methods were used for data-driven spatial modeling to produce a spatially

detailed map of plant-available topsoil P content in Estonia. Simple kriging performed the

best (R2 = 0.82, MAE = 19.8, and RMSE = 31.3) for agricultural lands where sampling den-

sity was high and spatially well distributed; however, due to the method’s reliance on

spatial autocorrelation, it is unable to predict values in areas with high spatial variance in

soil or land-use types, and in areas with missing sampling points or with low sampling

density where the active lag distance is exceeded.

The results of the bagging model were slightly inferior (R2 = 0.56, MAE = 30.4, and

RMSE = 47.3) to the kriging model. However, bagging models are not directly related to

certain sampling locations nor sampling density but incorporate supplementary spatial

variables (e.g., latitude, longitude, elevation land-use classes, and soil properties); there-

fore, these models are more capable of predicting available P values in agricultural lands

distant to any sampling point. Nevertheless, neither bagging nor kriging models, which

Agronomy 2023, 13, 1183 20 of 24

heavily rely on an extensive database of measured topsoil P content from agricultural

soils, could be used to predict topsoil P values outside of agricultural lands, where the

processes of P mineralization and immobilization, often biological in nature, are generally

different than in agricultural ecosystems where P is added as a fertilizer and crop rotation

affects the topsoil P status.

In large non-agricultural areas (wetlands and forests, covering 61% of Estonian land

area), topsoil plant-available P content was infrequently sampled with poor spatial cover-

age. To overcome the spatial soil P data scarcity, a two-parameter soil–land use ordination

matrix-based model was developed. The median topsoil P value of this ordination ap-

proach did not explain the spatial variance (R2 = 0.20, MAE = 36.7, and RMSE = 56.5) in

data-abundant agricultural areas; however, soil type alone explained about 15% of P var-

iance and land-use class contributed less (13%), thus making it sufficient for gap filling in

these spatially more homogeneous areas, especially for peat soils.

The final hybrid topsoil plant-available P map (R2 = 0.74, MAE = 18.8, and RMSE = 36.4)

was produced by combining arable-land polygons with the median value of the sampling

data, the bagging-based predicted values for agricultural lands with missing sampling data,

and the two-parameter soil–land use ordination matrix-based model values for forests and

wetlands.

The predicted topsoil plant-available P map can be used to assess soil P fertilization

need, estimate legacy P, and optimize P fertilizer application, taking into consideration

soil properties to reduce P loss and eutrophication of waterbodies.

Supplementary Materials: The following supporting information can be downloaded at

https://www.mdpi.com/article/10.3390/agronomy13051183/s1. Figure S1: Topsoil plant-available P

sample points taken from agricultural lands in 2005–2020. The count (n) and median P value (mg/kg)

are presented for every year; Table S1: Variables used in bagging model. P is the predicted variable,

and 32 variables were used to predict the P values. Index 1 after sand, silt, clay, rock, soc, bd, k, and

awc indicates the upper soil layer (topsoil). Different statistics (mean, median, and standard devia-

tion) of the topographic variables (TRI, TWI, LS, and SLP) were included in the bagging model.

References [63–68] are cited in Supplementary Materials.

Author Contributions: Conceptualization, A.K. (Ain Kull), P.P., and A.K. (Anne Kull); methodology,

A.K. (Ain Kull) and A.K. (Anne Kull); formal analysis, A.K. (Anne Kull), A.K. (Ain Kull), and T.K.; re-

sources, A.K. (Ain Kull) and P.P.; data curation, A.K. (Anne Kull), A.K. (Ain Kull), and T.K.; writing—

original draft preparation, A.K. (Anne Kull) and A.K. (Ain Kull); writing—review and editing, A.K. (Ain

Kull), A.K. (Anne Kull), T.K., and P.P.; visualization, A.K. (Anne Kull) and A.K. (Ain Kull); funding ac-

quisition, A.K. (Ain Kull) and P.P. All authors have read and agreed to the published version of the man-

uscript.

Funding: This research was supported by the Estonian Research Council research grants PRG352,

PRG1167, and SLTOM19384 (WaterJPI-JC-2018_13); the Estonian Environmental Investment Centre

grants SLOOM12006 and SLOOM14103; the Estonian State Forest Management Centre grant

LLTOM17250; and the Ministry of Rural Affairs.

Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.

Data Availability Statement: The high-resolution plant-available topsoil phosphorous hybrid map

generated during the study is available at the following site: https://datadoi.ee/handle/33/405?lo-

cale-attribute=en (accessed on 13 January 2023).

Acknowledgments: The authors would like to thank Elsa Putku for her valuable comments and

discussions during manuscript preparation process.

Conflicts of Interest: The authors declare no conflicts of interest. The funders had no role in the

design of the study; in the collection, analyses, or interpretation of data; in the writing of the manu-

script, or in the decision to publish the results.

Agronomy 2023, 13, 1183 21 of 24

References

1. Harrison, A.F. Soil Organic Phosphorous: A Review of World Literature; CAB international: Singapore, 1987.

2. Bardgett, R. The Biology of Soil: A Community and Ecosystem Approach; Oxford University Press: New York, NY, USA, 2005.

3. van der Wal, A.; de Boer, W.; Lubbers, I.M.; van Veen, J.A. Concentration and vertical distribution of total soil phosphorus in

relation to time of abandonment of arable fields. Nutr. Cycl. Agroecosystems 2007, 79, 73–79. https://doi.org/10.1007/s10705-007-

9097-3. Available online: https://agris.fao.org/agris-search/search.do?recordID=US201300792922 (accessed on 13 January 2023).

4. Batjes, N.H. Global distribution of soil phosphorus retention potential. ISRIC World Soil Inf. 2011, 6.

5. Penn, C.; Camberato, J. A Critical Review on Soil Chemical Processes that Control How Soil pH Affects Phosphorus Availability

to Plants. Agriculture 2019, 9, 120. https://doi.org/10.3390/agriculture9060120. Available online: https://www.mdpi.com/2077-

0472/9/6/120 (accessed on 13 January 2023).

6. Brady, N.C.; Weil, R.R. Soil organic matter. The nature and properties of soils. Prentice Hall: Upper Saddle River, NJ, USA, 1999.

7. Tarafdar, J.C.; Claassen, N. Organic phosphorus compounds as a phosphorus source for higher plants through the activity of

phosphatases produced by plant roots and microorganisms. Biol Fertil. Soils 1988, 5, 308–312. https://doi.org/10.1007/BF00262137.

Available online: https://www.researchgate.net/profile/J-Tarafdar/publication/226638701_Organic_Phosphorus_Com-

pounds_as_a_Phosphorus_Source_for_Higher_Plants_through_the_Activity_of_Phosphatases_Pro-

duced_by_Plant_Roots_and_Microorganisms/links/5816002208ae90acb240e93f/Organic-Phosphorus-Compounds-as-a-Phos-

phorus-Source-for-Higher-Plants-through-the-Activity-of-Phosphatases-Produced-by-Plant-Roots-and-Microorganisms.pdf

(accessed on 13 January 2023).

8. Beckett, P.H.T.; White, R.E. Studies ON THE PHOSPHATE POTENTIALS OF SOILS: PART III: THE POOL OF LABILE INOR-

GANIC PHOSPHATE. Plant Soil 1964, 21, 253–282. https://doi.org/10.1007/BF01377744. Available online:

https://www.jstor.org/stable/42932073 (accessed on 13 January 2023).

9. Hinsinger, P. Bioavailability of soil inorganic P in the rhizosphere as affected by root-induced chemical changes: A review. Plant

Soil 2001, 237, 173–195. https://doi.org/10.1023/A:1013351617532. Available online: https://www.researchgate.net/pro-

file/Philippe-Hinsinger/publication/225852665_Bioavailability_of_soil_inorganic_P_in_the_rhizosphere_as_affected_by_root-

induced_chemical_changes_A_review/links/53d751c60cf228d363eaebbd/Bioavailability-of-soil-inorganic-P-in-the-rhizo-

sphere-as-affected-by-root-induced-chemical-changes-A-review.pdf (accessed on 13 January 2023).

10. Lambers, H.; Raven, J.A.; Shaver, G.R.; Smith, S.E. Plant nutrient-acquisition strategies change with soil age. Trends Ecol. Evol.

2008, 23, 95–103. https://doi.org/10.1016/j.tree.2007.10.008. Available online: https://www.sciencedirect.com/science/arti-

cle/pii/S0169534707003576 (accessed on 13 January 2023).

11. Nelson, L.; Cade-Menun, B.J.; Walker, I.J.; Sanborn, P. Soil Phosphorus Dynamics Across a Holocene Chronosequence of Aeo-

lian Sand Dunes in a Hypermaritime Environment on Calvert Island, BC, Canada. Front. For. Glob. Chang. 2020, 3, 83.

https://doi.org/10.3389/ffgc.2020.00083. Available online: https://doaj.org/article/2fbe4a8913d14de3a470466da86e84c5 (accessed

on 13 January 2023).

12. Tian, J.; Dong, G.; Karthikeyan, R.; Li, L.; Harmel, R. Phosphorus Dynamics in Long-Term Flooded, Drained, and Reflooded

Soils. Water 2017, 9, 531. https://doi.org/10.3390/w9070531.

13. Roger, A.; Libohova, Z.; Rossier, N.; Joost, S.; Maltas, A.; Frossard, E.; Sinaj, S. Spatial variability of soil phosphorus in the

Fribourg canton, Switzerland. Geoderma 2014, 217, 26–36. https://doi.org/10.1016/j.geoderma.2013.11.001. Available online:

https://dx.doi.org/10.1016/j.geoderma.2013.11.001 (accessed on 13 January 2023).

14. Mondal, B.P.; Sekhon, B.S.; Sadhukhan, R.; Singh, R.K.; Hasanain, M.; Mridha, N.; Das, B.; Dhyay, A.C.; Banerjee, K. Spatial

variability assessment of soil available phosphorus using geostatistical approach. Indian J. Agric. Sci. 2020, 90, 1170–1175.

https://doi.org/10.56093/ijas.v90i6.104795. Available online: https://www.researchgate.net/publication/344326265_Spatial_vari-

ability_assessment_of_soil_available_phosphorus_using_geostatistical_approach (accessed on 13 January 2023).

15. Scull, P.; Franklin, J.; Chadwick, O.A.; McArthur, D. Predictive soil mapping: A review. Prog. Phys. Geogr. 2003, 27, 171–197.

https://doi.org/10.1191/0309133303pp366ra. Available online: https://jour-

nals.sagepub.com/doi/full/10.1191/0309133303pp366ra (accessed on 13 January 2023).

16. Rivero, R.G.; Grunwald, S.; Bruland, G.L. Incorporation of spectral data into multivariate geostatistical models to map soil

phosphorus variability in a Florida wetland. Geoderma 2007, 140, 428–443. https://doi.org/10.1016/j.geoderma.2007.04.026. Avail-

able online: https://www.sciencedirect.com/science/article/pii/S001670610700122X (accessed on 13 January 2023).

17. Ballabio, C.; Lugato, E.; Fernández-Ugalde, O.; Orgiazzi, A.; Jones, A.; Borrelli, P.; Montanarella, L.; Panagos, P. Mapping LU-

CAS topsoil chemical properties at European scale using Gaussian process regression. Geoderma 2019, 355, 113912.

https://doi.org/10.1016/j.geoderma.2019.113912. Available online: https://www.sciencedirect.com/science/arti-

cle/pii/S0016706119304768 (accessed on 13 January 2023).

18. Matos-Moreira, M.; Lemercier, B.; Dupas, R.; Michot, D.; Viaud, V.; Akkal-Corfini, N.; Louis, B.; Gascuel-Odoux, C. High-reso-

lution mapping of soil phosphorus concentration in agricultural landscapes with readily available or detailed survey data. Eur.

J. Soil Sci. 2017, 68, 281–294. https://doi.org/10.1111/ejss.12420. Available online: https://onlineli-

brary.wiley.com/doi/abs/10.1111/ejss.12420 (accessed on 13 January 2023).

19. Viscarra Rossel, R.A.; Bui, E.N. A new detailed map of total phosphorus stocks in Australian soil. Sci. Total Environ. 2016, 542,

1040–1049. https://doi.org/10.1016/j.scitotenv.2015.09.119. Available online: https://dx.doi.org/10.1016/j.scitotenv.2015.09.119

(accessed on 13 January 2023).

Agronomy 2023, 13, 1183 22 of 24

20. Hosseini, M.; Rajabi Agereh, S.; Khaledian, Y.; Jafarzadeh Zoghalchali, H.; Brevik, E.C.; Movahedi Naeini, S.A.R. Comparison

of multiple statistical techniques to predict soil phosphorus. Appl. Soil Ecol. 2017, 114, 123–131. https://doi.org/10.1016/j.ap-

soil.2017.02.011. Available online: https://www.sciencedirect.com/science/article/pii/S0929139316304085 (accessed on 13 Janu-

ary 2023).

21. Esfandiarpour-Boroujeni, I.; Shahini-Shamsabadi, M.; Shirani, H.; Mosleh, Z.; Bagheri-Bodaghabadi, M.; Salehi, M.H. Assess-

ment of different digital soil mapping methods for prediction of soil classes in the Shahrekord plain, Central Iran. Catena 2020,

193, 104648. https://doi.org/10.1016/j.catena.2020.104648. Available online: https://www.sciencedirect.com/science/arti-

cle/pii/S0341816220301983 (accessed on 13 January 2023).

22. Bogrekci, I.; Lee, W.S. Spectral Phosphorus Mapping using Diffuse Reflectance of Soils and Grass. Biosyst. Eng. 2005, 91, 305–

312. https://doi.org/10.1016/j.biosystemseng.2005.04.015. Available online: https://www.sciencedirect.com/science/arti-

cle/pii/S1537511005000814 (accessed on 13 January 2023).

23. Yang, X.; Post, W.M.; Thornton, P.E.; Jain, A. The distribution of soil phosphorus for global biogeochemical modeling. Biogeo-

sciences 2013, 10, 2525–2537. https://doi.org/10.5194/bg-10-2525-2013. Available online:

https://search.proquest.com/docview/1349834631 (accessed on 13 January 2023).

24. Iatrou, M.; Papadopoulos, A.; Papadopoulos, F.; Dichala, O.; Psoma, P.; Bountla, A. Determination of Soil Available Phosphorus

using the Olsen and Mehlich 3 Methods for Greek Soils Having Variable Amounts of Calcium Carbonate. Commun. Soil Sci.

Plant Anal. 2014, 45, 2207–2214. https://doi.org/10.1080/00103624.2014.911304. Available online: https://agris.fao.org/agris-

search/search.do?recordID=US201500196377 (accessed on 13 January 2023).

25. Szara, E.; Sosulski, T.; Szymańska, M.; Szyszkowska, K. Usefulness of Mehlich-3 test in the monitoring of phosphorus dispersion

from Polish arable soils. Environ. Monit. Assess. 2018, 190, 1–10. https://doi.org/10.1007/s10661-018-6685-4. Available online:

https://link.springer.com/article/10.1007/s10661-018-6685-4 (accessed on 13 January 2023).

26. Vares, A.; Jõgiste, K.; Kull, E. Early growth of some deciduous tree species on abandoned agricultural lands in Estonia. Balt. For.

2001, 7, 52–58.

27. Laas, A.; Kull, A. Application Of GIS For Soil Erosion And Nutrient Loss Modelling In A Small River Catchment. WIT Trans.

Ecol. Environ. 2003, 67, 525-534. https://doi.org/10.2495/SPD030491. Available online:

https://search.proquest.com/docview/2261516459 (accessed on 13 January 2023).

28. Vares, A.; Uri, V.; Tullus, H.; Kanal, A. Height growth of four fast-growing deciduous tree species on former agricultural lands

in Estonia. Balt. For. 2003, 9, 2–8.

29. Kull, A.; Kull, A.; Jaagus, J.; Kuusemets, V.; Mander, U. The effects of fluctuating climatic conditions and weather events on

nutrient dynamics in a narrow mosaic riparian peatland. Boreal Environ. Res. 2008, 13, 243–263. Available online:

https://helda.helsinki.fi/bitstream/handle/10138/234736/ber13-3-243.pdf?sequence=1 (accessed on 13 January 2023).

30. Kund, M.; Vares, A.; Sims, A.; Tullus, H.; Uri, V. Early growth and development of silver birch (Betula pendula Roth.) planta-

tions on abandoned agricultural land. Eur. J. Forest Res. 2010, 129, 679–688. https://doi.org/10.1007/s10342-010-0369-0. Available

online: https://agris.fao.org/agris-search/search.do?recordID=US201301854605 (accessed on 13 January 2023).

31. Varik, M.; Aosaar, J.; Ostonen, I.; Lõhmus, K.; Uri, V. Carbon and nitrogen accumulation in belowground tree biomass in a

chronosequence of silver birch stands. For. Ecol. Manag. 2013, 302, 62–70. https://doi.org/10.1016/j.foreco.2013.03.033. Available

online: https://dx.doi.org/10.1016/j.foreco.2013.03.033 (accessed on 13 January 2023).

32. Lutter, R.; Tullus, T.; Tullus, A.; Kanal, A.; Tullus, H. Forest Ecosystem Recovery in 15-Year-Old Hybrid Aspen (Populus Trem-

ula l. × p. Tremuloides Michx.) Plantations on a Reclaimed Oil Shale Quarry. Oil Shale (Tallinn Est. 1984) 2017, 34, 368–389.

https://doi.org/10.3176/oil.2017.4.05. Available online: https://search.proquest.com/docview/1973739673 (accessed on 13 Janu-

ary 2023).

33. Lutter, R.; Tullus, A.; Kanal, A.; Tullus, T.; Vares, A.; Tullus, H. Growth development and plant–soil relations in midterm silver

birch (Betula pendula Roth) plantations on previous agricultural lands in hemiboreal Estonia. Eur. J. Forest Res. 2015, 134, 653–

667. https://doi.org/10.1007/s10342-015-0879-x. Available online: https://link.springer.com/article/10.1007/s10342-015-0879-x

(accessed on 13 January 2023).

34. Uri, V.; Kukumägi, M.; Aosaar, J.; Varik, M.; Becker, H.; Soosaar, K.; Morozov, G.; Ligi, K.; Padari, A.; Ostonen, I.; Karoles, K.

Carbon budgets in fertile grey alder (Alnus incana (L.) Moench.) stands of different ages. For. Ecol. Manag. 2017, 396, 55–67.

https://doi.org/10.1016/j.foreco.2017.04.004. Available online: https://dx.doi.org/10.1016/j.foreco.2017.04.004 (accessed on 13 Jan-

uary 2023).

35. Uri, V.; Kukumägi, M.; Aosaar, J.; Varik, M.; Becker, H.; Morozov, G.; Karoles, K. Ecosystems carbon budgets of differently aged

downy birch stands growing on well-drained peatlands. For. Ecol. Manag. 2017, 399, 82–93.

https://doi.org/10.1016/j.foreco.2017.05.023. Available online: https://dx.doi.org/10.1016/j.foreco.2017.05.023 (accessed on 13 Jan-

uary 2023).

36. Uri, V.; Aosaar, J.; Varik, M.; Becker, H.; Kukumägi, M.; Ligi, K.; Pärn, L.; Kanal, A. Biomass resource and environmental effects

of Norway spruce (Picea abies) stump harvesting: An Estonian case study. For. Ecol. Manag. 2015, 335, 207–215.

https://doi.org/10.1016/j.foreco.2014.10.003. Available online: https://dx.doi.org/10.1016/j.foreco.2014.10.003 (accessed on 13 Jan-

uary 2023).

Agronomy 2023, 13, 1183 23 of 24

37. Uri, V.; Kukumägi, M.; Aosaar, J.; Varik, M.; Becker, H.; Aun, K.; Krasnova, A.; Morozov, G.; Ostonen, I.; Mander, Ü .; Lõhmus,

K.; Rosenvald, K.; Kriiska, K.; Soosaar, K. The carbon balance of a six-year-old Scots pine (Pinus sylvestris L.) ecosystem esti-

mated by different methods. For. Ecol. Manag. 2019, 433, 248–262. https://doi.org/10.1016/j.foreco.2018.11.012. Available online:

https://dx.doi.org/10.1016/j.foreco.2018.11.012 (accessed on 13 January 2023).

38. Varik, M.; Kukumägi, M.; Aosaar, J.; Becker, H.; Ostonen, I.; Lõhmus, K.; Uri, V. Carbon budgets in fertile silver birch (Betula

pendula Roth) chronosequence stands. Ecol. Eng. 2015, 77, 284–296. https://doi.org/10.1016/j.ecoleng.2015.01.041. Available

online: https://dx.doi.org/10.1016/j.ecoleng.2015.01.041 (accessed on 13 January 2023).

39. Aosaar, J.; Mander, Ü .; Varik, M.; Becker, H.; Morozov, G.; Maddison, M.; Uri, V. Biomass production and nitrogen balance of

naturally afforested silver birch (Betula pendula Roth.) stand in Estonia. Silva Fenn. 2016, 50. https://doi.org/10.14214/sf.1628.

Available online: https://www.silvafennica.fi/article/1628 (accessed on 13 January 2023).

40. Aosaar, J.; Drenkhan, T.; Adamson, K.; Aun, K.; Becker, H.; Buht, M.; Drenkhan, R.; Fjodorov, M.; Jürimaa, K.; Morozov, G.;

Pihlak, L.; Piiskop, K.; Riit, T.; Varik, M.; Väär, R.; Uri, M.; Uri, V. The effect of stump harvesting on tree growth and the infection

of root rot in young Norway spruce stands in hemiboreal Estonia. For. Ecol. Manag. 2020, 475, 118425.

https://doi.org/10.1016/j.foreco.2020.118425. Available online: https://dx.doi.org/10.1016/j.foreco.2020.118425 (accessed on 13

January 2023).

41. Aun, K.; Kukumägi, M.; Varik, M.; Becker, H.; Aosaar, J.; Uri, M.; Morozov, G.; Buht, M.; Uri, V. Short-term effect of thinning

on the carbon budget of young and middle-aged Scots pine (Pinus sylvestris L.) stands. For. Ecol. Manag. 2021, 492, 119241.

https://doi.org/10.1016/j.foreco.2021.119241. Available online: https://dx.doi.org/10.1016/j.foreco.2021.119241 (accessed on 13

January 2023).

42. Becker, H.; Aosaar, J.; Varik, M.; Morozov, G.; Kanal, A.; Uri, V. The effect of Norway spruce stump harvesting on net nitrogen

mineralization and nutrient leaching. For. Ecol. Manag. 2016, 377, 150–160. https://doi.org/10.1016/j.foreco.2016.07.005. Available

online: https://dx.doi.org/10.1016/j.foreco.2016.07.005 (accessed on 13 January 2023).

43. Paal, J.; Jürjendal, I.; Suija, A.; Kull, A. Impact of drainage on vegetation of transitional mires in Estonia. Mires Peat 2016, 18, 1–

19. https://doi.org/10.19189/MaP.2015.OMB.183. Available online: http://mires-and-peat.net/media/map18/map_18_02.pdf (ac-

cessed on 13 January 2023).

44. Tullus, T.; Lutter, R.; Randlane, T.; Saag, A.; Tullus, A.; Oja, E.; Degtjarenko, P.; Pärtel, M.; Tullus, H. The effect of stand age on

biodiversity in a 130-year chronosequence of Populus tremula stands. For. Ecol. Manag. 2022, 504, 119833.

https://doi.org/10.1016/j.foreco.2021.119833. Available online: https://dx.doi.org/10.1016/j.foreco.2021.119833 (accessed on 13

January 2023).

45. Kull, A. Buffer zones to limit and mitigate harmful effects of long-term anthropogenic influence to maintain ecological func-

tionality of bogs, stage II. 2016, 183. (In Estonian). Available online: https://4ce0b57b-a630-4e1d-8a88-

d1e1a7a51b96.filesusr.com/ugd/6b6658_446958f4118b44a2a68812820c31119b.pdf (accessed on 13 January 2023).

46. Asi, E.; Timmusk, T. Greenhouse gas emissions inventory studies for the national reporting on forest litter and soil. 2018, 69. (In

Estonian). Available online: file:///C:/Users/Administrator/Downloads/kik_metsanduse_projekt_12654_loppa-

ruanne_koos_lisaga_0-3.pdf (accessed on 13 January 2023).

47. Kulhánek, M.; Balík, J.; Černý, J.; Vaněk, V. Evaluation of phosphorus mobility in soil using different extraction methods. Plant

Soil Environ. 2009, 55, 267–272. https://doi.org/10.17221/43/2009-PSE. Available online: https://www.agriculturejour-

nals.cz/pdfs/pse/2009/07/01.pdf (accessed on 13 January 2023).

48. Wolf, A.M.; Baker, D.E. Comparisons of soil test phosphorus by Olsen, Bray P1, Mehlich I and Mehlich III methods. Commun.

Soil Sci. Plant Anal. 1985, 16, 467–484. https://doi.org/10.1080/00103628509367620. Available online: https://agris.fao.org/agris-

search/search.do?recordID=US8614205 (accessed on 13 January 2023).

49. Kmoch, A.; Kanal, A.; Astover, A.; Kull, A.; Virro, H.; Helm, A.; Pärtel, M.; Ostonen, I.; Uuemaa, E. EstSoil-EH: A high-resolution

eco-hydrological modelling parameters dataset for Estonia. Earth Syst. Sci. Data 2021, 13, 83–97. https://doi.org/10.5194/essd-13-

83-2021. Available online: https://essd.copernicus.org/articles/13/83/2021/ (accessed on 13 January 2023).

50. Kachinsky, N. Soil Physics. 1965; 320 p. (In Russian)

51. Abd-Elsamad, A.; Abdelmoein, N.M.; Mahmoud, A.H.; Rostom, M.; Hassan, S.M.; Gazni, R. Evaluation and comparison of

ordinary kriging and inverse distance weighting methods for prediction of spatial variability of some soil chemical parameters.

Res. J. Biol. Sci. 2009, 4, 93–102. Available online: https://www.researchgate.net/publication/266732920 (accessed on 13 January

2023).

52. Sahu, B.; Ghosh, A.K. Seema Deterministic and geostatistical models for predicting soil organic carbon in a 60 ha farm on In-

ceptisol in Varanasi, India. Geoderma Reg. 2021, 26, e00413. https://doi.org/10.1016/j.geodrs.2021.e00413. Available online:

https://www.sciencedirect.com/science/article/pii/S2352009421000584 (accessed on 13 January 2023).

53. Xue, W.; Pangara, C.; Aung, A.M.; Yu, S.; Tabucanon, A.S.; Hong, B.; Kurniawan, T.A. Spatial changes of nutrients and metallic

contaminants in topsoil with multi-geostatistical approaches in a large-size watershed. Sci. Total Environ. 2022, 824, 153888.

https://doi.org/10.1016/j.scitotenv.2022.153888. Available online: https://dx.doi.org/10.1016/j.scitotenv.2022.153888 (accessed on

13 January 2023).

54. Peet, R.K. Ordination as a tool for analyzing complex data sets. Vegetatio 1980, 42, 171–174. https://doi.org/10.1007/BF00048883.

Available online: https://citeseerx.ist.psu.edu/docu-

ment?repid=rep1&type=pdf&doi=182c33c89a9878641e1de783a1bffb26637f1326 (accessed on 13 January 2023).

Agronomy 2023, 13, 1183 24 of 24

55. Ter Braak, C.J.F.; Prentice, I.C. A Theory of Gradient Analysis. In Advances in Ecological Research; Elsevier Science & Technology:

Amsterdam, The Netherlands, 2004; Volume 34, pp. 235–282.

56. Breiman, L. Bagging predictors. Mach Learn. 1996, 24, 123–140. https://doi.org/10.1023/A:1018054314350.

57. Therneau, T.M.; Atkinson, E.J. An Introduction to Recursive Partitioning Using the RPART Routines. Introd. Recursive Partit.

Using RPART Routines 2022, 61. Available online: https://stat.ethz.ch/R-manual/R-patched/library/rpart/doc/longintro.pdf (ac-

cessed on 13 January 2023).

58. Petersell, V.; Ressar, H.; Carlsson, M.; Mõttus, V.; Enel, M.; Mardla, A.; Täht, K. Geochemical atlas of Estonian agricultural soil.

Eest. Geol. Sver. Geol. Undersökning. Tallinn Upps. Seletuskiri 1997, 75, 171 p. Available online: https://www.di-

gar.ee/viewer/et/nlib-digar:331323/291375/page/109. (In Estonian)

59. Werner, F.; Mueller, C.W.; Thieme, J.; Gianoncelli, A.; Rivard, C.; Höschen, C.; Prietzel, J. Micro-scale heterogeneity of soil

phosphorus depends on soil substrate and depth. Sci. Rep. 2017, 7, 3203–3209. https://doi.org/10.1038/s41598-017-03537-8. Avail-

able online: https://www.ncbi.nlm.nih.gov/pubmed/28600571 (accessed on 13 January 2023).

60. van Leeuwen, J.P.; Saby, N.P.A.; Jones, A.; Louwagie, G.; Micheli, E.; Rutgers, M.; Schulte, R.P.O.; Spiegel, H.; Toth, G.; Creamer,

R.E. Gap assessment in current soil monitoring networks across Europe for measuring soil functions. ERL 2017, 12, 124007.

https://doi.org/10.1088/1748-9326/aa9c5c. Available online: https://iopscience.iop.org/article/10.1088/1748-9326/aa9c5c (ac-

cessed on 13 January 2023).

61. Schillaci, C.; Saia, S.; Lipani, A.; Perego, A.; Zaccone, C.; Acutis, M. Validating the regional estimates of changes in soil organic

carbon by using the data from paired-sites: The case study of Mediterranean arable lands. Carbon Balance Manag. 2021, 16, 19.

https://doi.org/10.1186/s13021-021-00182-7. Available online: https://link.springer.com/article/10.1186/s13021-021-00182-7 (ac-

cessed on 13 January 2023).

62. Keshavarzi, A.; Omran, E.E.; Bateni, S.M.; Pradhan, B.; Vasu, D.; Bagherzadeh, A. Modeling of available soil phosphorus (ASP)

using multi-objective group method of data handling. Model. Earth Syst. Environ. 2016, 2, 1–9. https://doi.org/10.1007/s40808-

016-0216-5. Available online: https://link.springer.com/article/10.1007/s40808-016-0216-5 (accessed on 13 January 2023).

63. Zhang, Y.; Schaap, M.G. Weighted recalibration of the Rosetta pedotransfer model with improved estimates of hydraulic pa-

rameter distributions and summary statistics (Rosetta3). J. Hydrol. 2017, 547, 39–53. https://doi.org/10.1016/j.jhydrol.2017.01.00.

Available online: https://www.sciencedirect.com/science/article/pii/S0022169417300057?casa_to-

ken=Ej7anM0xjK4AAAAA:5r0KM6e6cHxb1FI60ASBz5XveBbA_r-ATIJxgGDs8kXTvUxK1kTbfAf5l9sntbmn6d_KodQSTg (ac-

cessed on 13 January 2023).

64. Dipak, S.; Abhijit, H. Physical and chemical methods in soil analysis. New Age Int. Ltd. New Delhi 2005, 193 p.

65. Tóth, B.; Weynants, M.; Pásztor, L.; Hengl, T. 3D soil hydraulic database of Europe at 250 m resolution. Hydrol Process 2017, 31,

2662–2666. https://doi.org/10.1002/hyp.11203. Available online: https://onlinelibrary.wiley.com/doi/10.1002/hyp.11203 (ac-

cessed on 13 January 2023).

66. Riley, S.J.; DeGloria, S.D.; Elliot, R. Index that quantifies topographic heterogeneity. Intermt. J. Sci. 1999, 5, 23–27. Available

online: http://download.osgeo.org/qgis/doc/reference-docs/Terrain_Ruggedness_Index.pdf (accessed on 13 January 2023).

67. Beven, K.J.; Kirkby, M.J. A physically based, variable contributing area model of basin hydrology/Un modèle à base physique

de zone d’appel variable de l’hydrologie du bassin versant. Hydrol. Sci. J. 1979, 24, 43–69.

https://doi.org/10.1080/02626667909491834. Available online:

https://www.tandfonline.com/doi/abs/10.1080/02626667909491834 (accessed on 13 January 2023).

68. Moore, I.D.; Grayson, R.B.; Ladson, A.R. Digital terrain modelling: A review of hydrological, geomorphological, and biological

applications. Hydrol Process 1991, 5, 3–30. https://doi.org/10.1002/hyp.3360050103. Available online: https://onlineli-

brary.wiley.com/doi/10.1002/hyp.3360050103 (accessed on 13 January 2023).

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual au-

thor(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to

people or property resulting from any ideas, methods, instructions or products referred to in the content.