Articles | Volume 20, issue 13
Research article
05 Jul 2023
Research article |  | 05 Jul 2023

Optimizing the carbonic anhydrase temperature response and stomatal conductance of carbonyl sulfide leaf uptake in the Simple Biosphere model (SiB4)

Ara Cho, Linda M. J. Kooijmans, Kukka-Maaria Kohonen, Richard Wehr, and Maarten C. Krol

Carbonyl sulfide (COS) is a useful tracer to estimate gross primary production (GPP) because it shares part of the uptake pathway with CO2. COS is taken up in plants through hydrolysis, catalyzed by the enzyme carbonic anhydrase (CA), but is not released. The Simple Biosphere model version 4 (SiB4) simulates COS leaf uptake using a conductance approach. SiB4 applies the temperature response of the RuBisCo enzyme (used for photosynthesis) to simulate the COS leaf uptake, but the CA enzyme might respond differently to temperature. We introduce a new temperature response function for CA in SiB4, based on enzyme kinetics with an optimum temperature. Moreover, we determine Ball–Woodrow–Berry (BWB) model parameters for stomatal conductance (gs) using observation-based estimates of COS flux, GPP, and gs along with meteorological measurements in an evergreen needleleaf forest (ENF) and deciduous broadleaf forest (DBF). We find that CA has optimum temperatures of 20 C (ENF) and 36 C (DBF), which is lower than that of RuBisCo (45 C), suggesting that canopy temperature changes can critically affect CA's catalyzation activity. Optimized values for the BWB offset parameter are similar to the original value (0.010 ± 0.003 mol m−2 s−1), and optimized values for the BWB slope parameter (ENF: 16.4, DBF: 11.4) are higher than the original value (9.0) at both sites. The optimization reduces prior errors on all parameters by more than 50 % at both stations. We apply the optimized gi and gs parameters in SiB4 site simulations, thereby improving the timing and peak of COS assimilation. In addition, we show that SiB4 underestimates the leaf humidity stress under conditions where high vapor pressure deficit (VPD) should limit gs in the afternoon, thereby overestimating gs. Furthermore, global COS biosphere sinks with optimized parameters show smaller COS uptake in regions where the air temperature is over 25 C, mostly in the tropics, and larger uptake in regions where the temperature is below 25 C. This change corresponds with reported deficiencies in the global COS fluxes, such as missing sinks at high latitudes and required sources in the tropics. Using our optimization and additional observations of COS uptake over various climate and plant types, we expect further improvements in global COS biosphere flux estimates.

1 Introduction

The leaf assimilation of the atmospheric trace gas carbonyl sulfide (COS) has been suggested as a proxy to overcome the limitations of estimating photosynthetic carbon dioxide (CO2) assimilation (Whelan et al., 2018). Observations of the net ecosystem exchange (NEE) of CO2 include both gross primary production (GPP) and ecosystem respiration, and those two individual components cannot be directly observed during daytime. COS follows the same diffusional pathway into leaves through plant stomata as CO2. COS is then destroyed through hydrolysis catalyzed by the enzyme carbonic anhydrase (CA) and is assumed not to be produced by any process within leaves (Protoschill-Krebs et al., 1996; Stimler et al., 2010). The CA chemistry is not light dependent (Protoschill-Krebs et al., 1996), in contrast to photosynthetic CO2 fixation, which requires light. Therefore, when the CA activity is accurately quantified, measurements of COS uptake can provide information on stomatal conductance (Kooijmans et al., 2017).

Atmospheric COS mole fractions vary around 500 parts per trillion (ppt) and are primarily influenced by biosphere uptake, ocean emissions, and anthropogenic emissions (Kettle et al., 2002). Depending on the environmental conditions, soils can act as a COS source or sink (Maseyk et al., 2014; Whelan et al., 2016). Recent studies have found that a source is missing in the tropical region (Berry et al., 2013; Glatthor et al., 2015; Kuai et al., 2015; Ma et al., 2021). Moreover, Berry et al. (2013) and Hu et al. (2021) showed that a sink is missing, or a source is overestimated at higher latitudes. These findings ask for careful evaluation of all sources and sinks, including the biosphere.

Biosphere models, such as the Simple Biosphere model, version 4 (SiB4) (Berry et al., 2013; Kooijmans et al., 2021) and the Organizing Carbon and Hydrology In Dynamic Ecosystems model (ORCHIDEE; Launois et al., 2015; Maignan et al., 2021; Remaud et al., 2022; Abadie et al., 2022) have been used to estimate ecosystem exchange of COS quantitatively. The SiB4 COS biosphere exchange was recently assessed against observations by Kooijmans et al. (2021). They stressed the need to account for spatial and temporal variations in atmospheric COS mole fractions, which largely reduce SiB4 COS biosphere uptake in the tropics (although observations to confirm this influence are lacking). The calculated reduction in the tropics was not large enough to explain the gap in the COS budget. Kooijmans et al. (2021) and Vesala et al. (2022) also found that SiB4 COS biosphere flux simulations were low compared to observations in the boreal region, consistent with the underestimations found by Ma et al. (2021). Our study follows one of the recommendations in Kooijmans et al. (2021) by focusing on the parameterization of the temperature dependence of the CA enzyme activity to improve simulations of the vegetation COS uptake in SiB4.

In SiB4, the COS assimilation is described as a series of resistances (i.e., inverse conductances) at the leaf boundary layer (gb), the stomatal pores (gs), and the leaves' interior (gi). The gb and gs of COS are scaled relative to conductances for water vapor or CO2 with diffusivity ratios and a calibration factor. For gi, previous studies found that both the CA enzyme activity (Badger and Price, 1994) and mesophyll conductance (Evans et al., 1994) scale with the maximum velocity of carboxylation by the enzyme RuBisCo (Vmax, rub). Therefore, the COS internal conductance in SiB4 is scaled to Vmax, rub through a single calibration factor α based on laboratory leaf gas exchange measurements (Stimler et al., 2010, 2011; Berry et al., 2013). However, the enzymatic control of COS and CO2 assimilation differs. COS molecules are hydrolyzed by the enzyme CA in the mesophyll cells (Protoschill-Kreb et al., 1996). In contrast, photosynthesis is further controlled by the enzyme RuBisCo. Thus, CO2 has a different point of uptake compared to COS. The enzyme activity depends on the enzyme abundance and is related to environmental parameters such as temperature and pH (Michaelis and Menten, 1913). In particular, the CA enzyme does not require light to catalyze COS hydrolysis, whereas the RuBisCo enzyme does require light (Stimler et al., 2010). Different temperature responses of RuBisCo and CA were reported by Boyd et al. (2015) with the C4 plant Setaria viridis. They measured that Vmax, rub increased with temperature in the range 10 to 40 C, whereas the CA activity decreased above 30 C. Currently, however, there is limited information about the temperature response function of CA.

Several studies found that the leaf relative uptake ratio (LRU; which is proportional to the ratio of COS and CO2 deposition velocities) varies with temperature under conditions where light was not limiting photosynthesis (Cochavi et al., 2021; Stimler et al., 2010; Sun et al., 2018; Kooijmans et al., 2019). More specifically, the LRU decreased with increasing temperatures above 15 C, indicating that COS uptake has a lower optimum temperature than CO2 uptake, possibly driven by different temperature responses of the CA and RuBisCo enzymes. Therefore, to accurately simulate the relation between COS and CO2 exchange in leaves, it is necessary to use separate temperature response equations for the internal conductance to CO2 and COS.

Besides uncertainties in gi, uncertainties in gs can also affect the accuracy of simulated COS assimilation. A common approach for simulating gs is the semi-empirical Ball–Woodrow–Berry (BWB) model (e.g., Ball et al., 1987; Ball, 1988; Collatz et al., 1992). This model is also applied in SiB4 and utilizes a set of related variables (e.g., photosynthesis, relative humidity, and CO2 concentration at the leaf surface) and two empirical constants. One of the constants (b1) describes the slope of the relation between gs and GPP. The other constant (b0) represents the residual gs in the dark. The current implementation of the BWB model in SiB4 has only one pair of b1 and b0 values for C3 plants and only one pair for C4 plants, whereas the BWB constants should ideally be prescribed for each plant functional type (PFT) separately to obtain accurate gs (Miner et al., 2017). To constrain b0 requires information on nighttime gs. However, obtaining gs estimates from nighttime water vapor flux measurements in the field is highly uncertain due to observational constraints (Papale et al., 2006; Wehr et al., 2017; Wehr and Saleska, 2021). As an alternative, nighttime COS uptake was previously reported (White et al., 2010; Belviso et al., 2013; Commane et al., 2013, 2015; Berkelhammer et al., 2014; Billesbach et al., 2014; Wehr et al., 2017; Kooijmans et al., 2017), and when the soil uptake is properly accounted for, this flux could provide information on stomatal opening. Several multi-year measurement datasets of CO2 and COS biosphere and soil fluxes are now available (Commane et al., 2015; Wehr et al., 2017; Vesala et al., 2022). Multi-year datasets make it possible to distinguish valid signals from noise and to use COS to provide information on gs and constrain the BWB model parameters.

This research aims to optimize the temperature response of CA and BWB model parameters to better estimate COS assimilation in the SiB4 model. To do so, we will use eddy covariance (EC) measurements of the COS leaf flux, GPP derived from NEE, and gs derived from the EC COS flux. The optimization will be based on observations from two PFTs: a boreal evergreen needleleaf forest (ENF) at Hyytiälä, Finland, and a temperate deciduous broadleaf forest (DBF) at Harvard Forest, USA. The optimized parameters will be applied in a global simulation of the SiB4 biosphere model to evaluate the effects on the global COS biosphere sink.

2 Methodology

2.1 Modeling COS leaf uptake

2.1.1 SiB4 biosphere model

The SiB4 model is a prognostic land surface model that calculates the COS flux as described in Berry et al. (2013). The main application of the model is to estimate land–atmosphere exchange of carbon, energy, and water budgets (Sellers et al., 1986; Sato et al., 1989). SiB4 has a time step of 10 min and operates on a spatial resolution of 0.5× 0.5. Unlike the previous SiB3 model, which relies on satellite information to specify the time-varying phenological leaf state, version 4 fully simulates the terrestrial carbon cycle using a process-based model (Haynes et al., 2019).

