Gridded maps of wetlands dynamics over mid-low latitudes for 1980–2020 based on TOPMODEL

The conceptual flow chart of the process is provided in Fig. 1. We used seven reanalysis SM data (Table 2) masked with soil temperature (ST) and soil freeze/thaw status to calculate water table depth, i.e. the input of TOPMODEL, given the obvious disagreements between the input datasets. The diagnostic algorithms based on TOPMODEL were used following Stocker et al. (ref. ²⁰) and Xi et al. (ref. ²⁵), where the optimized parameters were calibrated with long-term maximum wetland areas from four observation-based wetland datasets (Table 1). Details about these datasets and computational processing are shown as follows.

Reanalysis soil moisture datasets

Seven long-term reanalysis SM datasets used in this study include NCEP-DOE (National Centers for Environmental Prediction-the Department of Energy)²⁶, MERRA-Land (the Modern-Era Retrospective Analysis for Research and Applications)²⁷, MERRA-2 (ref. ²⁸), GLDAS-Noah v2.0 (the Global Land Data Assimilation System)²⁹, GLDAS-Noah v2.1 (ref. ²⁹), ERA5 (European Environment Agency)^30,31, and ERA5-Land^30,31. Key characteristics of the seven SM data are listed in Table 2. The datasets differ by their spatial and temporal resolutions, the time-period they cover, as well as the definition of the soil layers. More details are provided for each dataset below.

• NCEP-DOE

  NCEP-DOE is an updated version of the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) Reanalysis 1 project, which uses a state-of-the-art analysis/forecast system to perform data assimilation with past data from 1948 to the present³². NCEP-DOE features the newer physics and observed SM forcing and also eliminates several previous errors, such as oceanic albedo and snowmelt term during the entire period, and snow cover analysis error from 1974 to 1994 (ref. ²⁶). With a spatial resolution of about 210 km, there are two vertical soil layers in NCEP-DOE for both SM and ST: 0–0.1 and 0.1–2 m.

• MERRA-Land and MERRA-2

  MERRA-Land soil moisture is generated by driving the Goddard Earth Observing System model version 5.7.2 (GEOS-5.7.2) with meteorological forcing from the MERRA reanalysis product²⁷. The precipitation forcing in MERRA-Land merges MERRA precipitation with a gauge-based data product from the National Oceanic and Atmospheric Administration (NOAA) Climate Prediction Center, and the Catchment land surface model used in MERRA-Land is updated to the “Fortuna-2.5” version²⁷. MERRA-2 intends to replace the original MERRA reanalysis and ingests important new data types²⁸. The Catchment model in MERRA-2 has been updated with rainfall interception and snow model parameters of MERRA-Land, and the precipitation correction is a refined version of MERRA-Land. For MERRA-Land and MERRA-2, there is only one layer for SM from the surface to the bedrock, with “depth-to-rock” depending on local conditions. ST is computed on six vertical soil layers: 0–0.10, 0.10–0.29, 0.29–0.68, 0.68–1.44, 1.44–2.95, and 2.95–12.95 m.

• ERA5 and ERA5-Land

  ERA5 is the fifth generation ECMWF (European Centre for Medium-Range Weather Forecasts) reanalysis of global climate and weather, replacing ERA-Interim^30,31. Based on a decade of developments in model dynamics and data assimilation, there is a significantly enhanced horizontal resolution (31 km), temporal resolution (hourly) and uncertainty estimation. ERA5 covers 1979–2020 and continues to be updated in near-real-time. ERA5-Land is produced with a finer horizontal resolution of 9 km by running the land component of the ERA5 climate reanalysis but without data assimilation. By March of 2021, the ERA5-Land outputs are only available since 1981. SM and ST are computed on four vertical soil layers (0–0.07, 0.07–0.28, 0.28–1, and 1–2.89 m) for both ERA5 and ERA5-Land.

• GLDAS-Noah v2.0 and GLDAS-Noah v2.1

GLDAS is a global, moderate-resolution (0.25° × 0.25°) offline terrestrial modeling system developed by NASA Goddard Space Flight Center (GSFC) and the NOAA National Centers for Environmental Prediction²⁹, thus similar to ERA5. To produce optimal fields of land surface variables in near-real-time, it incorporates satellite- and ground-based observations. GLDAS-Noah drives the Noah land surface model and has two components: one forced with the Princeton meteorological forcing data (i.e. GLDAS-Noah v2.0) and the other forced with a combination of model and observation (i.e. GLDAS-Noah v2.1). GLDAS-Noah v2.0 covers the period 1948–2014, while GLDAS-Noah v2.1 is available from 2000 to the present. There are four vertical layers in the Noah land surface model for both ST and SM: 0–0.1, 0.1–0.4, 0.4–1, and 1–2 m.

