Evaluation of carbonyl sulfide biosphere exchange in the Simple Biosphere Model (SiB4)

. The uptake of carbonyl sulfide (COS) by terrestrial plants is linked to photosynthetic uptake of CO 2 as these gases partly share the same uptake pathway. Applying COS as a photosynthesis tracer in models requires an accurate representation of biosphere COS fluxes, but these models have not been extensively evaluated against field observations of COS fluxes. In this paper, the COS flux as simulated by the Simple Biosphere Model, version 4 (SiB4) is updated with the latest mechanistic 35 insights and evaluated with site observations from different biomes: one evergreen needleleaf forest, two deciduous broadleaf forests, three grasslands, and two crop fields spread over Europe and North America. We improved SiB4 in several ways to 2 improve its representation of COS. To account for the effect of atmospheric COS mole fractions on COS biosphere uptake, we replaced the fixed atmospheric COS mole fraction boundary condition originally used in SiB4 with spatially and temporally 40 varying COS mole fraction fields. Seasonal amplitudes of COS mole fractions are ~50-200 ppt at the investigated sites with a minimum mole fraction in the late growing season. Incorporating seasonal variability into the model reduces COS uptake rates in the late growing season, allowing better agreement with observations. We also replaced the empirical soil COS uptake model in SiB4 with a mechanistic model that represents both uptake and production of COS in soils, which improves the match with observations over agricultural fields and fertilized grassland soils. The improved version of SiB4 was capable of 45 simulating the diurnal and seasonal variation of COS fluxes in the boreal, temperate and Mediterranean region. Nonetheless, the daytime vegetation COS flux is underestimated on average by 8 ± 27 %, albeit with large variability across sites. On a global scale, our model modifications decreased the modelled COS terrestrial biosphere sink from 922 Gg S yr -1 in the original SiB4 model to 753 Gg S yr -1 in the updated version. The largest decrease in fluxes was driven by lower atmospheric COS mole fractions over regions with high productivity, which highlights the importance of accounting for variations in atmospheric 50 COS mole fractions. The change to a different soil model, on the other hand, had a relatively small effect on the global biosphere COS sink. The secondary role of the modelled soil component in the global COS budget supports the use of COS as a global photosynthesis tracer. A more accurate representation of COS uptake in SiB4 should allow for improved application of atmospheric COS as a tracer of local to global-scale terrestrial photosynthesis. et al., 2007) were used to optimize the sources and sinks of COS. We used global 2D surface layer fields of COS 315 mole fractions resulting from these optimized sources and sinks, as they were optimally consistent with the available COS observations for the period 2016-2018. mole to lower COS uptake rates, but would in turn also lead to a smaller drop in COS mole fractions; this feedback is currently not 320 accounted for.

improve its representation of COS. To account for the effect of atmospheric COS mole fractions on COS biosphere uptake, we replaced the fixed atmospheric COS mole fraction boundary condition originally used in SiB4 with spatially and temporally 40 varying COS mole fraction fields. Seasonal amplitudes of COS mole fractions are ~50-200 ppt at the investigated sites with a minimum mole fraction in the late growing season. Incorporating seasonal variability into the model reduces COS uptake rates in the late growing season, allowing better agreement with observations. We also replaced the empirical soil COS uptake model in SiB4 with a mechanistic model that represents both uptake and production of COS in soils, which improves the match

Introduction 55
Carbonyl sulfide (COS) uptake by the terrestrial biosphere is the main sink of atmospheric COS (Whelan et al., 2018). COS uptake in plants is closely related to photosynthetic CO2 uptake because of its shared uptake pathway through plant stomata and, as a consequence, COS can be used to help constrain the terrestrial carbon and water cycles (Seibt et al., 2010;Stimler et al., 2010;Whelan et al., 2018). Key plant processes such as photosynthesis and transpiration are difficult to observe at scales larger than the leaf level because they are contained within the net CO2 flux and evapotranspiration and are not separable from 60 other fluxes. Constraints on these fluxes at larger spatial scales are therefore needed to improve terrestrial biosphere models to better simulate the responses of photosynthesis and stomatal gas exchange to a changing climate. Recently, COS has been shown to be valuable for understanding changes in plant uptake, e.g., the inhibition of photosynthesis during a heat wave (Wohlfahrt et al. 2018), the growth of the terrestrial gross primary production (GPP) during the twentieth century (Campbell et al., 2017), the regional-scale partitioning of NEE into GPP and ecosystem respiration (Hu et al., 2021) and changes in between 368 and 1845 Gg S yr -1 with a mean of 1084 Gg S yr -1 over 9 different studies as summarized in Table 1 (Kettle et al. 2002;Montzka et al. 2007;Suntharalingam et al. 2008;Berry et al. 2013;Launois et al. 2015a;Kuai et al. 2015;Wang et al. 2016;Maignan et al. 2020;Ma et al. 2021). These estimates were made through different approaches, such as scaling COS vegetation uptake to the net (NPP) or gross primary production (GPP), and more recently also from mechanistic implementations ( Table 1). The mechanistic implementations of COS vegetation uptake in the biosphere models yield a smaller 85 range of 688-775 Gg S yr -1 than when the COS vegetation uptake is scaled to the CO2 vegetation sink (Table 1). The global soil COS sink estimates range from 130 to 510 Gg S yr -1 , but with most estimates between 130 and 176 Gg S yr -1 . However, until now, land surface models have still not adopted the available mechanistic soil models from either Sun et al. (2015) or Ogée et al. (2016).