As each vegetation type has different physiological and phenological characteristics, SiB4 simulates photosynthesis in a heterogeneous land cover with different plant functional types (PFTs) per site or grid cell, each with separate fractions. These PFTs consist of nine natural vegetation classes and three specific crop types (maize, soybeans, and winter wheat), plus the separation of C3 and C4 plants in generic cropland and grassland. Besides responses of plant growth to temperature, humidity, radiation, and precipitation, the model accounts for environmental stress factors as a limitation to plant growth: the leaf humidity stress (FLH), the root-zone water stress (FRZ), and the canopy temperature stress (FT). Several variables (e.g., Vmax, rub) are prescribed according to phenological stages: leaf out, growth, maturity, senescence, and dormant stages. The leaf-out stage begins when the environmental conditions are suitable for photosynthesis to take place, and the growth stage is determined when the canopy is large enough to support photosynthesis. The maturity starts when the leaf amount is maintained. When plants experience stress and photosynthetic capacity is reduced, it is prescribed as senescence. In the dormant stage, plants do not have leaves in the canopy, or conditions are unsuitable for photosynthesis (Haynes et al., 2020).

2.1.2 Module for COS vegetation uptake in SiB4

SiB4 simulates COS vegetation assimilation as a combination of three conductances from the laminar boundary layer to the chloroplast (gb, gs, and gi) multiplied by the atmospheric COS mole fraction (Berry et al., 2013):

(1) F COS = C COS 1.94 g s + 1.56 g b + 1 g i - 1 ,

where FCOS is the COS vegetation assimilation in the canopy (pmol m−2 s−1), and CCOS is the COS mole fraction in the canopy air space (pmol mol−1). The factors 1.94 and 1.56 account for the smaller diffusivity of COS with respect to H2O through the boundary layer and stomatal pores, respectively (Seibt et al., 2010; Stimler et al., 2010). Note that gi includes all conductances downstream of the stomata, such as the mesophyll conductance. Within SiB4, the aerodynamic conductance is used to connect the mole fraction in the canopy air space to the atmosphere.

The stomatal conductance gs (mol m−2 s−1) in SiB4 is calculated by using the BWB model. This model relates gs and GPP as a function of environmental factors with two empirical constants b0 and b1:

(2) g s = b 1 GPP SiB 4 CO 2 s F LH + b 0 LAI F RZ ,

where GPPSiB4 (mol C m−2 s−1) is the canopy CO2 assimilation, CO2s (mol C mol air−1) is the CO2 mole fraction at the leaf surface, FLH (–) is the leaf humidity stress factor, LAI is the leaf area index (–), and FRZ is a non-dimensional term that accounts for root-zone water stress. FLH is related to relative humidity at the leaf surface and is calculated as a ratio of the water vapor mixing ratio at the leaf surface to the water vapor mixing ratio in the leaf internal space (Sellers et al., 1992). The value of FLH for ENF has a lower bound of 0.7, making ENF more resilient to humidity stress. However, Smith et al. (2020) found that with the 0.7 threshold in place, SiB4 did not accurately simulate the drought response for European ENF ecosystems. Therefore, we removed this lower bound in the optimization but will show the impact in a sensitivity study in Sect. 3.5.1.

The empirical constant b1 is the slope of the linear relationship between gs and GPPSiB4; and FLH, CO2s-1, and b0 (mol m−2 s−1) is the intercept indicating minimum gs (Ball et al., 1987; Ball, 1988). The choice for b1 significantly impacts simulated transpiration (Leuning et al., 1998; Lai et al., 2000; Bauerle et al., 2014) and is prescribed in SiB4 as 9.0 for C3 plants and 4.0 for C4 plants. The coefficient b0 is 0.01 mol m−2 s−1 for most PFTs but 0.04 mol m−2 s−1 for crops and C4 plants. The prescribed b0 term is converted from the leaf to the canopy scale by multiplying by LAI.

GPPSiB4 is explicitly calculated in SiB4 as the minimum of three assimilation rates limited by enzyme activity (wc), light (we), and carbon compound export (ws) (Haynes et al., 2020). The three rates are calculated by functions fc,e,s described in detail in Sellers et al. (1996a) depending on a canopy temperature (Tcan, K):


Here, pCO2i (Pa) is the internal partial pressure of CO2, pO2(T) (Pa) is the temperature response of partial pressure of O2, APAR (mol m−2 s−1) is the absorbed photosynthetically active radiation, and γ (Pa) is the CO2 photo-compensation point. Note that GPPSiB4 is used in SiB4 to calculate the COS leaf flux via gs, as described in Eq. (2) and evaluated independently from GPP calculated by the BWB model (GPPBWB), which will be introduced in Sect. 2.3.1.

The COS molecules that have diffused into the leaf mesophyll cells are hydrolyzed in a reaction catalyzed by the CA enzyme (gi). Since the enzyme activity and mesophyll conductance are analogous and the terminal COS concentration is assumed to be zero, SiB4 presumes that the two conductances can be combined and that the gi (mol m−2 s−1) scales with Vmax, rub at 298 K (mol m−2 s−1) as follows (Berry et al., 2013):

(6) g i ( PS ) = α V max , rub ( PS ) f T can SiB 4 F LC ( PS ) F RZ P P sfc T can T 0 .

Vmax, rub varies with phenological stage (PS) (see Table 3) and is scaled with a Tcan response function f(Tcan)SiB4 that prescribes the relative increase per 10 K increase (Q10) as 2.1 as follows:

(7) f T can SiB 4 = 2.1 0.1 T can - 298 .

In SiB4, the canopy temperature Tcan is calculated from the temperature above the canopy using the leaf surface energy balance (Sellers et al., 1996b), and Tcan is normally obtained from a meteorological analysis dataset. In this study, however, we use the air temperature measured above the canopy to obtain Tcan needed by SiB4. Likewise, we use the specific humidity measured above the canopy, which is used by SiB4 to calculate the leaf humidity stress factor FLH at leaf surface level, needed in Eq. (2).

Other modifying factors in Eq. (6) are the ratio of atmosphere pressure (P; hPa) to the surface pressure (Psfc= 1000 hPa) and the ratio of the temperature to the reference temperature (T0= 273.15 K). FLC (–) is the scaling factor from leaf to canopy, accounting for a fraction of absorbed photosynthetically active radiation (FPAR) and other factors such as light scattering and leaf projection. The calibration parameter α (–) was obtained from simultaneous measurements of COS and CO2 uptake (Stimler et al., 2010, 2012; Berry et al., 2013) and was estimated as 1400 for C3 and 8862 for C4 plants. These numbers were derived from a limited number of observations, so the values of α do not capture variability between plant species and seasons. Kooijmans et al. (2021) derived α from ecosystem observations of six sites throughout the growing season and found an average α of 1616 ± 562 (C3 plants). Here, the standard deviation indicates large variability over time and between sites. The impact of α on gi will be described in Sect. 2.1.3.

2.1.3 A new approach to describe gi

Each enzyme has its own kinetic characteristics, with activity generally increasing with temperature up to an optimum temperature and decreasing above this temperature. To derive a more realistic enzyme activity that also accounts for an optimum temperature, we propose a temperature response (f(Tcan)new) based on an Arrhenius-type equation that applies Michaelis–Menten kinetics. The Arrhenius equation has been used for Vmax, rub and maximum rate of photosynthetic electron transport to estimate GPP (e.g., Dreyer et al., 2001; Galmés et al., 2016). A similar model was previously used in COS soil models (Sun et al., 2015; Ogée et al., 2016). The equation is described as (Peterson et al., 2004; Daniel et al., 2010)

(8) f T can new = A T T can exp - Δ H a R T can 1 + exp - Δ H eq R 1 T can - 1 T eq .

Here, three variables for enzyme kinetics are included: ΔHa (J mol−1) is the activation free energy of the CA enzyme, ΔHeq (J mol−1) is the enthalpy change when the enzyme converts from an activated to inactivated state, and Teq (K) is the temperature at which activated and inactive enzymes' concentrations are equal (Daniel et al., 2010; Sun et al., 2015). The factor AT normalizes Eq. (8) such that, equivalent to Eq. (7), f(Tcan)new= 1 at Tcan= 298 K. We adopt AT as the value of f(Tcan)new-1 when T is equal to Teq. R is the universal gas constant (8.3145 J K−1 mol−1). Figure 1 shows that α and the three kinetic parameters have different effects on the temperature response of gi (Eq. 6). The calibration parameter α affects the strength of gi (Fig. 1a), and its accuracy is therefore crucial for accurate COS flux simulations. With ΔHa increasing (Fig. 1b), gi decreases (increases) for temperatures above (below) the optimal temperature. ΔHeq has the opposite effect, albeit with a different response to ΔHa (Fig. 1c). Both ΔHa and ΔHeq affect gi depending on the temperature range. Finally, Fig. 1d shows that Teq determines the optimum of the temperature response curve without having impact on the magnitude of gi.

Figure 1Calculated gi as a function of canopy temperature for parameters (a) α, (b) ΔHa, (c) ΔHeq, and (d) Teq from Eqs. (6) and (8). Each parameter is set to five different values given in the caption to investigate the response of gi to temperature. While the target parameter changes, the other variables are fixed as α= 1400, ΔHa= 40 kJ mol−1, ΔHeq= 100 kJ mol−1, and Teq= 295 K, which are initial values for Hyytiälä (see Sect. 2.3.2).


2.2 Observations

In optimizing the parameters gs and gi, we used the following variables obtained from observation to calculate COS leaf uptake (Eq. 1): the COS ecosystem flux, the COS soil flux, CCOS, temperature, and specific humidity, as well as GPP partitioned from NEE measurements. These data were collected and derived at Hyytiälä in Finland during 2013–2017 (Kooijmans et al., 2017; Sun et al., 2018; Vesala et al., 2022) and at Harvard Forest in the United States during 2012 and 2013 (Commane et al., 2015, 2016; Wehr et al., 2017). To validate the optimization results, we used the observation-based gs and gi (Sect. 2.2.2).

2.2.1 COS flux, GPP, and mixing ratio

We used canopy COS uptake derived from COS EC measurements for Hyytiälä (Kohonen et al., 2020; Vesala et al., 2022) and Harvard Forest (Wehr et al., 2017). The effect of storage in the canopy airspace was included by collocated COS profiles (Kooijmans et al., 2017; Kohonen et al., 2020).