Observation-based wetland/flooded area data

In terms of large uncertainties in current wetland datasets (Table 1) we selected four widely used and available satellite/satellite-based wetland/flooded area data including GIEMS-2 (ref. ¹⁴), RFW (the Regularly Flooded Wetland map)¹⁰, WAD2M (a global dataset of Wetland Area and Dynamics for Methane Modeling)³³, and G2017 (the pantropical wetland extent from an expert system model)⁹ for parameter calibration. Among them, GIEMS-2 and WAD2M include monthly wetland dynamics, while RFW and G2017 are static. The comparison of the four wetland datasets is shown in Supplementary Fig. 1; details on each data are provided below.

• GIEMS-2

  The GIEMS-1 is the first global estimate of monthly inundated areas, derived from passive microwave land surface emissivity³⁴. With a 0.25° × 0.25° resolution, GIEMS-1 documents a mean annual maximum inundated area of 9.5 Mkm² for 1993–2007 (including open water, wetlands, and rice paddies, but excluding large lakes), which shows good agreement with existing independent, static inventories as well as regional high-resolution synthetic aperture radar observations³⁴. Based on similar retrieval principles with GIEMS-1, GIEMS-2 is developed to less depend on ancillary data with an updated microwave emissivity, and correct a known overestimation over low vegetated areas from GIEMS-1 (ref. ¹⁴). The period is extended to 1992–2015 for GIEMS-2 and can be updated with the availability of observations. Globally, the mean annual maximum and long-term maximum inundated extent after removing the rice paddies using the Monthly Irrigated and Rainfed Crop Areas dataset (MIRCA2000)³⁵ for the period 1992–2015 are 6.7 and 10.9 million km² (hereafter Mkm²; sum of mean annual maximum or long-term maximum inundated extent for each grid cell) respectively. The rice paddies are removed here as they are not natural wetlands and cannot be simulated with TOPMODEL.

• RFW

  RFW is a static, high-resolution map (15 arc-sec) of regularly flooded wetlands, developed by overlapping flooded areas (permanent wetlands and flooded vegetation classes) for 2008–2012 from the ESA-CCI land cover map³⁶, mean annual maximum inundated areas (including wetlands, rivers, small lakes, and irrigated rice) for 1993–2004 from GIEMS-D15 global inundation extent (downscaled using GIEMS-1)³⁷, and long-term maximum surface water areas for 1984–2015 from JRC global surface water bodies product¹³. The large permanent lakes and reservoirs are distinguished using the HydroLAKES database³⁸. Globally, RFW covers 9.7% of the land surface area (~13.0 Mkm²) including wetlands, river channels, deltas, and flooded lake margins, but excluding large lakes¹⁰. Due to the mean annual maximum or long-term maximum inundation/surface water extent for 1984–2016 from the three wetland data is used, we treated RFW as the long-term maximum wetland extent in this study. Besides, given that GIEMS-D15 includes artificial rice paddies, we removed them with MIRCA2000 from RFW (~11.9 Mkm² after removing rice paddies).

• WAD2M

  WAD2M dataset used in this study is an improved version of the SWAMPS v3.2 from Jensen et al. (ref. ¹⁵), covering the years 2000 to 2018. With a spatial resolution of 25 km × 25 km, this data was used as input wetland area data of phase 2 of the Global Methane Budget³³. Given that the initial SWAMPS failed to detect wetlands lacking surface inundation and to differentiate between lakes, wetlands, and other surface water bodies, Zhang et al. (ref. ³³) modified it using a series of independent static wetland distribution data^7,9,39,40,41 in an attempt to include missing wetlands under dense canopies. Besides, they removed inland waters (lakes, rivers, and ponds) and rice agriculture with JRC and MIRCA2000, respectively. Globally, the mean annual maximum and long-term maximum wetland extent for the period 2000–2018 estimated by WAD2M are 8.1 Mkm² and 13.2 Mkm² (sum of mean annual maximum or long-term maximum inundated extent for each grid cell) respectively.

• G2017