90
The temporal and spatial variability of atmospheric COS mole fractions have a considerable influence on the COS biosphere uptake (Ma et al. 2021) because the COS plant uptake is governed by a first order kinetic process (Stimler et al., 2010); that is, COS uptake is linearly proportional to the atmospheric COS mole fraction. A typical seasonal amplitude of atmospheric COS mole fractions of ~100-200 parts per trillion (ppt) around an average of ~500 ppt affects the fluxes by ~20-40% even if stomatal conductance remained constant. In contrast to CO2, where a seasonal amplitude of ~6-7 ppm around ~410 ppm could 95 affect the fluxes only by ~2 %. Although some previous studies have considered the impact of variable COS mole fraction when the biosphere flux was introduced to an atmospheric transport model (Berry et al., 2013, Kuai et al., 2015, Wang et al., 2016, it has not yet been adopted as a standard approach (Maignan et al. 2021, Ma et al., 2021. Inverse modeling studies that account for all known sources and sinks of COS imply a missing source of COS in the Tropics 100 (Berry et al. 2013;Le Kuai et al. 2015;Ma et al. 2021). Ma et al. (2021) revealed considerable seasonal variations of the missing source. Yet the exact reason for this missing source has not been resolved. Although the missing source can be anthropogenic or from the tropical ocean (Launois et al. 2015b;Kuai et al. 2015;Lennartz et al. 2017Lennartz et al. , 2019, an overestimated tropical biospheric sink is also possible. Moreover, Ma et al. (2021) identified a missing sink at the higher northern latitudes that required uptake larger than in the inversion a priori model (i.e. SiB4). This missing sink could be explained by an 105 underestimated biosphere sink, and would be equivalent to a 6 % underestimation of the biosphere sink north of 30 o N (Ma et al., 2021).
A source of uncertainty for COS uptake by land surface models is that simulations have not been extensively compared against field observations because field measurements of ecosystem and soil fluxes are sparse. Yet, several research groups have 110 performed field observations of COS ecosystem fluxes in the last decade (Asaf et al. 2013;Maseyk et al. 2014;Commane et al. 2015;Kooijmans et al. 2017;Wehr et al. 2017;Yang et al. 2018;Spielmann et al. 2019;Berkelhammer et al. 2020;Vesala et al., 2021) with observations covering evergreen needleleaf forests, deciduous broadleaf forests, grasslands, and crop fields.
These experimental efforts now offer the possibility to compare model simulations of COS biosphere exchange against field observations from different biomes. 130 In this paper, we compare these field measurements with the latest version of the SiB model, version 4 (SiB4) and evaluate the calculated global COS biosphere flux. When compared to SiB3 (Berry et al. 2013), SiB4 simulates variable carbon pool allocation, prognostic phenology, land cover heterogeneity, and crop phenology (Haynes et al. 2019a). We evaluate seasonal and diurnal cycles of ecosystem COS fluxes and the representativeness of nighttime COS uptake, where the latter is important 135 for an accurate COS sink estimate. We furthermore update the SiB4 model with knowledge obtained on soil exchange of COS during the last decade by implementing the mechanistic soil model from Ogée et al. (2016) for COS soil uptake and production rates varying with biome after Meredith et al. (2018Meredith et al. ( , 2019. Furthermore, we replace the fixed atmospheric COS mole fraction of 500 pmol mol -1 (a nominal background tropospheric mole fraction) with spatially and temporally varying COS mole fraction fields obtained from an inversion with the TM5-4DVAR atmospheric transport model (Ma et al., 2021

SiB4 model
The Simple Biosphere Model version 4 (SiB4) is a mechanistic, prognostic land surface model that integrates heterogeneous 170 land cover, environmentally responsive prognostic phenology, dynamic carbon allocation, and cascading carbon pools from live biomass to surface litter to soil organic matter (Haynes et al. 2019a,b). SiB4 predicts vegetation and soil moisture states, land surface energy and water budgets, and the terrestrial carbon cycle. Rather than relying on satellite data to specify phenology as in SiB3, SiB4 simulates the terrestrial carbon cycle by using the carbon fluxes to determine the above and belowground biomass, which in turn feeds back to impact carbon assimilation and respiration . SiB4 175 predicts plant phenology, divided into different stages, allowing the change of photosynthetic activity over seasons through specified maximum RuBisCO velocities in each phenological stage. To classify land surface vegetation, SiB4 uses plant functional types (PFTs), which group plants according to their function and physical, physiological, and phenological characteristics. In addition to nine non-crop vegetation PFTs, SiB4 includes three specific crops (maize, soybeans, and winter wheat), and two generic crops (C3 and C4) following the crop phenology model developed by Lokupitiya et al. (2009). SiB4 180 includes land cover heterogeneity by simulating multiple PFTs per grid cell.

COS plant and soil uptake
COS plant uptake in the SiB4 model is based on the formulation of Berry et al. (2013) and is simulated as a series of conductances (gt) from the leaf boundary layer to the site of COS hydrolysis in the mesophyll cells. These conductances include the conductance from canopy air to the leaf surface, or leaf boundary layer conductance (gb), the stomatal conductance (gs), 185 and the internal conductance (gcos). The latter represents both the diffusion of COS to the mesophyll cells and the efficiency of the leaf mesophyll carbonic anhydrase (CA) to hydrolyze COS. This leads to the following equation for the COS uptake rate by vegetation: (1) where Fcos,veg is the COS vegetation uptake rate (pmol m -2 s -1 ) and Ca is the COS mole fraction in the canopy air space (pmol mol -1 ) calculated from the mixed layer COS mole fraction (standard 500 pmol mol -1 , but see Sect. 2.1.3.) taking into account 190 uptake of COS by soil and vegetation in the previous timestep. gs and gb are the respective stomatal and boundary layer conductances for water vapor (mol m -2 s -1 ) and are scaled to account for the different diffusivity rates of COS and H2O (Seibt et al., 2010;Stimler et al., 2010). The stomatal conductance gs is derived following the Ball-Berry photosynthesis-conductance model as modified by Collatz et al. (1992) and gb follows the formulations described by Sellers et al. (1996). The internal conductance gcos is assumed to scale with maximum carboxylation rate of RuBisCO, Vmax (µmol m -2 s -1 ) (Berry et al. 2013), 195 inspired by previous findings that both CA activity (Badger and Price, 1994), and mesophyll conductance (Evans et al., 1994) (Sellers et al., 1992): 8+9 = ",8+9:; 2.1 =.>(@ ABC D:EF.=) , (2) gcos is then described as function of VCOS that represents the CA activity: (3) where FLC is a factor scaling the flux from a single leaf to the canopy that considers the canopy profile of absorbed photosynthetically active radiation (Sellers et al., 1996), FRZ is the rootzone water potential, an empirical scaling factor that reduces the biochemical activity when little soil moisture is available (e.g. during extended periods of drought), p/p0sfc adjust 210 the fluxes for altitude, where p is atmospheric pressure (hPa) and p0sfc the reference surface pressure (1000 hPa), and Tcan/T0 scales the flux to a reference temperature at T0 = 273.15 K. A calibration term a was added to scale gcos to COS flux observations of controlled gas exchange measurements (Stimler et al., 2010(Stimler et al., , 2011, which resulted in a = 1200 for C3 and 13000 for C4 species (Berry et al. 2013). These numbers were later updated to a = 1400 and 8862 for C3 and C4 species, respectively after reanalysis of the gas exchange data. Berry et al. (2013) already noted that the a value did not constrain the 215 variability between plant species well, likely due to plant variability in CA activity and/or differences in mesophyll conductance. In Sect. 2.3 we explain how we use field measurements to explore whether we can refine a values for different plant functional types separately and to make it variable over time.
The enzyme CA is expressed in microbial communities in soils as well, leading to COS uptake by soils (e.g. Kesselmeier et 220 al., 1999;Meredith et al. 2019). In SiB4, COS uptake in soils (hereafter called "the Berry soil model") is coupled to heterotrophic CO2 respiration under the assumption that in more productive regions there would be more litter and surface soil carbon for respiration, and these richer carbon environments would have more CA as well (Yi et al., 2007). Additionally, COS soil uptake in the model is regulated by diffusion, controlled by soil porosity and the fraction of water filled pore space (Van Diest and Kesselmeier, 2008;Ogée et al., 2016;Sun et al., 2015;Whelan et al., 2016). Initial implementations of soil COS 225 uptake made calculations for the entire soil column, but subsequent model versions considered only uptake in the top 20 cm of the soil (Wang et al. 2016), thereby decreasing global soil uptake estimates from 355 (Berry et al. 2013) to 159 Gg S yr -1 (Wang et al. 2016). In the next section, we describe our update to the SiB4 model based on advances in our knowledge on COS soil exchange obtained during the last decade.

Mechanistic COS soil model
Field and laboratory experiments in the last decade showed that COS is not only taken up by soil but is also produced due to abiotic thermal degradation and photodegradation of soil organic matter and is especially enhanced in agricultural soils 235 (Maseyk et al. 2014;Whelan and Rhew 2015;Meredith et al. 2018;Kaisermann et al. 2018a). Besides COS soil production being enhanced in fertilized soils, COS uptake was shown to be diminished in fertilized soils (Kaisermann et al. 2018b). These effects of fertilization on soil COS exchange were initially not simulated in the SiB4 model.
New empirical soil models (Whelan et al., 2016) and mechanistic models (Ogée et al., 2016;Sun et al., 2015) were developed 240 during the last decade. The mechanistic models describe the uptake and production pathways together with COS diffusion in a soil column. Ogée et al. (2016) derived a simplified analytical solution assuming a soil column with uniform temperature, soil moisture, and porosity and steady state conditions for comparison against laboratory measurements. The model from Ogée et al. (2016), hereafter called "the Ogée soil model", was then used by several laboratory studies to study patterns in uptake and production of COS in soils (Meredith et al. 2018;Kaisermann et al. 2018a,b). Due to these efforts, there are now 245 reaction rate parameter values available for a range of biomes and land use types. Because these reaction rate values were derived by fitting the Ogée soil model on data from mesocosm experiments, they should be used in combination with this model to estimate ecosystem-scale soil COS fluxes. Also, compared to the COS soil model proposed by Sun et al. (2015), the steady state solution of the Ogée soil model is computationally inexpensive and therefore more suitable for implementation in SiB4 for global COS soil flux calculations. In the following paragraphs we describe the implementation of the Ogée soil model 250 in SiB4.
For field conditions (assuming a zero COS vertical gradient at the bottom of the soil layer and steady state), the COS soil flux (mol m -2 s -1 ) calculation simplifies to (Ogée et al., 2016): where k is the CA reaction rate (s -1 ), B (m 3 water m -3 air) the solubility of COS in water that relates to Henry's law constant 255 and depends on temperature, θ the soil water content (m 3 m -3 ), D the soil COS diffusivity (m 3 air m -1 soil s -1 ), Ca the COS mole fraction at the soil-air interface, z1 2 = D/(kBθ), and P the COS production rate (mol m -3 s -1 ) uniform over depth zp (here assumed to be 1.0 m). For details of the model calculations we refer readers to Ogée et al. (2016), here we provide the information specific for the implementation in SiB4. We assume Ca to be identical to the COS mole fraction in the canopy air space. While implementing and testing the model we recognized the strong dependence of the soil fluxes on soil porosity, choice of tortuosity 260 functions, and the SiB4 soil layer selected for temperature and soil moisture. For the calculation of D we used the SiB4 soil porosity (m 3 m -3 ; calculated from sand fractions following Lawrence and Slater (2008)) that accounts for the volume of ice in the soil. The simulated soil water content and soil temperature are taken from the top 5 cm soil layer, where most of the COS uptake takes place. D also depends on tortuosity functions that describe the tortuous movement through the air-or water-filled pore space. Several tortuosity functions are described in the literature and also Ogée et al. (2016) acknowledged that the response of the soil COS fluxes to soil moisture varied with the chosen tortuosity functions. We chose the tortuosity functions of Deepagoda et al. (2011) for air and Millington and Quirk (1961) for water, as these functions do not depend on a pore-size 270 distribution parameter, which facilitates its implementation in SiB4.
COS is taken up in soils through hydrolysis in soil water, where the main consumption is enzymatic, and thus depending on soil CA enzyme activity. Here, and following other studies (i.e. Ogée et al. 2016, Meredith et al. 2019, we expressed the CA reaction rate k relative to the uncatalyzed reaction rate (kuncat) at a reference temperature (Tref) and pH (here assumed constant 275 at 4.5): where fCA is the CA enhancement factor, kuncat varies with soil pH according to Elliott et al. (1989) (Table 2). Here, the fCA for agricultural soils is substantially smaller than that of other vegetated biome types, thereby including the reduced COS uptake in fertilized (agricultural) soils (Kaisermann et al. 2018b).

285
The COS production was defined by Ogée et al. (2016) as a temperature response function modulated by the soil redox potential. Meredith et al. (2018) also measured COS production at a temperature range between 10 and 40 °C for a range of biomes. Measurements were then fitted to an exponential model: We used this exponential temperature model and the biome averaged a and b (Table 2) in our calculation of P in SiB4. We assume here that the controlled laboratory measurements by Meredith et al. (2018; can be used to estimate soil fluxes 290 under field conditions. The higher value of a for agricultural soils (Table 2) allows for higher COS soil production in this soil type. Ideally, we would include production parameter values for wetland soils so that we take into account the typically large production that has been observed in wetland soils (Meredith et al., 2018;Whelan et al., 2013). However, SiB4 does not discriminate between oxic and anoxic (wetland) soils, which precluded the implementation of wetland-specific COS soil production. 295

Variable atmospheric COS mole fractions 305
The atmospheric COS mole fraction in the planetary boundary layer affects both the COS vegetation and soil flux calculations (Eq. (1) and (4)). In SiB4 a standard constant "place-holder" COS mole fraction of 500 pmol mol -1 is used. Ma et al. (2021) estimated that the global biosphere sink would decrease from 1053 Gg S yr -1 to 851 Gg S yr -1 if the fixed COS mole fraction were replaced with monthly mean fields that account for the drawdown of COS near the surface in the peak growing season.
We thus changed the prescribed COS mole fraction from a fixed value to one varying in space and time, including seasonal 310 and diurnal variability. To this end we used the surface COS mole fraction fields at a global resolution of 4° × 6° (latitude × longitude) at 3-hourly time steps as retrieved from an atmospheric transport inversion performed with TM5-4DVAR by Ma et al. (2021) using the chemistry transport model TM5 in which COS exchange was recently implemented. Measurements of atmospheric COS mole fractions at 14 sites from the National Oceanic and Atmospheric Administration (NOAA) flask network (Montzka et al., 2007) were used to optimize the sources and sinks of COS. We used global 2D surface layer fields of COS 315 mole fractions resulting from these optimized sources and sinks, as they were optimally consistent with the available COS flask observations for the period 2016-2018. We repeated the average over those years as input to the SiB4 mixed layer COS mole fraction for each year in the simulation (see global maps of monthly mean surface COS mole fractions in supplementary figure S13 of Ma et al. (2021)). In the inversion of Ma et al. (2021), the changing (e.g. lower) COS mole fractions would lead to lower COS uptake rates, but would in turn also lead to a smaller drop in COS mole fractions; this feedback is currently not 320 accounted for.

Simulations
We used meteorological data from the Modern-Era Retrospective Analysis for Research and Applications, version 2 330 (MERRA2), which are available from 1980 onwards (Gelaro et al., 2017) as meteorological forcing to SiB4. To ensure realistic rainfall, the convective and large-scale precipitation values were scaled such that the monthly total rainfall matches with the monthly precipitation in the Global Precipitation Climatology Project, Version 1.2 product (Huffman et al. 2001;Baker et al. 2010;Haynes et al. 2019a,b). Up to 10 PFTs per grid cell (at 0.5° × 0.5° resolution) are prescribed following PFT maps based on MODIS data (Lawrence and Chase, 2007). The soil characteristics such as sand fraction (used for the calculation of soil 335 porosity) are provided by the International Geosphere-Biosphere Programme (IGBP) Global Soil Data Task Group (2000).
We ran SiB4 from 2000 to 2020, and the simulations were preceded by a spinup iterating five times over the years 2000-2020 using an accelerated equilibrium approach (Haynes et al., 2019b) to initialize the carbon pools to reach steady state. CO2 mole fractions were held constant at 370 µmol mol -1 during spinup and simulations. Research is ongoing to implement an accurate 340 representation of the effect of CO2 fertilization in SiB4. We performed two sets of four simulations (global and site level) with the same driver data and settings, but with a different temporal resolution of the output: 1. For global simulations, we used monthly averaged output. Moreover, SiB4 simulates multiple PFTs per grid cell. These were averaged, weighted by the fraction of land area occupied by each PFT; 2. To compare SiB4 with site observations (listed in Table 3), we run the SiB4 model with 3-hourly output for only the grid cells (at 0.5° x 0.5° resolution) in which the sites are located. For comparison with 345 observations we selected the PFT that best represents the measurement site. For these site comparisons we use MERRA2 meteorological data (instead of local meteorological observations) to provide consistency in data collection, availability and application across sites, and for consistency with the global run.
We run SiB4 with four different configurations: 350 1) the original SiB4 model containing the standard COS mole fraction of 500 pmol mol -1 and the Berry soil model (SiB4_500_Berry); 2) the Ogée soil model and the standard COS mole fraction of 500 pmol mol -1 (SiB4_500_Ogee); 3) the Berry soil model and variable COS mole fractions (SiB4_var_Berry), and 4) the Ogée soil model and variable COS mole fractions (SiB4_var_Ogee). 355

Field observations
We use existing field observations for comparison with the SiB4 model simulations. Several studies have collected field data in the last two decades and we used those sites where continuous hourly measurements of ecosystem COS fluxes are available for at least a month. The site locations, some of their characteristics, and basic information on the observations are summarized in Table 3. The locations of the sites are indicated in Fig. S1. The measurements represent evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), maize (MAI), winter wheat (WWT), and C3 grasslands (C3-GRA), more specifically alpine grassland, prairie grassland, and savannah grassland. 365 All COS observations were made at FLUXNET, ICOS or AmeriFlux sites with the benefit that additional long-term measurements of CO2 and water exchange (Pastorello et al., 2020) are often available (see Table S1 for an overview), allowing  (Commane et al., 2015), but validated against COS flask measurements at the station, which are part of the NOAA flask measurement network (Montzka et al., 2007).
For further details about the site characteristics and measurement and processing procedures we refer to the original data 390 publications as reported in Table 3 and Table S1.
For evaluation of the model against observations we calculate the mean bias error (MBE; pmol m -2 s -1 ) and root mean square error (RMSE; pmol m -2 s -1 ) for monthly, daytime, and nighttime average fluxes.   Berry et al. (2013) used the calibration factor a to scale gcos to match the simulated COS vegetation flux with laboratory measurements. They noted that the a value did not constrain the variability between plant species well, likely due to plant 410 variability in CA activity and/or mesophyll conductance. Here, we derived aobs from COS field measurements; however, as aobs still contains several SiB4 simulated parameters, it is not strictly an observationally derived value. This analysis is meant to explore its variability over time and the necessity to define a values specific for different PFTs. We did not retain aobs for global simulations. This analysis requires that field measurements of ecosystem and soil fluxes are available. Under the assumption that both Vmax (and thus photosynthesis) and the soil flux are accurately simulated, the application of aobs would result in simulated COS vegetation fluxes that match with observations. 425

Results and discussion
First, we evaluate the SiB4 COS flux against observations (Sect. 3.1). The accuracy of the modeled ecosystem flux is controlled by several factors, such as the accuracy in leaf phenology, differences in accuracy of the daytime and nighttime COS vegetation flux, the accuracy of the soil flux (of both the Berry and Ogée soil models), and the sensitivity to atmospheric COS mole fractions. We discuss the role of each of these factors in the evaluation of SiB4 biosphere fluxes against observations. All 430 results are based on the standard a values of 1400 and 8862 for C3 and C4 species, respectively. We present COS fluxes relative to the atmosphere (i.e. negative values indicate uptake by the ecosystem). Next, we study the variability of aobs between different PFTs and across seasons (Sect. 3.2), to investigate remaining model-data mismatches in the COS vegetation flux that could potentially be solved by re-calibrating a. Finally, we present global estimates of the COS biospheric sink with different model configurations (Sect. 3.3). 435

Seasonal variability
Simulated COS ecosystem fluxes capture the seasonal variation of monthly-averaged observations (Fig. 1), with similar results for vegetation fluxes alone (Fig. S2). Specifically, simulated COS uptake peaked in summer, as was observed at the three sites that contain COS flux measurements across different seasons (Fig. 1a, d, e). At the other sites, COS ecosystem fluxes were 440 only measured during one part of the growing season. Therefore, we also used multi-year NEE, GPP, and latent heat flux (LE) from EC sites to evaluate the SiB4 phenology, which affects both CO2 and COS seasonality (Fig. S3-5).
Based on the NEE, GPP, and LE observations (Fig. S3-5), the start and end of the growing season are typically well captured by the SiB4 simulations. The timing and length of the growing season for grassland sites has been previously evaluated by 450 Haynes et al. (2019b) using remotely-sensed leaf area index and showed that SiB4 was capable of simulating growing season timing and variability across temperature and precipitation gradients. Also, the timing of maximum NEE and GPP, which differs by PFT and climatic regions, was well captured; e.g., simulated and observed NEE and GPP peak in spring at the savannah grassland site ES-LM1 and at the winter wheat site US-ARM. All other sites show an observed and simulated summer maximum carbon uptake. Only AT-NEU is an exception, with SiB4 predicting the peak net CO2 uptake too late into the 455 summer compared to the observations, which can be explained by grass cutting that was not included in SiB4. Crop harvesting was included in SiB4, but the exact timing was difficult to simulate due to local weather events and considerations other than   For the sites where COS fluxes were only measured in one part of the growing season we assume that the timing of seasonal patterns in COS assimilation were well captured since seasonal patterns in NEE, GPP, and LE are properly simulated (Fig. S3-5) and the model scales the CA activity with Vmax, and gs with GPP.

485
We generally found larger underestimations of the ecosystem COS exchange at the higher latitudes (FI-HYY, DK-SOR, AT-NEU, Fig. 1a-c), which is consistent with findings by Ma et al. (2021) who found a missing sink at the higher latitudes that required uptake larger than in the inversion a priori model (i.e. SiB4) in summer (their Fig. 4b). The model-observation biases that we see in the ecosystem COS fluxes are consistent with biases in GPP for some sites. For example, the underestimation of the COS ecosystem flux at DK-SOR, AT-NEU, and FI-HYY is consistent with underestimations of GPP (Fig. S4a-c), which 490 will be further discussed in Sect. 3.1.3.

Effects of varying atmospheric COS mole fractions
Modifying the COS mole fractions to vary spatially and temporally significantly improved the comparison with observations in North America, as seen from the orange (variable COS) and green (fixed COS) line in Fig. 1d-f. Generally, COS mole fractions are lower in the second half of the growing season (Fig. S6), leading to lower COS uptake in that period. When a 495 variable COS mole fraction was used, the MBE value in July-August improved from 9.0 to 2.0 pmol m -2 s -1 at US-HA1, from 7.2 to -0.9 pmol m -2 s -1 for US-IB2, and from 28.6 to 5.4 pmol m -2 s -1 in US-BO1. The influence of the COS mole fraction on the biosphere flux was largest at sites within or close to the Corn Belt in Midwestern US with strong biosphere COS uptake (see also Fig. 5) that therefore has the largest summer-time drop in COS mole fractions (Fig. S6d,e) or the lowest COS mole fraction in general (Fig. S6f). The large COS uptake by maize (corn) is confirmed by the observed COS fluxes reaching ~70 500 pmol m -2 s -1 at midday (Fig. S8). In this region, the lower COS mole fractions lead to reduced COS uptake but would in turn lead to a smaller drop in COS mole fractions. As COS uptake and COS mole fractions are interconnected, SiB4 should ideally be directly coupled to an atmospheric transport model. At other sites (Europe) the variable COS mole fractions did not improve the model-observation bias, but instead caused a 505 slightly larger underestimation by the model. The comparison of COS mole fractions from the TM5-4DVAR inversion against those observed at the measurement sites (Fig. S6), did not indicate that the COS mole fractions were consistently better simulated over North America than over Europe. These results imply that the underestimation of COS fluxes over Europe is not likely caused by an underestimation of the COS mole fractions. On the other hand, the COS mole fractions observed at the measurement sites are not as consistently calibrated as the NOAA measurement network. 510 day-and nighttime, defined as 10 -15 hr and 21 -03 hr local time, respectively. These day-and nighttime definitions exclude transitions between day and night (see diurnal cycles in Fig. S8 and S9). On average across all stations, simulated daytime uptake between April through October was lower than the observations by 1.9 ± 6.5 pmol m -2 s -1 (8 ± 27 %). Even though the average model-observation difference is small, there is substantial variability between sites. The underestimation of daytime 535 COS uptake of 19.3 pmol m -2 s -1 (34 %) at DK-SOR was exceptionally large, consistent with the underestimation of daytime GPP in the same period (13.1 µmol m -2 s -1 , 34 %). The COS measurements at DK-SOR were made in June 2016, a period that was warmer than average at this site. As a result, observed GPP was 25 % higher in June 2016 compared to the 1996-2018 average (Fig. S4b). However, SiB4 simulates only a 7 % higher GPP in June 2016 compared to the 1996-2018 average. At the same time, LE is overestimated (Fig. S5b). A SiB4 run with observed meteorology as driver input showed a similar GPP 540 anomaly as when the MERRA2 driver input was used (Fig. S7a). The fact that SiB4 is not able to capture the GPP anomaly is thus not due to the driver data used. These results point to an underestimation of the RuBisCo and CA enzyme activity, and thus gcos, rather than gs, as LE is not underestimated but instead overestimated. Also at AT-NEU and FI-HYY the underestimation of COS vegetation uptake was consistent with underestimations of simulated GPP against long term timeseries (Fig. S4a,c), with a 9 % underestimation of the COS vegetation flux and 13 % underestimation of GPP at FI-HYY in the 545 months June to August. At US-ARM we saw a switch from underestimation to overestimation of daytime vegetation COS uptake over the months April and May, which may be due to COS emissions from other components than the soil, possibly associated with senescing vegetation (Geng and Mu, 2005;Maseyk et al., 2014), which is currently not represented in SiB4 (nor in other models). Overall, we found large variability in model-observation biases between sites, but no clear distinctions emerge from different PFTs for daytime fluxes. 550 The simulated nighttime uptake was too small by an average of 2.1 ± 3.4 pmol m -2 s -1 (35 ± 57 %). Observed nighttime uptake 570 was on average 25 % of the daytime uptake across sites between May-September, with the largest uptake at AT-NEU (11.0 pmol m -2 s -1 ), ES-LM1 (6.9 pmol m -2 s -1 ), and FI-HYY (5.9 pmol m -2 s -1 ). The small flux values during nighttime make the model-observation comparison sensitive to the different correction and processing procedures that were used for the different datasets. Ecosystem fluxes were only storage-corrected for FI-HYY and US-HA1. Kooijmans et al. (2017) showed for FI-HYY that nighttime storage fluxes were on average ~1 pmol m -2 s -1 in summer. Additionally, some datasets are filtered based 575 on a friction velocity threshold, while others are not. Kooijmans et al. (2017) noted that filtering data based on the friction velocity might bias the data to higher nighttime COS uptake as the uptake can be expected to be limited by the COS gradient at the leaf boundary layer under low turbulence conditions. Given these differences between datasets, and the typically large random noise of COS flux measurements, the average underestimation may not be significant overall. Still, we found a substantial underestimation of the nighttime COS uptake at the C3-GRA sites AT-NEU and ES-LM1, and an overestimation 580 in summer at the DBF sites US-HA1 and DK-SOR. These biases might point to an inaccurate intercept of the Ball-Berry stomatal conductance model (g0; i.e. when GPP is (near-)zero) in SiB4, which is currently set to 10 mmol m -2 s -1 for all PFTs in SiB4, except for C4 grasslands and crop types (both C3 and C4) (40 mmol m -2 s -1 ). Observed g0 values at AT-NEU (10-65 mmol m -2 s -1 , Wohlfahrt, 2004) are mostly higher than the 10 mmol m -2 s -1 used in SiB4 and support the hypothesis that the values per PFT and showed that g0 was typically several times larger than the value of 10 mmol m -2 s -1 currently used in SiB4. We adopted the g0 values of Lombardozzi et al. (2017) in SiB4 to test the effect of a modified g0 setting on the nighttime COS 590 vegetation flux (see Table S2 and Fig. S10). Using these updated g0 values, the simulated nighttime COS uptake for C3-GRA improved at AT-NEU, but had larger biases for other sites and PFTs, especially DBF ( Fig S10). As the g0 values from Lombardozzi et al. (2017) did not consistently improve the nighttime COS uptake we did not adopt these as standard SiB4 settings. SOR (Fig. 3b, d). The observed high soil COS uptake in April at FI-HYY is possibly related to snow melt and thawing of the soil and neither model captures this effect on soil COS exchange.
Soil COS emissions were observed at ES-LM1 and US-ARM. US-ARM was an agricultural site where emissions may build up after the peak growing season in the period associated with senescence and harvest (Maseyk et al., 2014). The Berry soil 625 model did not simulate soil COS emissions (Fig. 3h). In contrast, the increase in COS emissions at the agricultural site US-ARM was simulated by the Ogée soil model, although the increase of the emissions started later than in the observations. The soil emissions of COS were not simulated at the C3-GRA site ES-LM1. However, the soil at ES-LM1 was fertilized (Weiner et al. 2018), as well as that AT-NEU (Spielmann et al. 2020), which make these sites more representative of agricultural soils rather than grassland soils. When ES-LM1 was simulated as an agricultural soil (the same code, but with different uptake and 630 production parameter values, see Table 2), the model showed COS emissions more consistent with observations (green line in Fig. 3g). Also, the simulated fluxes at the AT-NEU site became smaller and in better agreement with observations when the site was considered as an agricultural soil. The accuracy of simulations of soil COS emissions depends on the accuracy of the production parameter a. The standard deviation of the production parameter a (7.3) is relatively large for agricultural soils compared to other soil types (Table 2) and is an indication of the uncertainty of using a single production value in the SiB4 model. Reasons for this uncertainty can 645 be the local variability in soil characteristics like nitrogen content, which has been shown to correlate well with COS production rates (Kaisermann et al. 2018b). Moreover, soil moisture and soil temperature were important parameters in the calculation of the COS soil flux. In general, we found that the variability and absolute values of soil moisture and especially soil temperature were well captured by the SiB4 model. We found a MBE across all sites of 0.01 m 3 m -3 and 0.1 °C (RMSE 0.06 m 3 m -3 and 2.1 °C) for soil moisture and temperature, respectively, calculated over all years available from the EC sites. Also, Smith et al. 655 (2020) showed that SiB4 was capable of reproducing the drop in soil moisture as a result of a regional drought in Europe, albeit with a delay. We did not find consistent patterns in model-observation biases of the soil COS fluxes that were consistent with that of soil moisture or temperature (Fig. S11, S12). Still, the soil moisture observations at US-ARM show a sharper drop in spring than the simulations (Fig. S12h), which could explain why the simulations show a delayed onset of soil COS emissions. Moreover, the exact role of thermal and photo-production of COS remains uncertain, as well as the interaction with 660 soil organic matter and litter, and thereby limits the accuracy of soil COS production simulations (Maseyk et al. 2014;Whelan and Rhew 2015;Meredith et al. 2018;Kaisermann et al. 2018a).
Overall, changing from the Berry soil model to the Ogée soil model had a relatively small effect on monthly average ecosystem fluxes (see SiB4_500_Berry (blue) and SiB4_500_Ogee (green) in Fig. 1), except for agricultural sites, where the Berry soil 665 model lacked COS soil emissions that contribute to fluxes at those sites.

Calibration factor a
The calibration factor a was derived to scale gcos to match SiB4 COS plant assimilation with COS flux observations of laboratory leaf gas exchange measurements (Berry et al. 2013). The aobs values that we derived based on field measurements of COS ecosystem and soil fluxes, together with simulated gcos, gs and gb and VCOS are close to the value 1400 (Fig. 4), which 670 support the initial calibration by Berry et al. (2013) using laboratory leaf gas exchange measurements. At the same time, we found aobs to vary in time and between sites (Fig. 4), indicating that a single a value was not able to capture the variation of measured COS vegetation fluxes across sites and seasons. The average summer-time aobs (June-August) of 1616 ± 562 was 15 % higher than the current value of 1400. This was consistent with our findings that, on average, SiB4 underestimates COS biosphere fluxes with both a variable and constant COS mole fraction (Sect. 3.1.3). We did not find patterns in aobs that apply 675 to all PFTs in the same way that would have helped to update a in SiB4. However, for DBF and C3-GRA sites we observed that the 2-weekly average aobs typically goes down with increasing air temperature for temperatures above ~ 16 °C (Fig. S13).
This observation requires further investigation from hourly data points and will be further discussed in our recommendations

Global biospheric COS sink
The simulated global patterns in COS uptake were similar to that of GPP (Fig. S14), due to the modelled vegetation COS 695 uptake being coupled to GPP through the RuBisCO enzyme activity and stomatal conductance. Globally, the largest portion of COS uptake took place in tropical regions of South America, Africa, and Asia (Fig 5). In the Northern Hemisphere (NH), COS was mainly taken up during the summer months (Fig. 5b). The spatial distribution of COS uptake was also similar to that presented by Maignan et al. (2021) based on ORCHIDEE simulations. Using the original SiB4 model, i.e., the original Berry soil model and fixed 500 pmol mol -1 COS mole fractions, the global COS biosphere sink amounts to 922 ± 11 (mean ± SD) 700 Gg S yr -1 over the years 2000-2020 with no substantial trend. 146 Gg S yr -1 of the total COS biosphere sink was due to soil uptake ( Table 1). The change from the original Berry soil model to the Ogée soil model lowered the soil uptake in most regions globally (Fig. 6a, S15). The tropical soil COS uptake reduced from ~ 4-5 pmol m -2 s -1 to ~ 2-3 pmol m -2 s -1 . In the NH, the soil uptake is also reduced due to the contributions of COS production in agricultural soils. The global COS soil sink thereby reduced by 29% from 146 to 104 Gg S yr -1 when we changed from the original Berry soil model to the Ogée soil model, but 705 this represents only a 5 % reduction of the total COS biosphere sink. On the other hand, moving from a fixed to a spatially and temporally varying COS mole fraction caused an additional reduction of the global COS biosphere sink to 753 Gg S yr -1 (Fig.   6b, S16). This 15 % reduction relative to a simulation with a constant and spatially uniform 500 pmol mol -1 COS mole fraction illustrates the importance of accounting for varying COS mole fractions. The largest decrease in the global COS biosphere sink (169 Gg S yr -1 , i.e. from 922 to 753 Gg S yr -1 ) occurs in the Tropics (113 Gg S yr -1 for latitudes between -23.5 and +23.5 °N) as the high productivity in the Tropics leads to the largest COS uptake 725 and the largest decrease in COS mole fractions. This update is a significant contribution to closing the gap in the COS budget of ~ 432 Gg S yr -1 (Ma et al., 2021); however, it does not fully eliminate the missing source in the COS budget. The biosphere flux resulting from inverse modelling by Ma et al. (2021) indicates COS emissions in the Amazon (Fig. S17). While biosphere emissions over the Amazon are unrealistic (Glatthor et al., 2015), it reflects the large missing source in the Tropics (land and ocean) that we are unable to attribute to the biosphere; see a comparison of our SiB4 biosphere flux with the inverted biosphere 730 flux by Ma et al. (2021) in Fig. S17. A potential reason for unrealistic attribution of missing sources of COS is that there are no NOAA observations in the Tropics and its upwind regions to constrain the TM5 inversions. For these reasons, it is unlikely that the gap in the COS budget is solely caused by an overestimated tropical biosphere sink. Still, flux observations in the Tropics would have to confirm this.

735
In Sect. 3.1.3 we found on average an 8 ± 27 % underestimation of the daytime COS vegetation flux as simulated by SiB4. If we assume that the daytime uptake dominates the total COS uptake, and we correct the COS vegetation sink for the underestimation globally, then we find a vegetation sink of 717 ± 179 Gg S yr -1 instead of 664 Gg S yr -1 , and a total biosphere uptake of 806 ± 179 Gg S yr -1 instead of 753 Gg S yr -1 . Note, however, that this scaling is highly uncertain, because we found substantial variability between sites, and a large fraction of the uptake occurs in the Tropics, for which we cannot validate the 740 SiB4 model due to a lack of observations.

Recommendations for COS-specific future model development 765
We found model-observation biases that could be ascribed to different components of the model (depending on the site), such as the soil COS flux or vegetation COS uptake, where the latter was caused by underestimated enzyme activity that also links to GPP. If sufficient COS flux observations were available, these could help as an extra constraint to improve the model enzyme activity and thereby GPP. Such an approach would require a number of advancements in the understanding and implementation of COS biosphere exchange in SiB4. We have identified a number of ways to improve the COS flux 770 simulations in SiB4, which might also apply to mechanistic COS implementations in other biosphere models: 1. The simulation of COS uptake is strongly coupled to GPP through gs and Vmax (which is included in gcos), and therefore relies on the accuracy of these model parameters. However, several studies have shown that the ratio of COS to CO2 deposition velocities (in literature also called "leaf relative uptake") varies with temperature (Cochavi et al., 2021;Stimler et al., 2010) and humidity (Sun et al. 2018;Kooijmans et al. 2019), in addition to the better known variability 775 with light. The temperature response of the COS uptake is currently taken from Vmax and is scaled with an empirical temperature function (Eq. (2)) and an additional factor Tcan/T0 (Eq. (3)), where the latter increases the COS uptake at higher temperatures. However, the Tcan/T0 term has been added as a simple correction, but has not been empirically derived. We suggest a refined calibration of the internal conductance gcos such that it captures the true temperature variation of COS vegetation fluxes. The temperature dependence of the CA enzyme activity could be 780 determined from laboratory experiments, to be able to keep effects other than temperature (e.g. on mesophyll conductance) constant. Field observations could then be used to scale the laboratory-based calibration to ecosystem level and to different PFTs.

Deleted: 7
Deleted: Based on the analysis of the calibration factor aobs 785 (Sect. 3.2) we suggest a refined calibration of the internal conductance gcos such that it captures the true temperature variation of COS vegetation fluxes. The current calibration factor of 1400 is based on laboratory leaf gas exchange measurements of different plant species. However, our analysis 790 based on field observations of COS vegetation fluxes shows that an alpha of 1400 mostly underestimates the COS uptake, and suggest that the laboratory-based calibration did not fully resemble the COS vegetation uptake in the field. We suggest a re-calibration of gcos based on both laboratory and field observations of the COS 795 vegetation uptake. Extra attention should be given to the temperature response of COS uptake. The modeled COS uptake is strongly coupled to GPP. supports the use of atmospheric COS as a global-and regional-scale photosynthesis tracer. A larger effect on the global COS biosphere sink was found by changing the fixed COS mole fraction of 500 pmol mol -1 to values that vary spatially and temporally. The reduction in the COS sink strength is most pronounced in regions with large biomass such as the Tropics. This analysis highlights the importance of accounting for variations in atmospheric COS mole fractions, which was not yet adopted as standard practice. We make a number of recommendations for future improvements of the model, including re-calibration 855 of the COS model parameters. However, we are limited by site-and leaf level data coverage to be able to accurately constrain the model over different PFTs and seasons. More campaigns and long-term observations in underrepresented PFTs, biomes and soil types and more laboratory measurements such as for CA sensitivity in leaves and soils would be key to continued improvement of the model.

Code and data availability 860
The SiB4 code is available online at https://gitlab.com/kdhaynes/sib4_corral. SiB4 simulation output used in this study is available at https://doi.org/10.5281/zenodo.5084644. COS campaign data are downloaded from the original data publications as reported in Table 3. Long-term CO2 flux timeseries from FLUXNET, AmeriFlux or ICOS are downloaded from the references listed in Table S1.