Net exchanges of methane and carbon dioxide on the Qinghai-Tibetan Plateau from 1979 to 2100

The TEM is a process-based ecosystem model that simulates the biogeochemical cycles of C and N between terrestrial ecosystems and the atmosphere (Zhuang et al 2004, 2010). Within the TEM, two sub-models, namely the soil thermal model (STM) and the updated hydrological model (HM), are designed to simulate the soil thermal profile and hydrological processes with a daily time step, respectively. The STM is well documented for arctic regions and the Tibetan Plateau (e.g. Zhuang et al 2010, Jin et al 2013). The HM, inherited from the water balance model (WBM) (Vörösmarty et al 1989), has six layers for upland soils and a single box for wetland soils (Zhuang et al 2004). We assume the maximum wetland water table depth to be 30 cm following Zhuang et al (2004), so that soils are always saturated below 30 cm. These physical variables then drive the carbon/nitrogen dynamics module (CNDM), which uses spatially explicit information on climate, elevation, soil, vegetation and water availability, as well as soil- and vegetation-specific parameters, to estimate the carbon and nitrogen fluxes and pool sizes of terrestrial ecosystems (Zhuang et al 2004, 2010).

The methane dynamics module (MDM) was first coupled with the TEM by Zhuang et al (2004), to explicitly simulate the processes of CH₄ production (methanogenesis), CH₄ oxidation (methanotrophy) and CH₄ transport between the soil and the atmosphere. During the simulation, the water table depth estimated by the HM, the soil temperature profile estimated by the STM and the labile soil organic carbon estimated by the CNDM, as well as other soil and meteorological information, were fed into the MDM so that the whole TEM and MDM are fully coupled (figure S1). CH₄ production, a strictly anaerobic process, is modeled as a function of the labile carbon substrate, temperature, pH, and the redox potential in saturated soils. CH₄ oxidation, occurring in the unsaturated zone, depends on the soil CH₄ and O₂ concentration, temperature, moisture and the redox potential. The net CH₄ fluxes at the soil/water–atmosphere boundary are summations of different transport pathways (i.e. diffusion, ebullition and plant-aided emission), with a positive value indicating a CH₄ source to the atmosphere and a negative value for a CH₄ sink (Zhuang et al 2004).

Data used to calibrate the TEM were measured at the Luanhaizi wetland on the northeastern Tibetan Plateau (37° 35' N, 101° 20' E) from July 2011 to December 2013. This wetland is classified as an alpine Marsh, and has accumulated rich soil carbon of 24.5% on average for the top 30 cm soil layer due to slow decomposition. An eddy covariance measurement system was installed at a height of 2 m above the wetland surface, which recorded CH₄ fluxes every half hour (Yu et al 2013). A micro-meteorology station was set up adjacent to the eddy covariance tower to measure major environmental variables at half hour frequency, including air and soil temperature, total precipitation, downward shortwave irradiance and relative humidity. We parameterized the STM using the measured soil temperature at a 5 cm depth. Major parameters for the MDM were calibrated so that the simulations match the observed CH₄ fluxes. Due to a lack of high quality measurements, the water table depth was calibrated indirectly by fitting the final CH₄ fluxes, and compared to the observed precipitation as a reference. Parameters for the CNDM were obtained from our previous modeling studies on the Tibetan Plateau (Zhuang et al 2010, Jin et al 2015).

The mathematical structure of the TEM is fixed and can be expressed as below:

$$\begin{eqnarray}&&\hat{Y}=f(X| \theta )+e\end{eqnarray} \tag{ 1 }$$