G2017 (ref. ⁹) is a static, pantropical wetland and peatland extent map (covering 60°S–40°N) at 232 m × 232 m resolution, derived from a hybrid expert model system. With three biophysical indices related to wetland and peat formation (long-term water supply exceeding atmospheric water demand, annually or seasonally waterlogged soils, and geomorphological position where water is supplied and retained), G2017 identifies not only permanently and seasonally wetland areas, but also soil wetness and topographic conditions that favor waterlogging in the absence of flooding for the end of the 20^th century. Given the broad coverage of different types of wetlands, we also treated this map as long-term maximum wetland areas. This ‘pantropical’ data (60°S to 40°N) offers the advantage to include non-flooded wetland areas that are missed in satellite-based wetland products. However, note that not all detected wetlands or peatlands in G2017 have been observed. Rice agriculture was also removed with MIRCA2000 from G2017. The resulting wetland and peatland area for 60°S–40°N is 4.0 Mkm².

The TOPMODEL-based diagnostic model

TOPMODEL as improved by Stocker et al. (ref. ²⁰) and Xi et al. (ref. ²⁵) was used to calculate the inundated fraction from WTD at grid-scale in this study. Based on the assumptions that the local hydraulic gradient is approximated by the local topographic slope and the water table variations can be assimilated to a succession of steady states with uniform recharge, the classical TOPMODEL establishes an analytical relationship between the soil moisture deficit and the distributions of local topographic index within a catchment. At grid-scale, the analytical relationship can be represented as:

where CTI indicates the topographic index, defined as the log of the ratio of contributing area to the local slope. We used the CTI data at 500 m × 500 m resolution from Marthews et al. (ref. ²²), where lakes, reservoirs, mountain glaciers, and ice caps are removed using the Global Lakes and Wetlands Database⁷. The (overline{CT{I}_{x}}) indicates the average of CTI_i of all sub-grids (index i) within the grid cell x. M indicates a tunable parameter that describes the exponential decrease of soil transmissivity with depth²¹. Γ_i is the water table of the pixel i and (overline{{{Gamma }}_{{x}}}) is the mean water table of the grid x. When Γ_i is at the soil surface (i.e. Γ_i = 0), the threshold (CT{I}_{x}^{* }) above which all pixels are flooded for the grid x is derived:

The wetlands area is defined as the flooded areas (i.e. Γ ≤ 0), the flooded fraction in the grid x (f_x) being the percentage of pixels with CTI_i larger than a threshold (CT{I}_{x}^{* }):

with

To reduce the computational costs from the high-resolution CTI data for predicting long time series of wetland area, we used the asymmetric sigmoid function from Stocker et al. (ref. ²⁰) to fit the “empirical” relationship (widehat{{Psi }}) between (widehat{f}) and Γ:

where v_x, k_x, q_x are three parameters of the function. Given a value of parameter M, the three parameters can be derived with a sequence of Γ_x spanning a plausible range of values (−1 m to 2 m) and corresponding f_x from the initial TOPMODEL approach (Eq. (3)). Thus, the wetlands in our study are defined as the flooded area simulated by TOPMODEL. As for the range of parameter M, Stocker et al. (ref. ²⁰) used a global uniform value for M (M = 8) after testing simulated wetland fraction for a range of M (7, 8, 9). Nevertheless, given that distinct topography, soil types, and other intrinsic characteristics in different regions, we considered M as a tunable, spatially heterogeneous, and grid-specific parameter, with a range of 1–15 following Xi et al. (ref. ²⁵). Thus, for each grid cell x there are 15 choices for M, and then 15 sets of (v_x, k_x, q_x). The optimized parameter combination of (v_x, k_x, q_x) is determined by selecting minimum root-mean-square-error (RMSE) between simulated inundated fractions and observations:

where O_i and P_i are observed and simulated wetland fraction, respectively. n represents the time-series length for wetland extent. For simulations calibrated with RFW and G2017, the RMSE was computed with the long-term maximum (hereafter called MAX) monthly wetland area because the two data sets are static and only record the MAX wetland extent. While for simulations calibrated with GIEMS-2 and WAD2M which include temporal variations of wetland area, we calibrated the parameters with all months, mean seasonal cycle, yearly maximum, and MAX wetland area, but only showed the optimal simulations calibrated with MAX wetland area in this work to keep consistency with RFW and G2017. Besides, to provide more choices for users, we combined all of the four wetland datasets (i.e. the union of long-term maximum wetland extent) to generate a new wetland map (hereafter called MAX_all), and then used the MAX_all to calibrate the parameters to produce seven sets of global wetland extent products with seven soil moisture datasets. The simulations calibrated with yearly maximum wetland area from GIEMS-2 and WAD2M and long-term maximum wetland area from MAX_all are also provided in our resulting products.

