Global modelling of soil carbonyl sulfide exchanges

. Carbonyl sulﬁde (COS) is an atmospheric trace gas of interest for C cycle research because COS uptake by continental vegetation is strongly related to terrestrial gross primary productivity (GPP), the largest and most uncertain ﬂux in atmospheric CO 2 budgets. However, to use atmospheric COS as an additional tracer of GPP, an accurate quantiﬁcation of COS exchange by soils is also needed. At present, the atmospheric COS budget is unbalanced globally, with total COS ﬂux estimates from oxic and anoxic soils that vary between − 409 and − 89 GgS yr − 1 . This uncertainty hampers the use of atmospheric COS concentrations to constrain GPP estimates through atmospheric transport inversions. In this study we implemented a mechanistic soil COS model in the ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems) land surface model to simulate COS ﬂuxes in oxic and anoxic soils. Evaluation of the model against ﬂux measurements at seven sites yields a mean root mean square deviation of 1.6 pmol m − 2 s − 1 , instead of 2 pmol m − 2 s − 1 when using a previous empirical approach that links soil COS uptake to soil heterotrophic respiration. However, soil COS model evaluation is still limited by the scarcity of observation sites and long-term measurement periods, with all sites located in a latitudinal band between 39 and 62 ◦ N and no observations during winter-time in this study. The new model predicts that, globally and over the 2009–2016 period, oxic soils act as a net uptake of − 126 GgS yr − 1 and anoxic soils are a source of + 96 GgS yr − 1 , leading to a global net soil sink of only − 30 GgS yr − 1 , i.e.

