Carbonyl sulfide: comparing a mechanistic representation of the vegetation uptake in a land surface model and the leaf relative uptake approach

Land surface modellers need measurable proxies to constrain the quantity of carbon dioxide (CO2) assimilated by continental plants through photosynthesis, known as gross primary production (GPP). Carbonyl sulfide (COS), which is taken up by leaves through their stomates and then hydrolysed by photosynthetic enzymes, is a candidate GPP proxy. A former study with the ORCHIDEE land surface model used a fixed ratio of COS uptake to CO2 uptake normalised to respective ambient concentrations for each vegetation type (leaf relative uptake, LRU) to compute vegetation COS fluxes from GPP. The LRU approach is known to have limited accuracy since the LRU ratio changes with variables such as photosynthetically active radiation (PAR): while CO2 uptake slows under low light, COS uptake is not light limited. However, the LRU approach has been popular for COS–GPP proxy studies because of its ease of application and apparent low contribution to uncertainty for regional-scale applications. In this study we refined the COS–GPP relationship and implemented in ORCHIDEE a mechanistic model that describes COS uptake by continental vegetation. We compared the simulated COS fluxes against measured hourly COS fluxes at two sites and studied the model behaviour and links with environmental drivers. We performed simulations at a global scale, and we estimated the global COS uptake by vegetation to be −756 Gg S yr−1, in the middle range of former studies (−490 to −1335 Gg S yr−1). Based on monthly mean fluxes simulated by the mechanistic approach in ORCHIDEE, we derived new LRU values for the different vegetation types, ranging between 0.92 and 1.72, close to recently published averages for observed values of 1.21 for C4 and 1.68 for C3 plants. We transported the COS using the monthly vegetation COS fluxes derived from both the mechanistic and the LRU approaches, and we evaluated the simulated COS concentrations at NOAA sites. Although Published by Copernicus Publications on behalf of the European Geosciences Union. 2918 F. Maignan et al.: Modelling vegetation carbonyl sulfide uptake the mechanistic approach was more appropriate when comparing to high-temporal-resolution COS flux measurements, both approaches gave similar results when transporting with monthly COS fluxes and evaluating COS concentrations at stations. In our study, uncertainties between these two approaches are of secondary importance compared to the uncertainties in the COS global budget, which are currently a limiting factor to the potential of COS concentrations to constrain GPP simulated by land surface models on the global scale.