where $f$ is a conceptual function that represents all processes within the TEM, $\hat{Y}=({y}_{1},{y}_{2},\ldots {y}_{N})$ is a vector of the model outputs (e.g. time series of daily soil temperature or CH₄ fluxes), X is the model input data, $\theta =({\theta }_{1},\;{\theta }_{2},{\theta }_{3},\ldots ,{\theta }_{n})$ are $n$ unknown parameters to be calibrated, and $e=\{e{(\theta )}_{1},e({\theta }_{2}),\ldots ,e({\theta }_{N})\}$ are independently and identically distributed (i.i.d.) errors with zero mean and constant variance. The goal of model calibration with classical methods is thus to find a parameter set $\theta$ such that the predefined statistics of $e$ can be minimized. In contrast, Bayesian theory treats $\theta$ as random variables having a joint posterior probability density function (pdf). The posterior pdf of $\theta$ can be evolved from the prior distribution with observations $Y$ such that:

$$\begin{eqnarray}&&{p}_{{\rm{post}}}(\theta | Y)=\displaystyle \frac{L(Y| \theta ){p}_{{\rm{prior}}}(\theta )}{{\int }^{}L(Y| \theta ){p}_{{\rm{prior}}}(\theta )}\end{eqnarray} \tag{ 2 }$$

where $L(Y| \theta )$ is the likelihood function. Assuming a non-informative prior in the form of $p(\theta )\propto {\sigma }^{-1}$ and residuals which are i.i.d normal (Vrugt et al 2003), the likelihood of a specific parameter set $\theta ^{\prime}$ can be computed as:

$$\begin{eqnarray}&&L(\theta ^{\prime} | Y)={\rm exp}[-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}{(\displaystyle \frac{e{(\theta ^{\prime} )}_{i}}{\sigma })}^{2}]\end{eqnarray} \tag{ 3 }$$

The influence of $\sigma$ can be integrated out when $L(\theta ^{\prime} | Y)$ is plugged into equation (2), so that the posterior pdf of $\theta ^{\prime}$ is:

$$\begin{eqnarray}&&{p}_{{\rm{post}}}(\theta ^{\prime} | Y)\propto {[\displaystyle \sum _{i=1}^{N}e{{(\theta ^{\prime} )}_{i}}^{2}]}^{-0.5\times N}\end{eqnarray} \tag{ 4 }$$

For complex nonlinear system models like the TEM, however, it is impossible to obtain an explicit expression for ${p}_{{\rm{post}}}(\theta ^{\prime} | Y),$ making analytical optimization out of the question (Vrugt et al 2003). Alternatively, Markov chain Monte Carlo (MCMC) methods are well suited to solving these problems.

In this study, we implemented an adaptive MCMC sampler, named the shuffled complex evolution Metropolis algorithm (SCEM-UA) for global optimization and parameter uncertainty assessment. While the theoretical bases and computational implementation of the SCEM-UA can be found in (Vrugt et al 2003), we outline the key steps with an example of STM optimization below:

• (i)  
  Initialize parameter space. Select parameters to be calibrated, and assign prior range to each parameter (table 1).
• (ii)  
  Generate sample. Randomly sample $s$ sets of parameter combination $\{\vec{{\theta }_{1},\;}\vec{{\theta }_{2}},\ldots ,\vec{{\theta }_{s}}\}$ from the uniform prior distribution, where $\vec{{\theta }_{i}}$ is a vector of eight parameters for the STM in table 1.
• (iii)  
  Rank sample points. Compute the posterior density of each $\vec{{\theta }_{i}}$ using equation (3), and sort the $s$ points in 8-dimensions with decreasing order of posterior density. Store the $s$ points and their corresponding posterior in array $D$ with dimensions of $s\times 9$.
• (iv)  
  Initialize Markov chains. Assign the first $k$ elements of $D$ as the starting locations of $k\;$ sequences.
• (v)  
  Partition D into complexes. Partition $s$ rows of $D$ into $k$ complexes $\{{C}^{1},{C}^{2},\ldots {C}^{k}\},$ such that each complex contain $m=s/k$ points. More specifically, the $i-th$ complex is assigned every $k(j-1)+i$ rows of D, where $j=1,\ldots ,m$.