GPP at Hyytiälä has been obtained from NEE using multi-year parameter fits (Kolari et al., 2014; Kohonen et al., 2022). For Harvard Forest, we chose to use the GPP derived from the isotope spectrometer measurements because it is more accurate and reliable with frequent and rigorous calibrations (Wehr et al., 2016).

COS soil flux measurements were available for the 2016 growing season at Hyytiälä and for the 2012 and 2013 growing seasons at Harvard Forest. For the soil flux in other years at Hyytiälä, we applied the monthly average diurnal cycle of the soil flux from 2016 to the other years (2013–2015 and 2017). The seasonal and diurnal variation of the soil flux is small compared to the total ecosystem uptake of COS (Sun et al., 2018). Hence, the averaged value of 2016 can be safely used for other years.

Figure 2Monthly diurnal variation of COS fluxes in 2016 at Hyytiälä (a) and 2012–2013 at Harvard Forest (b, c). Lines are median values of each 3 h period, and the filled areas indicate the interquartile range (25 to 75 percentile). Black: COS ecosystem flux, blue: soil flux, red: vegetation flux estimated as ecosystem minus soil flux.


To convert the data frequency of observations to SiB4's 3 h time resolution, we calculated the median value of each variable in each 3 h interval and for each month. We only used data points when more than three data points were present and when all variables required for the optimization were available. Figure 2 shows the resulting average diurnal cycle per month for COS ecosystem, soil, and vegetation fluxes (ecosystem flux minus soil flux). Note that positive fluxes indicate uptake. Again, we note that we use the averaged soil flux at Hyytiälä because its variability is much smaller than the leaf flux.

2.2.2 Conductances gs and gi

Observation-based gs was derived from sensible heat flux and evapotranspiration measurements using the flux gradient (FG) equations (Baldocchi et al., 1991; Wehr and Saleska, 2015, 2021). A key step in the derivation of gs is the estimation of transpiration from evapotranspiration. At Harvard Forest, transpiration was estimated by an empirical equation established during times of minimal non-stomatal evaporation (i.e., a few days after rain, removing mornings with dew evaporation), as described in Wehr et al. (2017). At Hyytiälä, we simply restricted our analysis to periods of minimal non-stomatal evaporation by eliminating data when the dew point was equal to or greater than the air temperature or when the accumulated precipitation for the past 2 d was more than 0.01 mm.

The FG approach leads to significant uncertainties for nighttime data because the leaf-to-air water vapor gradient is too small under stable conditions (Wehr et al., 2017). We thus excluded nighttime gs when the values were smaller than 0.05 mol m−2 s−1. To reduce the effect of random noise on gs, we used an average diurnal cycle (based on 3 h medians) for each month.

Observation-based gi was extracted by rewriting Eq. (1) as follows:

(9) g i = C COS F COS - 1.94 g s - 1.56 g b - 1 .

Here, we used the observation-based gs from the FG equation as discussed above and filtered observations of CCOS and FCOS. Additionally, we used simulated gb from SiB4, as we do not have observed gb available and as the value of gb only has a minor effect on FCOS, which will be further discussed in Sect. 3.1. Although outliers of observed gs, CCOS, and FCOS were removed already, a significant number of outliers in gi appeared because of error propagation. To avoid excessive noise, we only retained gi values in the interquartile range (25–75 percentile) of 3 h for each month.

2.3 Optimization

2.3.1 Procedure

In the optimization steps, we minimized a quadratic cost function J(x) based on Bayes' theorem (Tarantola and Vallette, 1982; Enting et al., 1993):

(10) J x = x - x a 2 2 σ a 2 + y - H x 2 2 σ y 2 .

Here, x represents the state, xa the prior settings of the state, and σa the error assigned to the parameters. In the second term, y represents the observations and H(x) the model evaluation using the state x. The error σy represents the observational error. The details of σa and σy will be described in Sect. 2.3.2 and 2.3.3, respectively.

To optimize the gs and gi parameters, we intend to use the information from GPP and COS leaf uptake measurements sequentially. Thus, we propose a two-step approach in combination with an iterative minimization of the cost functions, as outlined in Fig. 3. In the first step, we optimally estimate gs parameter b1 by minimizing J(x) which sums GPP differences between estimation (H(b1) in Eq. 10) and observation.

Figure 3Flow chart of the procedure to optimize COS leaf uptake's parameters. The procedure has two steps: (1) optimize b1 by minimizing deviations between GPPBWB and observations; and (2) optimize b0, α, and Teq by minimizing deviations between modeled and observed COS uptake. Variables highlighted in gray are from observations, and the other variables are estimated from SiB4.


We select GPP for the first step optimization rather than gs, because derived GPP from NEE has been evaluated more frequently than observation-based gs. We use only positive GPPobs values (uptake) because our target parameter b1 in the first step cannot be optimized when GPP is zero. Here, we do not use GPPSiB4 because SiB4 does not apply the BWB model for GPP calculation as described in Eqs. (3)–(5). For this reason, we cannot optimize BWB parameters with GPPSiB4. Instead, we estimated GPP by rewriting the BWB model using an observation-based gs (Sect. 2.2.2), modeled RH at the leaf surface (FLH), and simulated CO2s from SiB4. Hereinafter, the estimated GPP by the BWB model is called GPPBWB:

(11) GPP BB = g s - b 0 LAI F RZ CO 2 s b 1 F LH .

In the second optimization step, we optimize the b0 and gi parameters (α in Eq. 6 and Teq in Eq. 8). These parameters are optimized by minimizing the differences between calculated and observed FCOS. FCOS is calculated with three conductances using Eq. (1). Specifically, gs is estimated with Eq. (2) using the optimized b1 from step 1. Here, we used GPPSiB4 to satisfy our aim of optimizing the SiB4 model parameters. Note that GPPBWB from Eq. (11) cannot be used here because it would make the estimated gs equal to observation-based gs. Based on sensitivity studies in Appendix A, we decided to select α, b0, and Teq as target parameters and to fix ΔHeq and ΔHa at 100 and 40 kJ mol−1, respectively.

In the optimization procedure, we specifically exploit the fact that the nighttime COS flux carries information about nighttime gs through the parameter b0. The alternative, i.e., optimizing b0 already in step 1, would ignore the information of nighttime gs brought by COS flux observations. Consequently however, we have to iterate the procedure several times to reach convergence. Figure 3 specifies which observations and observation-based quantities are used in each step (gs, GPPobs, FCOS, obs, CCOS highlighted as grey) and which variables are simulated by SiB4 (e.g., Tcan,FLH, GPPSiB4, CO2s, gb).

We applied the simplicial homology global optimization (SHGO) from the SciPy python library to minimize the cost functions. SHGO is appropriate for solving non-continuous, non-convex, and non-smooth functions (Endres et al., 2018). SHGO also allows the definition of a valid parameter range, as will be discussed in Sect. 2.3.2 and in Appendix A.

The Vmax, rub was found to vary over the phenological stage and per PFT (Woodward et al., 1995; Wolf et al., 2006; Kattge et al., 2009; Walker et al., 2014), which also affects the calibration factor α. Therefore, we optimized α for each PFT and each phenological stage. In contrast, b0, b1, and Teq were only separately determined for the different PFTs, assuming local characteristics for each PFT. We did not include Vmax, rub in the state variables because this would require SiB4 CO2 simulations. These simulations need several parameters, like carbon cycle pools, which are difficult to estimate. Therefore, we focus this research on estimating Vmax, CA by optimizing gi-related parameters.

2.3.2 Initial parameters and prior errors

The first term in the cost function (Eq. 10) ties the values of the parameters to realistic values. We additionally confined the parameter values within realistic physical ranges using the SHGO algorithm. Initial parameters and prior errors were chosen based on thresholds outlined in Appendix A, and they will be compared with optimized results in Sect. 3.3. The variation in the resulting cost function shows distinct differences between Hyytiälä and Harvard Forest, which reinforces our strategy to optimize parameters for each station separately.

2.3.3 Observation errors

To quantify the observational errors σy, we first calculated the 3 h average coefficient of variation (CV) relative to the mean of the observed COS vegetation flux in each phenological stage and observed GPP for the entire growing season. Figure 4 shows the results of observational errors. The GPP error is applied in step 1, and the COS leaf uptake error is used in step 2 in the optimization. We multiplied the CV with the mean in each phenological stage. Here, we classify the error of the COS leaf uptake in each phenological stage because we optimized α in each stage. In Fig. 4, we found that the errors differ slightly per phenological stage. At Hyytiälä, the errors are larger in the growth stage compared to the maturity stage, possibly due to the unstable weather conditions in growth stage. The COS leaf uptake error is larger at both stations during nighttime than during daytime. A potential reason can be the relatively higher uncertainty in the EC method during stable nighttime conditions.

Figure 4Distribution of the observation error in the growth (a) and maturity stages (b) for COS leaf uptake and all stages for GPP (c). The upper panel is the errors for Hyytiälä (HYYT) and the lower panel is for Harvard Forest (HVFM).


2.4 SiB4 simulations

We utilized several simulated variables from SiB4 in our optimization. Specifically, calculated GPPSiB4, gb, Tcan, and Vmax, rub and functions FLH, FRZ, and FLC were used to calculate COS leaf uptake. In addition, LAI and CO2s were used to estimate GPPBWB. Furthermore, we introduced the new temperature function (f(Tcan)new) in the gi calculation (Sect. 2.2.3) to calculate COS leaf uptake and excluded PPsfc-1 and TcanT0-1 from Eq. (3) due to minor impacts of these factors for these ecosystems.

To simulate the vegetation assimilation FCOS at the two stations, we used the Modern-Era Retrospective Analysis for Research and Application, version 2 (MERRA-2) (Gelaro et al., 2017) as meteorological driver data. Only air temperature and leaf-specific humidity were taken from observations. To initialize the carbon pools, we spun up the model to equilibrate the pools. The spin-up was performed from 2000 to 2010 with 10 iterations. We used observed CCOS, and ambient CO2 mole fractions were prescribed at 370 ppm.