Finally, to avoid unrealistically high wetland fraction output from the function, the simulated maximum wetland fraction f_x is constrained by the observed MAX wetland area with a parameter ({f}_{x}^{max}) (Eq. (6)), which is different from Stocker et al. (ref. ²⁰). The determination of ({f}_{x}^{max}) is analyzed in the supplemental material in detail (Supplementary Text 1). Once the value of (v_x, k_x, q_x) are determined, the wetland fraction f_x can be directly derived from the monthly water table Γ_x according to Eqs. (4) and (6).

Calculation of water table depth

Water table depth is not computed by land surface models, given their coarse soil vertical discretization. We thus used the saturation deficit of soil moisture (θ_SD) as a surrogate of water table depth, θ_SD being defined as an index consisting of saturated volumetric water content and the “actual” soil depth modified by soil freeze/thaw status:

Subscript l represents the l^th soil layer, l₀ is the number of layers above the first frozen soil layer counted from the top (l = 1 at the soil surface), θ_l is the monthly volumetric water content in the l^th soil layer (m³ m⁻³), (Delta {z}_{l}) is the thickness of the l^th soil layer (m), θ_S is the saturated volumetric water content (in m³ m⁻³ units, uniform over depth).

As formulated in Eq. (7), ({z}_{{l}_{0}}) is the thickness of all soil layers (or depth to bedrock) when there is no frozen soil layer. If there exists at least one frozen layer, ({z}_{{l}_{0}}) is set to the depth of the uppermost frozen soil layer. We excluded the frozen soil layers here given that some important wetland processes such as methane production and transport are insignificant when the soils are frozen. In high latitudes, the presence of frozen soil layers may lead to an overestimation of the wetland fraction due to relatively large θ_SD values even if there is little liquid soil water above the uppermost frozen soil layer. Hence, we used monthly soil temperature (ST) at 70 cm, the Global Record of Daily Landscape Freeze/Thaw Status data⁴², and the Köppen climate classification system⁴³ to refine the frozen mask. When the monthly mean ST at 70 cm is below 0 °C, or soil freezing days are more than 5 in a month, or the grid is classified as the Hot desert (BWh) in the Köppen climate classification system, the wetland fraction for the grid is set to zero. However, it should be noted that the algorithm using the ST at 70 cm could omit some unfrozen soil layers above 70 cm, which could lead to bias in estimation of methane emissions from these unfrozen layers. We provided the global wetland maps in our resulting products, but the potential uncertainties in wetland estimation due to the omitted unfrozen layers should be considered, particularly at high latitudes. We used seven reanalysis SM products to compute θ_SD to provide the uncertainty in SM input (Table 2). All data are re-interpolated to 0.25° × 0.25° resolution.

Evaluation against wetland calibration data and independent satellite products

Although we calibrated parameters of the TOPMODEL-based diagnostic model with the observation-based wetland data, to what extent the simulations can reproduce the spatial patterns and temporal dynamics of the calibration wetland data must be evaluated. For spatial patterns, we calculated the RMSE of wetland area between our simulations and corresponding wetland calibration data following Eq. (5), and evaluated the spatial patterns of simulated wetland extent in two wetland hotspots including Amazon basin and Western Siberia lowlands with three independent global/regional water products. For Amazon basin, we used the global surface water dataset from JRC¹³ (optical satellite images) and the wetland map produced using mosaics of Japanese Earth Resources Satellite (JERS-1) L-band SAR imagery from Hess et al. (ref. ⁴⁴, hereafter H2015). For West Siberian lowlands, we used JRC and the Boreal–Arctic Wetland and Lake Dataset (BAWLD, only covers the north of ~55°N) produced using an expert assessment and extrapolated using random forest modelling from climate, topography, soils, permafrost conditions, vegetation, wetlands, and surface water extents and dynamics⁴⁵. For temporal dynamics, since we only used the static wetland area (long-term maximum) from all of the four observation-based wetland products to calibrate parameters, the simulated temporal dynamics can be evaluated with the two dynamic wetland products (GIEMS-2 and WAD2M). Besides, we also used the terrestrial water storage (TWS) from the Gravity Recovery and Climate Experiment (GRACE), which retrieves relative change in TWS from the monthly anomalies of the Earth’s gravity field for 2003–2016 measured by the twin GRACE satellites^46,47 to evaluate the simulated temporal dynamics.

Source: Ecology - nature.com