• (vi)  
  Evolve each complex and sequence. Use the sequence evolution Metropolis algorithm to evolve each sequence and complex (Vrugt et al 2003).
• (vii)  
  Shuffle complexes. Combine each point in all complexes back with $D$ in order of decreasing posterior density.
• (viii)  
  Apply stop rule. Stop when convergence criteria are satisfied; otherwise, go to step (v). A maximum of 100 000 runs for the STM and 200 000 runs for the MDM was imposed to override the stopping rule for computational considerations.

We selected $s=500$ and $k=10$ following Vrugt et al (2003) for complex model optimization. Given the complexity of our coupled-TEM, each run cost at least 1 minute, including reaching model equilibrium, model spin-up and transient simulations for the period 2011–2013. In this case, parallel implementation of the SCEM-UA was required to obtain computational efficiency. We developed a parallel R version of the SCEM-UA using the Open Message Passing Interface (MPI; Gabriel et al 2004) and the Rmpi package (Yu et al 2002) on Purdue University's Conte computer clusters. This parallel SCEM-UA program would evoke a master node that controlled the workflow and message communication among 10 slave nodes. Each slave node was responsible for the computation of posterior density (i.e. step (iii)) and for the evolution of a specific group of complexes and sequences (i.e. step (vi)). Compared to the parallel implementation of the SCEM-UA in Vrugt et al (2006), our program does not require the target model (e.g. the TEM in this study) to be written in the MPI structure, but only treats it as a black box to be executed using a system call in R. The source code for this parallel R version of the SCEM-UA is available upon request.

To make spatiotemporal estimates of the CO₂ and CH₄ fluxes on the QTP using the TEM, the model was run at a spatial resolution of 8 × 8 km and with a daily time step from 1979 to 2100. Static data including vegetation types, elevation, and soil texture were the same as those in Jin et al (2013). Soil pH data was derived from the China Dataset of Soil Properties for Land Surface Modeling by Wei et al (2013). The wet soil extent from Papa et al (2010) was used to determine the distribution of wetlands (CH₄ source) and uplands (CH₄ sink) within each pixel (figure 1(b)). For this temporally dynamical and 25 × 25 km resolution data set, we interpolated the maximum fractional inundation over 8 × 8 km using the nearest neighbor approach. It should be noted that the actual wetland distribution was less continuous than was shown in figure 1(b). The inland water body was excluded based on the IGBP DISCover Database (Loveland et al 2000). Seasonal flooding of the wetland was indirectly represented by the fluctuation of the water table depth below and above the soil surface in our model. Climate input data for the historical period of 1979–2011, including radiation, precipitation, air temperature and vapor pressure, were interpolated from the latest global meteorological reanalysis product ERA-interim (0.75° grid) published by the European Centre for Medium-Range Weather Forecasts (ECMWF) (Dee et al 2011). To reduce bias in air temperature caused by complex topography, we used elevation as a covariate during interpolation according to Frauenfeld et al (2005). Future climate scenarios from 2006 to 2100 were generated under four representative concentration pathways (RCPs) with six global climate models (GCMs) within Coupled Model Inter-comparison Project phase 5 (CMIP5), including: CESM1-CAM5, GFDL-CM3, GISS-E2-H, HadGEM2-ES, MPI-ESM-Mr and MRI-CGCM3 (figures S3, S4). A total of 24 simulations (4 scenarios × 6 GCM datasets) were processed to construct future projections. These selected GCMs, compared with many other candidates, in general had smaller biases in surface temperature and total precipitation across the Tibetan Plateau (Su et al 2013). A detailed description of these CMIP5 GCMs and the data processing method are provided in methods S1.

Zoom In Zoom Out Reset image size

