Jointly resolving the sign and magnitude of the CO 2 and CH 4 combined carbon-climate (C-climate) feedback is critical for quantifying the C-climate potential of wetland ecosystems in the Earth system. A modeling study shows that although the amount of CH 4 emissions was small relative to CO 2 uptake (∼10% of the amount of carbon uptakes on a molar basis), the CH 4 emission (in CO 2 equivalent GWP almost entirely offset the net CO 2 uptake over the 21st century under the no-policy climate change scenario, X. Zhu et al., 2013). In response to a future warming climate, both the CH 4 source and the CO 2 sink strengthened, but the combined CH 4 and CO 2 GWP response did not show a significant trend. In a similar investigation, Lawrence et al. (2015) use the Community Land Model to show that drier soil conditions accelerate organic matter decomposition with increased CO 2 emissions, but strongly suppress growth in CH 4 emissions, and the land-to-atmosphere CO 2 -equivalent C flux (namely the total CO 2 and CH 4 flux weighed by the radiative forcing of CH 4 relative to CO 2 ) decreases by more than 50%. While these studies provide invaluable insights on the potential role of CH 4 and CO 2 fluxes, model structures overwhelmingly rely on empirical parametrizations of CH 4 and CO 2 sensitivities to climate and soil states, which are ultimately key for accurately predicting CH 4 and CO 2 responses to a changing climate. Quantitative knowledge on the sign-and corresponding magnitude-of the C-climate feedback potential is critical for understanding the role of high latitude ecosystems in amplifying or dampening the atmospheric greenhouse gas burden in the coming decades.