To estimate the global impact of our findings, we performed a global SiB4 simulation from 2016 to 2018 to evaluate the influence of the new parameters on the monthly COS biosphere fluxes which are averaged for 3 years. The atmospheric COS mixing ratio CCOS were taken from optimizations using the TM5 chemical transport model (Ma et al., 2021; Kooijmans et al., 2021). We used 3 h CCOS averaged over 2016 to 2018 by Kooijmans et al. (2021). As we found that all target parameters differ between ENF and DBF (Appendix A), the application of the optimized parameters to other PFTs will likely be incorrect. Hence, we applied the optimized parameters only to ENF and DBF and used the standard values of SiB4 for the other PFTs. However, to confirm the f(Tcan)new effect on COS leaf uptake, we applied f(Tcan)new to all PFTs with averaged optimum Teq from the two stations (303 K) and fixed ΔHeq (100 kJ mol−1) and ΔHa (40 kJ mol−1) as described in Appendix A. The soil flux is estimated following Ogée et al. (2016) as implemented by Kooijmans et al. (2021).

To examine the humidity stress impact in SiB4, we performed a simulation with and without the lower threshold for FLH of 0.7 for ENF (see Sect. 2.1.2). Additionally, we replaced the RH at leaf level calculated by SiB4 by RH measured above the canopy. Results will be shown in Sect. 3.5.1. To account for the optimized humidity impact on the global COS leaf uptake, we simulated the global COS leaf uptake without the 0.7 threshold of FLH for ENF.

2.5 Error reduction and statistics

To determine the uncertainty in the optimized model parameters, we employed a Monte Carlo optimization procedure as described in detail in Appendix B. In short, 100 optimizations were performed. In each optimization, we perturbed the state with random Gaussian noise on the state and the observations (Chevallier et al., 2007; Bosman and Krol, 2023), according to the errors in the state and observations (Fig. 4). Posterior error statistics will be reported in Table 3.

Additionally, we quantified the performance of the optimization by calculating the root mean square errors (RMSEs), mean bias errors (MBEs), and the chi-square metric (χ2). The χ2 metric quantifies the average deviation from the observations, expressed in σy units. Thus, χ2= 1 signals that, on average, the model fits the observation within 1σ indicating a realistic error setting.

3 Results and discussion

3.1 Impact of each conductance

Figure 5 investigates which conductance contributes most to the total conductance (gt). In these plots, all conductances are prior values before optimization. gs and gi were derived from observations (Sect. 2.2.2). We find that gt is determined mainly by gi and gs. During daytime, gi is the lowest conductance in almost all months at Hyytiälä but is comparable to gs in Harvard Forest. The value of gb is the highest and hence has the smallest impact on gt. Therefore, to improve the accuracy of COS leaf uptake simulation effectively, parameters of gs and gi are evaluated and optimized, and gb is kept to its standard value.

Figure 5Monthly median value of diurnal conductances (black: gt, red: gs, orange: gi, blue: gb) at Hyytiälä (HYYT) and Harvard Forest (HVFM). gb is estimated by SiB4, gs and gi are calculated based on observations as described in Sect. 2.2.2. Negative values are not displayed. The total conductance gt is calculated from gb, gs, and gi according to Eq. (1).


3.2 Optimization performance

We obtained optimized parameters after five iterations. By design, the optimized results reduced the deviations between model and observation of GPP and COS leaf uptake. This improvement is quantified by statistical indexes in Tables 1 and 2, respectively. GPPBWB is improved slightly compared to the prior (Table 1), with RMSEs reduction from 4.08 to 3.84 µmol m−2 s−1 at Hyytiälä and 8.35 to 7.89 µmol m−2 s−1 at Harvard Forest. MBEs are decreased from 1.48 to 0.61 µmol m−2 s−1 at Hyytiälä and from 4.97 to 3.46 µmol m−2 s−1 at Harvard Forest (Table 1). The χ2 was reduced by about 0.09 at Hyytiälä and by 0.08 at Harvard Forest. The improvement in GPPBWB reflects the effect of optimizing b1 and b0 in the BWB model.

Table 1RMSE, MBE, and χ2 for the estimation of GPPBWB in daytime using prior stomata parameters (pri) and posterior parameters (post).​​​​​​​

Download Print Version | Download XLSX

Table 2Same as Table 1, but for COS leaf uptake, as applied in the original f(Tcan)SiB4 simulation with original conductance parameters (org), with the new temperature response function f(Tcan)new using the initial gs and gi parameters (pri), and with optimized parameters (post). Here, state parameters relevant for gi are α and Teq. The state parameters relevant for gs are b0 and b1.

Download Print Version | Download XLSX

The posterior result of COS leaf uptake (“post” in Table 2) shows a slight improvement compared to the original-state variables with f(Tcan)new in RMSE (from 7.67 to 5.73 pmol m−2 s−1 at Hyytiälä and from 10.45 to 9.54 pmol m−2 s−1 at Harvard Forest) but significantly improved MBE (from 5.10 to 0.01 pmol m−2 s−1 at Hyytiälä and from 2.53 to 1.99 pmol m−2 s−1 at Harvard Forest, see “Post” in Table 2). The large RMSE reflects the typically large random noise of COS flux observations (Kooijmans et al., 2016; Kohonen et al., 2020). However, χ2 drops by 0.41 at Hyytiälä and by 0.38 at Harvard Forest, confirming that the optimization properly reduced the mismatch between observations and the model within the error statistics. Figure 6 compares the optimized COS leaf uptake to the original SiB4 simulation in scatter plots. Where the original simulation with f(Tcan)SiB4 and previous-state variables was often underestimating the observations, the optimized results resemble the observations over a larger range of the data.

Figure 6Scatter plots between observed and estimated COS leaf uptake from original parameters with f(Tcan)SiB4 (left) and optimized parameters with f(Tcan)new (right) for two stations. The colors represent the density of data.


3.3 Optimized parameters

The optimized parameter values with posterior errors are listed in Table 3. The optimized SiB4 parameters differ between the stations, likely because the dominant PFT and the climate conditions differ between Hyytiälä and Harvard Forest. For instance, the optimum temperature is smaller at Hyytiälä (19.85 K) than in Harvard Forest (35.85 K), which are slightly smaller than Teq. Thus, the optimum temperature reflects the temperature dependence of the enzyme and its adaptation to temperature (Lee et al., 2007). This indicates that regional temperature information is important for correctly estimating gi globally. The optimum temperature can be compared with other observations. For instance, Burnell and Hatch (1988) observed increasing CA activity with maize grown in a temperate temperature range from 20 to 30 C, relative to a temperature of 17 C. Thus, we can assume the optimal temperature lies above 17 C. Another study by Boyd et al. (2015) observed the C4 plant Setaria viridis with a temperature of 28 C/18 C day/night, and a reduced CA activity is suggested at temperatures above 25 C. This optimum temperature falls between our values derived for Hyytiälä and Harvard Forest.

Table 3Original (org) and optimized (post) state vectors for Hyytiälä and Harvard Forest in different phenological stages as defined by SiB4. Values of posterior in parentheses indicate posterior errors. The definition of the prior values is outlined in Appendix A, and the error reduction is described in Appendix B.

Download Print Version | Download XLSX

The α, which is the enzyme activity of CA relative to the Vmax, rub, is reduced from the default value of 1400 to 1316 (in growth) and 1331 (in maturity) at Hyytiälä. At Harvard Forest, α values are larger than the original values in SiB4 for leaf-out (1780), growth (1740), and maturity (2224) phenological stages. Here it should be noted that the change of α should be interpreted in combination with the new temperature function f(Tcan)new of gi. Since we only optimize Teq for two PFTs (with identical and fixed values for ΔHa and ΔHeq) and Teq only shifts f(Tcan)new (Fig. 1), the magnitude of gi is primarily determined by parameter α. The different values of α derived for different phenological stages will be discussed in Sect. 3.4.

The optimized results of the BWB model parameters b0 are similar to the original values used in SiB4, but b1 values are mostly higher. The parameter values b0 for Hyytiälä (0.013 mol m−2 s−1) and Harvard Forest (0.007 mol m−2 s−1) are slightly changed compared to the initial value (0.010 mol m−2 s−1). For the optimized BWB model parameter b1, the empirical slope between gs and GPP, we find a considerable increase at Hyytiälä (16.38) and a slight increase in Harvard Forest (11.43), compared to the prescribed SIB4 value of 9.0. Our optimized values are larger than the values presented in a review paper for the evergreen gymnosperm tree which showed b1= 6.8 and are similar to b1= 8.7 for the deciduous angiosperm tree (Miner et al., 2017). As will be discussed in Sect. 3.5, the higher slope at Hyytiälä is possibly related to an incomplete separation of observed transpiration rates from the latent heat flux.

Concerning the estimated errors in b0, b1, and Teq, we find that errors have been reduced significantly compared to the prior error range. This indicates that the available data constrain these parameters well. Only the α parameters of Harvard Forest are less well constrained. Also, the skill of the optimization to independently optimize the parameters is high, as quantified by the posterior covariances that are presented in Appendix B.

3.4 Optimized temperature response

The optimized parameters show significant improvement in temperature response of the COS leaf uptake. Figure 7 presents the temperature dependency of gi and COS leaf uptake from the original and optimized simulations output and observations. As stated before, the original f(Tcan)SiB4 describes the CA enzyme activity as an exponentially increasing response to temperature, which does not resemble the observations. The optimized gi and COS leaf uptake follow the temperature dependence of the observation more closely than the original f(Tcan)SiB4. In Harvard Forest, an underestimated bias is shown at a lower temperature under 10 C, mostly corresponding to nighttime. This underestimate is related to the uncertainty in nighttime gS and the small data volume at low temperatures (details in Sect. 2.2.2).

Figure 7Temperature dependency on gi (a) and COS leaf uptake (b) at Hyytiälä (HYYT; left) and Harvard Forest (HVFM; right). The lines are medians, and the filled area represents the 25 to 75 percentiles of each temperature range with 3 C intervals. Black: data based on observations; blue: previous parameters with f(Tcan)SiB4; red: optimized parameters with f(Tcan)new and gs parameters of the BWB model; and orange: same as the red line, but now α is prescribed with original value (1400) and not optimized. The numbers indicate number of observations in each temperature bin.