To highlight the spatial pattern of the climate induced changes in GHG fluxes, we calculated the spatial difference between baseline simulations using ECMWF data in the 2000s and future predictions with four RCP scenarios in the 2090s. A statistical test of the difference in means was performed, and regions with statistical insignificance (α = 0.05) were masked out. To quantify the GHG effect of the continuous CH₄ and CO₂ fluxes from the QTP for the study period of 1979–2010, we calculated the net radiative forcing impact according to Frolking et al (2006). The lifetime of an individual net input of CH₄ into the background atmosphere was represented by a first-order decay function such that:

$$\begin{eqnarray}&&r{{\rm{CH}}}_{4}(t)={r}_{0}\times \mathrm{exp}(-t/{\tau }_{{\rm{CH4}}})\end{eqnarray} \tag{ 5 }$$

where ${r}_{0}$ is the initial CH₄ perturbation, and ${\tau }_{{\rm{CH}}4}$ is the adjusted time for CH₄ (∼12 y). A linear superposition of 5 different first-order decay pools was used to describe the more complicated behavior of CO₂ in the atmosphere:

$$\begin{eqnarray}&&r{{\rm{CO}}}_{2}(t)\sim \displaystyle \sum _{i=0}^{4}\;{f}_{i}{\rm exp}(-t/{\tau }_{i})\end{eqnarray} \tag{ 6 }$$

where the fraction of each pool, ${f}_{i},$ and the pool specific adjustment time, ${\tau }_{i},$ were set to be 26%, 24%, 19%, 14%, 18% and 3.4, 21, 71, 421, 10⁸ y, respectively (Frolking et al 2006). The instantaneous total radiative forcing from individual gas contributions since the reference time (here the year 1979) is given by:

$$\begin{eqnarray}&&{\rm{RF}}(t)=\displaystyle \sum _{i=0}^{5}({\xi }_{i}{A}_{i}{f}_{i}\cdot \displaystyle {\int }_{o}^{t}{{\rm{\Phi }}}_{i}(s){\rm exp}(-(t-s)/{\tau }_{i}){\rm{d}}s)\end{eqnarray} \tag{ 7 }$$

in which $i\ =\ 0-4$ is CO₂ and $i\ =\ 5$ is CH₄, and ${\xi }_{i}$ is a multiplier for indirect effects (1.3 for CH₄ and 1.0 for CO₂), ${A}_{i}$ is the GHG radiative efficiency (1.3 × 10⁻¹³ W m⁻² kg⁻¹ for CH₄ and 0.0198 × 10⁻¹³ W m⁻² kg⁻¹ for CO₂), ${f}_{i}$ is the fractional multiplier (1 for CH₄ and see equation (6) for CO₂ values), and ${{\rm{\Phi }}}_{i}(s)$ is the net flux of GHG $i$ at time $s$ relative to the reference year of 1979. The integral term is thus the cumulative flux of gas $i$ at time $t$ since the reference start point (i.e. the year 1979) after partial to complete decay in the atmosphere. The numerical integration was applied with an annual time step. It should be noted, however, that our calculation here was only for the goal of comparing the relative contributions of CO₂ and CH₄ fluxes since 1979 to the radiative forcing, rather than to give an accurate estimation of the absolute values of the GHG effect.

By applying the SCEM-UA method, initial parameter ranges evolved into narrower posterior intervals (table 1). Both the STM and MDM outputs were able to reproduce the seasonal dynamics of the observed daily soil temperature and CH₄ fluxes (figure 2). The adjusted R² and RMSE for the soil temperature simulated with the optimal STM are 0.95 and 1.88 °C, respectively. The model performance is comparable with other modeling results (Wania et al 2009, Jin et al 2013, Zhu et al 2014), with only one underestimation during the summer of 2013. The MDM also performed well (R² = 0.82, RMSE = 18.41 mg CH₄ m⁻² day⁻¹) in capturing the annual magnitude, cycling and a small peak of CH₄ burst in spring. Considering the uncertainties in the CH₄ flux measurements (Yu et al 2013) and the CH₄ model structures (Wania et al 2013), our model performance was comparable with similar studies (Wania et al 2009, Lu and Zhuang 2012, Zhu et al 2014). Parameter uncertainty was well constrained (figures 2(a), (c)), indicating that a global search near the optimal space would produce many parameter sets to allow model simulation matching the observations. The water table depth followed the observed daily precipitation pattern (figure 2(b)). A severe drought was detected for the 2012 summer, which was also reflected by the distinctively low CH₄ emission during that year.