Introduction
Humanity has to face the urgency of climate change if it hopes to limit adverse future impacts (Allen et al., 2018;IPCC, 2019a, b). In order to make reliable predictions of future climate, scientists have built powerful numerical Earth system models (ESMs), where they continuously integrate gained knowledge on a multitude of climate-related and climate-interacting processes. The carbon cycle is at the heart of the present global warming, caused by anthropogenic CO 2 emissions (Ciais et al., 2013). In the global carbon budget, the land component shows the largest uncertainty (Le Quéré et al., 2018;Bloom et al., 2016). Land surface models (LSMs) struggle to accurately represent the large spatial and temporal variability of the CO 2 gross and net fluxes (Anav et al., 2015). CO 2 is first assimilated through plant photosynthesis, before being respired by the ecosystem. The quantity of assimilated carbon is called gross primary productivity (GPP). All other carbon fluxes and stocks derive from this first gross assimilation flux. To help reduce uncertainties in the estimated GPP, LSMs can benefit from knowledge obtained through local eddy covariance measurements of the net ecosystem-atmosphere CO 2 exchange (Friend et al., 2007;Kuppel et al., 2014).
GPP proxies are also used, such as solar-induced fluorescence (Norton et al., 2019;Bacour et al., 2019), isotopic composition of atmospheric CO 2 (δ 18 O; Farquhar et al., 1993;Welp et al., 2011; δ 13 C: Peters et al., 2018), and carbonyl sulfide (COS) atmospheric concentrations (Hilton et al., 2015). Using atmospheric COS measurements as a tracer for terrestrial photosynthesis was first suggested by Sandoval-Soto et al. (2005) and Montzka et al. (2007), and Campbell et al. (2008) provided quantitative evidence using airborne observations of COS and CO 2 concentrations and an atmospheric transport model. COS is an atmospheric trace gas that has a molecular structure very similar to CO 2 and is likewise taken up by plants through stomates. COS is hydrolysed within the leaf, with this reaction being catalysed by the enzyme carbonic anhydrase (CA). This reaction is lightindependent (Protoschill-Krebs et al., 1996;Goldan et al., 1998) and, because of the high catalytic efficiency of this enzyme (Ogawa et al., 2013;Ogée et al., 2016;Protoschill-Krebs et al., 1996), COS hydrolysis inside the leaf seems therefore to be limited by COS supply driven by changes in stomatal conductance (Goldan et al., 1988;Sandoval-Soto et al., 2005;Seibt et al., 2010;Stimler et al., 2010). Leaves' uptake of COS and CO 2 is thus very similar, but leaves do not produce COS (Protoschill-Krebs et al., 1996;Notni et al., 2007), whereas they emit CO 2 through respiration. That is why vegetation COS fluxes could be used as a proxy for GPP. It is however to be noted that Gimeno et al. (2017) reported COS emissions by bryophytes during daytime.
The approach generally adopted to constrain GPP with COS relies on the determination of a leaf relative uptake (LRU), which is the ratio of COS to CO 2 uptake normalised by their atmospheric concentrations (Sandoval-Soto et al., 2005): where F COS is the flux of COS uptake (pmol COS m −2 s −1 ), GPP is the gross flux of CO 2 assimilation (µmol CO 2 m −2 s −1 ), [COS] a is the atmospheric COS mixing ratio (pmol COS mol −1 , ppt), and [CO 2 ] a is the atmospheric CO 2 mixing ratio (µmol CO 2 mol −1 , ppm). LRU can be estimated experimentally and then used as a scaling factor for estimating GPP, if F COS , [COS] a , and [CO 2 ] a are available. Measurements can be made at the leaf level using branch chambers (Seibt et al., 2010;Kooijmans et al., 2019); LRU can also be estimated at the ecosystem level: eddy-covariance flux towers measure the ecosystem total COS flux (Kohonen et al., 2020), and removing the soil contribution gives access to the vegetation part (Wehr et al., 2017). Soil can absorb and emit COS (Whelan et al., 2016;Kitz et al., 2020), with the magnitude of their flux being generally much lower than that of vegetation fluxes (Berkelhammer et al., 2014;Maseyk et al., 2014;Wehr et al., 2017;Whelan et al., 2018). Epiphytes (lichen, mosses) could also have a significant contribution to the ecosystem COS budget (Kuhn and Kesselmeier, 2000;Rastogi et al., 2018).
However, LRU does not appear constant under some environmental conditions. For example, the fixation of carbon from CO 2 relies on light-dependent reactions, unlike the uptake of COS by the CA enzyme, which is light-independent (Stimler et al., 2011). Because of these different responses of COS and CO 2 uptake in leaves, LRU varies with light conditions and decreases sharply with PAR increase (Stimler et al., 2010;Maseyk et al., 2014;Commane et al., 2015;Wehr et al., 2017;Yang et al., 2018). Consequently, LRU values are smaller at midday or in seasons with high incoming light (Kooijmans et al., 2019). Moreover, COS assimilation continues at night as stomatal conductance to gas transfer does not drop to zero, whereas CO 2 uptake by plants stops, leading to an infinite value of LRU. Note however that stomates mostly close at night, so the COS uptake at night is smaller than the COS uptake during the day. The diel (i.e. 24-hourly) variation in LRU with light may however be only of second-order importance as GPP is very low at low light, and Yang et al. (2018) found that considering sub-daily variations in LRU when computing daily mean GPP values had no importance. It has also been shown that LRU varies between plant species (Stimler et al., 2011), which is why different LRU values were estimated for different vegetation types (Seibt et al., 2010;Whelan et al., 2018). The variability of LRU with plant type and over a day and season (inferred by changes in light conditions) should therefore be carefully accounted for when COS concentrations or flux measurements are used to estimate GPP at the ecosystem and larger scales. We also have to acknowledge that there are still factors that are not accounted for if discrepancies between GPP and COS-based estimations are larger than their estimated respective uncertainties.
Before being able to use COS observations to constrain the simulated GPP, land surface models (LSMs) first need to have an accurate model to simulate vegetation COS fluxes. In a former study, Launois et al. (2015b) simply defined the COS uptake by vegetation as the CO 2 gross uptake simulated by LSMs, scaled with a constant LRU value for each large vegetation class. The goal of this study is to now simulate the uptake of atmospheric COS by continental vegetation in a more complex and realistic way using a mechanistic approach within an LSM and apply this model to evidence the shortcomings or pertinence of the LRU concept, depending on the studied scales.
i. We used the state-of-the art ORCHIDEE LSM (Krinner et al., 2015), and we implemented in it the vegetation COS uptake model of Berry et al. (2013) to simulate the COS fluxes absorbed at the leaf and canopy levels by the continental vegetation.
ii. We evaluated the simulated COS fluxes against measurements at two forest sites, namely the Harvard Forest, United States (Wehr et al., 2017), and Hyytiälä, Finland (Kooijmans et al., 2019;Kohonen et al., 2020;Sun et al., 2018a). We studied the high-frequency behaviour of the modelled conductances over the season and the dependency of the LRU on the environmental and structural conditions.
iii. We compared the simulated mechanistic COS fluxes at the global scale to former estimates; we studied LRU values estimated from monthly fluxes, which are pertinent for atmospheric studies, and compared them to monthly means of high-frequency LRU values.
iv. The mechanistic and LRU simulated COS fluxes were used with the atmospheric transport model LMDz , to provide atmospheric COS concentrations that were evaluated against measurements at sites of the NOAA network.  Kuppel et al., 2012). To compute the maximum photosynthetic capacity at leaf level, V max , the reference value is multiplied at a daily time step by the relative photosynthetic efficiency of leaves based on the mean leaf age following Ishida et al. (1999) (see Eq. A12 and Fig. A12 in Krinner et al., 2005). Leaves are very efficient when they are young and stay so till they approach their pre-defined leaf lifespan. The temperature dependence of the maximum photosynthetic capacity follows Medlyn et al. (2002) and Kattge and Knorr (2007). A water stress function varying between 0 and 1 depending on soil moisture and root profile (de Rosnay and Polcher, 1998) is applied on maximum photosynthetic capacity and conductances. The canopy is discretised in several layers of growing thickness, the number depending on the actual leaf area index (LAI). All the incoming light is considered to be diffuse, and no distinction is made between sun and shaded leaves. The light is attenuated through the canopy following a simple Beer-Lambert absorption law. The CO 2 assimilation, the stomatal conductance, and the intercellular CO 2 concentration C i are computed per LAI layer, provided LAI is higher than 0.01 and the monthly mean air temperature is higher than −4 • C. The CO 2 assimilation and the stomatal conductance are further summed up over all layers to compute GPP and the total conductance at canopy level. The scaling to the grid cell is made using means weighted by the plant functional type fractions. Phenology is fully prognostic with PFT-specific phenological models as described in Botta et al. (2000) and MacBean et al. (2015). ORCHIDEE can be run from the site scale to the global scale, coupled with an atmospheric general circulation model, or in off-line mode forced by meteorological fields. In this study, we prescribed the vegetation distribution for site simulations and used yearly PFT maps derived from the ESA Climate Change Initiative (CCI) land cover products for global simulations . The soil type is derived from the Zobler map (Zobler, 1986). To account for the CO 2 fertilisation effect, we considered global means of [CO 2 ] a with yearly varying values, as provided by the TRENDY model inter-comparison project . The impact of not taking into account the spatial and temporal variations in [CO 2 ] a on GPP has been studied in Lee et al. (2020); while this simplification has indeed no impact at a global yearly scale for GPP, this may be less true at site and seasonal scales. We used the recent ORCHIDEE version fine-tuned for the Climate Model Intercomparison Project (CMIP) 6 exercise (Boucher et al., 2020), forced by micro-meteorology fields at FLUXNET sites or by 2 • CRUNCEP reanalyses at the global scale (https://rda.ucar.edu/datasets/ds314.3/, last access: 19 April 2021).

The Berry model for plant COS uptake
In the ORCHIDEE LSM we implemented the mechanistic model of plant COS uptake based on Berry et al. (2013). In this model, COS follows a diffusive law from the atmosphere to the leaf interior, where it is consumed by CA in the chloroplasts. The uptake from the atmosphere is assumed to be unidirectional, reflecting the fact that COS is generally not produced by plants. The model distinguishes three conductances along the COS path between the atmosphere and the leaf interior: (1) the boundary layer conductance (g B_COS ) to gas transfer between the leaf surface and the atmosphere, (2) the stomatal conductance (g S_COS ), and (3) the internal conductance (g I_COS ). Internal conductance combines the mesophyll conductance and the CA activity into a single equivalent conductance. The stomatal and boundary layer conductances are associated with factors describing diffusion of COS relative to that of water vapour (1.94 and 1.56, respectively; Stimler et al., 2010). In the chloroplast, the COS hydrolysis is catalysed by the enzyme CA, following first-order kinetics. COS uptake depends on the amount of CA and its relative location to intercellular air spaces, which brings in the mesophyll conductance. These two factors have been shown to scale with the maximum reaction rate of the Rubisco enzyme, V max (µmol m −2 s −1 ) (Badger and Price, 1994;Evans et al., 1994). The mesophyll conductance and the first-rate constant are then regrouped into a single equivalent internal conductance, proportional to V max : (2) The parameter α takes two values depending on the plant photosynthetic pathway (C 3 or C 4 ). These values were determined experimentally by Berry et al. (2013), who estimated an α = 0.0012 for C 3 and an α = 0.013 for C 4 species. We thus have the final equation: where F COS is the flux of COS uptake (pmol COS m −2 s −1 ); [COS] a is the background atmospheric COS mixing ratio considered here to be a constant (500 ppt); g T_COS , g B_COS , g S_COS , and g I_COS are respectively the total, boundary layer, stomatal, and internal conductances to COS (mol COS m −2 s −1 ); and g B_W and g S_W are respectively the boundary layer and stomatal conductances to water vapour (mol H 2 O m −2 s −1 ). Note that in this work [COS] a is held constant when computing the COS fluxes, contrary to Berry et al. (2013) and Campbell et al. (2017), where [COS] a is dynamic and taken from the previous time step's PCTM (Parameterized Chemical Transport Model) value. The uncertainty introduced by this simplification is evaluated in the Discussion section. The vegetation COS flux and related conductances are computed for each LAI layer and then summed up to get total values at the canopy level. Unless specified otherwise, fluxes, conductances and LRU are further presented and discussed at the canopy level.

Minimal conductances
As plant CO 2 uptake only occurs under certain conditions such as with sufficient light, temperature, and water, CO 2 assimilation is not calculated in ORCHIDEE when these conditions are not fulfilled. Therefore, the stomatal conductance to CO 2 that is needed to obtain the stomatal conductance to COS is not always computed in ORCHIDEE. However, some studies have shown incomplete stomatal closure at night (Dawson et al., 2007;Lombardozzi et al., 2017;Kooijmans et al., 2019), leading to nighttime COS plant uptake (Berry et al., 2013;Kooijmans et al., 2017). Therefore, we had to define a minimal stomatal conductance to COS under these particular conditions when there is no CO 2 assimilation. The minimal conductance to CO 2 used in ORCHIDEE is based on the residual stomatal conductance if the irradiance approaches zero, represented as the g 0 offset in the stomatal conductance models (see Eq. 15 for C 3 and Eq. 25 for C 4 plants in Yin and Struik, 2009). In the absence of water stress, g 0 takes a constant value for C 3 (0.00625 mol CO 2 m −2 s −1 ) and C 4 (0.01875 mol CO 2 m −2 s −1 ) plants. This constant is multiplied by a water stress function to compute the minimal conductance. This minimal conductance to CO 2 was then applied under conditions when there is no CO 2 assimilation, multiplied by the ratio to convert the conductance to CO 2 into a conductance to COS. We thus model COS assimilation even at night, for all PFTs, and in winter for evergreen species, depending on water stress conditions.

Simulations protocol
All simulations were preceded by a "spin-up" phase to get to an equilibrium state where the considered carbon pools and fluxes are stable with no residual trends in the absence of any disturbances (climate, land use change, CO 2 atmospheric concentrations) (e.g. Wei et al., 2014). A few decades are enough to equilibrate above-ground biomass and GPP. As we transport not only COS, but also CO 2 (see Sect. 2.4 below), we need a longer spin-up where all carbon pools including those in the soil are stable and the net CO 2 fluxes oscillate around zero. Equilibrating the ecosystem photosynthesis with its respiration takes a long time as the slowest soil carbon pool has a residence time on the order of 1000 years. The ORCHIDEE model has a built-in spin-up procedure to accelerate the convergence towards this equilibrium state, using a pseudo-analytical iterative estimation of the targeted carbon pools, based on Lardy et al. (2011). For global simulations, we first performed a 340-year spin-up phase with non-varying pre-industrial atmospheric CO 2 concentration and vegetation map, cycling over the same 10 years of meteorological forcing files, where the final relative variation in the global slowest soil carbon pool was less than 5 %. Starting from this equilibrium state, a transient state simulation was then run applying climate change, land use change, and increasing CO 2 atmospheric concentrations, and COS and GPP fluxes were calculated from 1860 to 2017. We performed site simulations at the Harvard Forest (United States) and Hyytiälä (Finland) FLUXNET sites (see below). For the two sites, we first performed a spin-up simulation cycling over the available years of the FLUXNET forcing files, for around 340 years, using a constant atmospheric CO 2 concentration corresponding to the first year of the FLUXNET forcing file. We then performed the transient simulations over the available FLUXNET years, for each site, with a varying CO 2 atmospheric concentration.

Evaluation of vegetation COS fluxes at two FLUXNET sites
Vegetation COS fluxes can be measured using branch chambers or estimated using the difference between measurements of ecosystem and soil fluxes. Such measurements were available at the Hyytiälä (Finland) and Harvard Forest (United States) FLUXNET sites. The Hyytiälä site (61.85 • N, 24.29 • E) is a boreal evergreen needleleaf forest dominated by Scots pine (Pinus sylvestris). Branch measurements of COS fluxes were made in a Scots pine tree from March to July 2017 using gas-exchange chambers (Kooijmans et al., 2019); fluxes were derived from mole fraction changes when the chambers were closed once every hour. Measurements were made with an Aerodyne quantum cascade laser spectrometer (QCLS) and were calibrated against reference standards (Kooijmans et al., 2016). Fluxes from empty chambers were regularly measured to be able to correct for gas exchange by the chamber and tubing material (Kooijmans et al., 2019). We also used the Hyytiälä COS ecosystem fluxes (Kohonen et al., 2020); eddy covariance fluxes were measured during the years 2013-2017 at 23 m height, approximately 6 m above the canopy height. Flux data were processed, quality-screened, and gap-filled according to recommendations by Kohonen et al. (2020). Soil fluxes were also available for the year 2015 (Sun et al., 2018a). We thus derived the COS vegetation fluxes at the canopy scale for that year from the difference between ecosystem and soil fluxes. Soil fluxes were generally low compared to plant uptake. The Harvard Forest site (42.54 • N, 72.17 • W) is a temperate deciduous broadleaf forest with mainly red oak (Quercus rubra), red maple (Acer rubrum), and hemlock (Tsuga canadensis). Ecosystem COS eddy flux measurements were carried out from a tower from May to October, in 2012 and 2013, using an Aerodyne QCLS and calibrated using gas cylinders. They were further split into vegetation and soil components, using soil chamber CO 2 measurements and a sub-canopy flux-gradient approach (Wehr et al., 2017).
The simulated COS fluxes were evaluated against measurements using the root-mean-square deviation: where N is the number of considered observations, F Obs COS (n) is the nth observed COS flux, and F Mod COS (n) is the nth modelled COS flux, and the relative RMSD which is the RMSD divided by the mean value of observations.

F. Maignan et al.: Modelling vegetation carbonyl sulfide uptake
We also computed the bias, standard deviations, and correlation coefficient.
We used partial correlations to identify the main drivers of the modelled conductances. Given the high non-linearity of the equations linking the conductances to their predictors, we also used random forests (RFs) to simulate ORCHIDEE results, and we applied a permutation technique on these RF models to rank predictors (Breiman, 2001). RFs are well adapted for non-linear problems; they were for example used to rank variables of importance for soil COS fluxes in Spielman et al. (2020).

Global-scale flux estimates and comparisons with the LRU approach
We compared our estimate for plant COS uptake at global scale to former studies, with a focus on the LRU approach. We also applied the LRU approach to derive new estimates of global plant COS uptake for comparison, using a monthly climatology of our modelled GPP fluxes over the 2000-2009 period, a constant atmospheric concentration of 500 ppt for COS and global yearly values for CO 2 (from 368 ppm for the year 2000 to 386 ppm for the year 2009). We considered two sets of constant PFT-dependent LRU values. The first set (LRU_Seibt) was taken from Seibt et al. (2010), based on the observed LRU values displayed in their Table 3 (intermediate  column). The second set (LRU_Whelan) used constant values for C 3 (1.68) and C 4 (1.21) plants where the values are an average over different field and laboratory measurements as assembled by Whelan et al. (2018). Both sets are listed in Table 1. Reciprocally, we derived LRU values using Eq. (1) applied to the monthly climatology of our modelled COS and GPP fluxes over the 2000-2009 period; these will be further called LRU_MonthlyFluxes values. LRU_MonthlyFluxes values were computed for all strictly positive GPP values. For each PFT, we studied the spatio-temporal distribution of LRU_MonthlyFluxes values among grid cells where the PFT was present. We also compared these LRU_MonthlyFluxes values computed from a climatology of monthly fluxes to the climatology of monthly mean LRU values, directly computed from the original half-hourly LRU values and further called Monthly_LRU. Given the non-linearity of the problem, we expect LRU_MonthlyFluxes to be different from Monthly_LRU values. Considering that the objective of the LRU approach was to estimate COS fluxes from GPP using a constant value per PFT, the optimal LRU value for each PFT was obtained by linearly regressing monthly COS fluxes against monthly GPP fluxes multiplied by the ratio of the mean COS to CO 2 concentrations, with no offset. Thus with N the number of grid cell month simulated fluxes where the PFT is present in the monthly climatology. We compared this new set of optimal PFT-dependent LRU values against LRU_Seibt and LRU_Whelan.
We finally used the LRU_Opt values to re-compute the monthly mean COS fluxes from our modelled monthly mean GPP and compared with the mechanistic COS flux calculation. The differences, due to the non-linearity of the COS flux calculation, provide some information on the use of a simplified approach based on mean LRU values.

Simulations of COS concentrations and evaluation at NOAA air sampling sites
The vegetation COS fluxes, as well as all other sources and sinks of the global COS budget, based on their latest estimates, are transported with an atmospheric transport model, so that we are able to simulate 3D COS atmospheric concentrations and compare them to the NOAA surface measurements.

The atmospheric transport model LMDz
In order to simulate COS and CO 2 concentrations in the atmosphere, we used the version of the atmospheric component LMDz of the Institut Pierre Simon Laplace Coupled Model (IPSL-CM) (Dufresne et al., 2013), which has been contributing to the CMIP6 exercise. To reduce the computation time, we used its off-line mode: precomputed air mass fluxes provided by the full version of LMDz are used to transport the different tracers . This version is further called LMDz6 and is described in Remaud et al. (2018) and references therein for the transport of CO 2 . The horizontal winds are nudged towards ECMWF meteorological analyses (ERA-5, https://www.ecmwf.int/en/forecasts/datasets/ archive-datasets/reanalysis-datasets/era5) to realistically account for large-scale advection. The tropospheric OH oxidation of COS is calculated from OH monthly data that are produced from a first simulation done with the INCA tropospheric photochemistry scheme (Folberth et al., 2006;Hauglustaine et al., 2004Hauglustaine et al., , 2014. The photolysis reaction of COS in the stratosphere is not considered: the lifetime of COS in the stratosphere is 64 years (Barkley et al., 2008). The model is set up at a horizontal resolution of 3.8 • × 1.9 • (96 grid cells in longitude and latitude) with 39 hybrid sigmapressure levels reaching an altitude up to about 75 km, corresponding to a vertical resolution of about 200-300 m in the planetary boundary layer. The model time step is 30 min, and the output concentrations are 3-hourly averaged.

Atmospheric simulations: sampling methods and data processing
We ran the LMDz6 version of the atmospheric transport model described above for the years 2000 to 2009. The prescribed COS and CO 2 fluxes used as model inputs are presented in Tables 2 and 3. The GPP estimated by ORCHIDEE (148.1 Gt C yr −1 ) is in the high range among the model estimates (Anav et al., 2015), with a corresponding high respiration (145.7 Gt C yr −1 ) to ensure a realistic net ecosystem exchange (Friedlingstein et al., 2019). However, other high GPP estimates can be found in the literature such as Welp et al. (2011) that suggest a range of 150 to 175 based on δ 18 O data. Likewise, Joiner et al. (2018) have proposed a new GPP product, based on satellite data and calibrated on FLUXNET sites, with an estimate around 140 Gt C yr −1 for 2007. The fluxes are given as a lower boundary condition of the atmospheric transport model (LMDz), which then simulates the transport of COS and CO 2 by the atmospheric flow. The atmospheric COS seasonal variations are likely to be dominated by the seasonal exchange with the terrestrial vegetation, while the mean mole fractions result from all sources and sinks of COS, some of which are still largely unknown (e.g. ocean fluxes, Whelan et al., 2018). In this study, we only focus on the seasonal cycle and do not attempt to simulate the annual mean value; we thus started from a null initial state. The atmospheric transport is almost linear with respect to the fluxes: the linearity is a property of the atmospheric transport, though it is violated in LMDz because of the presence of slope limiters in the advection scheme. Overall, since all the other LMDz components are linear, LMDz transport is generally considered linear with fluxes . Relying on this relationship, we first transported each flux separately, and then we added all the simulated concentrations in the end, for each species. For all COS and CO 2 observations, the model output was sampled at the nearest grid point and vertical level to each station and was extracted at the exact hour when each flask sample had been taken. For each station, the curvefitting procedure developed by the NOAA Climate Monitoring and Diagnostic Laboratory (NOAA/CMDL) (Thoning, 1989) was applied to modelled and observed COS and CO 2 time series to extract a smooth detrended seasonal cycle. We first fitted a function including a second-order polynomial term and four harmonic terms, and then we applied a lowpass filter with either 80 or 667 d as short-term and long-term cut-off values, respectively, to the residuals. The detrended seasonal cycle is defined as the smooth curve (full function plus short-term residuals) minus the trend curve (polynomial plus long-term residuals).  (Seibt et al., 2010;Whelan et al., 2018) * A bug has been discovered in the parameterisation of direct COS emissions in the NEMO PISCES ocean model: the hydrolysis rate was 3 times too low, resulting in an artificial build-up of COS in seawaters. As a correction, we divided the total amount of oceanic COS fluxes within a year by 3, assuming that the bug does not affect the spatial pattern of direct emissions of COS. We used the NOAA/GML measurements of both CO 2 and COS at 10 sites located in both hemispheres, listed in Table 4. The samples have been collected as pair flasks one to five times a month since 2000 and are then analysed in the NOAA/GML's Boulder laboratories with gas chromatography and mass spectrometry detection. The measurements are retained only if the difference between the pair flasks is less than 6.3 ppt for COS. These COS measurements can be downloaded from the ftp site ftp://ftp.cmdl.noaa.gov/hats/ carbonyl_sulfide/ (last access: 19 April 2021). The CO 2 atmospheric measurements come from the NOAA's Glob-alView Plus Observation Package (ObsPack; Cooperative Global Atmospheric Data Integration Project, 2018).

Evaluation metrics
To evaluate and compare the performances of the mechanistic and LRU approaches at different NOAA surface sites, we used the normalised standard deviation (NSD) and the Pearson correlation coefficient (R). NSD is calculated as the ratio between the standard deviation of the simulated concentrations and the observed concentrations at the NOAA surface sites. NSD and R values closer to 1 indicate a better accuracy of the model. COS assimilation is at a minimum at night (between 20:00 and 04:00 local solar time) for observed and simulated fluxes (Fig. 1a). During night, uptake of modelled COS flux is around −8 pmol m −2 s −1 while field observations vary between −5 and 0 pmol m −2 s −1 . In the morning, both simulated and observed uptakes increase. However, while the simulation shows a maximum assimilation of −38 pmol m −2 s −1 at noon, the maximum assimilation for observations is reached at 10:00 with a flux of −49 pmol m −2 s −1 . Observed fluxes thus have a greater daily amplitude than simulated fluxes and are a little ahead of the simulation, but this shift does not seem significant given the large variability of observations, as represented by the 1 standard deviation in

Seasonal cycle
The simulated weekly seasonal vegetation COS uptake roughly follows the same trend as the observed one (r = 0.53, Fig. 1b). COS uptake increases in spring when the vegetation growing season starts and decreases in autumn at the end of the forest activity period. Simulated and observed fluxes also take similar values over the 2 years. There are however differences: in 2013 the start of the season is simulated about 2 weeks too late in May instead of late April, and measured fluxes peak in May-June and August-September, while the modelled fluxes peak in July. We notice that the amplitude of observed COS flux variations is larger than the one of modelled fluxes. Kohonen et al. (2020) have quantified the relative uncertainty of weekly-averaged ecosystem COS fluxes at 40 %, which is coherent with the large standard deviation computed for field data (Fig. 1b). RMSD for the seasonal cycle is 7.0 pmol m −2 s −1 , and the relative RMSD is 41 %. The bias is low (−0.3 pmol m −2 s −1 ), and the standard deviations are similar: 6.6 pmol m −2 s −1 for the observations and 7.7 pmol m −2 s −1 for the simulated fluxes. At the Hyytiälä site in the year 2015 (Fig. B1b), the RMSD for the seasonal cycle is 2.4 pmol m −2 s −1 , and the relative RMSD is 25 %; the bias is low too (0.2 pmol m −2 s −1 ) and the standard deviations are also close: 3.6 pmol m −2 s −1 for the observations and 3.5 pmol m −2 s −1 for the simulated fluxes. The correlation coefficient is 0.78. We selected an arbitrary PAR threshold of 50 µmol m −2 s −1 to split between daytime and nighttime fluxes. We see that the modelled nighttime flux varies across the growing season, with a maximum uptake of −10 pmol m −2 s −1 reached in July and a lower absorption in the enclosing colder months. This seasonal variation can be explained by the seasonal change in LAI and the conductance dependency on T air , which increases in summer.

Nighttime fluxes
The observed nighttime fluxes are of the same magnitude but present an opposite seasonal cycle with lower uptake at the summer peak, albeit variations are within the 1 standard deviation represented in Fig. 1a. The modelled nighttime fluxes account from 22 % of the total COS uptake at the peak of the growing season to 45 % in April at the very beginning. The observed ones exhibit slightly lower values, between 14 % and 37 %. At Hyytiälä, the modelled nighttime ratio is also slightly higher (between 30 % and 34 %) than the observed one (between 20 % and 25 %, Fig. B2). These ratios are in line with other studies: Maseyk et al. (2014) reported a ratio of 29 ± 5 % over a wheat field in Oklahoma, and Sun et al. (2018c) reported one of 23 % for the San Joaquin Freshwater Marsh site in California. The results may vary given the definitions adopted for nighttime and daytime periods.

Modelled conductances
To investigate the importance of each conductance in vegetation COS uptake, we compared the three simulated conductances: leaf boundary layer, stomatal, and internal, studying their variability and their drivers at the diel and seasonal scales. The boundary layer conductance to COS is higher than the two other conductances by a median factor larger than 25 (see Table A1 for more detailed statistics). As a high conductance value is equivalent to a low resistance to COS transfer, we focused only on the stomatal (g S_COS ) and internal (g I_COS ) conductances, which are the two most limiting factors to plant COS uptake. Figure 3 presents the mean diel cycles of the simulated total, stomatal, and internal conductances for each season, computed over 2012 at Harvard Forest and 2017 at Hyytiälä. For practicality, we shifted the month of December before the month of January of the same year to compute the winter mean. The seasonal variations are similar at both sites. The conductances, as well as the amplitude of their diurnal cycle, increase from winter to summer and decline in autumn. Harvard Forest is predominantly a deciduous forest, and winter values of the conductances are zero at this site as there are no leaves in that season. Hyytiälä on the other hand is an evergreen pine forest, such that daytime stomatal conductance in winter does not become zero. The stomatal conductance peaks between 09:00 and 13:00, depending on site and season, while the internal conductance peaks later in the afternoon. The total conductance is in general limited by the internal conductance. The stomatal conductance is limiting roughly between 18:00 and 06:00 from spring to autumn at Harvard and only in June-July-August roughly between 21:00 and 09:00 at Hyytiälä.
These results are consistent with the results obtained at branch level by Kooijmans et al. (2019), who found that the COS flux was limited by the internal conductance in the early season and later during daytime, while the effect of the stomatal conductance was larger at night. For the Harvard Forest site, Wehr et al. (2017) computed the stomatal conductance using both a water flux method and a COS flux method and obtained a close agreement between two different methods; the mesophyll conductance is modelled using an experimental temperature response, and the biochemical conductance, representing CA activity, is modelled using a simple parameter (0.055 mol m −2 s −1 ); both scale with LAI to get canopy estimates. Wehr et al. (2017) found similar maximum values around 0.27 mol m −2 s −1 during daytime, from May to October, for the stomatal conductance and for the biochemical conductance (their Fig. 4); adding the slightly larger mesophyll conductance (peaking around 1.0 mol m −2 s −1 ) to the biochemical conductance would thus also lead to a more limiting role of the internal conductance (peaking around 0.21 mol m −2 s −1 ) during daytime, albeit not as strong as for the modelled one (peaking around 0.13 mol m −2 s −1 ); the simulated stomatal conductance exhibits minimum and maximum values similar to the observation-based ones but peaks more sharply in the morning.
To better understand the conductance behaviour, we studied the relative importance of their drivers. These include environmental variables directly or indirectly involved in their modelling: air surface temperature (T air ), photosynthetically active radiation (PAR), vapour pressure deficit (VPD), and soil moisture (SM), as well as LAI, as leaf-level conductances are summed over LAI layers to provide canopy-level conductances. Partial correlations are computed for all halfhourly values of the variables associated with LRU values between 0 and 8 and are provided in Table A2. We also used half-hourly ORCHIDEE outputs associated with LRU values between 0 and 8 to train random forest models for conductances at the two sites, taking into account the same five predictors. A random predictor was also added to check that the variable importance was correctly estimated. All RF models have an accuracy of at least 96 %. Figures B3 and B4 present the relative ranking of the five predictors for the two conductances and the two sites. The ranking is different between the two methods (partial correlation versus RF), but they agree that at both sites the main driver for the internal conductance is air temperature and the main driver for the stomatal conductance is PAR.
As expected, g I_COS mainly depends on T air . This is explained by the fact that g I_COS is proportional to V max , which represents the Rubisco activity for CO 2 ; V max is assumed to be a measure for the mesophyll diffusion and for the CA activity for COS, which are the components of the internal conductance (Berry et al., 2013). V max depends on T air , considered here to be a proxy of the leaf temperature (Yin and Struik, 2009). This strong link explains why g I_COS is more limiting in winter, as T air is low with thus lower enzyme activities, and, as soon as T air rises in spring, g I_COS becomes less limiting, especially at night. PAR is the most important variable for the stomatal conductance at the two sites. Due to how g S_COS is simulated according to Yin and Struik (2009), there is a linear relationship with the CO 2 assimilation, which depends mainly on PAR.

LRU variability
LRU decreases as a function of PAR, as initially observed by Stimler et al. (2010). Kooijmans et al. (2019) made measurements in two branch chambers installed at the top of the canopy in two Scots pine trees in Hyytiälä. They plotted the response of LRU to light, as quantified by PAR. To compare the ORCHIDEE model behaviour to these field data, we determined an LRU using our modelled COS and GPP fluxes, considering a constant atmospheric concentration of 500 ppt for COS and global yearly values for CO 2 .
LRU increases with low PAR values for both branch chambers and for the model and converges towards a constant value for high PAR values (Fig. 4). This demonstrates that assuming a constant value for LRU, and not considering an increase in LRU under low-light conditions, will result in erroneous estimation of COS fluxes. The increasing LRU can be explained by the light dependence of the photosynthesis reaction contrary to the CA activity that is light-independent. Consequently, CO 2 fluxes tend to zero when PAR decreases while COS is still taken up in the dark, leading in theory to infinite values of LRU. The drop of LRU when PAR increases is however much sharper in the model than in the observations. It is to be noted that here we compare LRU values estimated from measurements at the branch level to modelled LRU estimated at canopy level. We conducted a similar modelling study considering only the top-of-canopy level and the associated COS and GPP fluxes, yielding similar results (not shown). This can be linked to the fact that the version of OR-CHIDEE we use considers all the incoming light to be diffuse To focus on LRU behaviour when PAR decreases, we plotted LRU response to PAR for PAR < 1500 µmol m −2 s −1 . and does not distinguish between sun and shaded leaves. We thus have similar LRU values at all canopy levels.
Following the model developed in Seibt et al. (2010, their Eq. 8), the LRU explicitly depends on only two variables: the g S_COS -to-g I_COS ratio and the ratio of the CO 2 intracellular concentration, C i , to [CO 2 ] a (equally named C a ) ratio. The modelled daily mean values for the C i /C a ratio computed at the two sites vary between 0.68 and 1.00 (Fig. B5). These variations are in agreement with Prentice et al. (2014), who state that the C i /C a ratio is pretty stable with only ±30 % variations. These values are in the upper part of the range reported in Seibt et al. (2010, their Table 2); following their Fig. 3, for a given C i /C a ratio a larger g S_COS -to-g I_COS ratio implies a lower LRU, consistent with our results.
We also performed a predictor ranking for LRU, as was done previously with conductances. The predictors rank similarly for the two sites. As shown in Fig. B6, the main factors explaining the variability of the simulated LRU at a halfhourly time step are PAR, T air , and LAI.

Comparison of plant COS uptake sink estimates
The mechanistic approach simulated in the ORCHIDEE model gives a plant COS uptake of −756 Gg S yr −1 over the 2000-2009 period. COS fluxes are the strongest in South America, Central Africa, and Southeast Asia (Fig. 5), as ex-pected as these regions are also the most productive ones for GPP. The more recent studies (Montzka et al., 2007;Suntharalingam et al., 2008;Berry et al., 2013;Launois et al., 2015b) show a higher global plant sink than the one initially found by Kettle et al. (2002) (Table 5). Kettle et al. (2002) used an LRU-like approach, based on net primary productivity (NPP) and on the normalised difference vegetation index (NDVI) temporal evolution, and already acknowledged their estimate was assumed to be a lower-bound one. Estimates from plant chambers and atmospheric measurements (Sandoval et al., 2005;Montzka et al., 2007;Campbell et al., 2008) confirmed that the COS plant sink should be 2-fold to 5-fold larger than estimated in Kettle et al. (2002). Suntharalingam et al. (2008) also found a low estimate of −490 Gg S yr −1 , using 3D modelling of COS atmospheric concentrations, constrained by surface site observations. We note that our estimate is similar to the −738 Gg S yr −1 found by Berry et al. (2013), which was implemented in the Simple Biosphere (SiB) 3 LSM. The reason for this similarity can be that, on top of using the same mechanistic model for vegetation COS uptake, the leaf photosynthesis and stomatal conductance in both LSMs are derived from the same classical models from Farquhar et al. (1980), Collatz et al. (1992), and Ball et al. (1987). Launois et al. (2015b) adopted an LRU approach, using constant LRU values for large MODIS vegetation classes, adapted from Seibt et al. (2010). Based on these values and a set of global GPP estimates from three LSMs (ORCHIDEE, LPJ, CLM4), the authors derived the corresponding global vegetation COS uptakes reported in Table 5

Dynamics of simulated LRU values
The PFT distributions of the LRU values, both those computed using Eq. (1) applied to the monthly climatology of mechanistic COS and GPP fluxes over the 2000-2009 pe-  riod (LRU_MonthlyFluxes) and the climatological monthly means computed directly from the original half-hourly values (Monthly_LRU), do not support the idea of a constant PFT-dependent LRU value (Fig. 6). The distributions are usually not Gaussian; nor are they all unimodal, as is the case for PFT 12 C 3 agriculture for C 4 PFTs (PFT 11 C 4 grass and PFT 13 C 4 agriculture) exhibit a large spread. The median values are represented by vertical red bars in Fig. 6 and listed in Table 1. The optimal values (LRU_Opt) obtained by linearly regressing monthly COS fluxes against monthly GPP fluxes multiplied by the ratio of the mean COS to CO 2 concentrations (see Fig. C1) are represented by vertical green bars and also listed in Table 1. They are usually higher than the median values, with a mean difference of 12.1 %. Using either monthly means or yearly means of fluxes gives very similar optimal LRU values, the mean difference being only −0.2 %.
The LRU values from monthly fluxes (LRU_MonthlyFluxes) tend to be lower than the monthly means of the LRU computed at a half-hourly time step (Monthly_LRU). This is visible in Fig. 6 where the blue distributions yield larger LRU values and in the bi-dimensional histogram of LRU_MonthlyFluxes against Monthly_LRU (Fig. C2). The bias is −0.2 and the correlation is 0.67. This shows that LRU is scale-dependent. The values to be considered should be coherent with their usage. For example, the optimal values we computed are lower than values estimated from measurements, but they are adapted to make the link with atmospheric COS studies.
LRU_Opt values are much smaller than LRU_Seibt values for all PFTs, roughly by a factor of 2. They are closer to the LRU_Whelan values, being smaller for all C 3 PFTs except the tropical broadleaved evergreen forests and higher for C 4 PFTs (Table 1). In the LRU_Opt set, the most productive PFTs (tropical forests and C 4 crops) have the highest values around 1.7, while the less productive PFTs (boreal forests and grasses) have the lowest values around 0.9. To the contrary, in the LRU_Seibt set, temperate broadleaved forests have the highest values (3.6) while needleleaf forests have the smallest value around 1.9.
Another way to understand the distribution of LRU values is to look directly at the scatter plots of monthly COS fluxes against GPP fluxes, multiplied by the ratio of COS to CO 2 concentrations (Fig. C1). For most PFTs, it is in fact obvious that the relationship shows non-linear features, disagreeing with the classical linear LRU model. Based on these findings, we fitted a simple exponential model as with two parameters a and b. However, given the large spread of the data around the model, the Akaike criterion is always favourable to the LRU linear model, so we will not investigate further with this exponential model. More specific research is needed here in order to bridge this data gap. Still, it is important to note that the larger COS fluxes will in gen- Figure 6. Distributions of the LRU values computed from the mechanistic approach over the 2000-2009 period. Each subplot represents one of the 14 vegetated PFTs used in ORCHIDEE, considering all grid cells where the PFT is present. The x axis represents the LRU value between 0 and 3, with 0.1 bins. The y axis represents the occurrences. For each PFT, the black distribution is computed using a monthly climatology of simulated COS and GPP fluxes (LRU_MonthlyFluxes), the blue distribution is computed using the monthly climatology of LRU values estimated at the original half-hourly time step (Monthly_LRU), the red vertical bar represents the median LRU value for LRU_MonthlyFluxes, and the green vertical bar represents the LRU optimal value that minimises the error between plant COS uptakes estimated at a monthly time step by the mechanistic approach and the LRU approach, for all pixels of the considered PFT (see names and abbreviations in Table 1). eral be underestimated using a linear LRU approach. It also appears that in certain PFTs (4, 5, 7) small COS fluxes will be underestimated.
We computed mean annual vegetation COS fluxes using our modelled GPP and this new LRU_Opt set of values and compared them to the mechanistic COS fluxes (Fig. 7a).
The maps of differences between the mechanistic and LRU_Opt-based COS fluxes (Fig. 7b), and relative differences (Fig. 7c), provide evidence for the spatial errors introduced by considering a constant LRU value. The differ-ences are always lower than 4 pmol m −2 s −1 in absolute values and are mainly positive, with the main exception over the Amazon region where the mechanistic approach shows a larger uptake than the linear LRU approach. The difference between the global estimates of the two approaches is less than 2 %; we could still improve the linear regression determining the LRU optimal value by weighting grid-cell fluxes with the corresponding surface of the PFT.
We also compared the mean seasonal cycles of the COS vegetation flux over the 2000-2009 period, for the mecha- nistic approach and the LRU_Opt-based approach, for each PFT (Fig. C3). The seasonal cycles are very similar; for PFT 13 C 4 agriculture, the LRU_Opt-based cycle is slightly in advance compared to the mechanistic cycle.

Simulating atmospheric COS concentration at surface stations
We transported the global COS and CO 2 fluxes (i.e. the ones obtained from the ORCHIDEE model plus the additional components of each cycle, listed in Tables 2 and 3) with the LMDz6 atmospheric transport model as described in Sect. 2.4.2. We analysed COS concentrations derived from simulated COS fluxes obtained with the mechanistic and LRU approaches with regards to observed COS concentrations from the NOAA at a few selected sites. Figure 8 shows the detrended temporal evolution of CO 2 and COS concentrations for the mechanistic and LRU approaches at Utqiaġvik (UTK, Alaska) and Mauna Loa (MLO, Hawaii). The MLO site samples air masses coming from all over the Northern Hemisphere (Conway et al., 1994). CO 2 seasonal amplitude at UTK reflects the contributions of surface fluxes from high-latitude ecosystems (Peylin et al., 1999), but also from regions further south due to atmospheric transport (Parazoo et al., 2011;Graven et al., 2013). These two stations have been used to detect large-scale changes in ecosystem functioning (Graven et al., 2013;. In spite of their importance, LMDz driven by the ORCHIDEE vegetation fluxes has difficulties in representing their seasonal cycles. For instance, at MLO, the simulated seasonal amplitude of CO 2 is overestimated and precedes the observations by 1 month. For COS, the simulated concentrations match relatively well with the observed seasonal variations and seem to be more in phase with the observations than for CO 2 . Such a feature could indicate that the phase issues with CO 2 are not primarily driven by GPP issues but by the other CO 2 flux components. The mechanistic model and its LRU optimal equivalent better reproduce the observed 1-month lag between the COS and the CO 2 simulation at MLO (i.e. the minimum COS lags the one of CO 2 ) than the other LRU approaches with values from Whelan et al. (2018) and Seibt et al. (2010). The simulations differ more in the amplitude than in the phase of their seasonal cycles. The mechanistic approach simulates an amplitude lower than the LRU ones. At MLO for example, the lower amplitude of the mechanistic model is in better agreement with the observations. At UTK, its seasonal amplitude is also lower but is now underestimated. The COS concentration at this station from the mechanistic approach varies between +30 and −50 ppt while it varies between +50 ppt (+37) and −71 ppt (−50) for the simulation based on Seibt et al. (2010)  . This is a direct consequence of lower COS fluxes with the mechanistic model compared to the fluxes based on the Seibt and Whelan LRU approaches. At both the MLO and UTK sites, the difference between the mechanistic model and its LRU optimal equivalent after being transported is lower than 8 ppt, within the range of the observation uncertainty. Table 6 presents the NSDs and Pearson correlation coefficients between simulated and observed COS concentrations for the mechanistic and LRU approaches. We see that the simulation with Seibt et al. (2010) intermediate LRU values overestimates the seasonal standard deviation and has the lowest accuracy for most stations. It is difficult to tell whether the mechanistic model is better than the LRU approach based on Whelan values. While the mechanistic approach captures known features of the temporal dynamics of the COS-to-CO 2 flux ratio, it underestimates the simulated concentrations at Alert (ALT, Canada) and Utqiaġvik (UTK, United States). It should be noted that, due to other sources of errors (in particular transport and oceanic emissions), the comparison presented here should be taken as a sensitivity study of the COS seasonal cycle to the vegetation scheme rather than a complete validation of one approach.  The mechanistic model links vegetation COS uptake and GPP fluxes through the stomatal conductance model, which includes the minimal conductance as an offset, and the common use of the carboxylation rate of Rubisco, V max , in the internal conductance formulation for COS, and in the Rubiscolimited rate of assimilation for CO 2 . The downside is the introduction of the somewhat uncertain α parameter that relates the COS internal conductance to V max . Using COS flux measurements to optimise the parameters of the stomatal and internal conductances would thus in principle benefit the simulated GPP. This optimisation may be done based on appropriate data assimilation techniques; for example, Kuppel et al. (2012) optimised key parameters of the ORCHIDEE model related to several processes including photosynthesis (see their Table 2), by assimilating eddy-covariance flux data over multiple sites. The approach relies on a Bayesian framework where a cost function including uncertainties on observations, model, and parameters is minimised (Tarantola, 1987). The results obtained in this study pave the way for a similar approach using COS fluxes to optimise key parameters controlling GPP; they can be used to define an optimal set-up for the a priori errors and the error correlations in a Bayesian framework. We acknowledge however the scarcity of available measurements for the time being, with no samples for most biomes, a few sites with less than 1 year of data, and only Hyytiälä allowing for interannual variability studies.

First step: improving the mechanistic modelling of vegetation COS fluxes
Without any calibration, the chosen mechanistic model was able to reproduce observed vegetation COS fluxes at the Harvard Forest and Hyytiälä sites with relative RMSDs on the order of 40 %. Regarding conductances, differences are also seen between the diel cycles of simulated and observationbased conductances from Wehr et al. (2017). Diel variations in atmospheric [COS] a , not accounted for in our model, cannot explain these differences, as they would only affect F COS but not the conductances. These discrepancies advocate for the assimilation of COS fluxes to optimise the parameters related to the internal and stomatal conductances. In our modelling framework, the internal conductance is assumed to be the product of V max by the α parameter. This parameter has been calibrated by Berry et al. (2013) using gas exchange measurements of COS and CO 2 uptake (Stimler et al., 2010(Stimler et al., , 2012. As this α parameter seems much more uncertain compared to the relatively well-known V max , we should first try to optimise α keeping V max fixed.

Exploiting the alternative dominant role between stomatal and internal conductances
Without being perfect, the mechanistic model could reproduce some expected behaviours, such as the limiting role of the internal conductance in winter and then during daytime in the growing season, in relation to the control of CA activity and mesophyll diffusion by temperature, as also depicted in Kooijmans et al. (2019). Determining the limiting conductances to COS uptake depending on the time of day provides useful information, as it can be used to better target which model parameters to optimise, using data assimilation approaches. Thus, observations made in the morning and early afternoon could be used to better constrain the α parameter when the internal conductance limits COS fluxes, at least as modelled on the C 3 species of the two sites, and we could investigate whether the α parameter should be further quantified per PFT rather than simply per photosynthetic pathway. It is to be noted that for C 4 species, the internal conductance is larger than for C 3 species by a factor of 10, so that stomatal conductance is limiting, and it could be difficult and useless to try optimising internal conductance using the α parameter. We have to acknowledge the large uncertainty regarding the modelling of the internal conductance.
In parallel to optimising the parameters of the internal conductance, an improvement could thus also be to replace it by the two factors it represents, i.e. the mesophyll conductance and CA activity. A model for the mesophyll conductance is already implemented in ORCHIDEE, with a simple parameter depending on temperature through a multiplication by a modified Arrhenius function following Medlyn and al. (2002) and Yin and Struik (2009). The impact of mesophyll conductance on photosynthesis and water use efficiency is now more studied (e.g. Buckley and Warren, 2014), even if its modelling remains challenging too: the temperature response has notably been reported as highly variable between plant species (von Caemmerer and Evans, 2015), which would imply having PFT-dependent parameters. Regarding measurements, 13 C discrimination of the isotopic composition of CO 2 exchanges allows for an estimation of the mesophyll conductance (Stangl et al., 2019). Concerning CA activity, we could test the simple model using a constant value presented in Wehr et al. (2017). Measuring CA activity can be done at a coarse frequency, using different techniques (Henry, 1991).

Exploiting nighttime conductances
Recent studies have shown that nighttime field measurements of stomatal conductances often exhibit larger values than the ones used in models (Caird et al., 2007;Phillips et al., 2010). In the ORCHIDEE model, minimum stomatal conductances to CO 2 , g 0 , take two different values: 6.25 mmol m −2 s −1 for C 3 species and 18.75 mmol m −2 s −1 for C 4 species. However, Lombardozzi et al. (2017), using data from literature, found that observed nighttime conductances to CO 2 range from 0 to 450 mmol m −2 s −1 with an overall mean value of 78 mmol m −2 s −1 . Moreover, they defined a mean value for each PFT (see Table A3) while the ORCHIDEE model uses one value for all C 3 species and another one for all C 4 species. Using higher nighttime stomatal conductances in models has the impact of increasing plant transpiration and reducing available soil moisture, which alters water and carbon budgets, especially in semi-arid regions (Lombardozzi et al., 2017). Lower VPD values at night, which could limit the impact of higher nighttime stomatal conductances, follow an increasing trend however (Sadok and Jagadish, 2020). A better representation of these minimal conductances in the model could then improve the constraint of gas exchange between the atmosphere and the terrestrial biosphere. It is to be noted that Barnard and Bauerle (2013) found, based on sensitivity analyses, that g 0 was the parameter having the largest influence on their modelled transpiration estimates. They also stress that g 0 should maybe be seen as an asymptotic minimal value, rather than an offset. During nighttime, the stomatal conductance limits COS uptake. In the model, the nocturnal stomatal conductance to COS is calculated from the above-mentioned minimum stomatal conductance values. For now, the absolute vegetation COS fluxes at night are slightly overestimated compared to observed fluxes (updated Fig. 1a for Harvard and Fig. B1a for Hyytiälä), thus hinting to overestimated nighttime stomatal conductances. Therefore, nighttime observations of COS fluxes could be used to optimise the minimum stomatal conductance values for each PFT.
We thus see that COS fluxes could be used, through standard data assimilation techniques, to optimise the model parameters related to conductances, thus contributing to the improvement of the GPP. However, many more COS flux measurements are needed over a large variety of biomes, first to assert the validity of the mechanistic COS model at global scale and second to be assimilated in order to improve simulated conductances and GPP estimates.

The mechanistic versus LRU approach
The mechanistic model is able to reproduce the hightemporal-frequency LRU variations observed at sites. It is thus legitimate to consider this approach to be more accurate than the classical linear LRU approach that uses a timeconstant LRU value per PFT to estimate COS fluxes from GPP. Furthermore we have shown that computing LRU values using Eq. (1) applied to monthly mean fluxes yields values lower than computing monthly means of high-frequency LRU values (Fig. 6). This may explain why the LRU values we have thus estimated from monthly mean fluxes show generally lower values than the ones derived from measurements, although these cover a large range from 0.7 to 6.2 (Seibt et al., 2010;Whelan et al., 2018). More recently, Spielman et al. (2019) estimated LRU values from ecosys-tem and soil measurements: 0.89 for an agricultural soybean field, 1.02 for a temperate C 3 grassland, 2.22 for a temperate beech forest, and 2.27 for a Mediterranean savanna ecosystem; our corresponding PFTs respectively give 1.37 (C 3 crops), 1.18 (Temp C 3 grass), 1.31 (TempBroSum), and 1.06 (TempBroEver), with thus higher estimates for herbaceous plants and lower ones for trees. It is difficult to say whether in situ and laboratory measurements are too sparse and not representative enough of the variability of plants and environmental conditions across the globe to have a reasonable confidence in their derived mean or median LRU values, or whether we can use these LRU values to falsify the modelled COS and/or GPP fluxes. We may also add that LRU values derived from measurements performed in leaf chamber measurements, which are well ventilated and thus associated with large leaf boundary layer conductances, may not be representative of the real-world transfer processes, where the boundary layer conductances vary with wind speed, temporally and within canopy depth (Wohlfahrt et al., 2012).
Without any calibration, the mechanistic approach performs similarly to LRU approaches based on monthly mean fluxes, when COS is transported using all known COS fluxes as inputs, and COS concentrations are evaluated at stations of the NOAA network. We now have a much finer representation of the COS fluxes as, at every time step, the model integrates the plant's response to environmental conditions in the calculation of the internal and stomatal conductances, unlike in the LRU approach which uses constant values for each PFT.
In order to quantify the first-order uncertainty on F COS related to the fact that we have used a constant [COS] a in our implementation of the Berry model, we computed an alternative F COS' , using the LRU approach based on a climatology of hemispheric monthly means of COS atmospheric concentrations (Montzka et al., 2007), the optimal LRU we derived in this study (given in Table 1), average yearly values for CO 2 atmospheric concentrations, and a climatological seasonal cycle of simulated monthly GPP per PFT. Over the 2000-2009 period, the mean difference between the mean seasonal COS fluxes computed with this method (F COS' ) and the ones simulated with the mechanistic model (F COS ) amounts to −7.9 % over the Northern Hemisphere. As expected, the seasonal amplitude of COS fluxes is dampened as [COS] a decreases with vegetation growth. We thus have to improve our methodology to consider a varying [COS] a as was done in Berry et al. (2013), either inside the ORCHIDEE model or as post-processing. This requires devising some trade-off between the high-frequency time step of ORCHIDEE and the cost of running the transport model. However, it is to be noted that there is no impact on the derived LRU values as the LRU does not depend on the considered [COS] a , as long as the same one is considered for the computation of the COS fluxes in the mechanistic model (Eq. 3) and for the computation of the LRU (Eq. 1) (i.e. whether fixed or varying monthly).
However, there is currently a larger uncertainty on other COS fluxes in the global COS budget, which have an important impact on simulated COS concentrations (Ma et al., 2020) and their relative seasonal changes. For example, if we use another estimation of the direct oceanic fluxes (Lennartz et al., 2017), which shows a seasonal cycle whose amplitude is comparable to the one from the vegetation in high latitudes, this results in an overestimated seasonal cycle at all sites, with the mechanistic approach having the most realistic seasonal amplitude (see Appendix D1 and Fig. D1). An additional sensitivity test was performed to assess the impact of indirect oceanic emissions via DMS oxidation on simulated seasonal cycles as the importance of these fluxes in the global COS budget is still debated . Whereas the impact on northern sites is negligible, the removal of indirect oceanic emissions via the DMS of Kettle et al. (2002) decreases the seasonal amplitude of southern sites (CGO and SPO) in the same proportion in all experiments (see Appendix D2 and Table D2). Transport errors also add uncertainties on the simulated concentrations, especially at elevated continental sites (Remaud et al., 2018). Plus, given the present discrepancies between the GPP estimates of different land surface models, it can be argued that using a mechanistic model instead of an LRU approach when comparing COS concentrations seems to be of a second-order importance Hilton et al., 2017). We nevertheless note in this study that we found an uncertainty on the global vegetation COS uptake of 40 % when considering three different LSMs (Launois et al., 2015b), to be compared to an uncertainty of 70 % when considering three LRU datasets.
Setting aside the uncertainty for the moment, how could we use atmospheric COS concentrations to constrain GPP? A first optimisation was performed with the ORCHIDEE model in Launois et al. (2015b), who optimised a single scaling parameter applied on the vegetation COS fluxes simulated with the LRU approach, thus equivalent to a scaling factor applied on the GPP or the LRU. They assimilated the atmospheric COS concentrations measured at the NOAA air sampling stations, using the LMDz transport model ) and a Bayesian framework as in Kuppel et al. (2012). The optimisation reduced in absolute value the estimated global vegetation COS uptake from −1335 to −708 Gg S yr −1 , more in line with this work's estimate based on a mechanistic modelling of vegetation COS uptake. A mid-term perspective is to go beyond a single scaling parameter and to optimise a set of ORCHIDEE parameters using both atmospheric COS and CO 2 data. Such an approach has been used in several studies with CO 2 data only (e.g. Rayner et al., 2005;Peylin et al., 2016). However, compared to CO 2 , the spatial coverage of COS surface observations is still too sparse to accurately constrain the GPP and therefore OR-CHIDEE parameters (Ma et al., 2020). There is some hope that new satellite retrievals of COS column content, such as with the IASI (Infrared Atmospheric Sounder Interferome-ter) instrument, could have enough accuracy to better constrain the surface fluxes (Serio et al., 2020).

Conclusions and outlooks
We have implemented the mechanistic model of Berry et al. (2013) inside the ORCHIDEE land surface model for COS uptake by the continental vegetation. Modelled COS fluxes were compared at site scale against measurements at the Harvard temperate deciduous broadleaf forest (USA) and at the Hyytiälä Scots pine forest (Finland), yielding relative RMSDs of around 40 % at both diel and seasonal scales. We found that the mechanistic model yields a lower and thus more limiting internal conductance compared to former works (Seibt et al., 2010;Wehr et al., 2017). The next step is to perform a sensitivity analysis (Morris, 1991;Sobol, 2001) and to optimise the most sensitive parameters related to the modelled fluxes and conductances, to get a better agreement with observations.
Our global estimate of COS uptake by continental vegetation of −756 Gg S yr −1 is in the lower range of former studies. An important finding is that the LRU computed from monthly values of the COS and GPP fluxes yields values lower than monthly means of high-frequency LRU values. This has consequences for atmospheric studies where COS concentrations integrate influences from fluxes at large spatial and temporal scales.
Using appropriate LRU values, we transported the monthly mean COS fluxes from the mechanistic and LRU approaches using the LMDz6 model. The evaluation of the modelled COS atmospheric concentrations against observations at stations of the NOAA network yields comparable results for both approaches.
As a general conclusion and for the moment, we can say that the mechanistic model is particularly valuable when studying small timescales or spatial scales using COS fluxes, while for global analyses using COS concentrations, both the mechanistic and LRU approaches give similar results. The fact that the global COS budget has so many components with a large uncertainty  limits the use of COS concentrations as a constraint for GPP in land surface models on the global scale, for the present time.
A further development will be to refine the estimation for COS soil fluxes and to implement a mechanistic model for soil COS fluxes inside ORCHIDEE (Ogée et al., 2016;Sun et al., 2015). Having both the vegetation and soil contributions, we will also be able to assimilate ecosystem COS fluxes to optimise COS-related parameters such as α in the internal conductance formulation from the Berry et al. (2013) model for vegetation uptake, and those related to the stomatal conductance (Wehr et al., 2017;Berkelhammer et al., 2020). We will also later look at the complementary constraints on GPP brought by COS and solar-induced fluorescence, another GPP proxy (Bacour et al., 2019;Whelan et al., 2020).
Appendix A: Additional tables related to conductances  Figure B1. (a) Mean diel cycle of observed vegetation COS flux derived from ecosystem COS flux (Kohonen et al., 2020) and soil COS flux (Sun et al., 2018a), and modelled COS vegetation flux in July-September 2015, at Hyytälä, using an atmospheric convention where an uptake of COS by the ecosystem is negative. The shaded areas above and below each curve represent 1 standard deviation of the considered half-hourly values over the July-September period. (b) Mean seasonal cycle of simulated and observed weekly average vegetation COS flux in 2015, at Hyytälä. The shaded areas above and below each curve represent 1 standard deviation of the daily means within the considered week.        We further tested the impact of the indirect COS fluxes through DMS on the simulated concentrations at NOAA sites. To do that, we compared the atmospheric concentrations given with and without prescribing indirect oceanic fluxes through DMS using the Launois et al. (2015a) oceanic fluxes. In our case, the removal of the DMS oceanic emissions decreases the seasonal amplitude at SPO and CGO but has very few impacts at other sites. We also performed the same experiment using the Lennartz et al. (2017) fluxes and reported no impact of DMS indirect fluxes on simulated concentrations at NOAA sites. Competing interests. The authors declare that they have no conflict of interest.