Figure 8Monthly diurnal cycle of COS leaf uptake (a), gs (b), and gi (c) at Hyytiälä (HYYT) from 2012 to 2016. Data include observations (black dots), the original SiB4 model with f(Tcan)SiB4 (solid blue line), and SiB4 with optimized gs and gi parameters with f(Tcan)new (solid red line). The filled area corresponds to the 25–75 percentile of the data in each 3 h interval of each month.


In the upper panel of Fig. 7, we see the different roles of α and f(Tcan)new in the improvement of gi response to temperature as the red and orange lines. Without the α correction applied (posterior with α= 1400; orange line), the optimized gi resembles the fluctuations in the observations, but there remains a bias in the amplitude. In contrast, when the optimized value of α is included (red line), the amplitude of gi is improved. Compared to the optimization that excluded α, the MBE is reduced from 0.006 to 0.003 mol m−2 s−1. Due to the different optimized α values in each phenological stage, the improvement of the red line shows the appropriate temperature responses. For instance, in Harvard Forest, α in leaf out and growth (1798 and 1740) mostly corresponds to lower temperatures. At these stages, the impact on gi is smaller because gi is smaller than that at high temperature. At high temperatures, there are more significant corrections of gi, which correspond to the maturity stage value of α (2224).

The temperature responses of gi and COS leaf uptake now show an optimum temperature. As can be observed in Fig. 7a, the optimum temperature of the observation-based gi is seen as 293 K (19.85 C) at Hyytiälä. At Harvard Forest, the optimum temperature is 309 K (35.85 C) but falls outside the observation range. The optimum curve affects the COS leaf uptake at both stations with changing peak temperatures (Fig. 7b). The accuracies of gi and FCOS are improved significantly at temperatures both below and above around the optimum temperature at both sites.

3.5 Application in SiB4

3.5.1 Monthly diurnal variation

Figures 8 and 9 display the SiB4 simulation results obtained with the original and optimized parameterizations compared to observations for Hyytiälä and Harvard Forest, respectively. As a result of the optimization, the monthly diurnal variation of the optimized COS vegetation flux, gs, and gi are closer to observations than the original SiB4 simulations. The observed COS leaf uptake and gs show diurnal and seasonal fluctuations at both measurement sites, with the highest values around midday and in summer. For gi, we observe a weak diurnal cycle throughout the year and higher daytime maximum values in summer, driven by the temperature dependence of CA.

Figure 9Same as Fig. 8 but for Harvard Forest (HVFM).


At Hyytiälä, COS leaf uptake in the original SiB4 model was underestimated during daytime in all months. The fluxes increased too slowly in the morning for all months (Fig. 8a). These issues are solved by optimizing the BWB model parameters and temperature response function. In the case of gs (Fig. 8b), the original SiB4 simulation showed the correct timing of the increase and decrease of gs in the morning and afternoon but underestimated the peak daytime values. The optimized model now better resembles the daytime gs values.

However, the model still overestimates gs in the late afternoon of summer months at Hyytiälä. We speculate that one of the reasons lies in an inaccurate humidity, or humidity stress in SiB4. Figure 10 shows a diurnal cycle of gs simulations averaged from April to August with different choices on how the humidity stress factor is treated (Sect. 2.5). When the default 0.7 threshold of humidity stress (FLH) in ENF is applied in SiB4, gs is overestimated in the afternoon at Hyytiälä (dotted blue line). When we removed the minimum threshold of FLH for ENF, gs simulations during midday are improved (note that the threshold was only implemented for ENF, not for DBF, and thus the blue dotted line is not visible for Harvard Forest in Fig. 10).

Figure 10Average diurnal cycle of gs at Hyytiälä (HYYT) and Harvard Forest (HVFM) from April to August. Data include observations (solid black line), SiB4 simulation with optimized parameters with minimum bounds of FLH (dotted blue line) and without bounds (dotted red line), and with observed RH in air (dashed orange line). Note that the blue and red lines overlap for HVFM.


However, SiB4 still tends to overestimate gs in the morning and late afternoon. The overestimated FLH in SiB4 can result from three factors: (1) an overestimated water vapor flux in the boundary layer to the leaf surface, (2) an underestimated boundary conductance, or (3) an underestimated leaf surface temperature. Since we do not have observations of the leaf surface temperature, we have confirmed that the estimated canopy temperature has a tight 1 : 1 relationship with the observed air temperature. We speculate that the main reason for the overestimated FLH is the uncertain water vapor flux. When we base the gs calculation on the observed RH above the canopy, the diurnal cycle is better simulated (dashed orange line in Fig. 10). The overestimated water vapor pressure implies that SiB4 tends to underestimate the humidity stress in the late afternoon when converting observed specific humidity above the canopy to humidity at leaf surface level. We suggest evaluating the boundary conductance (point 2 above) with observations.

The optimized model still underestimates gs at Hyytiälä in April, September, and October (Fig. 8b). This might indicate that we did not properly separate stomatal transpiration rates from the observed latent heat flux. The simulated mean ratios of evaporation to evapotranspiration in these 3 months are 66 %, 60 %, and 95 %, respectively, and these values are higher compared to the other months (43 % to 53 %). Thus, we speculate that the observed evapotranspiration does not solely represent stomatal transpiration in these months due to larger evaporation rates, leading to overestimated gs in the observations.

Figure 8c shows that the optimized gi often resembles the observed daytime gi better than the original SiB4 simulation. Only in April is the optimized gi overestimated at Hyytiälä. Again, this can likely be explained by the underestimated gs, which is used to derive observation-based gi (see Sect. 2.2.2).

Figure 11Monthly COS sink and averaged temperature in the mixed layer (Tm) over global (a) and specific regions (north boreal: b, north temperate: c, tropics: d) with the original (blue line) and the optimized (red line) SiB4 model. The grey bars represent the differences in COS sink between the original and the optimized model (right axis). The yellow bars are averaged temperatures in the mixed layer.


At Harvard Forest, the optimized SiB4 model generally simulates the magnitude of the COS leaf uptake well (Fig. 9a). The model overestimates the COS leaf flux only in the afternoon during the summer months. However, gs values are generally overestimated, and SiB4 simulates two peaks during daytime. This indicates that humidity stress is only briefly occurring at midday in SiB4. However, in reality, the humidity stress likely remains a limiting factor in the afternoon under conditions with high vapor pressure deficit (VPD) (Kooijmans et al., 2019). Observations show that gs typically peaks in the early morning and decreases in the afternoon due to higher afternoon VPD. Figure 10 shows that, similar to the Hyytiälä simulation, the afternoon decrease in gs at Harvard Forest is better simulated when we use the RH observed above the canopy. In Fig. 9c, the optimized gi during the daytime agrees well with the observation-based gi, except for several drops or peaks in July and October, likely caused by observational errors or uncertainty of the observed gs.

3.5.2 Global application

Figure 11 shows the SIB4 calculated changes in the monthly COS biosphere flux after applying the optimized temperature function and stomatal parameters. The global COS sink remains almost preserved (original: 701 Gg S yr−1, optimized: 704 Gg S yr−1), but the regional budgets change significantly. For example, the optimized model estimates larger COS uptake for all seasons in boreal and temperate regions and smaller uptake in the tropics. These changes are explained by the new temperature function of gi in Fig. 7. However, since the new temperature function is based on only two observation sites in the boreal and temperate regions, the calculated uptakes need more verifications with observations obtained in other areas and in different climate conditions, such as the tropics.

The higher uptake at high latitudes and lower uptake at the tropics are nevertheless consistent with inverse modeling results presented in previous studies (Ma et al., 2021; Hu et al., 2021) and would help towards closing the COS budget. Still however, the temperature response function and BWB parameters are now based on measurements of only two sites in only two biomes. With more measurements over different vegetation types, these parameters could also be optimized for a wider range of ecosystems.

4 Conclusion

To simulate more accurate COS leaf uptake in the SiB4 model, we have proposed a new temperature function f(Tcan)new for the CA enzyme and have optimized gs and gi parameters using observations in ENF (Hyytiälä) and DBF (Harvard Forest) systems. The optimized model reduced the MBE from 5.10 to 0.01 pmol m−2 s−1 at Hyytiälä and from 2.53 to 1.99 pmol m−2 s−1 at Harvard Forest. Furthermore, χ2 decreases by about 0.41 and 0.38, respectively.

The new function now considers an optimum temperature for enzyme activity, contrary to the initial temperature function used in SiB4 where an exponential increase of the temperature function was adopted from the RuBisCo enzyme activity. The new temperature function is characterized by an optimum temperature of 293 K (19.85 C) (Hyytiälä) and 309 K (35.85 C) (Harvard Forest). The new temperature response increases gi, and thereby the COS flux when the temperature is below the optimum temperature (mostly at high latitudes), and decreases the COS uptake at higher temperatures. (e.g., close to the Equator). Globally, these modifications help to close gaps in COS budget that were identified in earlier studies. In this study, we have interpreted the decreasing gi at higher temperatures as an optimum enzyme activity with the widely applied assumption that there are no COS emissions in leaves. However, COS emissions have recently been reported at high temperatures (Maseyk et al., 2014; Commane et al., 2016; Gimeno et al., 2017). To determine reasons for reducing COS leaf flux and internal conductance at high temperatures, it will be necessary to analyze the possibility that leaf emissions exist in observations in the future.

We have optimized the BWB model parameters for which we took advantage of the characteristics that the nighttime COS flux informs about nighttime gs and thus the parameter b0. The improved correspondence between model and observations shows that COS observations can help to constrain the relation between gs and GPP better. In addition, we showed that SiB4 underestimates the leaf humidity stress under conditions where high VPD should limit gs in the afternoon. This can be improved with more accurate relative humidity values and removing the threshold of humidity stress that was implemented in SiB4 specifically for ENF.

The optimized parameters show different values depending on the PFT. Therefore, extending our approach with more observations in different climate zones and over different PFTs will help obtain accurate COS fluxes on a global scale. This approach would reduce the uncertainty in the global COS budget and provide additional constraints on GPP.

Appendix A: State variable error settings