Zoom In Zoom Out Reset image size

Due to a limited number of available field studies of CH₄ fluxes on the QTP, our model validation was done by comparing values reported in the literature to our simulations for the nearest 8 × 8 km pixel (table 2). Model estimates reasonably match most field measurements with respect to the mean and range, except for an apparent underestimation of the exceptionally high CH₄ emission rate from the littoral zone in Chen et al (2009). Considering the high variation among field measurements from different wetland types and vegetation covers, the model simulations were able to fit the mean (R² = 0.87, RMSE = 67 for all data and R² = 0.96, RMSE = 21 mg CH₄ m⁻² day⁻¹ excluding Chen et al (2009)). The validation gave us confidence for model extrapolation and regional CH₄ budget estimation.

Our estimated annual regional mean net primary production (NPP) (145–170 gC m⁻² y⁻¹) increased slightly for the historical period, and fluctuated under future scenarios (figure 3(a)). NPP increased almost 40% under RCP6.0 ($241\pm 22.3gC$ m⁻² y⁻¹) and 50% under RCP8.5 ($252\pm 22.4gC$ m⁻² y⁻¹) at the end of the 21st century compared to the 2000s. RCP2.6 and RCP4.5 showed a similar trend in NPP under higher emission scenarios before the 2050s, but leveled off thereafter at around 210 and 230 gC m⁻² y⁻¹, respectively. Simulated regional mean NEP showed high variation around 0 gC m⁻² y⁻¹ for the period 1979–2011, and decreased continually under future projections except for RCP8.5 (figure 3(b)). The NEP density was lowest under RCP4.5 ($-33.1\pm 9.5gC$ m⁻² y⁻¹), followed by RCP6.0 ($-30.3\pm 11.4gC$ m⁻² y⁻¹) and RCP2.6 ($-26.8\pm 11.5gC$ m⁻² y⁻¹) for the 2090s. Under RCP8.5, the NEP density decreased to the 2060s but reversed to $-19.5\pm 16.1gC$ m⁻² y⁻¹, possibly because heterogeneous soil respiration (R_H; figure S5) was stimulated more than NPP under this warmest and wettest scenario.

Zoom In Zoom Out Reset image size

Simulated annual mean CH₄ emissions from potential wetland areas increased gradually from 6.3 gCH₄ m⁻² y⁻¹ in 1979 to 8.5 in the 2050s under different scenarios, but diverged noticeably thereafter (figure 3(c)). The highest emissions for the 2090s from RCP8.5 projections ($12.7\pm 1.9$ gCH₄ m⁻² y⁻¹) increased 70% relative to the beginning of the 21st century, followed by 46%, 30% and 16% under RCP6.0, RCP4.5 and RCP2.6, respectively. The relative change percentages of CH₄ emissions were in general higher than those of NPP, implying that CH₄ production was more favored than photosynthesis under future climate conditions. Annual mean CH₄ uptake density (increased from 0.13 in the early 1980s to 0.17 gCH₄ m⁻² y⁻¹ around the 2050s, and ended up between 0.16 to 0.23 gCH₄ m⁻² y⁻¹ in the 2090s) was much smaller than that of emissions, while the inter-annual variations were similar under the four RCP scenarios (figure 3(d)).

