A soil diffusion-reaction model for surface COS ﬂux: COSSM v1

. Soil exchange of carbonyl sulﬁde (COS) is the second largest COS ﬂux in terrestrial ecosystems. A novel application of COS is the separation of gross primary productivity (GPP) from concomitant respiration. This method requires that soil COS exchange is relatively small and can be well quantiﬁed. Existing models for soil COS ﬂux have incorporated empirical temperature and moisture functions derived from laboratory experiments, but not explicitly resolved diffusion in the 5 soil column. We developed a mechanistic diffusion-reaction model for soil COS exchange that accounts for COS uptake and production, relates source-sink terms to environmental variables, and has an option to enable surface litter layers. We evaluated the model with ﬁeld data from a wheat ﬁeld (Southern Great Plains (SGP), OK, USA) and an oak woodland (Stunt Ranch Reserve, CA, USA). The model was able to reproduce all observed features of soil COS exchange such as diurnal varia- 10 tions and sink-source transitions. We found that soil COS uptake is strongly diffusion controlled, and limited by low COS concentrations in the soil if there is COS uptake in the litter layer. The model provides novel insights into the balance between soil COS uptake and production: a higher COS production capacity was required despite lower COS emissions during the growing season compared to the post-senescence period at SGP,


Introduction
Carbonyl sulfide (COS) is an atmospheric trace gas, with average concentration around 480 pmol mol −1 (Montzka et al., 2007). COS uptake in leaves is coupled with CO 2 uptake through stomatal 20 diffusion and hydrolysis mediated by carbonic anhydrase (CA) (Sandoval-Soto et al., 2005;Stimler et al., 2010;Seibt et al., 2010). Therefore plant COS flux measurements can be used to partition net carbon exchange into gross primary productivity (GPP) and respiration using the quantitative relationship between leaf COS and CO 2 uptake (Campbell et al., 2008;Asaf et al., 2013). Since the COS flux observed at ecosystem level and above is the sum of plant and soil fluxes, soil COS 25 fluxes must be well quantified for COS-based flux partitioning. For many biomes where soils have active COS reactions, neglecting soil COS flux would result in significant biases of estimated GPP (Whelan et al., 2015). To enable COS-based flux partitioning, particularly at larger scales, a soil diffusion-reaction model is needed to generate estimates of soil COS exchange from soil parameters and environmental variables. 30 Laboratory studies have shown that soil COS uptake is a function of soil temperature and moisture, microbial CA enzyme activity, and ambient COS concentration (Kesselmeier et al., 1999;Van Diest and Kesselmeier, 2008). These empirical relationships were used by Kettle et al. (2002) to model global soil COS flux from climatological soil data. Recently, Berry et al. (2013) simulated global soil COS flux in the Simple Biosphere Model (SiB3) with the assumption that soil COS flux is 35 a function of heterotrophic respiration and soil water-filled pore space fraction (WFPS). Here the effect of WFPS on diffusion was represented with an empirical function. Diffusion is not explicitly resolved in any existing model. The effect of a litter layer is also not considered. A litter layer may either act as a barrier for COS diffusion to the soil, or participate in COS exchange. Plant litter has been found to contribute to the surface COS flux in forests (Kesselmeier and Hubert, 2002), with 40 magnitudes comparable to that of soil uptake (Berkelhammer et al., 2014).
Existing models have also not considered the possibility of COS production in oxic, upland soils, based on the assumption that they usually behave as a COS sink (e.g. Kesselmeier et al., 1999;Yi et al., 2007;Van Diest and Kesselmeier, 2008;Liu et al., 2010a). However, soil COS fluxes varied between uptake and strong emissions in a wheat field depending on soil environmental variables 45 (Maseyk et al., 2014). The presence of a compensation point for COS uptake also indicates that uptake and production can co-exist in soils (Kesselmeier et al., 1999;Conrad and Meuser, 2000;Liu et al., 2010a).
Here we develop a diffusion-reaction model for soil COS flux, with four major advantages compared to previous empirical approaches: 50 (1) Diffusion is explicitly resolved, providing a more realistic concentration-dependent uptake in the soil column.
(2) The model accounts for flux activity in surface litter layers.
(3) The model includes COS production terms.
(4) Litter and soil COS uptake and production terms can be constrained or optimized using field 55 data. We evaluate the model using field data from a wheat field in the Southern Great Plains (Billings, OK, USA) and an oak woodland in the Santa Monica Mountains in southern California (Stunt Ranch Reserve, CA, USA).