To evaluate the impact of the various parameters in f(Tcan)new in the optimization as state variables, we implemented a sensitivity test of a total cost function combined with cost1 and cost2, excluding the background term in the cost function equation (Eq. 10). Figure A1 shows the shape of the cost function when one parameter is varied within an acceptable range while the other parameters are fixed (Daniel et al., 2010; Sun et al., 2015) (details in Sect. 2.3.2). Based on the shape of the cost function, we used a pragmatic approach to select realistic parameter ranges. Variable values that push the cost function beyond 3.45 (Hyytiälä) and 5.14 (Harvard Forest) were considered outside the allowed physical range (red lines in Fig. A2). These thresholds are determined by the cost function value assuming that the modeled H(x) is the 75-percentile value of observation in 3 h observation in each month. Variables α, b0, b1, and Teq (Fig. A1a, b, c, and f) have more significant impacts on the cost function than ΔHa (Fig. A1d) and ΔHeq (Fig. A1e). Overall, costs in Harvard Forest are higher than at Hyytiälä, likely because DBF has larger diurnal and seasonal variations in the observed fluxes than ENF. We set the optimization range as an initial value ±1.5 state error to apply SHGO algorithm.

Figure A1Cost function values plotted against the value of the state vector elements at Hyytiälä (solid line) and Harvard Forest (dotted line). The red lines indicate a criteria cost calculated by H(x) as the 75-percentile value of every 3 h observation in each month. While the target parameter changes, the other variables are fixed as α= 1400 (Hyytiälä), 2000 (Harvard Forest), ΔHa= 40 kJ mol−1, ΔHeq= 100 kJ mol−1, Teq= 295 K (Hyytiälä), 310 K (Harvard Forest), b0= 0.02 (Hyytiälä), 0.01 (Harvard Forest), and b1= 17 (Hyytiälä), and 12 (Harvard Forest). These values were based on the value where the cost reached a minimum.


Figure A2Contour diagram of the cost function value as a function of Teq and (a) ΔHeq and (b) ΔHaat Hyytiälä (left) and Harvard Forest (right). While the target parameter changes, the other variables are fixed as Fig. A1.


Figure A2 shows contour diagrams of the cost function as a function of Teq and other parameters of f(Tcan)new. The gradient is the cost function indicates the relative importance of each parameter. ΔHeq does not interact with Teq, but ΔHa is inversely proportional to Teq to minimize the cost. The cost function is most sensitive to variations in Teq, and therefore we decided to fix ΔHeq and ΔHa at 100 and 40 kJ mol−1, respectively, and to base our optimization on the state variables α, b0, b1, and Teq.

Appendix B: Posterior uncertainties

To evaluate the ability of constrain the parameters, we performed an ensemble optimization with 100 different members. In each optimization, noise was added to the parameters (εa) and to the observation (εy). Random perturbations were drawn from a normal distribution with zero mean and standard deviations σa for the state parameters (x) and σy for the observations (y). The new cost function of an individual optimization thus becomes

(B1) J x = x - x a + ε a 2 2 σ a 2 + y + ε y - H x 2 2 σ y 2 .

We optimized each ensemble with the same y (observationally derived GPP and COS leaf uptake) and x but added noise to each ensemble member (Chevallier et al., 2007). Subsequently, we calculated the posterior uncertainty as the 1 standard deviation of the posterior distribution of the optimized parameters.

Figure B1 shows the prior and posterior distribution of the parameters at the two stations. All posterior parameters show considerable reductions of variations (error), with optimized values that are listed in the main text in Table 3.

Additionally, we calculated a correlation matrix between the posterior-state parameters at the two stations, which is shown in Fig. B2. Overall, each parameter does not interact significantly (covariances < 0.7).

Figure B1Error reduction of state variables in two stations (Hyytiälä (HYYT) and Harvard Forest (HVFM)). The red lines represent median values, and the boxes represent errors. Column “pri” shows the initial value and state error. Column “post” represents the mean of the optimized-state variables and the corresponding standard deviation. αLO, αG, and αM indicate α in each phenological stage (leaf out, growth, and maturity, respectively).


Figure B2Covariance matrix for all state variables at Hyytiälä (HYYT) and Harvard Forest (HVFM).


Code availability

The SiB4 code is available online at (SiB4 project members, 2020). TM5-4DVar model codes are available on the TM5-4DVAR website (; Krol et al., 2013).

Data availability

Observation data are downloaded from the original publications as mentioned in Sect. 2.2.

Author contributions

AC, LMJK, and MCK devised the study. AC optimized ecosystem parameters and analyzed the results with the consultation of LMJK and MCK. KMK and RW provided observation data and site-specific insights. AC wrote the paper, and all the authors provided comments.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank everyone that contributed to the collection of data for the Hyytiälä and Harvard Forest sites. The ecosystem dataset by eddy covariance at Hyytiälä was supported by ICOS Finland (319871) and the Atmosphere and Climate Competence Center (ACCC) Flagship (Vesala et al., 2022). The soil dataset at Hyytiälä was collected from Sun et al. (2018). Data from Harvard Forest are supported by the AmeriFlux Management Project (Wehr et al., 2017; Commane et al., 2015).

This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. We acknowledge computing resources from the Netherlands Organization for Scientific Research (NWO; grant no. NWO-2021.010).

Financial support

This research has been supported by the European Research Council, H2020 European Research Council (grant no. 742798).

Review statement

This paper was edited by Christopher Still and reviewed by Georg Wohlfahrt and one anonymous referee.


Abadie, C., Maignan, F., Remaud, M., Ogée, J., Campbell, J. E., Whelan, M. E., Kitz, F., Spielmann, F. M., Wohlfahrt, G., Wehr, R., Sun, W., Raoult, N., Seibt, U., Hauglustaine, D., Lennartz, S. T., Belviso, S., Montagne, D., and Peylin, P.: Global modelling of soil carbonyl sulfide exchanges, Biogeosciences, 19, 2427–2463,, 2022. 

Badger, M. R. and Price, G. D.: The Role of Carbonic Anhydrase in Photosynthesis. Annu. Rev. Plant Phys., 45, 369–392,, 1994. 

Baldocchi, D. D., Luxmoore, R. J., and Hatfield, J. L.: Discerning the forest from the trees: an essay on scaling canopy stomatal conductance, Agr. Forest Meteorol., 54, 197–226,, 1991. 

Ball J. T.: An analysis of stomatal conductance, Doctoral dissertation, Stanford University, (last access: 21 November 2022), 1988. 

Ball, J. T., Woodrow, I. E., and Berry, J. A.: A model predicting stomatal conductance and its contribution to the control of photosynthesis under different environmental conditions, in: Progress in Photosynthesis Research, edited by: Biggins, J., Springer, Dordrecht, 221–224,, 1987. 

Bauerle, W. L., Daniels, A. B., and Barnard, D. M.: Carbon and water flux responses to physiology by environment interactions: a sensitivity analysis of variation in climate on photosynthetic and stomatal parameters, Clim. Dynam., 42, 2539– 2554,, 2014. 

Belviso, S., Schmidt, M., Yver, C., Ramonet, M., Gros, V., and Launois, T.: Strong similarities between night-time deposition velocities of carbonyl sulphide and molecular hydrogen inferred from semi-continuous atmospheric observations in Gif-sur-Yvette, Paris region, Tellus B, 65, 20719,, 2013. 

Berkelhammer, M., Asaf, D., Still, C., Montzka, S., Noone, D., Gupta, M., Provencal, R., Chen, H., and Yakir, D.: Constraining surface carbon fluxes using in situ measurements of carbonyl sulfide and carbon dioxide, Global Biogeochem. Cy., 28, 161–179,, 2014. 

Berry, J., Wolf, A., Campbell, J. E., Baker, I., Blake, N., Blake, D., Denning, A. S., Kawa, S. R., Montzka, S. A., Seibt, U., Stimler, K., Yakir, D., and Zhu, Z.: A coupled model of the global cycles of carbonyl sulfide and CO2: A possible new window on the carbon cycle, J. Geophys. Res.-Biogeo., 118, 842–852,, 2013. 

Billesbach, D. P., Berry, J. A., Seibt, U., Maseyk, K., Torn, M. S., Fischer, M. L., Abu-Naser, M.​​​​​​​, and Campbell, J. E.: Growing season eddy covariance measurements of carbonyl sulfide and CO2 fluxes: COS and CO2 relationships in Southern Great Plains winter wheat, Agr. Forest Meteorol., 184, 48–55,, 2014. 

Bosman, P. J. M. and Krol, M. C.: ICLASS 1.1, a variational Inverse modelling framework for the Chemistry Land-surface Atmosphere Soil Slab model: description, validation, and application, Geosci. Model Dev., 16, 47–74,, 2023. 

Boyd, R. A., Gandin, A., and Cousins, A. B.: Temperature Responses of C4 Photosynthesis: Biochemical Analysis of Rubisco, Phosphoenolpyruvate Carboxylase, and Carbonic Anhydrase in Setaria viridis, Plant Physiol., 169, 1850–1861,, 2015. 

Burnell, J. N. and Hatch, M. D.: Low bundle sheath carbonic anhydrase is apparently essential for effective C4 pathway operation, Plant Physiol., 86, 1252–1256,, 1988. 

Chevallier, F., Bréon, F. M., and Rayner, P. J.: Contribution of the Orbiting Carbon Observatory to the estimation of CO2 sources and sinks: Theoretical study in a variational data assimilation framework, J. Geophys. Res., 112, D09307,, 2007. 

Cochavi, A., Amer, M., Stern, R., Tatarinov, F., Migliavacca, M., and Yakir, D.: Differential responses to two heatwave intensities in a Mediterranean citrus orchard are identified by combining measurements of fluorescence and carbonyl sulfide (COS) and CO2 uptake, New Phytol., 230, 1394–1406,, 2021. 

Collatz, G. J., Ribas-Carbo, M., and Berry, J. A.: Coupled photosynthesis-stomatal conductance model for leaves of C4 plants, Funct. Plant Biol., 19, 519–538,, 1992. 

Commane, R., Herndon, S. C., Zahniser, M. S., Lerner, B. M., McManus, J. B., Munger, J. W., Nelson, D. D., and Wofsy, S. C.: Carbonyl sulfide in the planetary boundary layer: Coastal and continental influences, J. Geophys. Res.-Atmos., 118, 8001–8009,, 2013. 