Merging in situ measurements of both CO 2 and CH 4 with dynamic ecosystem models provides a unique opportunity for resolving the potential sign and magnitude of the C-climate feedbacks of wetland ecosystems. Eddy covariance (EC) flux towers provide direct and near-continuous ecosystem-scale CO 2 and CH 4 flux measurements without disturbing the soil or vegetation (Aubinet et al., 2012;Knox et al., 2019). The rapidly growing number of CO 2 and CH 4 flux tower measurements presents new opportunities for estimating fluxes at the ecosystem, regional, and global scale; furthermore, the integration of these measurements into biosphere model structures and parameters can amount to a critical step toward determining the joint sensitivities of CO 2 and CH 4 fluxes to climate variability. In turn, terrestrial biosphere models provide the necessary means to mechanistically represent the ecosystem responses to climate variability. Bayesian model-data fusion (MDF) using EC data has become a valuable tool for optimizing model C states, fluxes and process parameters and provides quantitative insights on terrestrial C cycling (Bloom & Williams, 2015;Famiglietti et al., 2020;Keenan et al., 2013;Smallman et al., 2017;Stettz et al., 2021;Williams et al., 2005;Yang et al., 2021). Model-data integration efforts also provide the C cycle mechanistic insights necessary for explicitly resolving the sensitivity of C cycling to climate and anthropogenic forcings (Bloom et al., 2020a(Bloom et al., , 2020b;;Stettz et al., 2021;Yin et al., 2020).

MA ET AL.

10.1029/2022GB007524 3 of 18

In this study we investigate the potential role and magnitude of wetland CO 2 and CH 4 flux sensitivities to climate, and the sign (±) of the combined CH 4 and CO 2 climate sensitivity in terms of their combined GWP (henceforth CO 2 eC , which denotes the GWP-weighed sum of CH 4 and CO 2 fluxes). To investigate the CO 2 , CH 4 and CO 2 eC flux responses to climate variability, we assimilate CO 2 and CH 4 EC flux tower data into a mechanistic C cycle model using the CARDAMOM Bayesian data-model integration framework (Bloom et al., 2016;Quetin et al., 2020) to constrain model parameters and states at seven EC sites in Alaska. Based on the observation-informed model analyses at each site, we probabilistically evaluate two hypotheses on the combined CO 2 and CH 4 responses to changes in climate: H1: Climate change induces a positive CO 2 eC response: the combined response of CO 2 and CH 4 fluxes to a change in climate-weighed by the GWP-is greater than 0. H2: Climate change induces a negative CO 2 eC response: the combined response of CO 2 and CH 4 fluxes to a change in climate-weighed by the GWP-is less than 0.

These hypotheses imply that wetland ecosystem responses to climate changes amount to positive (H1) or negative (H2) C-climate feedback. The two hypotheses are illustrated in Figure 1: positive land-to-atmosphere flux responses for both CH 4 and CO 2 to climate change amount to an overall positive CO 2 eC response to climate (positive CC feedback), and negative land-to-atmosphere flux responses for both CH 4 and CO 2 amount to an overall negative response to climate (negative CC feedback). The combination of a negative CH 4 flux response and a positive CO 2 response, or vice versa, could amount to either a positive or negative CO 2 eC response depending upon the relative contributions of CH 4 and CO 2 weighted by their GWP. Ultimately, an increase or decrease in the combined CO 2 eC response of CO 2 and CH 4 to climate amounts to positive or negative C-climate feedback.

To probabilistically distinguish between the two proposed hypotheses, we quantify the EC-constrained CARDAMOM CO 2 and CH 4 flux sensitivities to individual climate variables, namely temperature and precipitation, and to hypothetical shifts in present-day climate. In Section 2, we introduce the CARDAMOM MDF approach, associated forcing data and observational constraints, and the finite difference approach for resolving CH 4 and CO 2 flux sensitivities to climate. In Section 3, we present and discuss our results, and our conclusion is in Section 4.

To test the hypotheses presented in Section 1, we use the CARDAMOM MDF framework to constrain terrestrial C cycle parameters and initial model states based on in situ observations of CH 4 and CO 2 fluxes from EC towers. We introduce the structure of the process-based model used in CARDAMOM in Section 2.1; we describe the site level forcings and observational constraints in Section 2.2; we describe Bayesian MDF framework in Section 2.3; we then use the data-constrained model to determine the sensitivity of CH 4 and CO 2 fluxes to the climate in Section 2.4.

CARDAMOM is a model-data fusion (MDF) system that constrains model parameters and initial states in a process-oriented model-the Data Assimilation Linked Ecosystem Carbon model (DALEC)-and a Bayesian model-data integration algorithm using adaptive Metropolis-Hastings Markov Chain Monte Carlo (MHMCMC) (Haario et al., 2001). The DALEC model structure is uniquely suited to model-data integration, where unknown model parameters and initial states are statistically optimized to ensure minimal mismatches between ecosystem

, where z = 28 is a scalar to account for the 100-year scale GWP of CH 4 (Myhre et al., 2013). The slope of the equal feedback line is 1/28, and the slope of the zero C-climate feedback line is -1/28. The unit for both gases is gram molecular. The solid purple line denotes the "zero-feedback" line (Δ CO 2 eC-CO 2 = -Δ CO 2 eC-CH 4 ), where the area above and below the solid line represents the CH 4 and CO 2 flux responses to climate as resulting in a positive (hypothesis one; H1) or negative (hypothesis two; H2) C-climate feedback. The dashed gray line represents an "equal CH 4 and CO 2 feedback line," where

10.1029/2022GB007524

4 of 18 observations and corresponding model states. In contrast to more complex process-based terrestrial biosphere models, the relatively parsimonious DALEC C cycle model structures employed in MDF analyses (Bloom & Williams, 2015;Famiglietti et al., 2020;Norton et al., 2023;Quetin et al., 2020;Williams et al., 2005;Yang et al., 2021) allows for joint optimization of highly uncertain process parameters (such as initial carbon and water states, photosynthetic parameters, C allocation, turnover rates and their dependencies on soil moisture) based on the available observations of carbon and water variables. While typically incorporation of complexity in models is beneficial for process representation, Famiglietti et al. (2020) show that model complexity advantages are fundamentally limited when insufficient data are available to support parameter inference.

The mechanistic knowledge gained from observations can then be used to diagnose the climatic sensitivity of land-atmosphere C fluxes and associated process controls. The DALEC model-embedded CARDAMOM MDF system has been applied to a range of spatial scales with a suite of EC and satellite data sets to (a) represent C cycles, (b) optimize critical parameters and states, and (c) infer C cycle-associated climate sensitivities (Bloom & Williams, 2015;Famiglietti et al., 2020;Norton et al., 2023;Quetin et al., 2020;Williams et al., 2005;Yang et al., 2021).

The Data-Assimilation linked Ecosystem Carbon model 2a (DALEC2a) consists of a mechanistic representation of carbon and soil water states and fluxes (Bloom et al., 2020a(Bloom et al., , 2020b;;Williams et al., 2005). The six C pools-C content of foliage, labile C in plants, woody stem and coarse root, fine roots, litter, and soil organic matter-link processes such as photosynthesis, autotrophic respiration, phenology, allocation, turnover rates, fire-removed C (Bloom & Williams, 2015;Famiglietti et al., 2020;Quetin et al., 2020;Williams et al., 2005;Yang et al., 2021). The hydrological balance is defined as the sum of precipitation inputs (P) and evapotranspiration (ET) and runoff (R) outputs. In turn, the plant-available H 2 O limits gross primary productivity through the conservation of the inherent water-use efficiency (Beer et al., 2009), where ET is calculated as a function of gross primary production (GPP) and atmospheric vapor pressure deficit (Bloom et al., 2020a(Bloom et al., , 2020b)). Another feature of DALEC is that parameterizations per pixel/site are not based on plant functional types (PFTs). Rather, each site has its own PFT characterized by parameters that are constrained by observations. For the sake of brevity, we refer the detailed descriptions of DALEC2a by Williams et al. (2005), Bloom et al. (2020aBloom et al. ( , 2020b)), Yang et al. (2021), and references therein.

Here, we extend the DALEC model structure to include a moisture and temperature-sensitive Joint CH 4 + CO 2 Respiration (JCR) scheme (Figure 2) and henceforth refer to the extended model as DALEC-JCR. We illustrate 2020), Bloom et al. (2016). P is precipitation, ET is evapotranspiration, R is runoff, GPP is gross primary production of CO 2 , Ra is autotrophic respiration, and Rh is heterotrophic respiration. Rh CO 2 is simulated from anaerobic and aerobic soils, while Rh CH 4 emission is calculated from anaerobic soil only. The volumetric fraction of anaerobic respiration is determined by soil moisture and site-specific parameters constrained by observations. The shape of the soil moisture-respiration response curve varies across sites and is determined by site-specific observations (Figure S2 in Supporting Information S1).

MA ET AL.

10.1029/2022GB007524

5 of 18 the main features of the JCR module in the main text and explain the details in Supporting Information S1 (Text S1 and S2).

A particular challenge in Land Surface Models is the vast uncertainty in the representations of soil moisture scalar and parameterizations used to resolve respiration fluxes (Exbrayat et al., 2013). Evidence from experimental and modeling studies points to a better performance using unimodal respiration-soil moisture response curve, rather than a linear relationship (Cox, 2001;Exbrayat et al., 2013;Sitch et al., 2003). A less well-characterized mechanism in models is how soil properties and vegetation types change the shape of optimum soil moisture response curves, which is found to be important for predicting the response of soil carbon to future climate scenarios (Moyano et al., 2012).

To address this challenge within a MDF framework, DALEC-JCR resolves site-specific data-constrained parameters to characterize the shape of the soil moisture-respiration curve (Text S2 in Supporting Information S1). While CH 4 production occurs as long as anaerobic conditions occur at a microscopic scale, water is not evenly distributed within the spatial domain of the target site, not least due to topographic heterogeneity, along with variability in other edaphic and biotic factors. Therefore, we expect anaerobic respiration to increase gradually-rather than step-wise-with increasing soil moisture. To account for this in DALEC-JCR, we diagnostically determine the dynamics of the aerobic fraction of the soil column based on an empirical continuous functional relationship between soil moisture and aerobic soil fraction (Text S1 in Supporting Information S1).

In the DALEC JCR module, CO 2 is respired both anaerobically and aerobically-from oxic and anoxic soil conditions, respectively-while CH 4 is only respired anaerobically. The JCR module calculates aerobic and anaerobic respiration from the soil column separately, both as a function of C turnover rate, soil temperature, soil moisture, and volumetric fraction of aerobic/anaerobic soil:

where C litter and C SOM are the litter and soil organic matter C pool size; k litter and k SOM are the basal decay rate of C pools, their value are independent of respiration type; f W is the soil moisture scalar on aerobic respiration rate; f Wc is the soil moisture scalar when the soil is saturated; the anaerobic soil is saturated (soil moisture always equals one) and thus the soil moisture scalar is given a constant value of f Wc in Equation 2, the value of f Wc is determined at each individual site and constrained by observations in our data-model fusion framework; ae is the aerobic fraction of the vertical soil column; 1 -ae is the anaerobic fraction of the soil. Equations S1-S4 in Supporting Information S1 describe how f W and ae are calculated; f T is the soil temperature scalar on respiration rate:

where Q 10 is the factor by which respiration rate increases with a 10°C increase in temperature. T is the mean air temperature of the current time step, T mean is the multi-year mean air temperature at the region.

The heterotrophic respiration terms in the form of CO 2 ( ℎCO 2 ) and CH 4 ( ℎCH 4 ) are then calculated as:

where CH 4 is the fraction of CH 4 in anaerobic respiration:

where CH 4 is the potential ratio of anaerobically mineralized C released as CH 4 ; 10 CH 4 is the factor by which CH 4 production rate increases with a 10°C increase in temperature, on top of the temperature sensitivity encountered in Equations 1 and 2 (f T ); The reason we put a 10 CH 4 on top of the general respiration temperature sensitivity term here is that studies have found higher temperature sensitivity in methane production than CO 2 respiration across microbial to ecosystem scales (Yvon-Durocher et al., 2014). T is the mean air temperature of the current time step, T mean is the multi-years mean air temperature at the region.

MA ET AL.

6 of 18

The parsimonious DALEC-JCR model structure (Figure 2) does not represent a range of soil hydrological and physical states, such as soil energy dynamics, permafrost thaw, snow cover, snowmelt infiltration, and soil consumption of CH 4 . To determine if the DALEC-JCR model, with data-informed parameters, can simulate realistic observable CO 2 and CH 4 fluxes (i.e., NEE and CH 4 ), we compare their temporal variability to fluxes from a more structurally complex and comprehensive ecosystem model (Terrestrial ECOsystem 2.0, TECO 2.0, Weng & Luo, 2008), at a well-calibrated bog site in northern Minnesota (Ma et al., 2017). TECO 2.0 is a process-based ecosystem model that has six submodules: canopy, vegetation dynamics, soil water, soil energy (including snow cover, thermal dynamics, and frozen depth), soil carbon/nitrogen, and aerobic/anaerobic respirations (CO 2 and CH 4 ). TECO 2.0 explicitly represents the transient and vertical dynamics of soil moisture, soil temperature, liquid fraction, frozen depth, CH 4 production, oxidation, and major transport pathways (diffusion, ebullition, and plant-mediated-transport); a bucket model is used to estimate water table level, determined by soil moisture change. The aerobic and anaerobic zones are separated at the water table (Granberg et al., 1999;Wania et al., 2010;Q. Zhu et al., 2014). A detailed description of TECO is available in Weng and Luo (2008) and Ma et al. (2017Ma et al. ( , 2023)). The DALEC-JCR emulation of TECO 2.0 is described in the manuscript supplement (Text S3 in Supporting Information S1); We find DALEC-JCR performs favorably against the TECO 2.0 variability in CO 2 and CH 4 fluxes ( 2CO 2 = 0.92 ; 2 CH 4 = 0.8 , Figure S12 in Supporting Information S1), which corroborates the overall DALEC-JCR model structure.

We configure DALEC-JCR forcings and observation constraints at seven high-latitude eddy tower sites (Table 1 and Figure S3 in Supporting Information S1). Five sites are low-lying vegetation wet tundra (sedge grasses, moss, lichen, US-Bes, US-Beo, US-Brw, US-Atq, US-Ivo), one site is fen covered by shrubs, grasses, and moss (US-ICs), and one site is boreal bog covered by trees, shrubs, grasses, and moss (US-BZB). Following Bloom and Williams (2015), we provide a continuous monthly meteorological forcing data set for DALEC-JCR using 0.5° resolution meteorological forcing (namely monthly temperature, precipitation, global incident shortwave radiation, vapor pressure deficit, and burned area) obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis Interim (ERA-interim, Berrisford et al., 2011) to drive DALEC-JCR model at each site. The model was run at a monthly timestep. Basic information about the EC sites is listed in Table 1. We use monthly averaged net CH 4 and CO 2 fluxes to constrain DALEC-JCR (see Section 2.3); months with less than 10 days of observations are excluded from our analysis.

We optimize JCR module parameters (Table S1 in Supporting Information S1)-along with other DALEC time-invariant model parameters and initial C and H 2 O states (in total 33 time-invariant parameters and 8 initial states)-within the CARDAMOM (CARbon Data Model fraMework) MDF system. The CARDAMOM model parameters are independently optimized at each one of the seven flux tower sites. The CARDAMOM framework uses Bayes' theorem to optimize the posterior probability of initial states and time-invariant process parameters (y), given observations O, p(y|O), as follows:

where p(y) is the prior probability distribution of y, and p(O|y) is proportional to the likelihood of y given O, L(y|O). At each tower site, the observation vector O consists

AmeriFlux site ID (references) Latitude/Longitude/Wetland characteristics Vegetation type ERA-interim temperature ERA-interim precipitation MAT (2001-2016, °C) IAV (2001-2016, °C) Trend (1970-2016, °C/yr) a MAP (2001-2016, mm) IAV (2001-2016, mm) Trend (1970-2016, mm/yr) a US-ICs (Euskirchen et al., 2017) 68.6/-149.3/wet tundra fen (alkaline) Wet sedge and dwarf shrubs -9.4 0.83 0.039 426 75 1.47 US-BZB (Euskirchen et al., 2014) 64.7/-148.3/thermalkarst bog (acidic) Lowland boreal forests -1.5 0.94 0.039 413 49 1.55 US-Bes (Davidson et al., 2016) 71.3/-156.6/wet tundra (acidic) Sedges, grasses, moss, lichens -9.5 0.96 0.038 244 34 1.47 US-Beo (Davidson et al., 2016) 71.3/-156.6/wet tundra (acidic) Sedges, grasses, moss, lichens -9.5 0.96 0.038 244 34 1.47 US-Brw (Davidson et al., 2016) 71.3/-156.6/wet tundra (acidic) Sedges, grasses, moss, lichens .

Table 1 Mean and Inter-Annual Variation (IAV) During 2001-2016, and the Trend of Temperature and Precipitation During 1970-2016 at the Eddy Covariance Sites MA ET AL. 10.1029/2022GB007524 7 of 18 of measurements including the NEE, and CH 4 . Assuming errors are uncorrelated, the overall likelihood of y given O can be expressed as

For NEE and CH 4 we derive the corresponding likelihood function L * (i.e., LNEE, and CH 4 ) as follows:

where o i and m i (y) correspond to the ith observation and the corresponding modeled quantity derived from control vector y, respectively; σ i accounts for the combined errors from the DALEC model structure, forcing drivers, and observations. Following (Bloom et al., 2020a(Bloom et al., , 2020b)), we retrieve the distribution of p(y|O) at each tower site by running four adaptive MHMCMC chains for 10 8 iterations. We used the Gelman-Rubin statistic (Gelman & Rubin, 1992) to check the convergence of sampling chains. The first half of the accepted parameters was discarded as the burn-in period, and the second half was used for posterior analysis. We evaluate modeled CO 2 and CH 4 fluxes against observed ones using the RMSE normalize by the standard deviation (NRMSE), and coefficient of determination (R 2 ).

To characterize the potential role of climate variability on CH 4 and CO 2 fluxes, we focus on warming and precipitation from a regional perspective to diagnostically resolve CO 2 /CH 4 sensitivities and the sign of their combined feedback (Figure 1 2). ERA Interim is used due to its consistency with model driver data, and is broadly consistent with gap-filled FLUXNET CH 4 meteorology data (Vuichard & Papale, 2015). We use the CRU-NCEP reanalysis data to quantify the scenario S5 climate variable trends, the CRU-NCEP data set covers a substantially longer time period . The CRU-NCEP centennial temperature variability (shown in Figure S4 in Supporting Information S1) reveals a steep increase of temperature beginning in 1970 (in contrast to the 1901-1970 time period); we therefore chose to use 1970-2016 trend to characterize scenario S5 climate trends (Table 2).

To quantitatively resolve H1 and H2, we denote the modeled C fluxes, F, as a function of time-dependent meteorological drivers (M) and time-invariant model parameters and initial states (y):

where the DALEC-JCR() operator denotes the DALEC-JCR model simulation driven by forcing M and parameters y.

Climate scenarios Short name Purpose Description Climate perturbation S1 Warmer Reflect year-to-year changes, based on the historical Inter Annual Variation (IAV) of 2001-2016 A hypothetical increase/decrease in temperatures relative to nominal IAV Standard deviation of mean annual temperature (MAT) (+ temp , unit °C/yr) S2 Wetter A hypothetical increase/decrease in precipitation relative to nominal IAV Standard deviation of mean annual precipitation (MAP) (+ prec , unit mm/yr) S3 Warmer + Wetter A hypothetical increase in temperatures and precipitation relative to nominal IAV. To test sensitivity of CO 2 and CH 4 to warming and wetting Standard deviation of MAT and MAP ( temp + prec ) S4 Warmer + Drier A hypothetical increase in temperatures and decrease in precipitation relative to nominal IAV. To test sensitivity of CO 2 and CH 4 to warming and drying Standard deviation of MAT and MAP ( temp -prec ) S5 a 1971-2016 trend (warmer and wetter at all sites) General trend A hypothetical increase in temperatures and precipitation following the 50 years trend. To test sensitivity of CO in Supporting Information S1. We take two steps to evaluate the responses of CO 2 and CH 4 fluxes to decadal-scale climate change:

We calculate the local gradient of F response curve against environment variables (M) by perturbing M with an arbitrarily small change to M, ∆m (e.g., ∆m = 1e -5 °C for temperature dimension, ∆m = 1e -8 mm/day for precipitation dimension), where:

To calculate the uncertainty of dF/dM, we use 4,000 samples of y accepted from the CARDAMOM Bayesian inversion (Section 2.3) and for each one we derive the corresponding dF/dM values from Equation 11. We derive ∕ (temperature) and ∕ (precipitation) using Equation 11. The combined effects of temperature and precipitation, ∕ can be expressed as:

where I represents the interactions between P and T anomalies on C fluxes. We tested the interactions between temperature and precipitation effects (shown in supplementary Text S4 in Supporting Information S1) and found it is more than an order of magnitude smaller than the individual effects; we therefore approximate combined temperature and precipitation effects as

We estimate the C flux response (ΔF) to a change in climate using the following equation:

where ΔM represents a climate anomaly (Table 2, either the 2001-2016 inter-annual variability (IAV) or the 1970-present trends of temperature and precipitation). The flux sensitivity to climate derived here represent the contemporary integrated sensitivity to climate variability, with the assumption that the integrated responses of C fluxes ∆F to climate changes ∆m can be linearly approximated as ∆ = ∆ ∕ . We discuss this assumption and its limitations in Section 3.3.

Using Equation 14, we calculate individual C flux responses to individual climate variables. We then estimate the joint CO 2 eC effect of CO 2 and CH 4 to temperature changes, ΔT, as follows:

where Δ CO ( ) denotes the CO 2 equivalent C flux response, and z is the GWP of CH 4

Similarly, we derive the joint CO 2 eC CO 2 and CH 4 response to precipitation changes as follows:

For Scenarios 1-4 (S1-S4), we derive ΔP and ΔT as a 1 perturbation to mean 2001-2016 precipitation, where 1-represents the 2001-2016 IAV of annually averaged P and T, respectively. For Scenario 5 ΔP and ΔT are the 1970-present trend (τ) of temperature and precipitation. The joint effects of CO 2 and CH 4 in a warmer and wetter scenario is then equal to:

and a warmer and drier scenario is MA ET AL.

The seasonal variability and IAV of EC net CO 2 and CH 4 at all sites are captured by DALEC-JCR: the normalized RMSE (NRMSE) ranges are 0.15-0.0.37 (CO 2 ) and 0.33-0.9 (CH 4 ), R 2 ranges 0.90-0.97 (CO 2 ) and 0.30-0.90 (CH 4 ) (Figure 3, Figures S5-S11 in Supporting Information S1). R 2 of modeled CO 2 and CH 4 fluxes are higher at sites with over 1 year of observations. The model captures cold season (September to May) CH 4 emissions, which is found to dominate the Arctic tundra CH 4 budget but is not currently well simulated in most terrestrial biosphere models (Zona et al., 2016). Modeled cold season CH 4 emissions account for 32%-66% of the annual budget across seven Alaskan sites, comparable to the estimations based on in situ observations from the same time periods (37%-64%, US-Bes, US-Beo, US-Brw, US-Atq, and US-Ivo, Zona et al., 2016): Notably, the model captures the sudden decrease of CH 4 at the end of the growing season and autumn (August and September) due to a decrease in soil moisture and thus lower anaerobic soil fraction at US-Beo, US-Bes, US-Brw (Figures S7-S9 in Supporting Information S1). In situ observations at those sites indicate lowest water table depth in the same period (Zona et al., 2016). Both the model and in situ measurements show a 1-month-delayed drop in soil moisture (August) in response to decreased precipitation (June-July), which could be due to the replenishment from thawing permafrost (Figures S7-S9 in Supporting Information S1).

The temperature sensitivity of heterotrophic CO 2 respiration (Q 10 Rh) ranges from 1.3 to 1.8 (mean of posterior distribution from the seven sites), and Q 10 CH 4 ranges from 1.1 to 2.2, and the range of CH 4 /CO 2 potential (rCH 4 ) is 0.1-0.3. These posterior parameter estimates generally fall within the empirical ranges we found in the literature (Riley et al., 2011) and are consistent across sites in similar latitudes and vegetation communities (Figure S13 in Supporting Information S1).

The regression between aerobic respiration and soil moisture follows a unimodal curve (the C2 response curve in Figure S2 in Supporting Information S1) at US-ICs, and a logarithmic curve (Figure S2 C1 in Supporting Information S1) at all the other sites. As shown in Figure S2 in Supporting Information S1, a unimodal response curve indicates that soil saturation suppresses aerobic CO 2 respiration, a logarithmic curve indicates no suppression on aerobic CO 2 respiration in saturated soils. At US-ICs, we find that aerobic CO 2 respiration peaks when soil moisture is at 70%, and saturated soils suppress the heterotrophic respiration rates by 60%. Notably, US-ICs is a fen ecosystem, characterized as a peat-forming wetland relying on groundwater input, with low to moderate soil alkalinity (pH = 5.5-6.9), while all the other sites (bog and wet tundra) have acidic soil (pH = 3.3-5.5) (Bourbonniere, 2009;Clymo et al., 1984). Our finding potentially indicates that redox potential (as a result of different pH) is a key determinant of how soil moisture controls aerobic respiration, as found in numerous studies (Bhanja & Wang, 2020;Brewer et al., 2018;Fiedler et al., 2007;Porter et al., 2004).

We find that warming and wetting alone increase CH 4 emissions across all sites. Temperature induces a +0.07 to +0.36 gCH 4 /m 2 /yr/σ T and precipitation induces a +0.45 to +1.69 g CH 4 /m 2 /yr/σ P response across all sites, where σ T and σ P represent a 1 climate perturbation in temperature and precipitation, respectively (the ranges encompass median values from the seven sites). Warming accelerates net CO 2 loss at all sites (+1.7 to +27.6 gCO 2 /m 2 / yr/σ T ), while precipitation can change net CO 2 in both directions (-203.2 to +18.3 gCO 2 /m 2 /yr/σ P ). We find warming alone increased CO 2 eC at all seven high-latitude sites: CH 4 component of CO 2 eC ( CO 2 eC-CH 4 ) plays a larger role than the CO 2 component of CO 2 eC ( CO 2 eC-CO 2 ) at two of the sites (US-BZB, a boreal wetland site near Fairbanks, Alaska, and US-BES, a tundra site near Utqiagvik), and CH 4 component of CO 2 eC players a smaller role than the CO 2 component at the other sites (Figure 4a, Figure S14a in Supporting Information S1). We find MA ET AL. 10.1029/2022GB007524 10 of 18 that the impacts of a 1 change in precipitation are on average five-fold larger than 1 change in temperature for CH 4 (Figures 4a and 4b).

In terms of the combined CH 4 + CO 2 responses ( CO 2 e-C ), a 1 change in precipitation induces a larger response than a 1 change in temperature (CO 2 equivalent responses are -168.3 to +65.6 gCO 2 eC/m 2 /yr/σ T and +8.3 to +31.4 gCO 2 eC/m 2 /yr/σ P , respectively) at all sites. Although net CH 4 fluxes are at least 10 times smaller than CO 2 , the impact of a 1 increase in temperature and precipitation (Scenario 3) on Δ CO 2 eC-CH 4 is larger than

Δ CO 2 eC-CO 2 at four sites (Δ CO 2 eC-CH 4 /Δ CO 2 eC-CO 2 = 1.2-2.8 fold, range of medians). However, at US-BZB, US-ICs, and US-Brw, we find Δ CO 2 eC-CH 4 < Δ CO 2 eC-CO 2 (Δ CO 2 eC-CH4 /Δ CO 2 eC-CO 2 = 0.22-0.38 fold, range of medians) (Figure 5a). MA ET AL. 10.1029/2022GB007524 11 of 18

We found that a warmer and wetter climate induces a positive CO 2 eC response (Δ CO 2 eC = 17.7-88.2 gCO 2 eC/ m2/yr) in wet tundra sites but a negative CO 2 eC in the fen and bog sites (Δ CO 2 eC = 15.8-159.7 gCO 2 eC/m 2 /yr), where heterotrophic respiration dramatically decreases due to increased soil moisture and is accompanied by increased photosynthesis (Figure S14 in Supporting Information S1). The bog site (US-BZB) is the only site with boreal forest, thus the ecosystem has a larger leaf biomass and litter pool to facilitate a quick response of C decomposition and uptake. The fen site (US-ICs) is also unique from the other wet tundra sites due to the presence of shrub species (larger plant biomass) and alkaline soil (resulting a different redox potential), which may drive a stronger response in both heterotrophic respiration and photosynthesis. de Vrese et al. (2021) predicted weak soil CO 2 respiration in the wet months of the year, which led to low soil CH 4 fluxes in permafrost regions under Shared Socioeconomic Pathway 5 and the Representative Concentration Pathway 8.5. We predict the same directional changes in soil CO 2 respiration due to increased soil moisture, but increased CH 4 due to warmer temperature and larger anaerobic fraction in the soil column. Shu et al. (2020) predict CH 4 increase by 30% and 64% at the See Table 2 for reference to Scenario 1 and Scenario 2. Each color represents one single site and the contour lines are two-dimensional kernel density derived from 4,000 random samples from posterior estimate constrained by eddy covariance net CO 2 and CH 4 data (50%, 75%, and 95% distributions) derived from CARDAMOM posterior states and parameters. Refer to Figure 1 for descriptions of lines and intersections.

Figure 5. Responses of wetland CO 2 and CH 4 fluxes to 1 changes in temperature and precipitation at high-latitude sites: (a) Scenario 3, warmer and wetter (+ temp , + prec ); and (b) Scenario 4, warmer and drier (+ temp ,prec ). See Table 2 for reference to Scenario 3 and Scenario 4. Each color represents one single site and the contour lines are two-dimensional kernel density derived from 4,000 random samples from posterior estimate constrained by eddy covariance net CO 2 and CH 4 data (50%, 75%, and 95% distributions) derived from CARDAMOM posterior states and parameters. Refer to Figure 1 for descriptions of lines and intersections.

12 of 18 end of century under RCP4.5 and RCP8.5 in CONUS wetlands, which are the same directional changes of CH 4 emission as predicted in our study. Lawrence et al. (2015) predicted that 10%-drier soil conditions across high-latitude regions will accelerate net CO 2 emissions but strongly suppress growth in CH 4 emissions, resulting in a negative C-climate feedback and a 50% lower CO 2 eC . While our results agree that drier soil (1 decrease in precipitation) decreases CH 4 emission (with 95% confidence), it affects CO 2 exchange in two directions depending on ecosystem types (Figures 4b and 5b). Specifically, we find that the drier soil conditions result in a negative C-climate feedback across wet tundra sites with a weakening of the gross CO 2 eC by 15%∼52% (median of posterior estimations); however, due to the substantial increase in heterotrophic respiration under drier conditions, the boreal forest site exerts a positive C-climate feedback with a 236% increase in the gross CO 2 eC (Figure S14 in Supporting Information S1). Our insights on the control of soil moisture on CH 4 emissions are broadly consistent with Watts et al. (2014), where annual summer CH 4 emission budgets were found to fluctuate by ±4% due to the wet/dry cycles; our 1 change in annual precipitation swings annual CH 4 emission by ±7% to ±27% (median of posterior estimations) across the seven Alaskan sites.

If the 1970-present warming and wetting trend continues, our results indicate all tundra sites will exhibit a positive gross We also find that ERCO 2 is predominantly driven by heterotrophic respiration flux response: the difference in ERCO 2 sensitivities is consistent with differences in underlying heterotrophic sensitivities to precipitation and temperature. We speculate that under current climate conditions, wet tundra respiration increases with increasing soil moisture (moisture response between A and B in Figure S2 in Supporting Information S1), while bog and fen respiration decrease with increasing soil moisture (moisture response between B and C2 in Figure S2 in Supporting Information S1). These contrasting moisture-respiration responses could be due to different soil inundation conditions, subsidence, soil pH, soil redox potential, vegetation types, microbial structures, diurnal or seasonal cycles of fluctuating water table in hummocks and hollows.

To determine if there are substantial differences in climatic sensitivity during the growing/non-growing season, we compared CH 4 /CO 2 responses to 2 for reference to Scenario 5. Each color represents one single site and the contour lines are two-dimensional kernel density derived from 4,000 random samples from posterior estimate constrained by eddy covariance net CO 2 and CH 4 data (50%, 75%, and 95% distributions) derived from CARDAMOM posterior states and parameters. Refer to Figure 1 for descriptions of lines and intersections. The CH 4 GWP metric (Myhre et al., 2013) has been extensively used across investigations to quantify the impact of CH 4 fluxes on the evolution of the Earth System. However, Neubauer and Megonigal (2015) have refined the characterization of the CH 4 radiative forcing impact relative to CO 2 , and advocate for a "sustained-flux GWP" (SGWP) metric-to account for the impact of sustained shift in fluxes over time; this amounts to a considerably higher greenhouse gas impact of CH 4 relative to CO 2 (CH 4 SGWP = 45), relative to the GWP assumed in our analysis (CH 4 GWP = 28). To test whether SGWP and GWP metrics lead to consistent or conflicting conclusions on the sign of the inferred C-climate feedback responses in our study, we re-derived the results presented in Figure 7 using a value of z = 45 in Equations 15 and 16. The fluxes summarized in Figure 7 ( CO 2 eC and component fluxes, namely ERCH 4 , ERCO 2 and GPP) derived using both SGWP and GWP are presented quantified in Table S2 in Supporting Information S1. While we find that the choice of SGWP and GWP has a substantial impact on both ERCH 4 and the overall CO 2 eC fluxes, we find that the sign of the CO 2 eC sensitivity to climate remains unchanged, and therefore our conclusions are not dependent on which CH 4 GWP metric is assumed. We use 100-yr scale GWP value throughout the main text and figures to be consistent with recent wetland CH 4 investigations (Jackson et al., 2020;Lawrence et al., 2015;Peltola et al., 2019;Saunois et al., 2016Saunois et al., , 2020;;Tian et al., 2016;Webster et al., 2018).

Overall, our results indicate that accurately resolving both (a) mean CH 4 and CO 2 flux rates across both tundra and boreal ecosystems and (b) their associated climate sensitivity will be critical for determining the sign and magnitude of the northern high latitude ecosystem C-climate feedbacks in the coming decades. Our data-model fusion approach demonstrates joint observational constraints on CO 2 and CH 4 , which is the key to understanding ecosystem C cycle responses to the climate in the coming decades. Future studies could use this approach to jointly constrain CO 2 and CH 4 fluxes at a larger number of EC sites with longer records of observations. We emphasize the importance of adding long-term ground and spaceborne observations with better spatial coverage in northern high latitude to constrain terrestrial-atmosphere C balance and the C-climate feedback. While DALEC-JCR does broadly capture the seasonal and IAV of both CH 4 and CO 2 fluxes at all sites, the model does not explicitly represent soil energy dynamics, permafrost thaw, snow cover, snowmelt infiltration, soil consumption of CH 4 , soil texture, soil pH, and soil redox potential, which are potentially key for accurately resolving C-H 2 O cycle sensitivities to climate. We anticipate that the addition of these processes, along with integrating higher temporal/spatial resolution observations will improve modeled seasonal/inter-annual variations of combined CO 2 -CH 4 fluxes and their climate sensitivity.

We infer CH 4 and CO 2 flux sensitivities to climate from the contemporary meteorological forcing and flux and state observations. This derivation assumes linear responses of carbon cycling to climate (Equations 14-16). While the linear sensitivities presented here provide a first-order estimate of flux sensitivity to contemporary climate variability, larger and/or sustained climate perturbations can potentially induce non-linear responses such as (a) non-linear functional responses to climate, for example, heterotrophic temperature and moisture sensitivities or vegetation functional responses (Norton et al., 2023), (b) cumulative legacy effects, and their propagation across the terrestrial carbon cycle states (e.g., productivity increases and subsequently lagged growth of soil organic C states), and (c) long-term processes unrepresented in the model-data integration analysis, for example, shrub expansion, permafrost thaw or nutrient cycling, and/or (d) a shift of ecosystem states beyond a climate threshold or tipping point (Au et al., 2023;Luo et al., 2011) We highlight that machine learning methods have been used to upscale site-level CO 2 /CH 4 emission to regional/global budgets (Peltola et al., 2019;Tramontana et al., 2016); in this context, our MDF approach can potentially be used to investigate the possibility of propagating parameter knowledge from one site to another, supervised by Bayesian probabilities, observations, and physical/biochemical laws from process-based models; this approach could be key to constrain regional/global CH 4 and CO 2 fluxes-and their associated climate sensitivities. We also note that further investigations on the seasonal sensitivities of CO 2 -CH 4 fluxes are also needed, given that climate changes may be unique to each season, and that the mechanisms affecting growing/non-growing season fluxes may differ (Natali et al., 2015;Treat et al., 2018aTreat et al., , 2018b;;Watts et al., 2014;Zona et al., 2016).

Finally, we highlight that a growing number of satellite greenhouse gas observations can help constrain regional/ global CO 2 -CH 4 budget and climate sensitivity estimations using our approach, such as wetland CH 4 emissions from inversions of GOSAT data (Ma et al., 2022a;Lu et al., 2021;J. D. Maasakkers et al., 2019;J. Maasakkers et al., 2020;Turner et al., 2015;Y. Zhang et al., 2021) and the GOSAT-derived net biosphere C exchange (NBE) data set (CMS-Flux;Liu et al., 2014Liu et al., , 2018)). Furthermore, integration of top-down greenhouse gas flux estimates into land biosphere data-model fusion analyses (e.g., Bloom et al., 2020aBloom et al., , 2020b;;Quetin et al., 2020;Yin et al., 2020) can be potentially extend the quantification of CH 4 and CO 2 flux climate sensitivities to regional or global scales.

