A globally robust relationship between water table decline, subsidence rate, and carbon release from peatlands

Systematic review

We searched relevant publications through Web of Science (all databases), Google Scholar, and the China National Knowledge Infrastructure Database between 1945 and March 2021 with the following combinations of keywords: (drain* OR lower* water table OR standing water depth OR ground water table level drawdown OR decline OR drought OR dry*) with (peatland* OR mire* OR fen OR bog OR swamp OR marsh*) with (soil respiration OR heterotrophic respiration OR microbial respiration OR soil CO₂ OR soil carbon decompos* OR soil carbon minerali* or peat subsidence). Using these search terms, we initially identified 2120 different publications. To reliably evaluate WT decline impacts on SR and peat subsidence-associated soil CO₂ emissions, the following further criteria were applied:

1) Only paired studies with pristine peatland (i.e., undrained, near-natural peatland without direct drainage history) as a control and pristine peatland with direct WT decline (due to drainage and land use or climate-induced drying) as a treatment were included by carefully checking the descriptions of field conditions from the publications. For the pristine peatlands, we included the peatland only if the peat soil had at least 30% dry organic matter, a peat depth of >40 cm¹, and did not have any direct drainage history². We acknowledge that few, if any, untouched and completely pristine peatlands currently exist, particularly in Europe.

2) WT decline in peatlands referred to only the WT depth lowered by drainage or climate-induced drying and/or additional management practices related to C or N input (e.g., manure/N fertilizers); treatments in which WT decline was combined with manipulated warming, elevated CO₂, N deposition, etc., were excluded, while individual treatments (i.e., peatlands affected by WT decline without additional warming, elevated CO₂, N deposition treatments, etc.) were included, as the primary objective of this study was to evaluate the responses of peatland C decomposition to WT decline.

3) Each individual study included SR or at least one of its components (HR and AR), and the measurement intervals were at least monthly. The in situ measurements of SR or its components (HR and AR) covered at least the growing or nongrowing season in temperate/boreal climate zones and the whole wet or dry season in (sub)tropical climate zones.

4) Both in situ and soil core/microcosm/mesocosm measurements of SR or its components (HR and AR) were included. SR and its components were exclusively measured using the chamber method. The results of the latter group were used to test the results of the former.

Finally, 386 paired in situ and 21 paired soil core incubation measurements of SR or its components (HR and AR) were extracted from 63 in situ studies and 9 soil core studies, respectively (see Supplementary Data A). Furthermore, to estimate HR emissions from global drained peatlands, the in situ measured paired peat subsidence rate (R_ps, cm yr^–1) and drainage duration (i.e., years since first drainage) and the proportion of peat subsidence rate attributed to oxidation (P_o, %) and drainage duration, as well as the soil (0–30 cm) organic C and bulk density in pristine peatlands, were extracted from peer-reviewed publications. In drained boreal and temperate peatlands, most studies measured the total subsidence (in meter) during a certain drainage period, therefore the average R_ps was calculated as the ratio of total subsidence and drainage years. It was assumed that the R_ps was faster at the beginning and lower at the end of drainage duration, so the average subsidence rate is the rate for the middle year of the drainage duration⁴¹. The remaining studies directly showed the in situ measured R_ps at the ith year of drainage. A similar procedure was applied for the P_o in the ith year of drainage. In sum, 230 paired R_ps–drainage duration observations and 49 paired P_o–drainage duration observations, as well as 76 SOC and 63 BD in pristine peatlands, were taken from 80, 25, 58, and 44 studies, respectively (see Supplementary Data B).

Data compilation