Commane, R., Meredith, L. K., Baker, I. T., Berry, J. A., Munger, J. W., Montzka, S. A., Templer, P. H., Juice, S. M., Zahniser, M. S., and Wofsy, S. C.: Seasonal fluxes of carbonyl sulfide in a midlatitude forest, P. Natl. Acad. Sci. USA, 112, 14162–14167,, 2015. 

Commane, R., Wofsy, S., and Wehr, R.: Fluxes of carbonyl sulfide at Harvard Forest EMS tower since 2010, Harvard Forest Data Archive: HF214 (v.10), Environmental Data Initiative,, 2016. 

Daniel, R. M., Peterson, M. E., Danson, M. J., Price, N. C., Kelly, S. M., Monk, C. R., Weinberg, C. S., Oudshoorn, M. L., and Lee, C. K.: The molecular basis of the effect of temperature on enzyme activity, Biochem. J., 425, 353–360,, 2010. 

Dreyer, E., Roux, X. L., Montpied, P., Daudet, F. A., and Masson, F.: Temperature response of leaf photosynthetic capacity in seedlings from seven temperate tree species, Tree Physiol., 21, 223–232,, 2001. 

Endres, S. C., Sandrock, C., and Focke, W. W.: A simplicial homology algorithm for Lipschitz optimisation, J. Global Optim., 72, 181–217,, 2018. 

Enting, I. G., Trudinger, C. M., Francey, R. J., and Granek, H.: Synthesis inversion of atmospheric CO2 using the GISS tracer transport model, Tech. Rep. 29, Division of Atmospheric Research Technical Paper, CSIRO, Australia, (last access: 21 November 2022), 1993. 

Evans, J. R., Caemmerer, S. V, Setchell, B. A., and Hudson, G. S.: The relationship between CO2 transfer conductance and leaf anatomy in transgenic tobacco with a reduced content of Rubisco, Funct. Plant Biol., 21, 475–495,, 1994. 

Galmés, J., Hermida-Carrera, C., Laanisto, L., and Niinemets Ü.: A compendium of temperature responses of Rubisco kinetic traits: variability among and within photosynthetic groups and impacts on photosynthesis modeling, J. Exp. Bot. 67, 5067–5091,, 2016. 

Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G.-K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S. D., Sienkiewicz, M., and Zhao, B.: The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2), J. Climate, 30, 5419–5454,, 2017. 

Gimeno, T. E., Ogée, J., Royles, J., Gibon, Y., West, J. B., Burlett, R., Jones, S. P., Sauze, J., Wohl, S., Benard, C., Genty, B., and Wingate, L.: Bryophyte gas-exchange dynamics along varying hydration status reveal a significant carbonyl sulphide (COS) sink in the dark and COS source in the light. New Phytol., 215, 965–976,, 2017. 

Glatthor, N., Höpfner, M., Baker, I. T., Berry, J., Campbell, J. E., Kawa, S. R., Krysztofiak, G., Leyser, A., Sinnhuber, B. M., Stiller, G. P., Stinecipher, J., and Von Clarmann, T.: Tropical sources and sinks of carbonyl sulfide observed from space, Geophys. Res. Lett., 42, 10082–10090,, 2015. 

Haynes, K., Baker, I., and Denning, S.: Simple Biosphere Model version 4.2 (SiB4) technical description, Mountain Scholar, 730 Colorado State University, Fort Collins, CO, USA, (last access: 21 November 2022), 2020. 

Haynes, K. D., Baker, I. T., Denning, A. S., Stöckli, R., Schaefer, K., Lokupitiya, E. Y., and Haynes, J. M.: Representing grasslands using dynamic prognostic phenology based on biological growth stages: 1. Implementation in the Simple Biosphere Model (SiB4), J. Adv. Model. Earth Syst., 11, 4423–4439,, 2019. 

Hu, L., Montzka, S. A., Kaushik, A., Andrews, A. E., Sweeney, C., Miller, J., Baker, I. T., Denning, S., Campbell, E., Shiga, Y. P., Tans, P., Siso, M. C., Crotwell, M., McKain, K., Thoning, K., Hall, B., Vimont, I., Elkins, J. W., Whelan, M. E., and Suntharalingam, P.: COS-derived GPP relationships with temperature and light help explain high-latitude atmospheric CO2 seasonal cycle amplification, P. Natl. Acad. Sci. USA, 118, e2103423118,, 2021. 

Kattge, J., Knorr, W., Raddatz, T., and Wirth, C.: Quantifying photosynthetic capacity and its relationship to leaf nitrogen content for global-scale terrestrial biosphere models, Glob. Change Biol., 15, 976–991,, 2009. 

Kettle, A. J., Kuhn, U., von Hobe, M., Kesselmeier, J., and Andreae, M. O.: Global budget of atmospheric carbonyl sulfide: Temporal and spatial variations of the dominant sources and sinks, J. Geophys. Res.-Atmos., 107, 4658,, 2002. 

Kohonen, K.-M., Kolari, P., Kooijmans, L. M. J., Chen, H., Seibt, U., Sun, W., and Mammarella, I.: Towards standardized processing of eddy covariance flux measurements of carbonyl sulfide, Atmos. Meas. Tech., 13, 3957–3975,, 2020. 

Kohonen, K.-M., Dewar, R., Tramontana, G., Mauranen, A., Kolari, P., Kooijmans, L. M. J., Papale, D., Vesala, T., and Mammarella, I.: Intercomparison of methods to estimate GPP based on CO2 and COS flux measurements, Biogeosciences Discuss. [preprint],, in review, 2022. 

Kolari, P., Chan, T., Porcar-Castell, A., Bäck, J., Nikinmaa, E., and Juurola, E.: Field and controlled environment measurements show strong seasonal acclimation in photosynthesis and respiration potential in boreal Scots pine, Front. Plant Sci., 5, 717,, 2014. 

Kooijmans, L. M. J., Uitslag, N. A. M., Zahniser, M. S., Nelson, D. D., Montzka, S. A., and Chen, H.: Continuous and high-precision atmospheric concentration measurements of COS, CO2, CO and H2O using a quantum cascade laser spectrometer (QCLS), Atmos. Meas. Tech., 9, 5293–5314,, 2016. 

Kooijmans, L. M. J., Maseyk, K., Seibt, U., Sun, W., Vesala, T., Mammarella, I., Kolari, P., Aalto, J., Franchin, A., Vecchi, R., Valli, G., and Chen, H.: Canopy uptake dominates nighttime carbonyl sulfide fluxes in a boreal forest, Atmos. Chem. Phys., 17, 11453–11465,, 2017. 

Kooijmans, L. M. J., Sun, W., Aalto, J., Erkkilä, K.-M. M., Maseyk, K., Seibt, U., Vesala, T., Mammarella, I., and Chen, H.: Influences of light and humidity on carbonyl sulfide-based estimates of photosynthesis, P. Natl. Acad. Sci. USA, 116, 2470–2475,, 2019. 

Kooijmans, L. M. J., Cho, A., Ma, J., Kaushik, A., Haynes, K. D., Baker, I., Luijkx, I. T., Groenink, M., Peters, W., Miller, J. B., Berry, J. A., Ogée, J., Meredith, L. K., Sun, W., Kohonen, K.-M., Vesala, T., Mammarella, I., Chen, H., Spielmann, F. M., Wohlfahrt, G., Berkelhammer, M., Whelan, M. E., Maseyk, K., Seibt, U., Commane, R., Wehr, R., and Krol, M.: Evaluation of carbonyl sulfide biosphere exchange in the Simple Biosphere Model (SiB4), Biogeosciences, 18, 6547–6565,, 2021. 

Krol, M. C., le Seger, P., and Segers, A. J.: TM5-4DVAR base model code, SourceForge [code], (last access: 21 November 2022), 2013. 

Kuai, L., Worden, J. R., Campbell, J. E., Kulawik, S. S., Li, K.-F. F., Lee, M., Weidner, R. J., Montzka, S. A., Moore, F. L., Berry, J. A., Baker, I., Denning, A. S., Bian, H., Bowman, K. W., Liu, J., and Yung, Y. L.: Estimate of carbonyl sulfide tropical oceanic surface fluxes using aura tropospheric emission spectrometer observations, J. Geophys. Res., 120, 11012–11023,, 2015. 

Lai, C. T., Katul, G., Oren, R., Ellsworth, D., and Schäfer, K.: Modeling CO2 and water vapor turbulent flux distributions within a forest canopy, J. Geophys. Res-Atmos., 105, 26333–26351,, 2000. 

Launois, T., Belviso, S., Bopp, L., Fichot, C. G., and Peylin, P.: A new model for the global biogeochemical cycle of carbonyl sulfide – Part 1: Assessment of direct marine emissions with an oceanic general circulation and biogeochemistry model, Atmos. Chem. Phys., 15, 2295–2312,, 2015. 

Lee, C. K., Daniel, R. M., Shepherd, C., Saul, D., Cary, S. C., Danson, M. J., Eisenthal, R., and Peterson, M. E.: Eurythermalism and the temperature dependence of enzyme activity, FASEB J. 21, 1934–1941,, 2007. 

Leuning, R., Dunin, F. X., and Wang, Y. -P.: A two-leaf model for canopy conductance, photosynthesis and partitioning of available energy. II. Comparison with measurements, Agr. Forest Meteorol., 91, 113–125,, 1998. 

Ma, J., Kooijmans, L. M. J., Cho, A., Montzka, S. A., Glatthor, N., Worden, J. R., Kuai, L., Atlas, E. L., and Krol, M. C.: Inverse modelling of carbonyl sulfide: implementation, evaluation and implications for the global budget, Atmos. Chem. Phys., 21, 3507–3529,, 2021. 

Maignan, F., Abadie, C., Remaud, M., Kooijmans, L. M. J., Kohonen, K.-M., Commane, R., Wehr, R., Campbell, J. E., Belviso, S., Montzka, S. A., Raoult, N., Seibt, U., Shiga, Y. P., Vuichard, N., Whelan, M. E., and Peylin, P.: Carbonyl sulfide: comparing a mechanistic representation of the vegetation uptake in a land surface model and the leaf relative uptake approach, Biogeosciences, 18, 2917–2955,, 2021. 

Maseyk, K., Berry, J. A., Billesbach, D., Campbell, J. E., Torn, M. S., Zahniser, M., and Seibt, U.: Sources and sinks of carbonyl sulfide in an agricultural field in the Southern Great Plains, P. Natl. Acad. Sci. USA, 111, 9064–9069,, 2014. 

