Investigating Wetland and Nonwetland Soil Methane Emissions and Sinks Across the Contiguous United States Using a Land Surface Model

We estimated the distribution of CH4 emissions and sinks from wetlands (including freshwater and coastal wetlands) and nonwetland (including wet and dry soils) with a newly developed vertically resolved soil CH4 model, integrated into a global land surface model (ISAM). We calibrated and tested this integrated model with CH4 observations at test sites in the Contiguous United States (CONUS). ISAM is applied across the CONUS to estimate CH4 emissions and sinks given both recent past observed climate and wetland extent, and future climate and wetland extent driven by two scenarios, RCP4.5 and RCP8.5. Estimated net CH4 emissions for the 2000s are 13.8 TgCH4 yr−1, mostly from wetland soils. Estimated net emissions under RCP4.5 and RCP8.5 are 30% and 64% higher, respectively, in the 2090s than in the 2000s due to (1) higher temperature and seasonal wetland extent (driven by higher precipitation in the climate scenarios), which increase modeled methanogenic activity more than methanotrophic activity in soils and (2) altered transport in the soil column and exchange with the atmosphere by modeled transport processes (diffusion, ebullition, and aerenchyma transport). Nonwetland soils emit CH4 (1.4 TgCH4 yr−1) in some areas and take up CH4 (−2.9 TgCH4 yr−1) in other areas, resulting in a net estimated sink for the 2000s; the net nonwetland soil sink increases by 15% and 46% by the 2090s under RCP4.5 and RCP8.5, respectively, mainly due to drier soil conditions, which enhances methanotrophic activity and oxidation of CH4 diffused into soil from a future atmosphere with higher CH4 concentration.


Introduction
Wetlands, including seasonal wetlands and wet soils, make up the largest natural source of methane (CH 4 ; Kirschke et al., 2013;Matthews & Fung, 1987), whereas dry soils are a CH 4 sink (Adamsen & King, 1993). The net exchange of CH 4 between these soils and atmosphere, as a consequence of methanogenic and methanotrophic activities in wetlands and soils, is highly variable in space and time (Topp & Pattey, 1997). The methanogenic process produces CH 4 under anaerobic conditions in wetlands and wet soils with a shallow water table depth (Angle et al., 2017). On the other hand, the methanotrophic process oxidizes CH 4 under aerobic conditions to form other compounds, such as formaldehyde, acts as a sink of CH 4 in dry soils. These processes are sensitive to soil conditions, including temperature, soil moisture, soil texture, soil pH, and water table depth (Cao et al., 1998).
Estimates for CH 4 emissions from wetlands at regional to global scales have received wide attention in the literature (Bloom et al., 2017;Potter et al., 2006;Saunois et al., 2016;Tian et al., 2010Tian et al., , 2012Xu et al., 2010). Several studies have applied different methods to estimate the CH 4 emission and/or sinks for the Contiguous United States (CONUS). Potter et al. (2006) applied an empirical method that links the net ecosystem production (NEP) calculated from the CASA model to the CH 4 emission by multiplying a uniform CH 4 :CO 2 conversion factor, estimating CONUS emissions to be 5.5 TgCH 4 yr −1 . Tian et al. (2010) developed a process-based global biogeochemical model, which coupled five core components (biophysics, plant physiology, soil biogeochemistry, dynamic vegetation, and land use management.), to estimate both CH 4 emissions and sinks from the soils, including wetlands. Xu et al. (2010) applied the same model to attribute the spatial and temporal variability of CONUS CH 4 emission to different environmental variables from 1979-2008 and found precipitation to be the most important factor. Tian et al. (2012) further synthesized the estimation of CH 4 emission and sink from multiple approaches for the whole North America and found a range of 3.3-9.1 TgCH 4 yr −1 for the CONUS region. The most recent global CH 4 budget study (Saunois et al., 2016) estimated CONUS wetland emissions of 8-22 TgCH 4 yr −1 , based on top-down and bottom-up analyses. Wetlands, including seasonal wetlands, act as a major natural emission for the CONUS, contributing 25-45% of the total CH 4 emissions estimated by the bottom-up models (Saunois et al., 2016).
Precipitation events have been found to introduce variability in the magnitude of CH 4 sources and sinks from nonwetland soils (Ni & Groffman, 2018;Saunois et al., 2016), altering water table depth along with oxic and anoxic conditions on relatively short timescales. However, previous publications (Potter et al., 2006;Saunois et al., 2016;Tian et al., 2010;Xu et al., 2010), with the exception of Tian et al. (2012), did not provide estimates of both CH 4 emissions and sinks from nonwetlands across the CONUS. Also, these studies accounted for the effect of major environmental conditions on CH 4 emission (mainly focusing on the modeling of the CH 4 biogeochemistry dynamics), but these studies did not include a detailed linkage between soil hydrology and soil biogeochemistry, which is necessary for the more accurate estimation of both soil CH 4 emission and CH 4 sinks. While the importance of methanotrophic and methanogenic conditions and their interactions with soil conditions (e.g., temperature, soil moisture, soil texture, soil pH, and water table depth) were discussed in recently published landscape-level studies (Angle et al., 2017;Kaiser et al., 2018), the effects of these interactions on CH 4 emissions and sinks were not studied at a regional scale by the previously published studies mentioned above.
The overall objective of this study is to advance our understanding of the spatial and temporal distribution of CH 4 emissions and sinks from the wetlands, seasonal wetlands and nonwetlands (including wet and dry soils) for historical and future environmental conditions across the CONUS region using a land surface model, the Integrated Science Assessment Model (ISAM). Here wetlands refer to both inland freshwater wetlands and coastal wetlands. It is important to note that CH 4 emissions of coastal wetlands with high salinity are lower than inland wetland emissions due to relatively higher sulfate availability and generally higher abundances of sulfate-reducing bacteria, which outcompete methanogenesis (Poffenbarger et al., 2011). ISAM estimates CH 4 emissions and sinks from both wetlands and nonwetland soils under historical , as well as future environmental conditions under two future scenarios-RCP4.5 and RCP8.5-over the period 2016-2100 across the CONUS. In addition, the model also estimates soil CH 4 production, consumption, and three major transport processes (by aerenchyma, ebullition, and diffusion), which result in the spatial variation of CH 4 emissions and sinks. This study provides a foundation for understanding the spatial and temporal distribution of soil CH 4 sources and sinks for interpretation of atmospheric CH 4 measurements and budget, as well as the feedback between soil CH 4 emissions and climate change.