To systematically evaluate the impacts of WT decline on SR in pristine peatlands and clarify the underlying mechanisms, we obtained data related to SR and its components (HR and AR) together with environmental variables such as the mean annual temperature [MAT], mean annual precipitation [MAP], peat depth [PD], WT depth [WTD], soil water content [SWC], soil temperature [T_s], soil redox potential [E_h], soil air oxygen level [O₂], soil bulk density [BD], soil pH [pH], soil organic carbon [SOC], soil total nitrogen [TN], soil total phosphorus [TP], soil ammonium [({{{{rm{NH}}}}}_{4}^{+})], soil nitrate [({{{{rm{NO}}}}}_{3}^{-})], soil dissolved organic carbon [DOC], microbial biomass carbon [MBC], microbial biomass nitrogen [MBN], dissolved total phosphorus [DTP], belowground biomass [BGB], iron [Fe³⁺, Fe²⁺] and sulfate [({{{{rm{SO}}}}}_{4}^{2-})] when possible. If available, other important information, such as geographic location (latitude, longitude), climate and WT decline driver and duration, intensity, peatland type, R_ps, P_o, nutrient type, inundated condition, microtopography, and plant functional types, was recorded. For WT decline intensity, net WT declines greater and less than 30 cm were defined as deep and shallow declines, respectively, according to the IPCC wetland report⁴². The abovementioned information about pristine peatlands and peatlands affected by WT decline is compiled in Supplementary Data A and B.

We subsequently extracted the mean ((bar{X})), standard deviation (SD) and replicates (n) from different publications. If studies reported standard error (SE) rather than SD, then SD was calculated by SE (sqrt{n}). If studies reported only the median, maximum, minimum, and 25th and 75th percentiles, then the mean and SD were estimated following the mathematical equations recommended by ref. ⁶⁰. If neither SD nor SE was reported, then the missing SD was estimated by multiplying the reported mean by the average coefficient of variation (CV) obtained from the remaining observations, resulting in both the mean and SD being reported⁶¹. The data were either obtained directly from tables and texts or extracted by digitizing graphs using Getdata Graph Digitizer software (version 2.26, Russia).

The final database consisted of 250 paired SR, 101 paired HR and 35 paired AR in situ observations. Only 35 paired observations simultaneously reported SR, HR, and AR. Twenty-one paired SR soil core incubation measurements were also collected to test the results of the in situ measurements. The dataset mainly originated from Europe, North America, and Southeast Asia, and most studies (>70%) were conducted in temperate and boreal peatlands in the Northern Hemisphere (Fig. 1a). Moreover, 230 paired R_ps–drainage duration observations and 49 paired P_o–drainage duration observations (Fig. 5a, b) and an additional 485 drainage year (Supplementary Fig. 9) observations classified by climate zone (i.e., boreal, temperate and tropical) and land use (i.e., agriculture, forestry, and grassland) were collected. A total of 76 SOC and 63 BD measurements from pristine peatlands categorized by climate zone (i.e., boreal, temperate, and tropical) were extracted to estimate R_ps by oxidation and associated soil HR from global pristine peatlands due to drainage activities (Supplementary Fig. 10 and Supplementary Data B). In this study, we were unable to estimate climate drying-induced net CO₂ emissions through soil HR, as the areas of pristine peatlands affected by climate drying currently remain unknown.

Meta-analysis

To assess the relative changes in SR and its components (HR and AR), as well as environmental variables (e.g., SOC, BD, T_s, etc.) due to WT decline, the log-transformed response ratio (RR) was used:⁶²

The results are presented as the percent change ((e^lnRR  – 1) × 100). The variance (v) of RR was estimated using the following equation:

where X_t and X_c indicate the means of the treatment and control, SD_t and SD_c indicate the SDs of the treatment and control and n_t and n_c indicate the numbers of replicates in the treatment and control, respectively.

However, in our study, approximately 60% of the WTD and E_h observations for the peatlands in pristine condition (control) and affected by WT decline (treatment) showed opposite signs; e.g., the pristine peatlands generally exhibited positive WTDs (higher than the peat surface) and negative E_h values, while those affected by WT decline exhibited negative WTDs (lower than the peat surface) and positive E_h values. Since it is impossible to calculate the logarithm of negative values, we introduced a new study index (net changes) for these two variables in our meta-analysis according to ref. ⁶³:

where X_t and X_c indicate the paired annual mean WTD and E_h for the treatment and control, respectively, and D indicates the difference between the treatment and control.

The SD and variance (v) of D were estimated using the following equation:

where SD_t and SD_c indicate the SD of the treatment and control and n_t and n_c indicate the number of replicates for the treatment and control, respectively.

The weighted mean RR or D was calculated by individual RR or D with bias-corrected 95% confidence intervals (CIs) using the rma.mv function in the metafor package in R software (R core team, 2019)⁶⁴, in which the variable “study” was regarded as a random effect to account for the dependence of observations derived from the same study. The impact of WT decline on a response variable was considered significant if the 95% CI did not overlap 0⁶⁵. Differences between subgroups (e.g., WT decline driver, climate zone, drainage duration) were considered significant if the 95% CIs did not overlap each other⁶⁵. The frequency distribution of RR was calculated to test variability among individual studies using the Gaussian function (i.e., normal distribution)⁶⁶.

Estimation of peat subsidence rate by oxidation and associated HR rate

Drainage has induced widespread peat subsidence and associated large CO₂ release through soil HR and consequently reduced the sustainable utilization of drained peatlands and contributed to global warming^11,12. In this study, we estimated the spatial patterns of R_ps by oxidation and associated soil HR from global drained peatlands. Using the 230 paired R_ps and drainage duration observations synthesized in this study, we first constructed empirical models between R_ps and drainage duration for drained peatlands categorized by climate zone (boreal, temperate and tropical climate) and land use (i.e., agriculture, forestry and grassland) (Fig. 5a, b). The values of R_ps for certain groups classified by climate zone and land use could be estimated by using the corresponding empirical models established in this study and reported drainage durations that were extracted from the literature. The empirical models categorized by climate zone and land use are listed below (Fig. 5a, b):

where R_ps indicates the peat subsidence rate (cm yr^–1), Dur is the drainage duration, and the numbers indicate coefficients for the established empirical models. Bor, Tem, and Tro indicate boreal, temperate, and tropical climate zones, respectively. Agr, For, and Gra represent agriculture, forestry, and grassland land uses, respectively. We note that it was not possible to further distinguish these models between boreal and temperate climate zones and among agriculture, forestry, or grassland land use in tropical climates, as there is currently a lack of sufficient measurements, which warrants more research.

However, the R_ps is triggered by a combination of processes such as physical compaction by heavy equipment or livestock trampling and shrinkage through contraction of organic fibers when drying, consolidation by loss of water from pores in the peat and oxidation owing to the breakdown of peat organic matter^10,11,12. Therefore, to reliably estimate the soil HR rate from R_ps due to oxidation, the proportion of R_ps attributed to oxidation (P_o, in %) should be considered¹². Using the 49 paired P_o and drainage duration observations synthesized in this study, we then constructed empirical models between P_o and drainage duration for drained peatlands that were also categorized by climate zone (boreal, temperate, and tropical climate) and land use (agriculture, forestry, and grassland) (Fig. 5c, d). Similarly, the P_o values of certain groups classified by climate zone and land use could be estimated by using the corresponding empirical models established in this study and reported drainage durations that were extracted from the literature. The empirical models categorized by climate zone and land use are shown below (Fig. 5c, d):

where P_o indicates the proportion of R_ps attributable to oxidation, Dur is the drainage duration, and the numbers indicate coefficients for the established empirical models. The abbreviations Bor, Tem, Tro, Agr, For, and Gra have been described previously. We note that the different land uses shared the same models across temperate and boreal climates and tropical climate due to a lack of sufficient global observations. This will also induce some uncertainties in our analysis.

Furthermore, the soil HR (F_HR, Mt C yr⁻¹) due to peat oxidation induced by drainage was estimated using the following equation according to ref. ¹¹:

where SOC (g kg^–1) and BD (g cm^–3) indicate the soil (0–30 cm) organic C concentration and bulk density of pristine peatlands, respectively; A (×10³ km²) indicates the drained peatland area; i indicates the climate zone (boreal, temperate or tropical); j indicates the land use (agriculture, forestry or grassland); and R_ps (cm yr^–1) and P_o (%) are described in Eqs. (6–11). Datasets of the SOC concentration and BD and R_ps due to oxidation were systematically reviewed and bootstrapped and categorized by climate zones and land uses (see Supplementary Fig. 10 and Supplementary Data B). Regarding the large uncertainties for areas of drained peatlands, we combined two previously published datasets (72, 61, 22, 37, 43, 26, 94, 109, and 39 × 10³ km² by ref. ¹⁸, and 37, 55, 4, 109, 63, 58, 96, 72, and 1 × 10³ km² by ref. ²⁰. for agriculture-, forestry- and grassland-drained peatlands in boreal, temperate and tropical climate zones, respectively) and obtained their mean values with 95% CIs (for details, see bootstrapping procedure in Data analysis). Uncertainties (i.e., 95% CI) in total HR (δF_HR) were propagated according to the Gaussian random error propagation principle as follows:

where δF_HR, δR_ps, δP_o, δSOC, δBD, and δA indicate the 95% CIs of total soil HR, R_ps, P_o, SOC, and BD and drained peatland area, respectively, and i and j indicate the climate zone (boreal, temperate, tropical) and land use (agriculture, forestry, or grassland), respectively.

To further estimate the total SR (F_SR, Mt C yr⁻¹) and its uncertainty (δF_SR) from global drained peatlands, the following equations were used:

where C_HR (%) indicates the mean relative contribution of HR to SR from simultaneously measured SR, HR, and AR from our meta-analysis (see Supplementary Fig. 11 and Supplementary Data A) and i and j indicate the climate zone (boreal, temperate, tropical) and land use (agriculture, forestry, or grassland), respectively. F_HR and δF_HR are given in Eqs. (12, 13). We note that the C_HR could be classified only by climate zone, as there is a lack of sufficient measurements of land use; that is, the different land uses under the same climate shared the same C_HR value, which may induce uncertainties in estimating the total SR from global drained peatlands.

Regarding the abovementioned lack of sufficient measurements for distinguishing between boreal and temperate drained peatlands, we also used another method to estimate the annual total HR and SR from global drained peatlands. Specifically, we obtained the mean values of R_ps by oxidation across boreal and temperate drained peatlands for each land use (i.e., climate zones were classified as boreal+temperate or tropical) (Supplementary Fig. 14 and Supplementary Table 1). The estimation process was the same as that previously described. The different estimation methods were likely to provide results with greater convergence.

Data analysis

Significant differences in observed variables were tested by performing nonparametric analysis. Specifically, tests with two independent samples (i.e., Mann–Whitney U test) were used for only two variables (i.e., to compare the contribution of HR to SR between pristine and drained peatlands), and tests with two or more independent samples (i.e., Kruskal–Wallis test and pairwise comparisons) were used if there were three or more variables (i.e., SOC, BD and R_ps due to oxidation in the boreal, temperate and tropical climate zones or different land uses). Linear or nonlinear regression analysis was performed to examine the relationships between the responses of SR and its components with environmental variables or the peat subsidence rate with drainage duration.

To reliably estimate the uncertainties in R_ps by oxidation, SOC, BD, drained peatland area, and relative contribution of HR to SR, bootstrap resampling with 10000 iterations was conducted using the boot package, and 95% CIs were calculated using the “basic” type. The ggplot 2 package in R software (R core team, 2019) was used for statistical analysis. Data were expressed as the means with their 95% CIs, and significance of the regression analyses was indicated at the level of p < 0.05 across the study.

Source: Ecology - nature.com