Michaelis, L. and Menten, M. L.: Kinetik der Invertinwirkung, Biochem. Z., 49, 333–369, 1913. 

Miner, G. L., Bauerle, W. L., and Baldocchi, D. S.: Estimating the sensitivity of stomatal conductance to photosynthesis: a review, Plant Cell Environ., 40, 1214–1238,, 2017. 

Ogée, J., Sauze, J., Kesselmeier, J., Genty, B., Van Diest, H., Launois, T., and Wingate, L.: A new mechanistic framework to predict OCS fluxes from soils, Biogeosciences, 13, 2221–2240,, 2016. 

Papale, D., Reichstein, M., Aubinet, M., Canfora, E., Bernhofer, C., Kutsch, W., Longdoz, B., Rambal, S., Valentini, R., Vesala, T., and Yakir, D.: Towards a standardized processing of Net Ecosystem Exchange measured with eddy covariance technique: algorithms and uncertainty estimation, Biogeosciences, 3, 571–583,, 2006. 

Peterson, M. E., Eisenthal, R., Danson, M. J., Spence, A., and Daniel, R. M.: A new intrinsic thermal parameter for enzymes reveals true temperature optima, J. Biol. Chem., 279, 20717–20722,, 2004. 

Protoschill-Krebs, G., Wilhelm, C., and Kesselmeier, J.: Consumption of carbonyl sulphide (COS) by higher plant carbonic anhydrase (CA), Atmos. Environ., 30, 3151–3156,, 1996. 

Remaud, M., Chevallier, F., Maignan, F., Belviso, S., Berchet, A., Parouffe, A., Abadie, C., Bacour, C., Lennartz, S., and Peylin, P.: Plant gross primary production, plant respiration and carbonyl sulfide emissions over the globe inferred by atmospheric inverse modelling, Atmos. Chem. Phys., 22, 2525–2552,, 2022. 

Sato, N., Sellers, P. J., Randall, D. A., Schneider, E. K., Shukla, J., Kinter III, J. L., Hou, Y.-T., and Albertazzi, E.: Effects of implementing the Simple Biosphere Model in a general circulation model, J. Atmos. Sci., 46, 2757–2782,<2757:EOITSB>2.0.CO;2, 1989. 

Seibt, U., Kesselmeier, J., Sandoval-Soto, L., Kuhn, U., and Berry, J. A.: A kinetic analysis of leaf uptake of COS and its relation to transpiration, photosynthesis and carbon isotope fractionation, Biogeosciences, 7, 333–341,, 2010. 

Sellers, P. J., Mintz, Y., Sud, Y. C., and Dalcher, A.: A Simple Biosphere Model (SiB) for use within general circulation models. J. Atmos. Sci, 43, 505–531,<0505:ASBMFU>2.0.CO;2, 1986. 

Sellers, P. J., Berry, J. A., Collatz, G. J., Field, C. B., and Hall, F. G.: Canopy reflectance, photosynthesis, and transpiration. III. A reanalysis using improved leaf models and a new canopy integration scheme, Remote Sens. Environ., 42, 187–216,, 1992. 

Sellers, P. J., Randall, D. A., Collatz, G. J., Berry, J. A., Field, C. B., Dazlich, D. A., Zhang, C., Collelo, G. D., and Bounoua, L.: A revised land surface parameterization (SiB2) for atmospheric GCMs. Part I: Model Formulation, J. Climate, 9, 676–705,<0676:ARLSPF>2.0.CO;2, 1996a. 

Sellers, P. J., Tucker, C. J., and Collatz, G. J.: A revised land surface parameterization (SiB2) for atmospheric GCMs. Part II: The generation of global fields of terrestrial biophysical parameters from satellite data, J. Climate, 9, 706–737,<0706:ARLSPF>2.0.CO;2, 1996b. 

SiB4 project members: Simple Biosphere Model Version 4.2, Git-lab [code], (last access: 21 November 2022), 2020. 

Smith, N. E., Kooijmans, L. M. J., Koren, G., van Schaik, E., van der Woude, A. M., Wanders, N., Ramonet, M., XuerefRemy, I., Siebicke, L., Manca, G., Brümmer, C., Baker, I. T., Haynes, K. D., Luijkx, I. T., and Peters, W.: Spring enhancement and summer reduction in carbon uptake during the 2018 drought in northwestern Europe, Philos. T. Roy. Soc. B, 375, 20190509,, 2020. 

Stimler, K., Montzka, S. A., Berry, J. A., Rudich, Y., and Yakir, D.: Relationships between carbonyl sulfide (COS) and CO2 during leaf gas exchange, New Phytol., 186, 869–878,, 2010. 

Stimler, K., Berry, J. A., Montzka, S. A., and Yakir, D.: Association between carbonyl sulfide uptake and 18delta during gas exchange in C3 and C4 leaves, Plant Physiol., 157, 509–517,, 2011. 

Stimler, K., Berry, J. A., and Yakir, D.: Effects of carbonyl sulfide and carbonic anhydrase on stomatal conductance, Plant Physiol., 158, 524–530,, 2012. 

Sun, W., Maseyk, K., Lett, C., and Seibt, U.: A soil diffusion–reaction model for surface COS flux: COSSM v1, Geosci. Model Dev., 8, 3055–3070,, 2015. 

Sun, W., Kooijmans, L. M. J., Maseyk, K., Chen, H., Mammarella, I., Vesala, T., Levula, J., Keskinen, H., and Seibt, U.: Soil fluxes of carbonyl sulfide (COS), carbon monoxide, and carbon dioxide in a boreal forest in southern Finland, Atmos. Chem. Phys., 18, 1363–1378,, 2018. 

Tarantola, A. and Valette, B.: Inverse problems = Quest for information, J. Geogr., 50, 159–170, 1982. 

Vesala, T., Kohonen, K.-M., Kooijmans, L. M. J., Praplan, A. P., Foltýnová, L., Kolari, P., Kulmala, M., Bäck, J., Nelson, D., Yakir, D., Zahniser, M., and Mammarella, I.: Long-term fluxes of carbonyl sulfide and their seasonality and interannual variability in a boreal forest, Atmos. Chem. Phys., 22, 2569–2584,, 2022. 

Walker, A. P., Beckerman, A. P., Gu, L., Kattge, J., Cernusak, L. A., Domingues, T. F., and Woodward, F. I.: The relationship of leaf photosynthetic traits – Vcmax and Jmax – to leaf nitrogen, leaf phosphorus, and specific leaf area: A meta-analysis and modeling study, Ecol. Evol., 4, 3218–3235,, 2014. 

Wehr, R. and Saleska, S. R.: An improved isotopic method for partitioning net ecosystem–atmosphere CO2 exchange, Agr. Forest Meteorol., 214–215, 515–531,, 2015. 

Wehr, R. and Saleska, S. R.: Calculating canopy stomatal conductance from eddy covariance measurements, in light of the energy budget closure problem, Biogeosciences, 18, 13–24,, 2021. 

Wehr, R., Munger, J. W., McManus, J. B., Nelson, D. D., Zahniser, M. S., Davidson, E. A., Wofsy, S. C., and Saleska, S. R.: Seasonality of temperate forest photosynthesis and daytime respiration, Nature, 534, 680–683,, 2016. 

Wehr, R., Commane, R., Munger, J. W., McManus, J. B., Nelson, D. D., Zahniser, M. S., Saleska, S. R., and Wofsy, S. C.: Dynamics of canopy stomatal conductance, transpiration, and evaporation in a temperate deciduous forest, validated by carbonyl sulfide uptake, Biogeosciences, 14, 389–401,, 2017. 

Whelan, M. E., Hilton, T. W., Berry, J. A., Berkelhammer, M., Desai, A. R., and Campbell, J. E.: Carbonyl sulfide exchange in soils for better estimates of ecosystem carbon uptake, Atmos. Chem. Phys., 16, 3711–3726,, 2016. 

Whelan, M. E., Lennartz, S. T., Gimeno, T. E., Wehr, R., Wohlfahrt, G., Wang, Y., Kooijmans, L. M. J., Hilton, T. W., Belviso, S., Peylin, P., Commane, R., Sun, W., Chen, H., Kuai, L., Mammarella, I., Maseyk, K., Berkelhammer, M., Li, K.-F., Yakir, D., Zumkehr, A., Katayama, Y., Ogée, J., Spielmann, F. M., Kitz, F., Rastogi, B., Kesselmeier, J., Marshall, J., Erkkilä, K.-M., Wingate, L., Meredith, L. K., He, W., Bunk, R., Launois, T., Vesala, T., Schmidt, J. A., Fichot, C. G., Seibt, U., Saleska, S., Saltzman, E. S., Montzka, S. A., Berry, J. A., and Campbell, J. E.: Reviews and syntheses: Carbonyl sulfide as a multi-scale tracer for carbon and water cycles, Biogeosciences, 15, 3625–3657,, 2018.  

White, M. L., Zhou, Y., Russo, R. S., Mao, H., Talbot, R., Varner, R. K., and Sive, B. C.: Carbonyl sulfide exchange in a temperate loblolly pine forest grown under ambient and elevated CO2, Atmos. Chem. Phys., 10, 547–561,, 2010. 

Wolf, A., Akshalov, K., Saliendra, N., Johnson, D. A., and Laca, E. A.: Inverse estimation of Vcmax, leaf area index, and the Ball–Berry parameter from carbon and energy fluxes, J. Geophys. Res-Atmos., 111, D08S08,, 2006. 

Woodward, F. I., Smith, T. M., and Emanuel, W. R.: A global land primary productivity and phytogeography model, Global Biogeochem. Cy., 9, 471–490,, 1995. 

Short summary
Carbonyl sulfide (COS) is a useful constraint for estimating photosynthesis. To simulate COS leaf flux better in the SiB4 model, we propose a novel temperature function for enzyme carbonic anhydrase (CA) activity and optimize conductances using observations. The optimal activity of CA occurs below 40 °C, and Ball–Woodrow–Berry parameters are slightly changed. These reduce/increase uptakes in the tropics/higher latitudes and contribute to resolving discrepancies in the COS global budget.
Final-revised paper