The diffusion equation with source-sink terms
We construct a two-phase 1D diffusion model with microbial COS source-sink terms in soil and litter layers to calculate surface COS fluxes. Diffusion of COS in the soil and litter is described by Fick's law in porous media,  For computational efficiency, diffusion in the aqueous phase is not prognostically evaluated but forced by the solubility equilibrium. By invoking Henry's law, the above equation assumes that the aqueous concentration is always in equilibrium with the gaseous concentration, ensuring mass balance. Approximating aqueous diffusive flux with this assumption does not significant bias the surface flux (see SI), since aqueous diffusivity of COS is several orders of magnitude smaller than 75 the gas-phase diffusivity (Ulshöfer et al., 1996).
Advective transport is not considered since the evaluation datasets were not affected by advection.
However, advection effects on surface COS flux may not be negligible when there is strong wind pumping causing pressure fluctuations (Massman et al., 1997). Advective transport, and prognostic aqueous processes (diffusion, dissolution and possibly ebullition) will be considered in future 80 development once relevant data are collected.
where N = 25. The thicknesses of control volumes are, The depths of the layer interfaces are defined with, The use of face-centered rather than node-centered control volumes generally gives better evaluation of diffusive fluxes across interfaces.

Discretization
We use the Crank-Nicolson method to discretize the diffusion-reaction equation, which ensures numerical stability when using large time steps (Crank and Nicolson, 1947). The method is implemented with second order accuracy in space and in time. The discretized finite difference equations are a sparse linear system that can be easily solved. 100 We first discretize spatially Eq. (1) on the defined grid by integrating in each control volume (Durran, 2010;Tang and Riley, 2014), where η i = (k H θ w + θ a ) at z = z i is a simplifying notation, J a→0 represents diffusive flux from the atmosphere to soil layer 0, and J i−1→i represents diffusive flux from layer i − 1 to layer i. The flux 105 term J can be evaluated with a central difference scheme, The diffusivity at the interface, D i−1/2 can be interpolated linearly, the interface locates at the center of two computational nodes. For the ease of denotation, we define as the conductance at the interface.

110
The flux at the topmost layer is thus, Again, D a→0 is the diffusivity at the soil-atmosphere boundary, and G a→0 is the conductance. We use the harmonic mean of soil diffusivity D 0 and air diffusivity D a for D a→0 , following the con-ductance model in Tang and Riley (2013), By rewriting Eq. (5) in a matrix form, we obtain an ordinary differential equation (ODE) system for time, where 120 C = (C 0 , C 1 , . . . , C N ) The above linear ODE system is then discretized at time step t = (n + 1/2)∆t, Therefore the evolution in time is, At each time step, the concentration profile is evolved to the next time step with the diffusion-reaction operator. Thus the model can simulate transient conditions in real time. The diffusivity of COS in soil media (D, m 2 s −1 ) is calculated following Moldrup et al. (2003), where θ sat (m 3 m −3 ) is the soil total porosity (θ sat = θ a +θ w ), D m (m 2 s −1 ) is the COS-air diffusiv-135 ity at T ref = 25 • C and 1 atm, n = 1.5 (Bird et al., 2002), and b is parameterized with water retention function (Wingate et al., 2008;Clapp and Hornberger, 1978). We assume that the same function holds for litter layer diffusivity since litter is also a porous medium like the soil.

COS solubility function
In solubility equilibrium, the chemical potentials of COS in gas phase (µ g ) and in aqueous solution (µ aq ) are equal, where p is the partial pressure in gas phase, x is the molar fraction in aqueous phase, φ is the fugacity coefficient, γ is the activity coefficient, and p is the standard pressure. Here, µ g and µ * aq are both chemical potentials under chosen standard conditions. Note that µ g depends only on temperature, whereas µ * aq depends on both temperature and pressure.