Spatial patterns of the decadal mean NME showed substantial variations over the QTP (figure 4). For the 2000s, net CH₄ emissions simulated using reanalysis data were similar to the average of the RCP scenarios (figures 4(a) and (b)), with high net CH₄ emissions occurred at Zoige wetland and in the southwest part of the QTP. Net CH₄ sinks based on the two maps were comparable in magnitude, but high CH₄ uptakes were found at the southern edge of the QTP in figure 4(a) rather than in the northeastern part in figure 4(b). At the end of the 21st century, the NME increased noticeably with respect to the magnitude over the QTP except for some southeastern parts (figures 1(c)–(f)). Among the models, RCP8.5 led to the strongest CH₄ emissions in wetlands and the highest CH₄ consumption in uplands. The difference in means between the RCP scenarios and the 2000s historical mean was statistically significant (α = 0.05) for most of the study area (figures S6(a)–(d)). Future climate change would strengthen the upland CH₄ sink in the inner QTP, and substantially stimulate wetland CH₄ emission in northeastern regions. For simulated NEP, most of the study region had a negative value (i.e. CO₂ sink), with the magnitude up to −30 gC m⁻² y⁻¹ during the 2000s (figures 5(a) and (b)). The patterns were comparable with contemporary estimation for alpine grasslands by Zhuang et al (2010) and Piao et al (2012). Strong CO₂ sinks occurred in the southeast of the QTP, where forest and shrubs were the dominant vegetation cover (figure 1(a)). CO₂ sources accounted for a very limited portion of the total study area. Spatial patterns of NEP evolved substantially under future scenarios (figures 5(c)–(f)). The QTP became a uniform carbon sink under RCP2.6, indicating disproportional climate impacts on different ecosystem types. Increases in NEP were profound under two median scenarios, with differences statistically significant (α = 0.05) for most of the grassland (figures S6(f), (g)). Among these, the highest sink (up to 50 gC m⁻² y⁻¹) occurred in the forest ecosystem under RCP4.5 and in meadow regions under RCP6.0. NEP under RCP8.5 was distinctively lower than results under other scenarios and not statistically different from the historical 2000s mean (figure S6(h)), suggesting similar changes in the magnitude of NPP and R_H for all ecosystem types.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our estimate of the total CH₄ emissions from natural wetland over the QTP was 0.95 Tg CH₄ y⁻¹ during the 2000s, which was within the estimation range of 0.13–2.47 Tg CH₄ y⁻¹ from several other studies (table 3). The high variation among different estimates was mainly due to the uncertainty in wetland area estimates. Total CH₄ consumption from those upland soils were 0.19 Tg CH₄ y⁻¹ (i.e. 20% of the regional CH₄ emissions) for the 2000s, and increased slower than emissions under future scenarios (table 4), indicating that the QTP was likely to be a stronger CH₄ source in the 21st century. The simulated regional NEP of $-10.22$ Tg C y⁻¹ for the period 2006–2011 was lower than those for other modeling studies covering similar regions (Zhuang et al 2010, Piao et al 2012), but higher than that of Yi et al (2014). The near neutral properties of the historical NEP (also see figure 3(b)) were consistent with the results of Fang et al (2010), who found that soil C stock in China's grassland did not show a significant change during the past two decades. Future NEP increased under all RCP scenarios (table 4). Counter-intuitively, higher increases in NEP occurred under RCP4.5 and RCP6.0 instead of the warmest and wettest RCP8.5, indicating that carbon accumulation was more likely favored with modest warming and wetting.

Instantaneous radiative forcing due to net CH₄ emission and net CO₂ sequestration since 1979 increased to a plateau during the 1990s, dropped to zero around the 2010 s, and decreased almost linearly afterwards (figure 6(a)). Therefore, climate change in the 21st century is likely to trigger negative feedback (cooling) in the climate system. Among these, the fastest decreased rate occurred under RCP4.5, indicating that this moderate climate warming scenario stimulated vegetation production much more than the methanogenesis process. A seemingly level-off trend after the 2080s was identified for the RCP8.5 scenario, most likely because of declining NEP (figure 3(b)). Cumulative radiative forcing was positive until the 2030s, and become negative in an accelerated manner (figure 6(b)). Thus, the cumulative GHG fluxes from the QTP will exert a slight warming effect on the climate system until the 2030s, but an increasingly stronger cooling effect thereafter. The tipping point was 50 years after the reference zero time point, roughly the time required for nearly complete removal of the initial CH₄ perturbation input. At the end of the 21st century, cumulative mean radiative forcing was between −0.25 and −0.35 W m⁻², depending on different RCP scenarios. Overall, our results show that given a sustained CH₄ source and CO₂ sink on the QTP, the net GHG warming effect will only peak after a few decades and will eventually contribute a more cooling effect to the climate system.