We find that the CH 4 (in CO 2 equivalent GWP) response to 1970-present climate change is 50% greater than CO 2 at tundra sites but 30% less at the boreal site. Contemporary CH 4 and CO 2 flux sensitivities indicate that high-latitude wetland ecosystems will amplify C-climate feedbacks in tundra but dampen them in boreal forests if the 1970-present climate change trend continues. Precipitation dominates CH 4 sensitivities to climate through changes in the soil moisture. The CO 2 sensitivities are predominantly temperature driven at the tundra sites but are dominated by precipitation (through Rh suppression) at the boreal site.

Our results agree with previous studies that although the amount of CH 4 emissions is about one magnitude smaller than CO 2 on a molar basis, their responses to climate change play an important role in C-climate feedback (X. Zhu et al., 2013;Zhuang et al., 2007). Over the permafrost region, Lawrence et al. (2015) predicted that 10%-drier soil conditions will accelerate net CO 2 emissions but strongly suppress growth in CH 4 emissions, resulting in a negative C-climate feedback and a 50% lower GWP. We find distinct responses of CO 2 emissions to soil moisture changes between wet tundra and forest bog sites, which requires further investigation across individual ecosystem types. Our results highlight the relative importance of both CO 2 and CH 4 biogeochemical sensitivities to climate, which ultimately need to be jointly quantified to accurately resolve the sign and magnitude of high-latitude ecosystem C-climate feedback in the coming decades.