150
For a dilute solution, as x → 0 and p/p → 0, we have φ → 1 and γ → 1. This is valid for COS in natural environments because its molar fraction is only ∼ 10 −10 in the atmosphere (Montzka et al., 2007) and ∼ 10 −12 in seawater (Ulshöfer et al., 1996;Von Hobe et al., 1999). Therefore, where k H,xp is the Henry's law constant defined as the ratio of aqueous phase molar fraction to gas 155 phase partial pressure. According to van 't Hoff equation, where ∆ sol H m is the standard partial molar enthalpy of dissolution. We assume that it is constant within the temperature range of ambient conditions. With respect to a reference state T 0 , integrating yields, Temperature dependence of k H,xp can thus be described by the following function, with A and B constants, The dimensionless Henry's law constant, k H,cc , is defined as the ratio of aqueous molar concentration 165 to gaseous molar concentration.

Soil fluxes
Because both COS uptake and production activities exist in soils, and the net COS flux exhibits a linear response to COS concentration (Kesselmeier et al., 1999;Conrad and Meuser, 2000), we 175 assume that uptake and production are separable terms.
Soil COS uptake (U ) is formulated based on Michaelis-Menten kinetics, with dependence on soil temperature and moisture: where K m = 1.9×10 −6 mol m −3 is the Michaelis constant for COS hydrolysis by CA (Ogawa et al., We assume that the temperature limitation function for soil COS uptake reflects the temperature dependence of enzyme activity, which is described as (Peterson et al., 2004;Daniel et al., 2010), where ∆ ‡ G cat (J mol −1 ) is the activation free energy of the transition state of CA enzyme in terms of COS reactions, ∆H eq (J mol −1 ) is the enthalpy change during conversion from activated to inactivated form of the enzyme, T eq (K) is the temperature at which the concentrations of the activated enzyme and the inactivated enzyme are equal, A T is the pre-exponential factor to normalize the equation so that max f (T ) = 1. The function has a temperature optimum, T opt , defined at which the 190 function reaches its maximum value. This temperature optimum is close to the parameter T eq , and always smaller than it (Peterson et al., 2004). Other enzymes in soil microbes may also contribute to COS uptake, for example, COSase, nitrogenase and CS 2 hydrolase (Ogawa et al., 2013). These enzymes are not considered here due to lack of information on their kinetics and distribution in the microbial community. They are also unlikely to be as ubiquitous as CA.

195
Typical fit values of ∆ ‡ G cat for a range of enzymes are between 51 and 88 kJ mol −1 , whereas ∆H eq varies from 86 to 826 kJ mol −1 (Daniel et al., 2010). We assume that ∆ ‡ G cat = 84.10 kJ mol −1 (20.1 kcal mol −1 ) for CA-catalyzed COS hydrolysis, based on calculations of CA-COS nu-cleophilic attack, the rate-determining step of COS hydrolysis, which used [(H 3 N) 3 Zn(OH)] + as a structural analog of the active site of CA (Schenk et al., 2004).

200
The temperature dependence (Eq. 22) has a similar mathematical form as that in Kesselmeier et al. (1999), but the temperature limitation of COS uptake is now linked to enzyme kinetics that can be determined in the laboratory. Using ∆ ‡ G cat = 84.10 kJ mol −1 , ∆H eq = 358.9 kJ mol −1 , and T eq = 293.14 K, we can approximate the temperature function ϕ(T ) in Kesselmeier et al. (1999) ( Fig. 3a). In this study, T eq will be an adjustable parameter to fit temperature dependence for different 205 soils (see section 2.4.2 and Table 2).
For the parameterization of moisture dependence, we use a simple bell-shaped function, described by the form of Rayleigh distribution function, where w opt (m 3 m −3 ) is the optimal water content for COS uptake, and A W is the normalization 210 factor such that max g(w) = 1. The parameter w opt can be adjusted to fit observed data (see section 2.4.2 and Table 2). The function has a wider range than that in Kesselmeier et al. (1999) (Fig. 3b), and does not approach zero at high soil moisture. It represents the dependence of microbial uptake per se instead of a combination of diffusion and microbial uptake, because diffusion has already been accounted for. We found no empirical reason to assume that microbial uptake of COS to be zero at 215 high soil water content, if it is not diffusion-limited.
Soil COS production is represented as an exponential function of temperature, where V SP,max (mol m −3 s −1 ) is the soil COS production capacity, k T (K −1 ) is the temperature dependence factor that is equivalent to Q 10 = 1.9, and T ref = 25 • C. The Q 10 value here is set as 220 the average value of the Q 10 associated with the two regression lines of COS flux vs. temperature in Maseyk et al. (2014). Lab incubation of various soils has found that soil COS emissions generally exhibit exponential responses with temperature (Whelan et al., 2015). Although COS may be produced by both soil microbes and roots (Maseyk et al., 2014), we use a single production term in Eq. (24) due to a lack of data on the individual temperature responses of the different components.

Litter fluxes and litter properties
Litter was present at one of the model evaluation sites (Stunt Ranch). We obtained an equation for litter COS fluxes based on an incubation experiment with litter at the Stunt Ranch site (Sun et al., 2015),

230
where V LU,max (mol m −3 s −1 ) and V LP,max (mol m −3 s −1 ) are the litter COS uptake and production capacities, k T (K −1 ) is the temperature dependence factor for COS production equivalent to Q 10 = 1.9, and k L (dimensionless) is the moisture limitation factor for COS uptake. The hyperbolic sine function is chosen to describe the exponential dependence on litter water content.
Litter porosity is usually much larger than that of soils. In the model, the porosity at the grid 235 point near the litter-soil interface is interpolated so as to prevent numerical instability caused by discontinuity (Fig. 4).
Temperatures of the litter layers were interpolated between soil and chamber air temperature by assuming a logarithmic temperature profile in the surface layer, according to mixing length theory (Prandtl, 1925). The roughness length was assumed to be 0.01 m (Stull, 1988

Site-specific parameters
Site-specific parameters are summarized in Table 2. We used site-specific values for soil porosity 255 that were constant with depth except in the top few centimeters when there is loose topsoil. Soil porosity is estimated as 0.50 m 3 m −3 at the SGP site from the maximum observed soil moisture (Fischer et al., 2007;Maseyk et al., 2014), but that in the top 2 cm is assumed to be 0.60 m 3 m −3 as the topsoil is looser due to agricultural activity. At the SR site, soil porosity was measured as 0.35 m 3 m −3 . We assume the field capacity (θ FC ) is 0.20 m 3 m −3 at the SR site based on soil textures 260 (Or and Wraith, 2002). For the parameter b in Eq. (13), we use 4.9 for SR (sandy loam) and 5.3 for SGP (silt loam) based on Clapp and Hornberger (1978).
Soil temperature profiles are modeled with the observations at 5 cm depth. The temperature signals are considered as a superposition of fast varying diurnal signals (T F ) and slowly varying synoptic scale signals (T S ) (Van Wijk and de Vries, 1963), where z T is the damping depth for diurnal temperature waves, ω = 2π day −1 is the angular frequency, ψ is the phase constant. z T is determined from, where α T (m 2 s −1 ) is the soil thermal diffusivity, calculated from soil mineral composition and water 270 content using empirical formulae in Peters-Lidard et al. (1998) and Johansen (1975). Soil thermal diffusivities are 2.5-5.0 × 10 −7 m 2 s −1 at the SGP site, and 6.8-8.1 × 10 −7 m 2 s −1 at the SR site, assuming typical mineral compositions. The thermal diffusivities are translated to average damping depths of 0.11 m (SGP) and 0.14 m (SR), respectively.
For the SR data, we observed that the temperature optimum for COS uptake is T opt ≈ 13 • C (Sun 275 et al., 2015). This is reproduced by setting T eq at 15 • C in Eq. (22). For the SGP soil, because COS flux did not show a temperature optimum for uptake during the observed temperature range (7-47 • C) in Maseyk et al. (2014), we set T eq = 10 • C (equivalent to T opt ≈ 7 • C) to get a monotonic relationship with temperature.
Soil moisture was measured at 5 cm depth at both sites, and additionally at 30 cm at SGP. We 280 generate soil moisture profiles for the SGP site by interpolating between the 5 and 30 cm data, and assuming constant soil moisture below 30 cm (cf. 1D simulations with gravity drainage boundary condition in Walker, 1999). For the SR site, we assume soil moisture at 1 m depth is at field capacity (θ FC , m 3 m −3 ), and interpolated smoothly between the values at the surface and 1 m depth (Fig. 4).
The conditions below 10 cm do not significantly affect surface COS flux according to Fig. 11 (see 285 also section 4.2). The optimum water content for soil COS uptake, w opt , is assumed to be 0.20 m 3 m −3 for SGP soils, corresponding to a threshold for the COS flux-temperature response (Maseyk et al., 2014). For SR soils, we set w opt = 0.14 m 3 m −3 , estimated from lab incubation data of the same soils (Whelan et al., 2015).
Since litter was present at the SR site, litter layers are included in model simulations when val-290 idating the model with the SR dataset. The 2 cm thick litter layer is represented in the top 6 grid points. Litter porosity was measured to be 0.94 m 3 m −3 . Litter water content (LWC) is an important variable controlling litter COS fluxes. Since LWC was not measured at SR, we test two scenarios: (1) fast decreasing LWC, and (2) slowly decreasing LWC following the precipitation event on the day before the measurements started. We also evaluate the surface COS flux in the absence of litter 295 fluxes.

Evaluation results
In the evaluation, soil COS uptake and production capacities (V SU,max and V SP,max ) are chosen to fit the general trend of observed data. The aim here is to demonstrate that the model can reproduce the data with the proper settings, and thus has the potential to be applied at large scales. Model 300 performance in the evaluation tests is summarized in Table 3.

Evaluation against data from the wheat field soil (SGP)
The simulated COS fluxes at the SGP site are generally in good agreement with the observations (Fig. 5), especially the diurnal variations and the transition from uptake to emissions. Model results and the data lie on a 1:1 line, albeit with large scatter at high emission fluxes (Fig. 5c). As shown in 305 the distribution of model residuals with respect to surface soil temperature and moisture (Fig. 5d), the model tends to underestimate the high emissions at high temperatures. In the growing season (before DOY 130), the model captures the high uptake instances, but tends to underestimate high emission instances at mid-day. This is mainly because some high emissions are not associated with high soil temperature. During the senescence and post-harvest period (after DOY 134), the model 310 reproduces the strong diurnal variations, except that the mid-day peak emissions are underestimated.
Examples for COS profiles during overall soil COS uptake and emissions show exponential decrease and increase of the COS concentration in soil air with depth, respectively (Fig. 7a).  (Table 3) given that there were only a few incubation experiments.
The observed COS emissions were higher in the senescence and post-harvest stages (after DOY 134) than during the growing season. However, when the temperature dependence is considered, the COS production capacity in the growing season (before DOY 130) must be higher than after harvest 325 to account for the high fluxes under relatively low soil temperatures (< 20 • C) in this period ( Fig. 5a and Table 3). We hypothesize that root production of COS may have peaked during the late growing season following increased root sulfate uptake during seed ripening (Zhao et al., 1999), and then decreased during the following senescence stage.
Surprisingly, to simulate the strong diurnal variability of COS emissions in the senescence and 330 post-harvesting stages, the high emissions need to be counteracted by continuing COS uptake, with the same uptake capacity as during the growing season. Without uptake activity, COS would accumulate in the soil column from high production in the daytime, and still exhibit high emissions at night. An additional contribution to daytime emissions could have come from photochemical production of COS, as observed for SGP soil in lab incubations (Whelan and Rhew, 2015). However, 335 this would only be a minor component since the maximum photochemical production rate translated to a per area basis was around 2 pmol m −2 s −1 , less than 10% of the observed diurnal variability. In addition, photochemical production would not have affected the fluxes measured inside the chamber because it was opaque. However, it may have affected the fluxes measured at the ecosystem scale (Billesbach et al., 2014;Maseyk et al., 2014). Hence, photochemical production will be included in 340 a future model version.

Evaluation against data from the oak woodland soil (SR)
The observation period at SR started with a rain event, and thus at high LWC. We set the starting LWC to 0.32 g g −1 to fit the observed high uptake fluxes, and simulated a rapid water loss (Fig. 8b) due to strong evaporation at typically warm and sunny conditions. The LWC was assumed to de-345 crease to 0.06 g g −1 following an exponential decay function. This assumption fits well with the observed fluxes, both for the diurnal variations and the overall trend (Fig. 6). As an alternative, we used a scenario with slowly decreasing LWC (Fig. 8b), or assuming no litter flux at all. In the first case, the uptake during the dry period was overestimated, whereas the no-litter flux scenario did not reproduce the general trend, even after soil uptake capacity was increased by 10-fold to compensate 350 for the missing litter uptake (Fig. 8a). The soil moisture change in this period was not strong enough to drive the overall trend from strong uptake to weak uptake.
The simulated COS profiles show that the presence of litter layers on top of the soil reduces the COS supply available to the soil, especially when there is strong uptake in the litter layers, thus limiting the contribution of soil uptake to the overall surface uptake (Fig. 7b). When LWC is low and 355 litter uptake is weak, this limiting effect is not significant.

Sensitivity tests: diffusion-controlled uptake
Laboratory studies have found that soil COS uptake has an optimum at 19-30% water-filled pore space fraction (WFPS) and approaches zero as WFPS becomes higher (Kesselmeier et al., 1999;360 Van Diest and Kesselmeier, 2008). We conducted a set of idealized simulations to find what factors control such dependence. Based on model results, the shape of the surface flux vs. soil moisture curve is the net result of two competing effects: with increasing soil moisture, the higher microbial activity is counteracted by stronger diffusion limitation. The resultant surface COS uptake has an optimum value of 20% WFPS, much lower than the optimum value of microbial uptake per se (Fig. 9). The 365 shape of the microbial activity curve is not important at high soil water content where the diffusion limitation is dominant.
We also show that for diffusion-controlled, concentration-dependent uptake, the surface flux does not increase linearly with soil COS uptake capacity (V SU,max ), but follows a logarithmic relationship (Fig. 10). Hence, if we can constrain the range of soil COS uptake capacity from laboratory experiments, the uncertainty in fluxes arising from this parameter would not be dominant. The non-linear increase also shows that empirical methods that assume a linear response of fluxes to soil COS uptake capacity would not give accurate results.

Advantages of a depth-resolved model
One of the main advantages of using a depth-resolved model is that it enables the analysis of changes 375 in uptake and production capacities over time, or between sites. For example, an unexpected finding was that the production capacity parameter at SGP needed to be decreased after the senescence stage (from 1 × 10 −10 to 3.3 × 10 −11 mol m −3 s −1 , Table 3). Surface COS emissions strongly increased after senescence, but were associated with much higher temperatures. The post-senescence decrease in production capacity indicates that soil COS production may be related to plant activity. 380 We also found that during the period of strong surface COS emissions after the harvest at SGP, the pronounced diurnal variations required the same soil COS uptake capacity as during the growing season, indicating that soil COS uptake does not directly depend on plant activity. Soil COS uptake could be partly abiotic (Liu et al., 2008(Liu et al., , 2010b or largely due to microbial activity (Kesselmeier et al., 1999), but even in the latter case may not respond quickly to changes in plant cover.

385
The uptake and production parameters that best fit at each site and during distinct phenological periods can be obtained by optimization procedures used with a particular soil dataset, another significant advantage of a depth-resolved model. Given reasonable initial guesses of the uptake and production parameters, the optimization runs iteratively over the data with the gradient descent or Newton's method until the minimal sum of square errors is attained. An example of data optimization 390 applied to the SR dataset is in Sun et al. (2015).
We found that soil COS uptake is largely determined by activity in the top 10 cm of the soil. For each layer, we calculated how much the surface uptake would increase as a result of increasing the COS uptake capacity (V SU,max ) in that layer by a factor of 10 (Fig. 11). Surface uptake is sensitive to V SU,max changes in the top 10 cm, even though the layer sizes are the smallest, but not sensitive 395 to changes below 10 cm. Thus, physical parameters in the top layers of the soil are most important for soil COS fluxes. On the other hand, it also means that information on sources and sinks of COS in the deeper soil cannot be obtained from surface flux measurements.

Which biome-specific parameters are needed for a global simulation?
The COS flux model can be integrated into more comprehensive land surface models (e.g. CLM4.5 400 and SiB3) to simulate global fluxes. Soil temperature and moisture are usually generated from prognostic equations in these models. Litter layers would need to be added as they are often not represented separately but embedded in soil layers as soil organic carbon (e.g. in CLM4.5). Several parameters may be biome-specific: (1) the uptake and production capacities of COS (V SU,max , V SP,max ); (2) the optimum temperature and moisture of soil COS uptake (T opt , w opt ); (3) the Q 10 parameter 405 of COS production from microbial, root, or abiotic sources; and (4) the function form of litter COS uptake.
The V max parameters should be controlled by microbial biomass and enzyme content in the soil, which in turn are related to soil physical state and organic carbon content. These parameters can be estimated by optimizing model results against field data. The T opt and w opt parameters can be 410 determined from laboratory experiments, and extended to soils in similar environments, considering that they are likely related to local climate regimes and soil hydrology. Root COS production as a function of temperature needs to be determined from laboratory experiments. More studies are also needed to constrain litter COS fluxes, and how they change with litter type, level of litter degradation, and litter moisture content and temperature. With more field observations and lab incubations, we 415 expect a database of these parameters in different biomes will be constructed for regional and global simulations of soil COS fluxes.

Conclusions
This study presents a mechanistic diffusion-reaction model coupling physics and biogeochemistry to simulate soil COS flux, as well as its evaluation with field data. The model explicitly accounts 420 for diffusion in the soil column, COS production, and COS exchange in the litter layer. The model reproduced well observed fluxes at two sites, and has enabled us to gain novel (and unexpected) insights such as the higher COS production capacity at SGP during the growing season despite lower soil COS emissions, and the continuing COS uptake capacity required to partly counteract the large COS emissions at SGP after harvest. We also demonstrate that diffusion must be considered to accu-425 rately simulate surface COS fluxes. Diffusion control on surface uptake is evident in its sensitivity to soil water content, and the sub-linear response of uptake flux to soil COS uptake capacity. For largescale simulations of soil COS fluxes, further lab and field studies are needed to establish a database of soil and litter COS uptake and production capacities and parameters for typical soil types across key biomes.
This study with data kH,cc = (T /K) exp α + β T /K mol mol −1 α = −20.00 from Elliott et al. (1989) β = 4050   VLU,max 1.68 × 10 −9 mol m −3 s −1 Litter COS uptake capacity Sun et al. (2015) VLP,max 1.33 × 10 −11 mol m −3 s −1 Litter COS production capacity Sun et al. (2015) VSU,max mol m −3 s −1 Soil COS uptake capacity VSP,max mol m −3 s −1 Soil COS production capacity Temperature (K) Henry's law constant (dimensionless) Wilhelm et al. (1977) De Bruyn et al. (1995 This study, fit with data from Elliott et al. (1989) pH < 9, no NaCl pH ≥ 9, no NaCl 0.5 M NaCl solution k H =Texp(4050.32/T−20.0007) Figure 2. COS solubility function used in this study (pink), obtained from regression with Elliott et al. (1989) data (pink). Solubility data measured under pH < 9, non-saline conditions (pink diamonds) were used for the regression. Also plotted are the COS solubility functions from Wilhelm et al. (1977) andDe Bruyn et al. (1995).  Note that in (c) and (d), the profiles also show the porosity in the litter layers and litter water content (g g −1 ) in red. The dashed lines denote litter-soil interfaces, and the light gray lines represent air-litter interfaces (same for Fig. 7). One computational node near litter-soil interface is interpolated for porosity in order to prevent numerical instability. Measured data of soil moisture at 5 and 30 cm for SGP, and 5 cm for SR, are shown as red crosses.    COS uptake activity is included in these simulations. Also shown are the two model variables that dominate the simulated flux response, the effective soil COS diffusivities (blue) and microbial COS uptake turnover rates (green) at the same soil temperatures. Note the y-axis of soil COS flux is reversed to show the uptake. [COS] atm =500 pptv Figure 11. Increase of surface COS uptake (in percentage) from 10-fold changes of VSU,max at each depth.
The black line shows the changes for an uptake-only case (VSU,max = 10 −8 mol m −3 s −1 ), while the red line shows the changes when both uptake and source activities are included (VSU,max = 10 −8 mol m −3 s −1 , VSP,max = 10 −11 mol m −3 s −1 , but the soil is still a sink). Soil conditions are the same as for Fig. 10.