Zoom In Zoom Out Reset image size

Overall, our model calibration results were satisfactory, but a few tips should be mentioned when applying the global optimization method to complex system models with high parameter dimensions. Global optimization for problems without explicit analytic expressions of the objective function is challenging, because the algorithm must avoid being trapped by several local optima, while maintaining robustness in the presence of parameter interaction and non-convexity of the response surface, and having high efficiency in searching high dimensional space (Duan 2013). When applying the SCEM-UA, a tradeoff between goodness of fit and computational cost is still a problem for users.

The speed of algorithm convergence highly depends on the model structure and parameter interactions. In our case, major parameters evolved to a narrow range after 100 000 total iterations (figure S2). In contrast, the MDM optimization converged much slower despite the smaller number of parameters to be optimized. We argue that this is mainly due to the multi-scalar function method used to simulate CH₄ production and oxidation. For instance, if scalars ${f}_{1}(x)$ and ${f}_{2}(x)$ are shaped with two parameters only, it is intuitive to imagine that many pairs of these two parameters can produce similar results, as long as they can adjust scalars ${f}_{1}(x)$ and ${f}_{2}(x)$ in the opposite direction. Given that the scalar approach is extensively implemented in current ecosystem models (Zhuang et al 2004, 2010, Wania et al 2013, Zhu et al 2014), convergence criteria (such as the Gelman and Rubin diagnostic suggested by Vrugt et al (2003)) can hardly be reached for most parameters when calibrating these models.

On the other hand, sub-optimal is usually sufficient in practice. Increasing the sampling size from parameter space can push the posterior of the objective function to the high-end value (figure S2), but is by no means a guarantee of better model performance. Most likely, many behavioral parameter combinations with a similar capability of reproducing the observations can be found by a search conducted in the feasible parameter space (Vrugt et al 2003). As shown in our example, the top 500 parameter sets (whose values can differ) for either the STM or MDM produced close goodness-of-fit and small variations in simulations (figure 2), indicating that the benefits from additional searching are marginal.

To quantify and compare the net GHG effect of a sustained CH₄ source and CO₂ sink on the QTP, the total radiative forcing rather than the more widely known metrics of the GWP was computed in our study according to Frolking et al (2006). As a tool originally designed for evaluating and implementing policies to control multi-GHGs, the GWP is defined as the time-integrated radiative forcing due to a pulse emission of a given component, relative to a pulse emission of an equal mass of CO₂ (IPCC, 2013). It has usually been integrated over a somewhat arbitrary time horizon of 20, 100 or 500 years, of which the choice of 100 years is most commonly adopted (e.g. CH₄ is 28 times CO₂ in terms of warming effect). Applying the GWP to biogeochemical cases could be problematic as GHG fluxes are often sustained and temporally dynamical from a long existing sink or source, even though many studies continue to use it because of its simplicity (e.g. Zhu et al 2013, Gatland et al 2014, Vanselow-Algan et al 2015). The method proposed by Frolking et al (2006) goes beyond the standard GWP approach such that (1) persistent emission or uptake of GHGs can be accounted for, and (2) the instantaneous radiative forcing for each year of simulation is quantified to allow a comparison of multiple gases in common units at any given time. Although the assumption of a near constant background atmosphere in Frolking's method is open to debate, and the decision on the appropriate pulse emission to consider (the new flux rate versus the change in the flux rate) could greatly change the results of the simulation, it appears a useful method in biogeochemical studies when accounting for the budget of multiple sustained GHG fluxes. Frolking's method has evolved with new findings in the atmospheric sciences. For example, the fraction of CO₂ remaining in the atmosphere after a pulse input can be divided into four components (Joos et al 2013) instead of the original five-pool setting, and the multiplier of the indirect effect for CH₄ is now 1.65 according to IPCC (2013). A test of these and other parameter changes is beyond the scope of this study.