Model Description
To accomplish our overall objective, we first couple ISAM (El-Masri et al., 2015;Jain et al., 2005;Song et al., 2013;Yang et al., 2009) to a newly developed CH 4 module containing parameterizations of soil CH 4 dynamics. Each model grid cell is of 0.5°× 0.5°spatial resolution and is subdivided into 24 plant functional types (PFTs), including 20 vegetation types, desert, polar desert, bare soil, and land ice, based on the fractional area of PFTs (Meiyappan & Jain, 2012). The model is run at a 1-hr temporal resolution.
ISAM includes three coupled components: (1) a biogeophysical component, (2) a dynamic vegetation component, and (3) a soil biogeochemistry component. The biogeophysical component consists of soil energy, soil hydrology (Barman et al., 2014a(Barman et al., , 2014b, and snow processes (Barman & Jain, 2016) and calculates energy and water fluxes (Barman & Jain, 2016;El-Masri et al., 2013). The hydrological model includes precipitation, dew formation, soil and canopy evaporation, canopy transpiration, surface, and subsurface runoff, groundwater discharge/recharge, and changes in water quantity (Barman et al., 2014b, Song et al., 2016. Soil water content is calculated based on Richard's equation and is redistributed with dynamic root growth and distribution in soil layers extended to a depth of 3.5 m globally, with the thickness of each soil layer increasing exponentially with the depth. The groundwater discharge/recharge process is implemented through a dynamic coupling between the bottom soil layers and an unconfined aquifer. The model accounts for the bedrock depth distribution based on the conterminous United States multilayer soil characteristics data set (CONUS-SOIL; Miller & White, 1998). The model calculates the thermal and hydrological properties for each soil layer as a function of soil liquid and ice water content, soil temperature, soil texture, and soil organic C and gravel content (Lawrence & Slater, 2008). To better estimate soil hydrological properties in fine-textured soils, the model also considers the impacts of root growth, which changes the soil structure and properties, notably by affecting soil hydraulic conductivity in fine-textured soils (Song et al., 2016).
The dynamic vegetation component describes the above-ground C and N cycle through parameterization of the physiological processes, including photosynthesis, carbon, and nutrient allocation in the plant, plant phenology, and plant rooting depth for both natural vegetation (El-Masri et al., 2015;Yang et al., 2009) and cropland (Song et al., , 2014. In addition to CH 4 dynamics, a 10-layer soil biogeochemistry component was recently incorporated in ISAM to estimate the vertical soil organic C (SOC) and N (SON) profiles. This vertically resolved soil biogeochemistry component includes (1) vertical transport of SOC (including bioturbation and cryoturbation processes) and (2) vertically resolved SOC decomposition (Yang et al., 2009) and related environmental modifiers. The vertical transport of SOC among the soil layers is represented by diffusion-advection processes (Elzein & Balesdent, 1995). The environmental conditions are also linked to the SOC decomposition at each soil layer. The module is calibrated and evaluated against multiple observations of SOC and soil Δ 14 C profile. The vertically resolved SOC module enables better representation of the SOC substrates for CH 4 production at different soil layers.
The maximum hydrologically active soil column depth used in ISAM was extended to 3.5 m (Barman et al., 2014a). The soil biogeochemistry component accounts for the processes governing the accumulation and formation of SOC and SON profiles include the vertical transport of SOC and SON, multilayer SOC decomposition, depth-dependent abiotic control on SOC decomposition and below-ground N processes (Yang et al., 2009). The model structure and parameters have been evaluated in previous literature (Barman et al., 2014a(Barman et al., , 2014bEl-Masri et al., 2015;Gahlot et al., 2017;Jain et al., 2013;Meiyappan et al., 2015) and the performance of ISAM under multiple scales have been extensively evaluated in intercomparison projects (FACE-MDS: Walker et al., 2014;NACP: Schaefer et al., 2012;MsTMIP: Huntzinger et al., 2017;Mao et al., 2015;Zscheischler et al., 2014;GCP: Ahlström et al., 2015;Jung et al., 2017;Zhao et al., 2016).

Model Extension to Include CH 4 Dynamics
The CH 4 module ( Figure 1) is coupled to ISAM through soil hydrology and soil biogeochemistry components to estimate soil CH 4 and O 2 sources and sinks, concentrations, and net emissions from the land to the atmosphere. In this study, we consider both inland freshwater wetlands and coastal wetlands. The CH 4 module shares the same 10 soil layers (to 3.5 m depth) as the soil biogeophysics and biogeochemistry modules, with the thickness of each soil layer increasing exponentially with depth (Barman et al., 2014a). Both O 2 and CH 4 transports in the soil column are represented in the model. O 2 has several orders of magnitude lower diffusivity in the water than in the air. Hence, soil saturation conditions below the water table tend to lead to anoxic zones, whereas drier conditions above the water table leading to oxic conditions. In anoxic conditions, organic matter is decomposed by methanogens that produce CH 4 . Once CH 4 is produced in anoxic soils, it can diffuse to the atmosphere. On its way to the atmosphere, CH 4 can be partly or completely lost due to oxidation by methanotrophic bacteria in the soil oxic zone, which uses CH 4 as their only source of energy and carbon (Segers, 1998). CH 4 from the atmosphere can also diffuse into the nonwetland soil column (Equation S1 .27 in the supporting information), and this diffused CH 4 can also be oxidized (See Text S1.4 for the detailed description). A detailed description of the module can be found in supporting information Text S1, and all the related equations are described in supporting information Table S2. All the abbreviations appearing in the paper main text and in the equations are described in Table S1. The CH 4 module developed here accounts for all the major processes affecting the aqueous and gaseous phase volumetric concentration of both CH 4 and O 2 in the soils.
The CH 4 module accounts for six specific processes: (1) soil CH 4 production in the soil column by methanogenic activities (P), (2) soil CH 4 released to the atmosphere by gas diffusion (FD ↑ ), (3) transport of CH 4 from the root zone of aerenchymatous plants to the atmosphere (A), (4) ebullition from wetlands to the atmosphere (E), (5) soil CH 4 loss in the soil column by CH 4 oxidation in the soil (O), and (6) atmospheric CH 4 uptake by soil due to gas diffusion (FD ↓ ). Here we report the net diffusion flux (FD = FD ↑ − FD ↓ ) calculated using the Equation S1.28. A positive FD represents a CH 4 emission to the atmosphere.
The soil production and losses of CH 4 are linked to the hydrological dynamics in the model. ISAM first calculates soil moisture with a hydrology module, and then the diffusion of the O 2 in the soil column to determine the O 2 concentration in soil. Production of CH 4 happens in anoxic soil where the soil is saturated, and O 2 concentration is low. Anoxic conditions also exist in unsaturated soil where anoxic microsites have limited access to O 2 due to the soil tortuosity (Angle et al., 2017). Losses of CH 4 are considered in the model through the CH 4 oxidation where O 2 concentration is high in dry soil layers.
Since ISAM soil hydrology dynamics are coupled with CH 4 production and oxidation, the variation of the water table calculated from ISAM affects the magnitude of CH 4 production and oxidation, and hence net CH 4 emission. Soil with a shallower water table contains more saturated soil for CH 4 production and less unsaturated soil for CH 4 oxidation and thus has a higher CH 4 emission than the soil with a deeper water table. This linkage in the model can result in changes in CH 4 emissions due to the variation of the water Figure 1. The schematic diagram illustrates the CH 4 dynamic processes represented in the ISAM CH 4 module. The module includes CH 4 production through methanogens in the anaerobic zone; CH 4 transport in the soil column and to the atmosphere through diffusion; CH 4 transport to the atmosphere through aerenchyma and ebullition; and CH 4 consumption through methanotrophs in the aerobic zone. The module treats wetlands and nonwetlands differently. For permanent wetlands or seasonal wetlands, the entire soil column is considered anaerobic, and ebullition is included as a direct pathway for CH 4 transport to the atmosphere. For nonwetlands, CH 4 oxidation in the aerobic zone consumes the CH 4 transported upward from the anaerobic zone below the water table, or CH 4 that may defuse downward from the atmosphere insufficiently to dry soils. The model couples all these functions with the dynamics of soil biogeophysics (energy and hydrology) and biogeochemistry (carbon and nitrogen cycles), which enables the estimation of soil CH 4 sources and sinks under a variety of climate and soil conditions across the CONUS. The aqueous and gaseous phase concentrations of CH 4 and O 2 for each soil layer are calculated using Equations S1.1a-S1.1d. The various phases of diffusion rates of CH 4 and O 2 at each soil layer (Equations S1.7-S1.11 and supporting information Text S1.1) are determined based on soil texture, the soil water table depth (z wt ), and soil moisture (Millington & Quirk, 1961;Moldrup et al., 2003) with boundary conditions given by Equation S1.12. The atmospheric CH 4 and O 2 concentrations are calculated using Equation S1.13. The sum of all sources, sinks and transportation fluxes of CH 4 and O 2 from soil to atmosphere (Equation S1.14), including P (Equations S1.15-S1.17), O (Equations S1.18-S1.20), E (Equations S1.21-S1.23), and A (Equations S1.24-S1.26) are added to the diffusion equation, Equation S1.1. The diffusion equation is then solved using the Crank-Nicholson finite differencing method. Finally, the sum of wetlands (N w ) and nonwetlands (N nw ) soil CH 4 emission EM (= N w + N nw ) equals to the sum of FD, E, and A fluxes, where FD (= FD w + FD nw ) is the sum of wetlands (FD w ) and nonwetlands (FD nw ), is described by Equations S1.28-S1.29.
ISAM accounts for wetland and nonwetland soil types within a grid cell (defined as a time-varying fraction of the grid cell), rather than assuming the entire grid cell to be either wetlands or nonwetlands. The wetland differs from the nonwetland in the calculation of soil hydrology and soil biogeochemistry. We have added SOC and SON pools for all wetland-specific biogeochemical processes in the soil biogeochemistry component. Aerenchyma development can happen for plants in seasonal anoxic soil and also in nonwetlands under seasonal wet soil with partial waterlogging conditions (Jackson & Armstrong, 1999). Therefore, ISAM includes aerenchyma transport not only for wetlands but also for nonwetlands. The model also accounts for the salt water intrusion for the coastal wetlands (hereafter referred to as coastal wetland effect), which reduces the methanogenesis. This effect is parameterized by calibrating the value of CH 4 /CO 2 production ratio parameter (r CH4:CO2 ). The effects of sea level rise on the future changes in coastal wetland extent are also considered using the parameterization scheme described in section 3.
Two types of wetlands (marsh and saltwater marsh) considered in ISAM are dominated by the herbaceous plants (more than 95%), a small amount of forest (i.e., mangrove, about 0.53%) (Krauss et al., 2011), and crop (i.e., rice, about 2.3%; Portmann et al., 2010) over the CONUS. In ISAM herbaceous wetland plants are C3 and C4 grasses, the mangrove forest is forest PFT (Krauss et al., 2011), and the crop is rice. We calibrated/evaluated the model for the herbaceous wetland plants using seven sites (see section 5.1). For rice, we applied a specific calibrated parameterization describing rice phenology (Lin et al., 2017) and use one rice site for model calibration. There is no other rice FLUXNET site in the CONUS that provides sufficient information to evaluate ISAM. In addition, there is no AMERIFLUX site for mangrove wetlands (Knox et al., 2019) in the CONUS. Therefore, we are not able to calibrate/evaluate ISAM for forest wetland type.
The coupled model was first calibrated and evaluated with site-level data, including CH 4 sources and sinks, which is described next.

Model Calibration and Evaluation
Measured CH 4 fluxes from 10 FLUXNET sites across the CONUS region were used to calibrate and evaluate model parameters (Table 2), including aerenchyma transport factor (F aere_s ), r CH4:CO2 , potential CH 4 oxidation rate (O PM ), Q 10 coefficient and scaling factor of ebullition threshold (F ebul_s ), and evaluate ISAM-calculated CH 4 fluxes at each site. Based on the FLUXNET network, we found site data for three different types of wetlands (i.e., marsh, rice, and salt marsh) and upland. Marsh and rice are the freshwater wetlands and salt marsh is a saline water wetland. Overall, we used five sites (two sites for marsh and one site each for the salt marsh, upland, and rice) for model calibration, and five sites (three sites for marsh and one site each for salt marsh and upland) for model evaluation. To evaluate modeled emissions for upland, we used the results from Curry (2007) for the Colorado State measurement site and for the regional scale to show that ISAM results are consistent with the Curry (2007) in section 5.2. While there are emissions data available for additional sites (e.g., Knox et al., 2019;Olefeldt et al., 2013;Turetsky et al., 2014), missing input data for ISAM at these sites limits their use for evaluation. The data is also available for the peatland sites (e.g., Iversen et al., 2018), but we can not use these sites data, because the model does not account for peatland soil dynamics. We note that while the US-Tw1 marsh calibration site and US-Myb marsh evaluation site are located less than 15 km away from each other, their soil texture, porosity, and water table depths are different  Table 1).
Site-specific observation data for seven climate variables (surface air temperature, specific humidity, precipitation rate, surface wind speed, surface air pressure, incoming shortwave radiation, and incoming longwave radiation) and water table depth were inputs into ISAM to drive ISAM's emission flux response. A brief description of each site is provided in Table 1.
Site-level ISAM simulations for calibration and evaluation were performed in two steps. First, the model was spun-up for each site with the original values of the five model parameters given in Table 2; these parameters    were then calibrated in the second step. For this step, atmospheric CO 2 and CH 4 concentrations, and N deposition were set at 1981 levels. Since climate forcing and water table depth from site observations are only available for a limited number of years (Table 1), we recycled the site-specific climate data and water table depth until the SOC profile reached a steady state. Second, we calibrated five different model parameters (Table 2) by iteratively running the model over the period 1981-2010 using time-dependent observed atmospheric CO 2 and CH 4 concentrations, N deposition and climate forcing data, as well as by varying the parameter values until we found a set of parameter values that minimized the root-mean-square errors (RMSEs) between the simulated and the observed daily net CH 4 flux. We used the FORTRAN implementation of the Feasible Sequential Quadratic Programming (FFSQP; Zhou et al., 1997) numerical optimization scheme to adjust the parameter values in each iteration. In cases where the optimization scheme was not able to find a set of parameters with an RMSE less than 0.5 (e.g., in cases where observations have a high year-to-year variation), we manually adjust the parameters via the trial and error method to minimize the RMSE. Next, we evaluated the model for two independent sites using the average values of each calibrated parameter for sites. The site-level evaluation is done using the refined Willmott's index (Willmott et al., 2012) method described in supporting information Text S4.

CH 4 Emissions Across the CONUS: Input Data, Model Experiments, and Results
The coupled ISAM model was used to estimate CH 4 emissions over historical conditions , and future scenarios (2016-2100) across the CONUS region using the datasets described next.

Input Data
1. Gridded climate forcing. For the historical simulation over the period 1981-2015, we used CRUNCEP climate forcing reanalysis data for the same seven near-surface climate variables we used for model calibration described in the previous section (Harris et al., 2014). For the future simulations under scenarios RCP4.5 and RCP8.5 over the period 2016-2100, we used the CESM1BGC fully coupled model results for climate forcing data (Ren et al., 2018). The CESM1BGC future atmospheric variables have biases as compared to the CRUNCEP reanalysis, which we addressed by adjusting the mean and standard deviation of the future datasets (see details in Text S2). 2. Atmospheric CO 2 and CH 4 concentrations and N depositions. Historical inputs for atmospheric CO 2 and CH 4 concentrations were taken from Le Quéré et al. (2018) and Dlugockenky (2018), and N deposition from Lamarque et al. (2011). The projection of these variables for future scenarios RCP4.5 and RCP8.5 was obtained from the output of the CESM1-BGC coupled simulation (Ren et al., 2018). 3. Gridded land cover (LC) or PFT data. We applied LC data for the year 2000 that was reconstructed from the Historical Database of the Global Environment (HYDE 3.1) by Meiyappan and Jain (2012). LC data did not change with time in our analysis, as our study focused on natural terrestrial CH 4 emissions. We assumed that for the wetland fraction with permanent surface water coverage, the dominating LC would be marshes with strong anoxic conditions. In the model, we treated grasses as marshes, if they exist along with wetlands in a grid cell, since the presence of marsh prevents other dense root PFTs to survive. For grid cells with seasonal wetlands and/or the seasonal nonwetland fractions, any PFT along with seasonal wetlands and/or seasonal nonwetlands could exist in a grid cell. 4. Gridded soil texture and topographic slope data. We forced ISAM with the newly available global soil texture dataset GSDE (Shangguan et al., 2014), which was generated by merging various regional and national soil databases and soil maps. The topographic slope data used is from the U.S. TOPO Map (http://nationalmap.gov/ustopo/index.html). 5. Gridded surface inundation data. The fractional area of permanent wetlands and seasonal or inundated wetlands within each 0.5°× 0.5°model grid cell was prescribed using the Surface WAter Microwave Product Series dataset (SWAMPS, Schroeder et al., 2015) for the period 1992-2012. This data determined the wetland extent within each model grid cell based on the fractional surface area of the grid saturated with water (FW; Figure 3a). The fraction of nonwetlands within each model grid cell was then determined by subtracting the FW from 1.00 (Figure 3b). To better account for the role of coastal wetlands in our estimation, we used the coastal wetland extent fraction data of Schuerch et al. (2018), which we combined with wetland extent data in each grid cell along the coastline.
Within a grid cell, both wetlands (including freshwater and coastal) and nonwetlands can coexist. Since the SWAMP data does not include open water, for example, river and lake, we did not calculate CH 4 emissions from these waters. Logistic functions fit to rainfall data extending from the most recent decade through the future scenarios were used to represent freshwater wetland area under the two future scenarios for 2016-2100. To project coastal wetland extent, we first calculated the wetland fraction using the same approach as for freshwater wetlands and then subtracted from this the fraction of the coastal wetlands lost due to the sea level rise provided by Schuerch et al. (2018). According to this data, the estimated loss of coastal wetlands due to sea level rise over the period 2010-2100 was 6% for the RCP4.5 and 19% for the RCP8.5. See Text S3 for details of this method.

Experimental Design of Model Results
ISAM was used to estimate the net CH 4 emissions (EM) equation (Equations S1.28a and S1.28b) over the CONUS for a historical period and two alternative future climate scenarios: RCP4.5 and RCP8.5. For regional simulation over the CONUS, we directly applied the fraction of each PFT based on the land use data (Meiyappan & Jain, 2012). For the mangrove forest we applied mangrove mask data of Giri et al. (2011), and for the rice crop we applied the rice mask data of Portmann et al. (2010).
Before the start of these two simulations, the model was first spun-up using a similar approach as used for site-level simulation: that is, forcing the model by repeating CRUNCEP climate forcing from 1981-2015 over 750 years with atmospheric CO 2 and CH 4 concentrations fixed at 1981 levels (337.7 ppm and 1,504.9 ppb respectively) and spatially resolved N deposition also fixed at 1981 levels.
Next, the model was run for the historical period 1981-2015 with the varying historical atmospheric CO 2 , CH 4 , N deposition, and climate forcing data discussed above. For the two future scenarios, RCP4.5 and RCP8.5, the model was forced over the period 2016-2100 with spatial and temporal varying climate forcing data, atmospheric CO 2 and CH 4 concentrations, and N deposition data (see section 4.1).
Model-estimated CH 4 emissions from wetlands and nonwetlands for the current and future environmental conditions over the CONUS are reported in section 5. Wetlands consist of permanent and seasonal wetlands (hereafter wetlands). Nonwetlands can be seasonal wet soil or seasonal dry soil (hereafter nonwetland soil). The spatial distribution of CH 4 emissions and sinks for each grid cell and the regional total values for the CONUS are reported for wetlands and nonwetlands (wet and dry soil emissions combined) averaged for the period 2001-2010 (hereafter referred to as 2000s) and the change for the averaged period 2091-2100 (hereafter referred to as 2090s) relative to the 2000s under two future scenarios, RCP4.5 and RCP8.5. Model results are also compared with other existing model results of historical wetland emissions.
To evaluate ISAM performance for the CONUS, ISAM-estimated wetland CH 4 emissions are compared with a recent product of the global wetland CH 4 emission model ensemble for use in atmospheric chemical transport models (WetCHARTs, Bloom et al., 2017), which we compiled for the period 2009-2010.

Model Calibration and Evaluation
Model calibration and evaluation results for EM for five calibrated sites (Figures 2a, 2c, 2f, 2h, and 2j), and for five evaluated sites (Figures 2b, 2d, 2e, 2g, and 2i) show a large variations in estimated seasonal and annual EM across all sites, because these sites are diverse (i.e., eight wetland sites and two upland sites) but with different PFT distribution. Therefore, CH 4 production, ebullition, and aerenchyma transport are also diverse across all sites. This is also reflected by the relatively large variations in calibrated parameter value for F aere_s , r CH4:CO2 , and F ebul_s across five calibrated sites, suggesting that the processes and soil inputs (e.g., soil texture), which these parameters relate to, are diverse across calibrated sites.
Nevertheless, ISAM results are able to capture the overall observation of the seasonal variations of EM at calibrated and evaluated sites. The agreement between model results and observations for most of the sites, with the exception of one salt marsh site (US-StJ) and both upland sites, is also supported by higher (lower) refined Willmott's index and R 2 values (Table 3). The sites with high Willmott's index and R 2 have a good fit for both the absolute values and seasonal variations between model simulations and observations, while sites with low mean error (ME) show a good fit of absolute values.

Global Biogeochemical Cycles
While the model estimated EM magnitudes for most of the sites compared reasonably well with observation data, model is not able to capture peak EM for few sites. For example, model is not able to capture the peak of CH 4 emission within a year at the US-Myb site, suggesting the underestimation of the EM during the period with peak plant productivity . The plausible cause of such discrepancy is that model is missing site-level input data, including the soil texture and plant community. In addition, most of the available site-level datasets do not provide observation information for CH 4 production, oxidation, and three different transport pathways, which makes the model validation task challenging.
The CH 4 fluxes show high seasonal variation governed by the aerenchyma transport and ebullition, which contribute the most to the total CH 4 emissions for various marsh sites (Figures 2a-2e). However, it is important to note that marsh sites are showing two different types of seasonal behavior based on the modeled emissions fluxes. While the seasonal peak in Tw1 and DPW are dominated by ebullition, the seasonal peak in the other three marsh sites is dominated by aerenchymous transport. Although published literature on these sites do not report any information about the contribution of different transport pathways of the CH 4 EM, a general linkage between plants' root tissues and a relatively larger root volume supports a larger aerenchyma transport at the time of higher productivity (Bhullar et al., 2013), which is also supported by ISAM results for Myb, WPT, and OWC sites, suggesting that plant roots contain more biomass at these sites compared to Tw1 and DPW sites. Thus, yield estimations being dominated by aerenchyma transport at these sites. In contrast, ISAM results for Tw1 and DPW suggest the production of CH 4 causing a buildup of soil CH 4 bubbles during the days after peak production. As these trapped CH 4 bubbles grow in size over time, the level of the soils where CH 4 bubbles residing in also rises slowly as well, and eventually trigger the ebullition, resulting in the release of CH 4 bubbles to the atmosphere in late summer and early fall. This is the reason that the peak emissions for Tw1 and DPW sites occur not at the time of peak CH 4 production, but at a  Figure 2f). While US-LA1 does not have a large impact from the tides, as reflected by low standard deviation from the observation, the calibrated site US-StJ results show a better Willmott's Index and R 2 but higher ME due to higher coastal tide effect, which ISAM does not account for. The impacts of the tide and salt water intrusion during summer leads to lower CH 4 emissions, which model is able to capture reasonably well. However, as Figure 2f shows, model is unable to capture late-season high CH 4 fluxes at this site, suggesting that causes of the higher CH 4 emission in later winter at US-StJ site cannot be solely explained by the changes in salinity and tides. These high emissions have not been discussed in the published literature. While our understating of changes in salt water marsh cycling due to changes in biochemical and environmental variables remain limited, a recent study (Huertas et al., 2019) has discussed two additional sources of CH 4 which may lead to higher CH 4 emissions: (1) lateral transport of nearby freshwater CH 4 through riverine overflow and (2) sedimentary methanogenesis. In an earlier study, the higher winter rates were hypothesized to be caused by the release of readily metabolized organic substrates as marsh plants mature and die (Bartlett et al., 1987).
Regarding two upland sites, one calibrated (US-Ho1) and another one evaluated (US-CRT) site, both sites have a low ME value, suggesting the model is able to capture the magnitude of CH 4 emission, but the low Willmott's Index and R 2 values describing a relatively poor representation of the seasonal variation of CH 4 emission. There could be several reasons for such low performance of the model. Since the CH 4 emission/sink from the upland ecosystem is sensitive to the soil hydrology , the seasonality of water table depth is critical for CH 4 emission, which depends on soil texture and porosity. For example, soils with higher clay percentage tend to create fewer oxygen microsites (Wagner, 2017) and to increase methanogenic activity by decreasing permeability, increasing soil stratification, and anaerobic conditions. However, higher soil clay percentage also adsorbs more substrate, which prevents methanogenesis, and thus reduce the CH 4 production (Segers, 1998). The net effect of soil clay contents on methanogenesis depends on the relative magnitude of these two opposite effects. However, two upland FLUXNET sites do not provide texture and porosity profile information. Therefore, we use the global dataset (GSDE) for these variables (Table 1). In addition, the emissions are impacted by the site history and the heterogeneity of the landscape nearby, as the site could be a drained swamp, which is the case for US-CRT (Chu et al., 2014). Overall, without considering the detailed site-specific information, the model is not able to capture observed seasonal variations in the soil moisture for these two upland sites. Therefore, the estimated Willmott's index and R 2 values for these two upland sites are low.
For the calibrated rice site (Us-Twt), the seasonal variation of CH 4 emission is primarily controlled by the specific rice phenology , thus producing several emission peaks during each year

Global Biogeochemical Cycles
growing season. The model accounts for six phenological stages of rice growth (Lin et al., 2017), which help to produce seasonal emissions that are consistent with the observation, although the model is not able to reproduce peak emissions during the autumn of 2013, which is similar to the results of the nearby marsh site (US-Tw1), that a drastic increase of ebullition after the growing season can explain this finding, while other carbon fluxes and environmental variables are not able to explain this large interannual variability of CH 4 emission .
CH 4 emissions/sinks from the upland ecosystems are also sensitive to the soil hydrology. The seasonality of water table depth (which depends on soil texture and porosity) drives the seasonality of CH 4 emissions. For example, soil with high clay percentage tends to adsorb more substrate and prevent methanogenesis, and thus decrease the CH 4 production. However, FLUXNET sites do not provide texture and porosity profile information. Therefore, we applied the global dataset (GSDE) for these input variables (Table 1).

Historical Soil CH 4 Emissions Across the CONUS, 2001-2010
ISAM-estimated net emission EM over the 2000s is 13.8 TgCH 4 yr −1 (Table 4), which is the sum of wetland total flux N w = 15.3 TgCH 4 yr −1 and nonwetland total flux N nw = −1.5 TgCH 4 yr −1 across the CONUS. Marsh and salt marsh wetlands contribute 60% and 40% of total modeled wetland emissions. Wetland CH 4 production P was 15.5 TgCH 4 yr −1 ( Figure S2), a slightly higher value than N w , suggesting almost all wetland production is emitted to the atmosphere via the three emission pathways, FD (39%), A (43%), and E (18%) over the decade. If the effects of salt water intrusion and sea level rise on coastal wetlands are not included, the ISAM-estimated EM and N w for the 2000s would be 20.7 and 22.2 TgCH 4 yr −1 , or about 49% and 44% higher respectively, a proof of the importance of including these effects.
Nonwetland grid cells with higher P nw in anaerobic microsites than O nw contributed 1.4 TgCH 4 yr 1 soil emissions, whereas grid cells with higher O nw than P nw contributed −2.9 TgCH 4 yr −1 soil sinks, and so overall the nonwetland is a net sink of CH 4 over the CONUS.
ISAM calculated the spatial pattern of wetland flux N w over the CONUS region ( Figure 4a) shows higher emissions in the eastern and southeastern coastal regions, Mississippi River Delta zone, and near Lake Superior. While increasing wetland extent (Figure 3a) is the major cause of higher emissions, a higher temperature is an additional factor contributing to higher emissions, because methanogens in warm and moist environments consume soil carbon more rapidly. Grid cells with the same wetland fraction but higher temperatures show higher N w ( Figure S3). This effect is stronger when the wetland extent is 0.4 or less.
The majority of nonwetlands in the midwest, northeast, and Mississippi River delta zone shown in Figure 4d are a net source of N nw for the 2000s with A nw and FD nw each contributed approximately 50% to the total transport of these emissions from the soil column to the atmosphere. Other nonwetland regions, especially the Rocky and Appalachian mountainous, are a net sink for N nw in Figure 4d.
Both P nw and O nw are important determinants for the overall pattern of N nw . Soil moisture of the top 30 cm is the dominant factor impacting P nw and O nw and determining whether the soil is a CH 4 source or sink ( Figure S4). For regions with soil moisture higher than 50% ( Figure S4), such as the U.S. Midwest and the coastal region of Northeastern CONUS (Box (a) in Figure S5), both P nw ( Figure S6d) and O nw ( Figure S6g) are higher, but P nw exceeds O nw because of colder temperatures ( Figure S4) and nonwetland soil is a Note. The contributions from the wetlands (subscript w) and nonwetlands (subscript nw) are reported separately.Positive values represent a CH 4 emission to the atmosphere. a A nw occurs when the soils are wet as also shown in the observation study (Jackson & Armstrong, 1999).

10.1029/2019GB006251
Global Biogeochemical Cycles source. In regions with soil moisture lower than 50%, such as the large mountainous region in the western CONUS (Box (b) in Figure S5), O nw ( Figure S6g) is many times larger than P nw ( Figure S6d). Thus, nonwetland soil is a sink for atmospheric CH 4 ( Figure S4).
Of the three transport mechanisms for wetland and nonwetland combined, we find A (7.2 TgCH 4 yr −1 ) to be the largest contributor to the transport of total CH 4 (i.e., wetland + nonwetland) to the atmosphere (Table 4). Wetland regions with high A are in the southern CONUS wherevascular plants with higher aerenchyma were prevalent, such as shrubs and forests. In addition, the higher temperature also increases the plant activity and thus A. The contribution of FD to EM is the second largest of the CONUS (3.8 TgCH 4 yr −1 ) and is the dominant transport mechanism emissions for seasonal wetlands and nonwetland regions (Figure 5b). FD is the only transport mechanism for soil CH 4 uptake (i.e., negative EM) in the nonwetland western CONUS and the southeastern mountainous region, because these regions have a deeper water table depth. Transport by E (2.8 TgCH 4 yr −1 ) is the smallest and is limited to wetland regions with high temperature, especially in the southeastern coastal region of the CONUS (Figure 5c).

Future Soil CH 4 Emissions Across the CONUS, 2011-2,100
In future climate scenarios, the increases in temperature and wetland extent lead to increases in EM. By the 2090s, the expansion of the wetland extent (from 0.45 to 0.48 Mkm 2 under RCP4.5 and 0.49 Mkm 2 under RCP8.5) caused by higher precipitation and sea level rise contributes to higher P (Figures 5k and 5q). Thus, EM is higher under two future climate scenarios. ISAM estimated EM under two scenarios ranges 17.9-22.6 TgCH 4 yr −1 (RCP4.5-RCP8.5) by the 2090s, which are about 32-87% higher than the 2000s EM (Table 4). Increase in P nw (8.7-17.4 TgCH 4 yr −1 ) is much higher than the increase in P w (4.5-9.8 TgCH 4 yr −1 ) by the 2090s ( Figure S2). However, O nw also increases at higher rates (9.9-20 TgCH 4 yr 1 ), resulting in the overall increase in EM under two scenarios. N w is the major contributor to the increase in EM under both scenarios (19.4-24.7 TgCH 4 yr −1 ) by the 2090s due to wetter and warmer conditions. Note that if the effect of salt water intrusion were not accounted for, these emissions would have been increased by 23-27%, and sea level rise by 2% and 4% under two scenarios.
There is no clear spatial pattern of N w change found, except for the higher increase values in N w in the northwestern and southeastern CONUS (Figures 4b and 4c) due to both the large increase in the wetland extent (Figures 3b and 3c) in these areas.
Of the three different transport pathways (A, E, and FD), the increase in A w along northern Great lakeshore and in the southern CONUS under both scenarios is the primary contributor to the larger N w (Table 4 and Figures 5j and 5p), mainly due to an increase in biomass as a result of both the warmer climate and the CO 2 fertilization effect. On the other hand, an increase in E (or E w ) under RCP4.5 by the 2090s is much smaller (19% relative to the 2000s) than under RCP8.5 (78% relative to the 2000s) (Table 4 and Figures 5i and 5o). E is higher mainly in productive wetlands (Table 4), because of higher CH 4 concentration due to higher methanogenesis through higher litter input, and thus greater bubble formation. The larger increase in the productive land under RCP8.5 by the 2090s is the result of a larger increase in the wetland area extent, especially the increase in fractional water near the eastern coastal region and the Mississippi Delta (Figures 3b and 3c). Warmer temperature decreases the CH 4 solubility in the soil water thus increases E, but this increase in CH 4 emissions is smaller compared to the increase in N w due to the expansion of the wetland extent. The increase of FD is relatively small as compared to A under both scenarios partly due to FD nw , that is the diffusion of atmospheric CH 4 to the soils (Table 4), and the small response of FD w for both future climate scenarios by the 2090s, because of the drier topsoil and lower CH 4 production.
The nonwetland N nw contributes a small sink to EM under both scenarios (RCP4.5-RCP8.5) by the 2090s (Table 4). The increase in CH 4 sinks is larger than the CH 4 sources by the 2090s. Overall, the net CH 4 sink from the nonwetland region becomes larger under both scenarios by the 2090s, but still is less than 10% of EM.
The spatial patterns of N nw under both scenarios show a transition from CH 4 source to sink in the northeastern CONUS (Figures 4e and 4f) compare to the pattern over the 2000s, because of the reduction of the soil moisture over this region under the future climate projections. As for the aerenchyma transport pathway (i.e., A nw ), our model results show that it is almost doubled under RCP4.5 and tripled under RCP8.5 (Table 4), despite the fact that the soil moisture is decreased under future scenarios. This is because aerenchyma transport is also positively correlated with the plant productivity and leaf area index (Wania et al., 2010) (see model Equations S1.24-S1.26). ISAM results show that plant productivity (and leaf area index) increases, and soil moisture decreases under future scenarios. Overall, the effect on A nw of increase productivity outweighs that of reduced soil moisture. The FD nw from the atmosphere to the soil column is also increased, mostly due to the increase in temperature.

Comparison of ISAM Results With Other Studies
Note that we compare ISAM-estimated N w with other model results, because most of the reported results are only for N w (Table 5). Tian et al. (2012) did provide an estimate of the CH4 sink across the CONUS of 1.3 TgCH4 yr−1, slightly lower than the ISAM estimates of 1.5 TgCH4 yr−1. Our model also reflects a good match to the previous publications of the hotspot regions for CH 4 emission (Bloom et al., 2017;Potter et al., 2006;Saunois et al. 2016;Tian et al., 2010), and these include southeastern CONUS, Mississippi delta, Western coast, and the Upper Midwest. Among these regions, the southeastern CONUS contributes the largest emission due to the higher temperature. All the hotspot regions contain a certain fraction of the saline wetland except the Upper Midwest.
In particular, ISAM-estimated spatial distribution for N w is consistent with the WetCHARTs product ( Figure S7), except for the Midwest region, where ISAM N w shows the seasonal wetland area contributed substantially to N w , while WetCHARTs product shows almost negligible emission from this region. In spite of the similarity between ISAM and WetCHARTs results, there are various differences between the two studies. First, while the annual mean values for the wetland extent for the 2000s used in ISAM are lower (0.45 Mkm 2 ) than WetCHARTs product (0.47 Mkm 2 ), the SWAMPs product used in this study shows a higher seasonal wetland extent areas, particularly in the Midwest region, compared to WetCHARTs product. This difference has led to higher wetland emissions over the Midwest in ISAM case. Second, ISAM omitted anaerobic oxidation in wetlands, especially in the coastal wetlands, resulting in the lower N w . Third, ISAM uses a higher Q 10 value of 2.2 (which we had calibrated with site-level data) compared to other studies (Riley et al., 2011;Wania et al., 2010), which use a Q 10 of 2. With slightly higher Q 10 , ISAM does estimate higher N w (Table S4). For example, ISAM estimates a higher N w in the warmer eastern coastal region, but a lower N w in the inland western CONUS ( Figure S7). Fourth, a vertically resolved soil biogeochemistry component of ISAM calculates soil heterotrophic respiration and thus P for each soil layer up to 3.5 m and thus estimates additional N w from deeper soil layers. For example, ISAM-estimated N w based on first 0.3 m soil (five soil layers) were about 2.5 TgCH 4 yr −1 lesser than estimated based on 3.5-m soil (Table S4). Fifth, ISAM also accounts for the effect of salt water intrusion on methanogenesis along the coastlines. This effect alone reduced the ISAM-estimated CH 4 emissions for the 2000s by 6.9 TgCH 4 yr −1 (the difference between with and without coastal wetland effect, Table 4). Finally, ISAM calibrated r CH4 : CO2 for the freshwater wetlands is about 0.06 compared to a published range of 0.0001-0.1 and median of 0.04, which are estimated based on the multiple freshwater wetland sites (Segers, 1998). When calculations were performed with r CH4:CO2 of 0.04 for freshwater, the ISAM estimated N w turned out to be about 14.8 TgCH 4 yr −1 . This result is consistent with the mean of the WetCHARTs product value (14.8 TgCH 4 yr −1 , Table 5). However, with the coastal wetland effect, ISAM would have an estimated lower N w than WetCHARTs product value.
We also compared ISAM estimated soil CH 4 uptake from the dry or upland soils (i.e., mainly FD for the dry soils) for a Colorado State measurement site studied by Curry (2007). ISAM estimated CH 4 uptake for this site is 0.3 gCH 4 m −2 yr −1 , which exactly matched with the Curry (2007) observation-based estimate. ISAM estimated spatial pattern of FD for the upland (Figure 5b) across the CONUS shows the highest uptake in the mountainous region and southwest, where soil moisture is low and soil pores are aerated well, which also matches well with the conclusions of the previous studies (Curry, 2007;Tian et al., 2016). Soil emperature is another factor determining the dry soil CH 4 uptake rate. For example, compared to the colder Midwest region, higher soil uptake is found in the upland soils between northeastern and southeastern areas due to their warmer temperature increasing methanotroph activity in dry soils and thus increasing the soil sink.
There is no other modeling study available in the literature on future CH 4 scenario analysis for the CONUS that we can compare with our model results directly. Zhang et al. (2017), however, performed simulation using a simple algorithm to calculate CH 4 emission from wetlands at a global total case and found an increase in global total CH 4 emissions by 20% under RCP4.5 and 70% under RCP8.5. Our modeling study also shows similar CH 4 emissions increasing trends under the two scenarios across the CONUS (Table 4).

Conclusions
A newly developed vertically resolved soil CH 4 dynamics module has been coupled to a land surface model, ISAM. After calibrating and evaluating ISAM using site-level soil CH 4 emissions data, the model is applied to calculate CH 4 emissions and sinks from wetlands and nonwetland soils across the CONUS over historical time and under two future scenarios.
ISAM-estimated EM across the CONUS for the 2000s is 13.8 Tg CH 4 yr −1 . The results show that wetlands, including seasonal wetlands, dominate soil CH 4 emissions (N w = 15.3 Tg CH 4 yr −1 ), mainly in the eastern and southeastern coastal regions, Mississippi River delta zone and the regions near Lake Superior.
Model-estimated future EM under RCP4.5 and RCP8.5 scenarios by the 2090s were about 30% and 64% higher than over the 2000s, mostly due to an increase in temperature and an increase in seasonal wetland extent due to higher precipitation. The overall increase in CH 4 emissions under two future scenarios suggests positive feedback in the climate system with increases in radiative forcing leading to higher CH 4 emissions.
In contrast to wetland emissions, nonwetland soils were estimated to be a small CH 4 sink (N nw = −1.5 Tg CH 4 yr −1 ) for the 2000s, and this sink increased by 17% and 58% by the 2090s under RCP4.5 and RCP8.5 scenarios. Overall, the magnitude of N nw is less than 10% of EM.
We find that the strength of different CH 4 transport mechanisms differs across wetlands and nonwetland soils, and they are sensitive to environmental conditions. Diffusion dominates in nonwetland soils with low soil moisture, whereas ebullition is more important in wetlands, and aerenchyma transport is important in both wetlands and nonwetlands. Under future scenarios, ISAM simulations show higher CH 4 emissions, which are mainly transported by ebullition due to the increase in wetland extent and by increased aerenchyma transport caused by aerenchyma growth due to the CO 2 fertilization effect.
This is the first study that estimates CH 4 emissions and sinks not only from wetlands (including freshwater and coastal wetlands) but also from the nonwetland soils (including wet and dry soils) under future climate scenarios across the CONUS. Applying this approach to soils worldwide, could provide the information needed to make quantitative estimates of the strength of the climate feedback between climate change and soil CH 4 emissions. To do this, there is a need for additional CH 4 emission modeling analyses from wetlands and nonwetland soils for other regions ranging from the arctic to the tropics.

Data Availability Statement
The estimates of annual mean wetlands and nonwetlands CH 4 emissions for the CONUS region using ISAM under three cases: (1) 2001-2010, (2) 2091-2,100 under RCP4.5, and (3) 2091-2,100 under RCP8.5 scenarios can be downloaded for free from the ISAM website (http://climate.atmos.uiuc.edu/Methane). The data files contain the following methane fluxes for both wetlands and nonwetlands soils: (1) CH 4 emission, (2) CH 4 production, (3) soil CH 4 oxidation, (4) CH 4 emissions through aerenchyma transport, (5) CH 4 emissions through ebullition, and (6) CH 4 emissions/sink through diffusion. The wetland extent described by the fractional area of surface water is also stored. We are sharing the data at 30 m spatial resolution at 0.5 × 0.5°resolution. All the data are stored in NetCDF files. The data can be accessed from the ISAM website (http:// climate.atmos.uiuc.edu/Methane_ShuEtAl/).