global scale but inferred form from field or lab experiments. 95 In this study, our goal is to provide and evaluate new global estimates of net soil COS exchange. To this end: 96 i.
We implemented an empirical-based and a mechanistic-based soil COS model in the ORCHIDEE 97 LSM; 98 ii. We evaluated the soil COS models at seven sites against in situ flux measurements; 99 iii. We estimated soil contributions to the COS budget at the global scale; 100 iv. We transported all COS sources and sinks using an atmospheric model and evaluated the 101 concentrations against measurements of the National Oceanic and Atmospheric Administration 102 (NOAA) air sampling network. 103 This uncatalyzed hydrolysis is quite low compared to the COS hydrolysis catalysed by soil microorganisms, which 240 is the main contribution of COS uptake by soils (Kesselmeier et al., 1999;Sauze et al., 2017;Meredith et al., 241 2018). The enzymatic reaction catalysed by CA follows Michaelis-Menten kinetics. The turnover rate (s -1 ) 242 and the Michaelis-Menten constant (mol m -3 ) of this reaction depend on temperature. The temperature 243 dependence of the ratio is expressed as (Ogée et al., 2016), 244 where , and are thermodynamic parameters, such as = 40 kJ mol -1 , = 200 kJ mol -1 and 246 = 660 J mol -1 K -1 . 247 The total COS consumption rate by soil (s -1 ) is described with respect to the uncatalyzed rate at = 298.15 K 248 and = 4.5 (Ogée et al., 2016) where is the CA enhancement factor, which characterizes the soil microbial community that can consume 251 COS. The CA enhancement factor depends on soil CA concentration, temperature, and pH. Ogée et al. (2016) 252 reported that its values range between 21 600 and 336 000, with a median value at 66 000. We adapted the values 253 of found in (Meredith et al., 2019) to have a CA enhancement factor that depends on ORCHIDEE biomes 254 (Appendix A, Table A1). 255 8 a mis en forme : Justifié Oxic soil COS production 256 Abiotic oxic soil COS production has been observed at high soil temperature (Maseyk et al., 2014;Whelan and 257 Rhew, 2015;Kitz et al., 2017Kitz et al., , 2020Spielmann et al., 2019Spielmann et al., , 2020. However, photodegradation has also been 258 proposed as an abiotic production mechanism in oxic soils (Whelan and Rhew, 2015;Kitz et al., 2017Kitz et al., , 2020. 259 Abiotic COS production is still not well understood but was assumed to originate from biotic precursors (Meredith 260 et al., 2018). 261 In Ogée et al. (2016), the production rate is described as independent of soil but depends on soil temperature 262 and redox potential. This dependence on soil redox potential enables us to consider the transition between oxic 263 and anoxic soils. However, because little information is available on soil redox potential at the global scale, its 264 influence cannot yet be represented in a spatially and temporally dynamic way in a land surface model such as 265 ORCHIDEE. Thus, we decided to use the production rate described in Whelan et al. (2016) that only depends on 266 soil temperature and land use type, 267 where is expressed in pmol g -1 min -1 , is soil temperature (°C) and and are parameters determined by 269 Whelan et al. (2016) for each land use type using the least-squares fitting approach. We adapted the values of 270 and given for four land use types to ORCHIDEE biomes (Appendix A Table A2). Values of and for deserts 271 could not be estimated by Whelan et al. (2016) because COS emission for this biome was not found to increase 272 with temperature. Figure 11 in Whelan et al. (2016) shows that COS emission from a desert soil is always near 273 zero for temperatures ranging from 10°C to 40°C. Moreover, COS emission from a desert soil is also found to be 274 near zero in Fig. 1 of Meredith et al. (2018). This could be explained by a lack of organic precursors to produce 275 COS (Whelan et al., 2016). Therefore, we considered that desert soils, which correspond to a specific non-276 vegetated PFT in ORCHIDEE, do not emit COS. For other ORCHIDEE biomes, COS production was estimated 277 using and for each PFT and the mean soil temperature over the top 9 cm. The unit of was converted from 278 pmol g -1 min -1 to mol m -3 s -1 (in equation 3) using soil bulk density information from the Harmonized World Soil Several studies have shown direct COS emissions by anoxic soils (Devai and DeLaune, 1997;de Mello and Hines, 283 1994;Whelan et al., 2013;Yi et al., 2007). This has been linked to a strong activity of sulfate reduction 284 metabolisms in highly reduced environments such as wetlands (Aneja et al., 1981;Kanda et al., 1992;Whelan et 285 al., 2013;Yi et al., 2007). A previous approach developed by Launois et al. (2015) was based on the representation 286 of seasonal methane emissions by Wania et al. (2010) in the LPJ-WHyME model to represent anoxic soils in 287 ORCHIDEE. The mean values of soil COS emissions from Whelan et al. (2013) were used to attribute to each 288 grid point a value of soil COS emission. In this approach by Launois et al. (2015), salt marshes were not represented 289 despite their strong COS emissions found in Whelan et al. (2013). Emissions from rice paddies were also neglected. 290 Thus, COS emissions from anoxic soils peaked in summer over the high latitudes, following methane production. 291 Because of the scarce knowledge on anoxic soil COS exchange, here we propose another approach to represent 292 the contribution of anoxic soils, which could be compared to the previous approach developed by Launois et al. 293 (2015). To represent the distribution of anoxic soils we selected the regularly flooded wetlands from the map 294 9 a mis en forme : Justifié developed by Tootchi et al. (2019), as represented in Fig. 1. The regularly flooded wetlands cover 9.7% of the 295 global land area, which is among the average values found in the literature ranging from 3% to 21% (Tootchi et 296 al., 2019). Then, in ORCHIDEE each pixel is either considered as anoxic following the wetland map distribution 297 from Tootchi et al. (2019), or as oxic for the rest of the land surfaces. The pixels defined as anoxic soils are 298 considered flooded through the entire year: the seasonal variations of the flooding, as happening during the 299 monsoon seasons, are consequently neglected. 300 On anoxic pixels, we represent anoxic soil COS flux with aThe production rate for anoxic soils is based on the 301 expression developed by Ogée et al. (2016), 302 (mol m -2 s -2 ) the reference production term, a reference soil temperature (K) and 10 the 304 multiplicative factor of the production rate for a 10 °C increase in soil temperature (unitless). As anoxic soil 305 production ranges from 10 to 300 pmol m 2 s -1 for salt marshes and is usually below 10 pmol m -2 s -1 for freshwater 306 wetlands (Whelan et al., 2018), the reference production term was set to 10 pmol m -2 s -1 . 307 All the variables and constants of the empirical and mechanistic models are presented in Appendix A Tables A3 308 and A4. 309

2.1.3
The atmospheric chemistry transport model LMDZ 310 To simulate the COS atmospheric distribution, we use an "offline" version of the Laboratoire de Météorologie 311 Dynamique General Circulation Model (GCM), LMDZ 6 (Hourdin et al., 2020), which has been used as the 312 atmospheric component in the IPSL Coupled Model for CMIP6. The LMDZ GCM has a spatial resolution 313 3.75°long.×1.9°lat. with 39 sigma-pressure layers extending from the surface to about 75 km, corresponding to a 314 vertical resolution of about 200-300 m in the planetary boundary layer, and a first level at 33 m above sea or 315 ground level. The model u and v wind components were nudged towards winds from ERA5 reanalysis with a 316 relaxation time of 2.5 hours to ensure realistic wind advection (Hourdin and Issartel, 2000;Hauglustaine et al., 317 2004). The ECMWF fields are provided every 6 hours and interpolated onto the LMDZ grid. This version has 318 been shown to reasonably represent the transport of passive tracers (Remaud et al., 2018). The off-line model uses 319 pre-computed mass-fluxes provided by this full LMDZ GCM version and only solves the continuity equation for 320 the tracers, which significantly reduces the computation time. In the following, we refer to this offline version as 321 LMDZ. The model time step is 30 minutes, and the output concentrations are 3-hourly averages. 322 The atmospheric COS oxidation is computed from pre-calculated OH monthly concentration fields produced from 323 a simulation of the INCA (Interaction with Chemistry and Aerosols) model (Folberth et al., 2006;Hauglustaine et 324 al., 2004Hauglustaine et 324 al., , 2014 coupled to LMDZ. The atmospheric OH oxidation of COS amounts to 100 GgS yr -1 in the model. 325 Similarly, the COS photolysis rates are also pre-calculated with the INCA model, which uses the Troposphere 326 Ultraviolet and Visible (TUV) radiation model (Madronich et al., 2003) adapted for the stratosphere (Terrenoire 327 et al.,in prep.). The temperature-dependent carbonyl sulfide absorption cross-sections from 186.1 nm to 296.3 nm 328 are taken from (Burkholder et al., 2019). The calculated photolysis rates are averaged over the period 2008-2018 329 and prescribed to LMDZ. Implemented in LMDZ, the COS photolysis in the stratosphere amounts to about 30 330 10 a mis en forme : Justifié GgS yr -1 , which of the same order of magnitude as previous estimates: 21 GgS yr −1 (71% of 30 GgS yr −1 ) by Chin 331 and Davis (1995), between 11 GgS yr −1 and 21 GgS yr −1 by Kettle et al. (2002) and between 16 GgS yr −1 and 40 332 GgS yr −1 by Ma et al. (2021). 333

2.2
Observation data sets 334

2.2.1
Description of the sites 335 The description of the studied sites is given in Table 1.  Table 1). The aboveground vegetation was removed one day 339 before the measurements if needed and the fluxes were derived from concentration measurements using a Quantum 340 Cascade Laser (see Kitz et al., 2020and Spielmann et al., 2020. At AT-NEU, DK-SOR, ES-LMA and IT-341 CRO, a Random Forest model was calibrated against the manual chamber measurements, and then used to simulate 342 half-hourly soil COS fluxes in Spielmann et al. (2019). We compared the ORCHIDEE half-hourly simulated fluxes 343 to half-hourly outputs of the Random Forest model. This enabled to study the diel cycle, and to compute daily 344 observations with no sampling bias for the study of the seasonal cycle. Soil COS fluxes for ET-JA were derived 345 by using the same training method than as the one used in Spielmann et al. (2019). 346 At FI-HYY, soil COS fluxes were measured using two automated soil chambers in 2015. These chambers were 347 connected to a quantum cascade laser spectrometer to calculate soil COS fluxes from concentration measurements 348 (see Sun et al. (2018) for more information on the experimental setup). Any vegetation was removed from the 349 chambers before the measurements. 350 At US-HA, soil COS fluxes in 2012 and 2013 were not directly measured but derived from from flux-profile 351 measurements, connected to CO2 soil chamber measurements and profiles.eddy covariance COS and CO2 352 measurements and soil chamber CO2 measurements conducted in 2012 and 2013. A sub-canopy flux gradient 353 approach was used to partition canopy uptake from soil COS fluxes. For more information on this approach and 354 its limitations, see Wehr et al. (2017). 355 In the study of soil COS fluxes, the difficulty of performing soil COS flux measurements must be acknowledged, 356 as well as the differences between experimental setups and methods to retrieve soil COS fluxes. These limitations 357 are illustrated in the set of observations selected here. Aboveground vegetation had to be removed at some sites to 358 not measure the plant contribution in addition to soil COS fluxes (Sun et al., 2018;Spielmann et al., 2019;Kitz et 359 al., 2020). Vegetation removal prior to the measurements might lead to artefacts in the observations. Some 360 components of the measuring system can also emit COS. In this case, a blank system is needed to apply a post-361 correction to the measured fluxes (Sun et al., 2018;Kitz et al., 2020). Litter was left in place at the measurement 362 sites. 363

COS concentrations at the NOAA/ESRL sites 364
The NOAA surface flask network provides long-term measurements of the COS mole fraction at 14 locations at 365 weekly to monthly frequencies from the year 2000 onwards. We use an extension of the data initially published in 366 Montzka et al. (2007). The data were collected as paired flasks analyzed using gas chromatography and mass 367 11 a mis en forme : Justifié spectrometry. The stations located in the northern Hemisphere sample air masses coming from the entire northern 368 hemisphere domain above 30 degrees. Among them, the sites LEF, NWR, HFM, WIS have a mostly continental 369 footprints (Remaud et al., 20212) while the sites SPO, CGO, PSA sample mainly oceanic air masses of the southern 370 hemisphere (Montzka et al., 2007). The locations of these sites are depicted in Appendix B, Fig. B1. 371

Spin-up phase 373
A "spin-up" phase was performed before each simulation, which enabled all carbon pools to stabilize and the net 374 biome production to oscillate around zero. Reaching the equilibrium state is accelerated in the ORCHIDEE LSM 375 thanks to a pseudo-analytical iterative estimation of the carbon pools, as described in Lardy et al. (2011)

2.3.3
Atmospheric simulations: sampling and data processing 386 We ran the LMDZ6 version of the atmospheric transport model described above for the years 2009 to 2016. We 387 started from a uniform initial condition and we remove the first year as it is considered to be part of the spin-up 388 period. The prescribed COS fluxes used as model inputs are presented in Table 2 2. The fluxes are given as a lower 389 boundary condition, called the surface, of the atmospheric transport model (LMDZ), which then simulates the 390 transport of COS by large-scale advection and sub-grid scale processes such as convection and boundary layer 391 turbulence. In this study, we only evaluate the sensitivity of the latitudinal gradient and seasonal cycle of COS 392 concentrations to the soil COS fluxes. The horizontal gradient aims at validating the latitudinal repartition of the 393 surface fluxes, while the seasonal cycle partly reflects the seasonal exchange with the terrestrial sink, which peaks 394 in spring/summer. This study does not aim at reproducing the mean value as the top-down COS budget is currently 395 unbalanced, with a source component missing (Whelan et al., 2018;Remaud et al., 20222022, and see Table Table  396 53). 397 For each COS observation, the 3D3DD simulated concentration fields were sampled at the nearest grid point to 398 the station and at the closest hour of the measurements. For each station, the curve fitting procedure developed by 399 the NOAA Climate Monitoring and Diagnostic Laboratory (NOAA/CMDL) (Thoning et al., 1989) was applied to 400 modelled and observed COS time series to extract a smooth detrended seasonal cycle. We first fitted a function 401 including a first-order polynomial term for the growth rate and two harmonic terms for seasonal variations. The 402 residuals (raw time series minus the smooth curve) were fitted using a lowpass filter with either 80 or 667 ddaysdd 403 12 a mis en forme : Justifié as short-term and long-term cut-off values. The detrended seasonal cycle is defined as the smooth curve (full 404 function plus short-term residuals) minus the trend curve (polynomial plus long-term residuals). Regarding The closer NSD and r values are to 1, the better the model accuracy is. 420 421

Data assimilation 422
One of the main difficulties with the implementation of a model is to define the parameter values that lead to the 423 most accurate representation of the processes in ORCHIDEE. Calibrating the model parameters is of interest as 424 Ogée et al. (2016) indicate that some of the model parameters such as and the production term parameters have 425 to be constrained by observations. Moreover, the default values for the soil COS model parameters used in this 426 study (Appendix A Tables A1 and A2) are determined by laboratory experiments (Ogée et al., 2016;Whelan et 427 al., 2016), that is why it is interesting to study how the values obtained by calibration against field observations 428 differ from these default values. Data assimilation (DA) aims at producing an optimal estimate by combining 429 observations and model outputs. In this study, we used data assimilationDA to find the model parameter values 430 that improve the fit between simulated and observed soil COS fluxes from the empirical and the mechanistic 431 models. We used the ORCHIDEE Data Assimilation A System (ORCHIDAS), which is based on a Bayesian 432 framework. ORCHIDAS has been described in detail in previous studies (Bastrikov et al., 2018;Kuppel et al., 433 2014;MacBean et al., 2018;Peylin et al., 2016;Raoult et al., 2021), so below we only briefly present the method. 434 Assuming that the observations and model outputs follow a Gaussian distribution, we aim at minimizing the 435 following cost function ( ) by optimizing the model parameters (Tarantola, 2005), 436 selection (Goldberg, 1989;Haupt and Haupt, 2004) and is described for ORCHIDAS in Bastrikov et al. (2018).

Sensitivity analyses 448
We conducted sensitivity analyses at two contrasting sites (FI-HYY and US-HA) to determine which model 449 parameters have the most influence on the simulated soil COS fluxes from the empirical and the mechanistic 450 models. Sensitivity analyses can help to identify the key parameters before aiming at calibrating these parameters. 451 Indeed, focusing on the key model parameters for calibration limits both the computational cost of optimization 452 that increases with the number of parameters and the risk of overfitting. 453 The Morris method (Morris, 1991;Campolongo et al., 2007) was used for the sensitivity analysis as it is relatively 454 time-efficient and enables ranking the parameters by importance. This qualitative method requires only a small 455 number of simulations, (p+1)n, with p the number of parameters and n the number of random trajectories generated 456 (here, n=10). 457 We selected a set of parameters for the Morris sensitivity analyses based on previous sensitivity analyses conducted 458 on soil parameters in ORCHIDEE (Dantec-Nédélec et al., 2017;Raoult et al., 2021;Mahmud et al., 2021). A 459 distinction is made between the soil COS model parameters called first-order parameters ( , and for the 460 mechanistic model and for the empirical model), and parameters called second-order parameters related to 461 soil hydrology, carbon uptake and allocation, phenology, conductance, or photosynthesis (18 parameters, see 462   Tables S3 and S4). The range of variation of the second-order parameters are described in previous studies using 463 ORCHIDEE (Dantec-Nédélec et al., 2017;Raoult et al., 2021;Mahmud et al., 2021). For the first-order 464 parameters, the range of variation is described in Yi et al. (2007) for (±1.08 pmol COS µmol -1 CO2) and in 465 Table 1 in Meredith et al. (2019) for . The ranges of variation for and parameters are not directly given in 466 the literature and were calculated based on information from the production parameters defined in Meredith et al. 467 (2018) (Text S1 and Table S5). The empirical model mainly differs from the mechanistic model with a stronger seasonal amplitude of soil COS 473 fluxes (34% higher), except at the sites where a net COS production is found with the mechanistic model in summer 474 14 a mis en forme : Justifié (ES-LMA and IT-CRO). At all sites, the empirical model shows that the simulated uptake increases in spring 475 reaching a maximum in summer, and decreases in autumn with a minimal uptake during winter. The strong COS 476 uptake in summer from the empirical model can be explained by the proportionality of soil COS uptake to 477 simulated soil respiration, which increases with the high temperatures in summer. In contrast, the mechanistic 478 model depicts almost no seasonality at all the sites where no net COS production is found over the year. As the 479 mechanistic model represents both soil COS uptake and production, the increase in COS production due to higher 480 temperature in summer compensates part of the COS uptake (Appendix C Figure C1). While the uptake from the 481 empirical model is often higher than the one computed with the mechanistic model in summer, soil COS uptake 482 in winter is stronger with the mechanistic representation. 483 The scarcity of field measurements at AT-NEU, ES-LMA, IT-CRO, DK-SOR and ET-JA does not allow an 484 evaluation of the simulated seasonality of COS fluxes. However, at US-HA, the absence of seasonality from May 485 to October in the observations is also found in the mechanistic model, while a maximum net soil COS uptake is 486 reached with the empirical model. 487 We found that the mechanistic model is in better agreement with the observations for 4 (IT-CRO, ET-JA, FI-HYY, 488 US-HA) out of the 7 sites (Table 3), with a mean of 1.58 pmol m -2 s -1 and 2.03 m -2 s -1 for the mechanistic and 489 empirical model, respectively. However, the mechanistic model struggles to reproduce soil COS fluxes at AT-490 NEU and ES-LMA, with an overestimation of soil COS uptake or an underestimation of soil COS production at 491 AT-NEU and a delay in the simulated net COS production at ES-LMA. We might suspect that the removal of 492 vegetation at these sites prior to the measurements could have artificially enhanced COS production in the 493 observations. Indeed, the removal of vegetation could change soil structure and increase the availability of soil 494 organic matter to degradation (Whelan et al., 2016). AT-NEU and ES-LMA are grassland sites for which soils are 495 expected to receive higher light intensity than forest soils. These sites also show a high mean soil temperature of 496 about 20°C during the measurement periods. Therefore, high soil temperature and light intensity on soil surface 497 could enhanced soil COS production as it was related to thermal or photo degradation of soil organic matter (Kitz 498 et al., 2017(Kitz 498 et al., , 2020Whelan et Rhew, 2015;Whelan et al., 2016Whelan et al., , 2018. This is not the case at FI-HYY, ET-JA or 499 DK-SOR, where soil temperature is much lower (mean value about 10°C at FI-HYY and 15°C at ET-JA and DK-500 SOR during the measurement periods) and the forested cover decreases the radiation level reaching the soil. Note 501 that herbaceous biomass is also likely to be higher in grasslands than in forests. Besides, AT-NEU and ES-LMA 502 are managed grassland sites with nitrogen inputs. Then, soil COS production could also be enhanced by a high 503 nitrogen content as suggested by several studies (Kaisermann et al., 2018;Kitz et al., 2020;Spielmann et al., 2020), 504 which is not represented in our models. The mechanistic model is able to represent a net COS production at IT-505 CRO but overestimates it. This might highlight the importance of adapting the production parameters ( and ) 506 in this model to adequately represent a net COS production. In this model, the net soil COS production is related 507 to an increase in soil temperature. However, it is to be noted that IT-CRO is an agricultural site with nitrogen 508 fertilization. Therefore, soil COS production in the observations could also be enhanced by nitrogen inputs. As 509 expected, the empirical model is unable to correctly simulate the direction of the observed positive soil COS 510 exchange rates at ES-LMA and IT-CRO. 511 15 a mis en forme : Justifié

3.1.2
Soil COS flux diel cycles 512 Figure 3 shows the comparison between the simulated and observed mean diel cycles over a month. The 513 observations show a minimum net soil COS uptake or a maximum net soil COS production reached between 11 514 am and 1 pm at AT-NEU, ES-LMA, IT-CRO and DK-SOR. A minimum net soil COS uptake is also observed at 515 US-HA but in the afternoon. At AT-NEU and ES-LMA, neither model is able to represent the observed diel cycle. 516 At these grassland sites, Spielmann et al. (2020) and Kitz et al. (2020) found that the daytime net COS emissions 517 were mainly related to high radiations reaching the soil surface, which impact is not represented in the soil COS 518 models. At IT-CRO and, DK-SOR and US-HA, the diel cycles simulated by the mechanistic model show patterns 519 similar to the observations with a peak in the middle of the day, but with an overestimation of the net soil COS 520 production and a delay in the peak at IT-CRO, and an overestimation of the net soil COS uptake at DK-SOR. The content that is too low to influence soil COS flux. In ORCHIDEE, the simulated range of temperature at US-HA 527 is larger than the one measured on site and temperature is the main driver of the decrease in net soil COS uptake 528 at this site (not shown). Therefore, the enhancement of soil COS production by soil temperature could be only 529 found in the simulated flux, or it could be totally compensated by soil COS uptake in the observations. Therefore, 530 the enhancement of soil COS production by soil temperature could be only found in the simulated flux. Another 531 possibility is that, or it could be totally compensated by soil COS uptake in the observations. The mismatch 532 between the model and the observations could be due to several factors including: i) an insufficient representation 533 of the vegetation complexity by the division in PFTs; ii) a poor calibration of the PFT-specific parameters ( , , 534 ); or iii) missing processes in the model, such as considering the effect of nitrogen content on soil COS fluxes. 535 As the mechanistic model includes PFT-specific parameters ( , , ), we can think that these parameters would 536 need to be calibrated to improve the model performance at the site-scale. The empirical model shows a maximum 537 soil COS uptake around 3 pm at ET-JA, FI-HYY, US-HA and IT-CRO, which is not found in the observations at 538 FI-HYY and is in contradiction with the observed diel variations at IT-CRO and ES-LMA. Considering all sites, 539 the mechanistic model leads to a smaller error between the simulations and the observations, with a mean RMSD 540 of 1.38 pmol m 2 s −1 against 1.87 pmol m 2 s −1 for the empirical model (Table 4) CRO and ES-LMA), soil COS uptake globally generally decreases with increasing soil water content, which 546 appears to be the main driver of soil COS fluxes. This behaviour can be explained by a decrease in COS diffusivity 547 through the soil matrix with increasing soil moisture, reducing soil COS availability for microorganism 548 consumption. Furthermore, an optimum soil water content for net soil COS uptake is found between 10% and 549 15%, which was also observed .%.%. This optimum soil moisture is also represented in Ogée et al. (2016) and was 550 16 a mis en forme : Justifié described in several field studies to be around 12% (Kesselmeier et al., 1999;Liu et al., 2010;van Diest and 551 Kesselmeier, 2008). ThiseTheThe optimum soil water content for soil COS uptake is related to a site-specific 552 temperature optimum, which is found between 13°C and 15°C at US-HA for example. Indeed, Ogée et al. (2016)  553 also describe a temperature Similarly, a temperature optimum was described in Ogée et al. (2016) and in empirical 554 studies with a nanan optimum valuevalue that also depends on the studied site (Kesselmeier et al., 1999;Liu et al., 555 2010;van Diest and Kesselmeier, 2008). At IT-CRO and ES-LMA where a strong net soil COS production is 556 simulated by the mechanistic model, the main driver of soil COS fluxes becomes soil temperature. At these sites, 557 the net soil COS production increases with soil temperature, due to the exponential response of soil COS 558 production term to soil temperature. The increase in soil COS production with soil temperature at IT-CRO and 559 ES-LMA is supported by the observations ( Figure S1). 560 Contrary to the mechanistic model, soil COS uptake computed with the empirical model is mainly driven by soil 561 temperature, with a soil COS uptake that increases with increasing soil temperature. This response of the empirical 562 model to soil temperature is due to its relation to soil respiration, which is enhanced by strong soil temperature. 563 Howeverer, this net increase in soil COS uptake with soil temperature at all sites is not found in the observations 564 ( Figure S1). HoweverIt can be noted that HoweverHowever, low soil moisture values were found to limit soil COS 565 uptake for the empirical model, as seen at ES-LMA for a soil water content below 8%. as it directly scales soil respiration to soil COS fluxes. The following parameters to which soil COS fluxes are the 578 most sensitive are the scalar on the active soil C pool content (soilC) and the temperature-dependency factor for 579 heterotrophic respiration (soil_Q10). Indeed, the soilC parameter determines the soil carbon active pool content, 580 which can be consumed by soil microorganisms during respiration, therefore impacting soil COS fluxes from the 581 empirical model. soil_Q10 impacts soil COS fluxes at both sites as it determines the response of soil heterotrophic 582 respiration to temperature, which is included in the proportionality of soil COS fluxes to the total soil respiration 583 in the empirical model. Similarly, one of the second order parameters, the minimum soil wetness to limit the 584 heterotrophic respiration (min_SWC_resp), has an impact on soil COS fluxes from the empirical model only. The 585 importance of min_SWC_resp for soil COS fluxes is found at US-HA but not at FI-HYY. This can be explained 586 by the difference in soil moisture between the two sites, with an annual mean of 16.2% at US-HA and reaching a 587 minimum of only 8.8%, against an annual mean of 17.5% with a minimum of 12.4% at FI-HYY. 588 17 a mis en forme : Justifié Contrary to the empirical model, soil COS fluxes computed with the mechanistic model are more sensitive to two 589 second-order parameters, the Van Genuchten water retention curve coefficient n (n) and the saturated volumetric 590 water content (θSAT). These two second-order parameters are strongly linked to soil hydrology and determine the 591 soil water content, which affects COS diffusion through the soil matrix and its uptake. The Van Genuchten 592 coefficients occur in the relationships linking hydraulic conductivity and diffusivity to soil water content (van 593 Genuchten, 1980). At both sites, the strong impact of the Van Genuchten water retention curve coefficient n on 594 soil COS fluxes simulated with the mechanistic model highlights the critical importance of soil architecture. Thus, 595 soil COS fluxes computed with the mechanistic model are expected to strongly vary according to the different soil 596 types. Then, the first-order parameters ( , and ) also influence soil COS fluxes from the mechanistic model. 597 However, the uptake parameter ( of PFT 15, boreal C3 grass) has the most influence on soil COS fluxes at FI-598 HYY, while it is the production-related parameter ( of PFT 6, temperate broadleaved summergreen forest) that 599 has the largest impact at US-HA. The stronger influence of the production parameter involved in the temperature 600 response at US-HA might be explained by the difference of temperature between the two sites, which ranges from 601 -10°C to 25°C at US-HA with an annual mean of 7.5°C in 2013, while only ranging from -5°C to 15°C with an 602 annual mean of 4.3°C at FI-HYY in 2015. Similar to the difference in the main driver of soil COS fluxes found in 603 Fig. 4, the most important first-order parameters to which soil COS fluxes are sensitive seem to differ between 604 uptake and production parameters depending on the site conditions. It is to be noted that at US-HA, the most 605 important production parameters are the ones of the dominant PFT at this site (PFT 6), which also correspond to 606 a stronger response of the production term to temperature than for PFT 10 (temperate C3 grass). However, at FI-607 HYY the most influential uptake parameter is for PFT 15 (boreal C3 grass) that only represents 20% of the PFTs 608 at this site while PFT 7 (boreal needleleaf evergreen forest) is the dominant PFT. This can be explained by the 609 range of variation that is assigned to of PFT 7 by Meredith at al. (2019), which is larger than the one of for 610 PFT 15 (9000 against 3100). 611 Finally, a set of parameters related to photosynthesis, conductance, phenology, hydrology, and carbon uptake has 612 an impact on soil COS fluxes computed with both the empirical and the mechanistic models at the two sites. The 613 specific leaf area (SLA), maximum rate of Rubisco activity-limited carboxylation at 25°C (Vcmax25), residual 614 stomatal conductance (g0) and minimum photosynthesis temperature (Tmin) have an impact on soil COS fluxes 615 as they also indirectly affect soil moisture through their influence on transpiration and stomatal opening. The 616 second-order parameters related to soil hydrology (a, Ks, Zroot, θWP, θFC, θR, θTransp_max) impact the soil 617 water availability, which affects soil respiration for the empirical model and soil COS diffusion and uptake in the 618 mechanistic model. For example, the parameter for root profile (Zroot) determines the density and depth of the 619 roots, and therefore how much water can be taken up by roots. At FI-HYY, the difference between prior and posterior soil COS fluxes from the empirical model seems to mainly 632 come from the change in soil_Q10 value (Appendix E, Figure E1). soil_Q10 value drops from 0.83 to 0.53, which 633 corresponds to a prior Q10 value of 2.29 versus a posterior value of 1.70, decreasing the heterotrophic respiration 634 response to soil temperature. Soil COS fluxes computed with the empirical model were found to be strongly 635 sensitive to soil_Q10 ( Figure 5). The posterior value of this parameter has nearly attained the lower bound of its 636 variation range. Since the range of variation represents the realistic values this parameter can take, we need to be 637 careful about the fact that this parameter is trying to take values close to, or potentially beyond, these meaningful 638 values. Furthermore, the optimization deviates the Q10 value at FI-HYY from the ones calculated in the 639 observations over the measurement period (3.0 for soil chamber 1 and 2.5 for soil chamber 2). We could assume 640 that should be defined as temperature-dependent for linking soil COS flux to soil respiration (Berkelhammer 641 et al., 2014;Sun et al., 2018), instead of being considered as a constant. Thus, the optimization of the empirical 642 model could in fact be aliasing the error of onto soil_Q10 because of the impossibility to account for the 643 temperature-dependence of soil COS to CO2 uptake ratio (Sun et al., 2018). At US-HA, the optimization also leads 644 to a decrease of soil_Q10 but to a lesser extent, the parameter remaining comfortably within its range of variation. 645 For the mechanistic model, the optimization reduces the enhancement factor value ( ) for PFT 10 at US-HA and 646 increases the value of the production parameter for the dominant PFT (PFT 6). This enhances the reduction in 647 net soil COS uptake, which was slightly overestimated with the prior model parametrization. At FI-HYY, the 648 optimized parameters show higher values of and of for PFT 15, and of both production parameters ( and 649 ) for the dominant PFT (PFT 7). This increase in both soil COS uptake and production after optimization could 650 correspond to an attempt to better simulate the larger range of variation found in the observations compared to the 651 modelled fluxes. 652 Finally, the optimization also affects hydrology-related parameters for both models. However, while it improves 653 the simulated water content compared to the observations for the mechanistic model at the two sites (RMSD 654 decreases by 28% at FI-HYY and 22% at US-HA), it leads to a degradation at FI-HYY for the empirical model 655 (RMSD increases by more than 3 timesnot shown). Since the empirical model is quite a simplistic model with few 656 parameters, it relies on parameters from different processes to help better fit the observationssometimes 657 degrading the fit to the other processes. The mechanistic model is able to both improve the fit to the COS 658 observations and soil moisture values implying its parameterization is more consistent. 659 This optimization experiment has been promising, highlighting how observations can be used to improve the 660 models. However, since we only optimized over two sites due to the scarcity of soil COS flux observations, for 661 the global scale simulations in the rest of this study, we will rely on the default parameter values of each 662 parameterization. 663 664 19 a mis en forme : Justifié

Soil COS fluxes 666
The spatial distribution of oxic soil COS fluxes shows a net soil COS uptake everywhere except in India, in the 667 Sahel region and some areas in the tropical zone, where net soil COS production is simulated ( Figure Figure 7a). 668 The strongest uptake rates are found in Western North and South America, and in China, with a mean maximum 669 uptake of -4.4 pmol COS m -2 s -1 over 2010-2019. The difference in magnitude between the maximum uptake value 670 and the maximum of production can be noticed, with a net production reaching 67.2 pmol COS m -2 s -1 in the Sahel 671 region. India and the Sahel region, where oxic soil COS production is concentrated, are represented in ORCHIDEE 672 by a high fraction of C3 and C4 crops ( Figure S3S3S3).S34). In the mechanistic model, crops are associated with 673 the lowest value due to overall lower fungal diversity and abundance in agricultural fields (Meredith et al.,674 2019), and the strongest response of oxic soil COS production to temperature as observed by Whelan et al. (2016). 675 Thus, these PFT-specific parameters combined with high temperature in the tropical region can explain the net 676 oxic soil COS production found in these regions. C3 crops are also dominant in China near the Yellow Sea ( Figure  677 S3S3S3).S3S4). However, the mean soil temperature in this region is about 15°C lower than the mean soil 678 temperature in India, leading to a lower enhancement of soil COS production. The highest atmospheric COS 679 concentration is also found in this region with about 800 ppt ( Figure S2S32S2 (Figure Figure 7b). This region is characterized by rice paddies, which were also 685 associated with strong COS production in previous studies (Zhang et al., 2004). 686 The total soil COS fluxes (oxic and anoxic) computed with the mechanistic model ( Figure Figure 7c) show a very 687 different spatial distribution than the one obtained with the empirical model ( Figure Figure 7d). Soil COS fluxes 688 from the empirical model are on the same order of magnitude for net COS uptake than the mechanistic model, 689 with a mean maximum uptake of -6.41 pmol COS m -2 s -1 . However, most soil COS uptakes simulated by the 690 empirical model is located in the tropical region, where soil respiration is strong due to high temperature. The 691 distribution and magnitude of soil COS flux from the empirical approach is similar to the one presented in 692 The difference of soil COS fluxes between the mechanistic model and the empirical model ranges from -4.1 pmol 699 COS m -2 s -1 to +68.0 pmol COS m -2 s -1 (Appendix D, Figure D1). Over western North and South America, northern 700 and southern Africa, western Asia, and eastern, northern and Central Asia, the net COS uptake from the 701 mechanistic model exceeds the uptake from the empirical model. On the contrary, soil COS uptake from the 702 empirical approach is higher than the net COS uptake simulated with the mechanistic model over Eastern North

Interhemispheric gradient 729
We transported total COS fluxes for the different configurations (i.e. including the soil fluxes but also other 730 components of the COS atmospheric budget, listed in Table 2 2) with the LMDZ6 atmospheric transport model as 731 described in Sect. 2.1.3. We analyzed COS concentrations derived from simulated COS fluxes obtained with the 732 mechanistic and two empirical approaches with regards to the COS concentrations observed at 14 NOAA sites 733 depicted in Appendix B, Fig. B1. Note that atmospheric mixing ratios of COS result from the transport of all COS 734 sources and sinks and that, due to other sources of errors (transport and errors in the other COS fluxes), the 735 comparison presented in the following should be taken as a sensitivity study of COS seasonal cycle and inter-736 hemispheric gradient to the soil exchange fluxes rather than a complete validation of one approach or the other. 737 identify vegetation as an underestimated sink in the high latitudes (Ma et al., 2021;Remaud et al., 2022), and the 753 missing source as being the tropical oceanic emissions as being the missing source (Berry et al., 2013;Launois et 754 al., 2015;Le Kuai et al. 2015;Ma et al. 2021;Remaud et al., 20222022;Davidson et al., 2021). The present anoxic 755 soil fluxes have little impact on the surface latitudinal distributions and therefore are unlikely to shed new light on 756 the tropical missing source. An explanation for the small impact is that they are located outside areas experiencing 757 deep convection events (e.g. the Indian monsoon domain) and thus the surface concentrations are less sensitive to 758 these fluxes. 759 Seasonal cycle at NOAA sites 760 Figure 10 shows the detrended temporal evolution of COS concentrations for the mechanistic and empirical 761 approaches at Alert (ALT, Canada) and Harvard Forest (HFM, USA). Because of the mean westerly flow, the 762 HFM site is influenced by continental regions to the west (Sweeney et al., 2015), and is more sensitive to the soil 763 fluxes than other mid-latitude sites located to the west of the ocean (MHD, THD), see Fig. 1 in Remaud et al. 764 (2022Remaud et al. 764 ( (2022. The ALT site samples air masses coming from high-latitude ecosystems (Peylin et al., 1999), but 765 also from regions further south due to atmospheric transport (Parazoo et al., 2011). The reader is referred to 766 Appendix B, Table B2 for the other sites. At both sites, the mechanistic approach tends to weaken the total seasonal 767 amplitude and increase the model-data mismatch. mismatch. At HFMAt ALT, the seasonal amplitude is marginally 768 decreased, while at HFM it is divided by two. At ALT, BRW and SUM, the too high atmospheric concentrations 769 and the too weak seasonal amplitude given by the mechanistic approach are consistent with an underestimated soil 770 absorption at sites ET-JA (Estonia) and FI-HYY (Finland) (see Figure 2). As for Harvard Forest, since the 771 mechanistic soil model shows overall good agreement with the observed soil fluxes (e.g. Figure 2Figure 2 from Kettle et al. (2002) and Launois et al. (2015). These results illustrate the necessity of well constraining both 782 the soil and vegetation fluxes in order to optimize the GPP with the help of atmospheric inverse modelling. 783 4 Discussion 784

Soil budget 785
According to the mechanistic approach of this study, the COS budget for oxic soil is a net sink of -126 GgS yr -1 786 over 2009-2016, which is close to the value of -130 GgS yr -1 found by Kettle et al. (2002) ( Table Table 35). This 787 net COS uptake by oxic soils is higher than the one found in SiB4 by Kooijmans et al. (2021) with -89 GgS yr -1 , 788 also based on the mechanistic model described in Ogée et al. (2016). In SiB4 and in ORCHIDEE, the mechanistic 789 model gives the lowest oxic soil COS net uptake comparedThe mechanistic model gives the lowest oxic soil COS 790 net uptake compared to all previous studies, which were using empirical approaches. This budget is also 41% 791 lower than the one found with the Berry empirical approach in this study, with an uptake of -214 GgS yr -1 . The 792 anoxic soil COS budget computed with the mechanistic approach is +96 GgS yr -1 , which is close to the budget 793 found by Launois et al. (2015) of +101 GgS yr -1 based on methane emissions. However, while COS emissions 794 from anoxic soils were only located in the northern latitudes in Launois et al. (2015), the COS production in this 795 study is also distributed in the tropical region. Thus, we can expect that despite similar budget values for anoxic 796 soils, the difference in flux distribution will impact the latitudinal gradient of COS fluxes. Finally, adding anoxic 797 soil COS budget to oxic soil COS budget results in a total soil COS budget of only -30 GgS yr -1 for the mechanistic 798 approach. 799 When computing the net total COS budget considering all sources and sinks of COS (Table 2), the net total 800 fromempirical we found that neglecting the potential COS production of oxic soils and COS emissions from anoxic 801 soils leads to an overestimation of COS sink or an underestimation of COS source to close the budget (-165 GgS 802 yr -1 ). On the contrary, the total COS budget computed with the mechanistic soil model is closed given the 803 uncertainties on each component (Table 2). However, despite a closed budget, the mismatch between the observed 804 and simulated latitudinal gradients of atmospheric COS concentration highlights errors in COS flux component 805 distributions ( Figure 9). 806 It When computing the net total COS budget considering all sources and sinks of COS, the net total from the 807 empirical approach is closer to zero (-(--35 GgS yr -1 ) than the net total from the mechanistic model (+model 808 (+model (+149 GgS yr -1 ). In the empirical approach, neglecting the potential COS production of oxic soils and 809 COS emissions from anoxic soils leads to aa small overestimation of COS sink or underestimation of COS source 810 to close the budget. On the contrary, the mechanistic approach leads to anann overestimation of COS source or an 811 underestimation of COS sink components. This positive net global budget could be due to an underestimation of 812 vegetation COS uptake in the northern hemisphere, participating in the underestimation of the COS concentration 813 drawdown (Figure 9), but the absence of anthropogenic emission seasonality could also play a role. The two net 814 totals obtained in this study are closer to closing the COS budget than the previous approach from Launois et al. 815 (2015). 816 23 a mis en forme : Justifié Despite a net COS budget closer to zero with the empirical model, it is also to be noted that the mechanistic model 817 better simulates the lack of seasonality in the soil COS flux at US-HA compared to the empirical model ( Figure  818 2). US-HA is represented by 80% of PFT6 (temperate broadleaved summergreen forest) and the absence of 819 seasonality by this PFT has also been reported at a mid-latitude site at Gif-sur-Yvette (Belviso et al., 2020). This 820 PFT is largely found in the temperate region such as in Europe and in the southern United-States. Moreover, NWR, 821 HFM and LEF stations are mainly influenced by COS exchanges from the PFT6. Therefore, the use of the 822 mechanistic model would be recommended to carry out new comparisons at these mid-latitude sites. 823

Variable atmospheric COS concentration versus constant atmospheric COS concentration 824
We studied the impacts of using a constant versus a variable atmospheric COS concentration on soil COS fluxes. 825 At the site-scale we found a distinction between the sites where soil COS production is strong (IT-CRO and ES-826 LMA) and the sites mainly showing a net soil COS uptake. The impact of using a constant atmospheric COS 827 concentration is lower at IT-CRO and ES-LMA because the atmospheric COS concentration does not directly 828 impact the soil COS production term but participates in the net soil COS flux. Our study shows that at the sites 829 where a net soil COS uptake is dominant, using a constant atmospheric COS concentration leads to an 830 underestimation oflower soil COS flux in winter and an overestimation higherof soil COS flux from spring to 831 autumn (not shown). Indeed, during the growing season, plant uptake decreases atmospheric COS concentration 832 ( Figure S1S21S1S1), which reduces COS availability for soil COS diffusion, whereas during winter, a higher 833 atmospheric COS concentration enhances COS diffusion into the soil. shown). We found that using a constant atmospheric COS concentration leads to an overestimation increase of net 844 soil COS uptake over the whole year in the southern latitudes and from June to February in the northern latitudes 845 (reaching 1.62 pmol m -2 s -1 ). This overestimation increase is highers over the regions with the lowest atmospheric 846 COS concentrations, which limits COS diffusion through the soil matrix. On the contrary when atmospheric COS 847 concentration is high in the northern latitudes between April and May, considering a constant atmospheric COS 848 concentration leads to an underestimation ofdecreases the net soil COS uptake. We notice that this underestimation 849 lower net soil COS uptake with a constant atmospheric COS concentration can be found as early as March over 850 Europe, where atmospheric COS concentration is higher in this region. In eastern Asia, where atmospheric COS 851 concentration is higher than 800 ppt, the underestimation decreaseo fin the net soil COS uptake can reach -2.34 852 pmol m -2 s -1 when considering a constant atmospheric COS concentration. 853 It is to be noted that the modelled COS concentrations we used have their own uncertainty, which is however 854 smaller than their difference with the fixed value (Remaud et al., 20212021). 20212022). 2021). 855 24 a mis en forme : Justifié

Foreseen improvements 856
The mechanistic representation of soil COS fluxes was found to be in better agreement with the observations at 857 field sites. However, there can be strong differences between the simulated fluxes and the observations at some 858 sites, especially at AT-NEU and ES-LMA. In the mechanistic approach, the influence of light on soil COS fluxes 859 is not considered. Several field studies have reported light-induced emissions in oxic soils (Kitz et al., 2017;860 Meredith et al., 2018;Spielmann et al., 2019;Whelan and Rhew, 2015), assumed to be related to the effect of light 861 on soil organic matter. Spielmann et al. (2019) related strong soil COS emissions during daytime to light at the 862 sites where direct solar radiations reached the surface, such as ES-LMA and AT-NEU. At these sites, the 863 mechanistic model was unable to represent the soil COS emission peak during daytime. The optimization we 864 performed showed that, as expected, adjusting the parameters to site observations improves the fit between the 865 simulated and observed fluxed. However, it is necessary to represent all important processes in the mechanistic 866 approach before calibrating the parameters. Thus, a next step in our modelling approach could be to include the 867 light influence on soil COS fluxes, which can be of major importance for the sites where radiations strongly affect 868 soil COS fluxes. Mellillo and Steudler (1989) Several studies also found that soil COS production could be related 869 to nitrogen content, which increases with nitrogen fertilizer application (Kaisermann et al., 2018;Meredith et al. 870 2018Meredith et al. 870 , 2019. At the sites where soil is enriched with nitrogen inputs, such as agricultural fields or managed and 871 fertilized grasslands and forests, the fertilization Then crop management practices might would also need to be 872 included when representing the dynamics of soil COS fluxes. However, the soil nitrogen content and soil microbial 873 nitrogen biomass vary not only with fertilization, but also with location. Then, in addition to indications on land 874 use, information on the total soil nitrogen content should be included in the model to consider nitrogen impact on 875 soil COS flux. In the soil COS models, the impact of snow cover is also not represented. Indeed, due to the scarcity 876 of soil COS flux observations in winter and with snow cover, its effect on soil COS flux could not be implemented 877 in soil COS models yet. However, Helmig et al. (2009) found that COS uptake was not zero when soil is covered 878 by snow at Niwot Ridge, Colorado. 879 Moreover, one difficulty with the study of soil COS fluxes arises from the scarcity of field measurements that 880 could be used for data assimilation. Therefore, more field measurements would help to build a larger field 881 observation database for model validation and calibration. Besides, the observation sites considered here are all 882 located in a small latitudinal range between 39°N and 62°N. Measurements in the tropics and in the Southern 883 hemisphere are needed. Especially, soil COS flux observations in Northern India could help to validate the net soil 884 COS production simulated in both SiB4 and ORCHIDEE. In the tropical rainforest, soil COS flux measurements 885 were performed at La Selva Biological Station in Costa Rica (personal communication). In particular, t In the 886 tropical rainforest, soil COS flux measurements were performed at La Selva Biological Station in Costa Rica (Sun 887 et al., 2014). When available, these measurements could allow a first comparison between the observed and 888 simulated soil COS flux in a tropical region. 889 Then, the characterization of the soil microbial community should also be addressed to improve the scaling of CA 890 content and activity, represented by the parameter (Meredith et al., 2019). 891 The mechanistic model from Ogée et al. (2016) has also recently been implemented in the LSM SiB4 (Kooijmans 892 et al., 2021). The implementation of the soil COS flux mechanistic model from Ogée et al. (2016) iIn SiB4 893 (Kooijmans et al., 2021), the simulated soil COS fluxes with the mechanistic model shows a seasonal cycle with 894 a maximum net soil COS uptake in summer for the sites without crops, while the fluxes computed in ORCHIDEE 895 25 a mis en forme : Justifié show almost no seasonality. The expression of the production term differs between the two models, which is 896 based on Meredith et al. (2018) in SiB4 and on Whelan et al. (2016) in ORCHIDEE. The observation sites that are 897 common to the two studies (FI-HYY, US-HA, AT-NEU and DK-SOR) are also represented by different fractions 898 of biomes between SiB4 and ORCHIDEE, which changes the parameterization to compute soil COS fluxes. 899 Finally, the parameter values for the enhancement factor for grass differ as the value for tropical grass is also 900 assigned to C3 and C4 grass in SiB4. Soil COS flux field data are mainly available in summer, therefore having 901 field measurements over a whole year could better inform the seasonality of observed soil COS fluxes to compare 902 to the simulations. 903 The optimization does not modify the respective seasonality of both soil COS models, with a seasonal cycle that 904 agrees with the one of soil respiration for the empirical model and a lack of seasonality for the mechanistic model. 905 The lack of observations in winter does not enable to validate or constrain soil COS fluxes in winter. Therefore, 906 having field observations over a whole year could help to determine if both models could be calibrated with a 907 constrain over the whole year instead of only during summer and autumn. Moreover, the optimized set of 908 parameters for the empirical models leads to a degradation of the simulated soil water content compared to the 909 observations at FI-HYY, while the optimized parameters of the mechanistic model improve the representation of 910 soil water content at US-HA and FI-HYY. Thus, the mechanistic approach is to be preferred over the empirical 911 model and should be selected for future COS studies in ORCHIDEE. 912 The sensitivity analyses showed the importance of the hydrology-related parameters in the computation of soil 913 COS fluxes with the mechanistic model. Thus, assuming an accurate representation of soil COS fluxes, soil COS 914 fluxes could have the potential to add a new constraint on hydrology-related parameters. 915 In this work, soil COS fluxes are computed in the top 9 cm, which assumes that soil COS uptake and production 916 depend on the conditions in the first soil layers. Indeed, soil COS uptake depends on diffusive supply of COS from 917 the atmosphere. However, since soil COS production does not depend on COS supply, deeper soil layers could 918 also contribute to soil COS production. A study by Yang et al. (2019) presents COS profile measurements in an 919 orchard, which shows a non-zero COS concentration in deeper soil layers, but no direct evidence for attributing it 920 to soil COS production. Thus, we could consider deeper soil layers in the future to study the impact on soil COS 921 fluxes compared to considering only the top soil layers. 922 The anoxic soil map of regularly flooded wetlands from Tootchi et al. (2019) enables to approximate the spatial 923 distribution of anoxic soil. However, in our approach, seasonality is only represented through soil temperature 924 seasonality. Anoxic soil temporal dynamic was initially included in the model described by Ogée et al. (2016) with 925 the soil redox potential but is not implemented in land surface models such as ORCHIDEE yet. We could also 926 refine our approach by distinguishing between the different types of wetlands and define a Pref value for each 927 wetland type instead of a global value of 10 pmol COS m -2 s -1 . . Then, a distinction could also be made for anoxic 928 soil COS fluxes from boreal peatlands, as Meredith et al. (2019) give a value of specific to this biome. 929 Moreover, indirect COS emissions from DMS oxidation in anoxic soils have been reported (Kettle et al., 2002;930 Watts, 2000) but are not represented in this study. Finally, the anoxic map used here represents 9.7% of the global 931 land area, but the distribution of anoxic soils can greatly vary depending on the study (between 3% and 21%, 932 Tootchi et al., 2019). Therefore, it would also be interesting to investigate the impact of anoxic soil coverage on 933 soil COS flux uncertainty. The soil COS budget at global scale over the 2009-2016 period is -30 GgS yr -1 , resulting from the contribution of 943 oxic soils that represent a net sink of -126 GgS yr -1 , and of anoxic soils that represent a source of +96 GgS yr -1 . It 944 is to be noted that the contribution from anoxic soils, while leading to a similar global budget to Launois et al. 945 (2015), has a different spatial distribution based on the repartition of regularly flooded wetlands from Tootchi et 946 al. (2019). This repartition seems more accurate as it also includes anoxic soil COS flux in the tropical region and 947 considers a larger variety of anoxic soils, such as salt marshes and rice paddies.