This is the first study to quantify both CH₄ emissions and net carbon exchanges on the QTP. However, the quantitative analysis is uncertain due to incomplete representation of physical and biogeochemical processes in the model (Bridgham et al 2013, Melton et al 2013, Bohn et al 2015), inaccurate model assumptions (Meng et al 2012), variations in the forcing climate data (Su et al 2013) and the extrapolation from site to regional scale. While progresses in model structures and mechanisms are usually slow and incremental, efforts to reducing data uncertainty are more feasible in a short period to improve model predictability.

First, seasonal and inter-annual dynamics of the wetland extent is critical to CH₄ modeling. Synthetic aperture radar (SAR) is currently the first choice to delineate wetland distribution, such as the monthly distribution of surface water extent with ∼25 km sampling intervals used in this study (Papa et al 2010). However, optical remote sensing is highly sensitive to cloud or vegetation cover. Alternative data from passive and active microwave systems that will penetrate cloud and vegetation cover is favored (Schroeder et al 2010), but is currently unavailable over our study area. Some contemporary methane models, such as VIC-TEM and SDGVM, are capable of outputting dynamical wetland extent as an internal product (Hopcroft et al 2011, Lu and Zhuang 2012), but they tend to substantially overestimate wetland area (Melton et al 2013). In order to capture the seasonal flooding area, our study combined a static max potential inundation map and simulated water table depth to represent the wetland fluctuation. Improving the modeling abilities to capture wetland distribution and extent, especially the seasonality of water table dynamics, should be a research priority in the future (Zhu et al 2011).

The second uncertainty in quantifying CH₄ emissions on a regional scale is from the spatial-scale extrapolation across highly heterogeneous but poorly mapped wetland complexes (Bridgham et al 2013). In this study, we only calibrated our model at one alpine wetland site by assuming that the remaining wetlands over the QTP have the same inert characteristics, but differ with climatic conditions. As a matter of fact, alpine wetlands on the QTP can be classified into peatlands, Marshes, and swamps, which have distinct characteristics in vegetation cover, hydrological processes, and soil history (Chen et al 2013). However, due to severe field experimental conditions, high quality observational data are mostly limited to the Lunhaizi wetland in Haibei Station (Hirota et al 2004, Yu et al 2013) and the Zoige wetland (Chen et al 2009, 2013). Field researchers should consider expanding field measurement footprints in the future, and data sharing with modelers is highly recommended.

Finally, much of the uncertainty in future projections is due to the poor agreement among CMIP5 GCMs under the RCP scenarios (figures S3, S4). The majority of the GCMs have cold biases of 1.1–2.5 °C in air temperature for the QTP, while they overestimate annual mean precipitation by 62%–183% (Su et al 2013). A multi-model ensemble approach, as was suggested in the IPCC AR5 report (IPCC, 2013), was used to configure the climatology uncertainty. While the statistical interpolation method used to generate fine resolution data is simple and computationally efficient (Wilby et al 1998), dynamical downscaling with regional climate models was more ideal to compensate for the shortage of coarse output from GCMs and to capture details of the complex surface properties of the QTP (Ji and Kang 2012). However, due to the high computational cost, dynamically downscaling climate data that cover representative GCMs under all RCP scenarios is generally not available to ecosystem modelers. A publicly accessible database like this would greatly benefit the research community in future studies.