Optimizing seasonally variable photosynthetic parameters based on joint carbon and water �ux constraints

Abstract


Introduction
Terrestrial ecosystems account for approximately 36% of global CO2 uptake via plant photosynthesis (Friedlingstein et al., 2020) and return about 60% of terrestrial precipitation to the atmosphere through evapotranspiration (Jasechko et al., 2013;Oki and Kanae, 2006).
Precise modeling of gross primary productivity (GPP) and evapotranspiration (ET) in terrestrial ecosystems is important in quantifying the intricate dynamics of global terrestrial carbon and water cycles under climate change (Jung et al., 2010;Peters et al., 2018).State-ofthe-art terrestrial biosphere models (TBMs), that integrate biogeochemical and biogeophysical processes in the soil-vegetation-atmosphere continuum, have been developed to simulate terrestrial water and carbon fluxes and to quantify their roles in global carbon and water cycles (Cao and Woodward, 1998;Sitch et al., 2003).Within the essential biophysical components involved in terrestrial ecosystem processes, plant stomata regulate the exchange of carbon and water between ecosystems and the atmosphere (Berry et al., 2010;Lin et al., 2015;Wolz et al., 2017).Hence, gaining a better understanding of how to model stomatal conductance (  ) is important for precise quantification of GPP and ET (Berry et al., 2010).
A common approach to estimating carbon and water fluxes is to integrate an enzyme-kinetic biochemical photosynthesis model developed by Farquhar et al. (1980) with a semi-empirical stomatal conductance model, known as Ball-Woodrow-Berry (BWB) (Ball et al., 1987).The Farquhar model quantifies how leaf photosynthesis responds to environmental variables such as incoming radiation, air temperature, CO2 concentration and to leaf physiological conditions such as leaf nitrogen content.The pivotal parameter in the model is the maximum carboxylation rate normalized to 25 °C (  25 ), representing the leaf photosynthetic capacity.Other model parameters, such as the maximum rate of electron transport (  ) and dark respiration (  ), can be inferred from   25 (Medlyn et al., 1999;Wullschleger, 1993).The BWB model describes stomata conductance   as: where  is the photosynthetic rate (mol m -2 s -1 ),  is the relative humidity at the leaf surface (%),   is the CO2 concentration at the leaf surface (µmol mol -1 ),  is the slope of linear regression between   and  • /  (unitless), and  0 is the intercept of the regression representing minimum   .The linkage between   and  in the BWB model shows the impact of  on   estimation depends on the environmental conditions and model parameters (i.e.

𝑉 𝑐𝑚𝑎𝑥 25
) that influence . serves as a key parameter in determining the rate of carbon and water fluxes between plants and the atmosphere.Therefore, in the photosynthesis-stomata coupling schemes,  and   25 are crucial parameters to depict carbon and water cycles within terrestrial ecosystems (Ainsworth and Rogers, 2007;Bauerle et al., 2014).
The prescribed values of  and   25 in TBMs can significantly affect the simulated transpiration (Barnard and Bauerle, 2013;Lai et al., 2000;Leuning et al., 1998) and GPP (Bauerle et al., 2014;Luo et al., 2001). and   25 are often set as constant values for each biome in TBMs.For example, the Community Land Model (CLM) version 4.0 assigned  = 6,

25
= 52 µmol m - 2 s -1 for boreal deciduous broadleaf forests (Oleson et al., 2010).However,  and   25 have been found to vary seasonally with environmental conditions (Miner and Bauerle, 2017;Misson et al., 2004), and the uncertainties in model parameterization can be comparable to model structure errors (Bonan et al., 2011;Cai et al., 2019;Ricciuto et al., 2018)  25 with only leaf-level measurements, which cannot be extended for other species or regional analysis.At the ecosystem level,   25 can be estimated from assimilation of suninduced chlorophyll fluorescence (SIF) data (He et al., 2017) and leaf chlorophyll content (LCC) (Croft et al., 2020), and the estimations of   within these two approaches have been shown to have good accuracy compared to ground-based data (Chen et al., 2022b).To estimate  at the ecosystem level, the model inversion method is often utilized by forcing models to regenerate flux-derived GPP (Ju et al., 2010;Ono et al., 2013), latent heat fluxes (Wolf et al., 2006;Xu and Baldocchi, 2003), or measured sap flow (Lai et al., 2000;Sala and Tenhunen, 1996).However,  estimates were only indirectly validated using GPP and ET simulations (Ju et al., 2010).It is unclear whether the estimated  from model inversions are accurate, without direct comparisons with  upscaled from gas exchange measurements.Furthermore, in the model inversions for  estimation,   25 was assumed to be constant throughout the growing season (Dang et al., 1998;Lai et al., 2000;Wolf et al., 2006), which may result in errors in GPP and ET simulations (Chen et al., 2022b;Chen et al., 2019).with gas exchange measurements.

As pivotal parameters in
2) evaluate the improvements in modeling GPP and ET with seasonally varied  and   25 .
3) investigate the sensitivities of  and   25 on modelled GPP and ET estimates.

Site profile
The Borden Forest Research Station is a mixed forest site located in southern Ontario (44°19' N, 79°56' W) near the southern end of Georgian Bay within the Great Lakes/St.Lawrence Forest region (Froelich et al., 2015).It lies in an ecotone with southern temperate forest species and northern boreal forest species, which extends across eastern North America between 44° and 47°N.The mean annual temperature at this site is 7.4 °C and mean annual total precipitation is 784 mm (Froelich et al., 2015).The mean canopy height is 22 m and the dominant species include red maple (Acer rubrum, 52.2% stem count), eastern white pine (Pinus strobes, 13.5% stem count), large-tooth aspen (Populus grandidentata, 7.7% stem count), trembling aspen (Populus tremuloides, 2.5% stem count), and white ash (Fraxinus americana, 7.1% stem count) (Lee et al., 1999;Teklemariam et al., 2009).The fetch of largely uninterrupted forest extends to distances of 1.5-4 km in southeastern and southwestern quadrants, and to 1 km in the northeastern direction (Froelich et al., 2015).The northwest quadrant contains a white pine plantation and was excluded from the footprint of eddy-covariance calculations at this site (Luo et al., 2018b).The soil type at this site is sandy loam (Gonsamo et al., 2015).

Field measurements
Leaf gas exchange measurements were conducted for four dominant C3 broadleaf species (red maple, large-tooth aspen, trembling aspen, and white ash) every 7-15 days during the growing seasons of 2014-2018, using a LI-6400XT portable photosynthesis system (LI-COR Inc., Lincoln, NE USA).Leaves from each species were selected from top-of-canopy branches directly accessed from the 44 m flux tower.The CO2 and   response curves were measured under an artificial saturated light source at photosynthetic photon flux density (PPFD) levels of 1800 µmol m -2 s -1 , and stepwise ambient CO2 in the chamber of 400,200,100,50,400,400,600,800,1000,1200,1500, 1800 µmol mol -1 .Prior to logging response curves, the system was warmed up and leaves acclimated in the chamber at PPFD of 1800 µmol m -2 s -1 , ambient relative humidity, a temperature of ~25 °C, and a CO2 concentration at 400 µmol mol -1 .The photosynthetic model parameter,   25 was calculated from the CO2 response curves using a curve-fitting tool developed by Kevin Tu (www.landflux.org),following Ethier and Livingston (2004) and then normalized to 25 °C using the Arrhenius equation for comparability amongst parameters in terrestrial biophysical models (Sharkey et al., 2007).The stomatal conductance model parameter, the Ball-Berry slope  , was fitted using ordinary least square regression between stomatal conductance   and the Ball index,  • /  recorded in gas exchange experiments.In the regression, only data when   >360 µmol mol -1 were used following Xu and Baldocchi (2003) and only  with intercept  0 <0.05 mol m -2 s -1 and Pearson's R 2 >0.80 were selected as measured  values (Miner et al., 2017).The

25
and  for each species on each field measurement day were upscaled to an ecosystem level by multiplying the relative proportion among the four species and then aggregated to monthly parameters.
Leaf area index (LAI) and canopy structural parameters were collected on the same days as leaf gas exchange experiments, along a 100 m line transect extending from the base of the tower in a direction from the north to the south.Effective LAI (  ) values were measured using the LAI-2000 plant canopy analyzer (LI-COR, Lincoln, NE, USA).The true LAI values were then calculated from   following the methods by Chen et al. (1997): where  is the ratio of woody area to total area (set as  = 0.17 accounting for the interception of radiation by branches and tree trunks that results in artificially high LAI values),   is the ratio of needle area to shoot area (set at 1 considering individual leaves of broadleaf species as foliage elements), and Ω  is the element clumping index, measured as 0.95 using the TRAC (Tracing Radiation and Architecture of Canopies) instrument (Chen and Cihlar, 1995).

Flux and meteorological measurements
The ecosystem fluxes were obtained based on measurements of turbulent fluctuations of wind velocity, temperature, water vapor and CO2, using a sonic anemometer (K-Type, Applied Technologies Inc., CO, USA) located at a height of 33 m coupled with a closed-path infrared gas analyzer (IRGA, model LI-6262, LI-COR Inc., Lincoln, NE, USA).Both the IGRA and the anemometer are operated at 10Hz and the IGRA is placed in a temperature-controlled hut at the base of the tower.The high frequency eddy covariance (EC) fluxes were post-processed and aggregated into half-hourly fluxes based on the methods in Froelich et al. (2015).
The net ecosystem exchange (NEE, µmol m -2 s -1 ) is calculated as the sum of   2 , the vertical CO2 fluxes, and   2 , the change of CO2 storage in the air below the sensors at 33 m.  2 (µmol m -2 s -1 ) was calculated from the WPL (Webb-Pearman-Leuning) corrected vertical turbulent flux of the CO2 mole fraction as , where ′ is the vertical wind velocity,  2 ′ is the atmospheric CO2 mole fraction,   is the density of air,   is the molecular weight of dry air, and the WPL eliminates the density effects due to water vapor fluctuations (Webb et al., 1980).Other terms such as horizontal turbulent fluxes, vertical and horizontal advections were assumed to be negligible (Staebler and Fitzjarrald, 2004).The CO2 storage change   2 was estimated as the vertical integral of change of CO2 concentration with time () as , where  is the height of sensors at 33 m.NEE data were filtered out when wind direction was >285° or <20° due to the unrepresentativeness of the vegetation in that direction.A change point detection method by Barr et al. (2013) was used to determine the friction velocity threshold ( * ) to exclude eddy covariance measurements with poor turbulence for computed fluxes.Gaps in the NEE data were then filled using the model of Barr et al. (2004), which partitions NEE between gross ecosystem productivity (GPP) and ecosystem respiration (RE) as NEE = GPP -RE.RE was assumed to be equal to NEE during nights or cold periods while during warm seasons or in data gaps, RE was empirically modeled based on air and soil temperature (Froelich et al., 2015).Half-hourly GPP was calculated as the sum of NEE and RE and was gap-filled using a Michaelis-Menten type model based on PAR (Barr et al., 2004;Froelich et al., 2015).
In parallel with the CO2 flux measurements, the latent heat flux (LE, W m -2 ) is also measured , where  is the latent heat of vaporization,   is the density of air,  ′  2  ′ ̅̅̅̅̅̅̅̅̅̅  is the WPL-corrected vertical turbulent flux of the water vapor mole fraction, and   ,   are water, air molar masses, respectively (Froelich et al., 2015;Teklemariam et al., 2009).
Several half-hourly auxiliary meteorological variables were also obtained at the site to initialize and drive the TBM.Incoming shortwave radiation (W m -2 ) on top of the canopy was measured at the top of the flux tower.Air temperature (°C) and relative humidity (%) were measured at two heights (33 and 41 m), but the 33 m data were used as model inputs.Wind speed (m s -1 ) and wind direction were measured by the sonic anemometer mounted at 44 m.Soil temperature (°C) and soil moisture (m 3 m -3 ) at the depth of 5-30 cm were selected for model initialization.
Precipitation data (mm h -1 ) were obtained from Environment Canada weather station at Egbert, which is the data source for meteorological data gap-filling at the site (Froelich et al., 2015).
All the flux and meteorological data during the growing season in 2014-2018 were used for parameter optimization, except for August in 2016 due to some problems with condensation in the sample tubing during that period.

Terrestrial ecosystem biophysical model
The Boreal Ecosystem Productivity Simulator (BEPS) is a two-leaf diagnostic enzyme-kinetic model originally developed for regional carbon uptake and the water cycle for boreal ecosystems (Chen et al., 1999;Liu et al., 2003;Liu et al., 1997).It has been revised and upgraded since the original daily version to simulate carbon and water cycles at an hourly time step at sites over various PFTs (Ju et al., 2006;Luo et al., 2018b;Luo et al., 2019), regionally (Gonsamo et al., 2013;Liu et al., 2003), and globally (Chen et al., 2019).Previous inter-model comparisons and site-level validations have substantiated the reliability of simulated carbon and water fluxes from BEPS (Amthor et al., 2001;Luo et al., 2018b;Potter et al., 2001), and also shown the validity of biophysical model hyperparameters inferred from BEPS (Chen et al., 2022b;He et al., 2019).The modules related to GPP and ET simulations in BEPS are described in the following subsections.

Leaf-level photosynthesis and transpiration
BEPS adopts the Farquhar biochemical model (Farquhar et al., 1980), calculating the instantaneous leaf photosynthetic rate as: where   is the photosynthetic rate (µmol m -2 s -1 ) under RuBP carboxylase/oxygenase limited conditions,   is the photosynthetic rate under electron transport limited conditions in the RuBP regeneration process, and   is the dark respiration rate.
where   is the maximum carboxylation rate (µmol m -2 s -1 ) calculated from preset   25 and a temperature dependent function (Sharkey et al., 2007),  is the electron transport rate (µmol m -2 s -1 ),   is the intercellular CO2 concentration (µmol mol -1 ), Γ is the CO2 compensation point without dark respiration (µmol mol -1 ), and  is a Rubisco enzyme kinetics function as  =   /(1 +   /  ).  ,   are Michaelis-Menten constants for CO2 (µmol mol -  =   ×   + 2.1  (7.) and   is obtained using the equation by Medlyn et al. (1999): BEPS follows the method by Baldocchi (1994) to solve  and   from the simultaneous equations in the Farquhar model and the BWB model until the energy balance converges in each time step.The leaf-level transpiration is then calculated using the Penman-Monteith equation in BEPS: where   is the instantaneous net radiation on the leaf surface (W m -2 ),  is the heat storage of the leaf,  is the density of air (kg m -3 ),   is the specific heat of air (J kg -1 °C-1 ),  is the vapor pressure deficit on the leaf surface (kPa), Δ is the rate of change of saturation specific humidity with air temperature (kPa °C-1 ),  is the latent heat of water (J kg -1 ),  is the psychrometric constant (kPa °C-1 ), and   is the leaf boundary layer resistance to water vapor (m s -1 ).

Two-leaf scheme upscaling from leaf to canopy
Based on the first order of the incoming solar radiation on leaves, the canopy can be separated into sunlit and shaded leaves groups (Chen et al., 1999;Norman, 1982), and such two-leaf scheme has been proven to be capable and reliable for GPP and ET simulations than the bigleaf approach (Chen et al., 1999;Luo et al., 2018a).Following the two-leaf scheme, leaves are separated into sunlit and shaded leaves groups based on incoming solar radiation orientation and canopy geometry.Sunlit leaves receive both direct and diffuse radiation while shaded leaves receive only diffuse radiation.Sunlit and shaded LAI (  and  ℎ ) are calculated in BEPS following the method described in Chen et al. (1999): where   and   are the instantaneous photosynthetic rate and the transpiration from a sunlit leaf,  ℎ and  ℎ are the instantaneous photosynthetic rate and the transpiration from a shaded leaf.
In addition to transpiration, BEPS also simulates evaporation from soil (  ) and from wet leaf surface (  ): , can be estimated.
To explore how GPP and ET affect the retrieval of  and   25 , three scenarios were tested: 1) Scenario 1 (carbon): only the photosynthesis constraint is considered in the parameter optimization, and the cost function is the modified mean squared error (MMSE) between modeled GPP from BEPS (  ) and observed GPP from fluxes (  ): 2) Scenario 2 (water): only the water constraint is considered in the parameter optimization, and the cost function is the MMSE between modeled ET from BEPS (  ) and observed LE from fluxes (  ): 3) Scenario 3 (carbon-water coupling): both photosynthetic and water constraints are considered in the parameter optimization, and the cost function is the sum of normalized  1 and  2 to eliminate the units of GPP and ET: where  is the total number of hourly simulations in a parameter optimization interval,      is the ordinal least square regression (OLS) slope between   and   ,      is the OLS slope between   and   .The MMSE in cost functions combines the mean squared error and the deviation of OLS slopes between modeled and observed variables from the identity line.By minimizing MMSE, modeled variables and observed variables can better fit the identity line, which yields the two optimal parameters,  and   25 , that can generate a model of higher accuracy.
In this study,  and   25 were estimated monthly.In the optimization process,  was set in range [1, 15],   was set in range [0, 80], the minimum stomatal conductance in the BWB model  0 was set to be smaller than 0.02, because forcing  0 to zero can impact estimates of  and induce errors (Falge et al., 1996;Miner et al., 2017;Van Wijk et al., 2000).To eliminate the stochastic fluctuations of estimated  and   25 from the global optima,  and   25 were the mean values from ten optimizations.In each round of optimizations, BEPS was run for 500 times (Snoek et al., 2012).The  values optimized in the three scenarios show great discrepancies. optimized from the carbon-water coupling scenario (Figure 1e) show the highest correlation with measured  (R 2 =0.70), while R 2 for  in the carbon (Figure 1a) and water scenarios (Figure 1c) are 0.00 and 0.52, respectively.These results indicate that by constraining ET only in Scenario 2, the seasonal variations and anomalies were optimized.Nonetheless, the magnitude of  was underestimated when only ET is constrained, as demonstrated by a regression that falls below the identity line (the slope of regression between optimized and measured  is 0.75).When only GPP is constrained in Scenario 1,  cannot be accurately optimized.However, the incorporation of both GPP and ET constraints in  show the better correlation with measured GPP (R 2 =0.94 and R 2 =0.85 for hourly and daily results, respectively), smaller derivation from the identical line of measured GPP (y=0.95x+0.020g m -2 h -1 and y=0.85x+1.466g m -2 d -1 for hourly and daily results, respectively), and lower RMSE of 0.137 g m -2 h -1 and 2.038 g m may not precisely simulate the stomata functions during the short periods at the beginning, peak, and end of the growing season.The tendency to overestimatee fluxes when they are low and underestimate them when they are high have been a common issue for TBMs (Schaefer et al., 2012;Winkler et al., 2019).Incorporating seasonally variable  and   25 may help address this issue.Figure 4a and Figure 4c show that GPP is not sensitive to  increment (20~75 g m -2 year -1 to +50% ) but sensitive to  decrement (-70~-300 g m -2 year -1 to -50% ), and the sensitivity depends on the   25 level.This indicates that though  increments result in higher   , the photosynthesis is limited by photosynthetic capacity   25 .Figure 4b and Figure 4d show that ET is sensitive to  and the sensitivity is influenced by   25 .ET is sensitive to  increment (100~160 mm year -1 to +50%  ) or decrement (-130~-220 mm year -1 to -50%  ).The sensitivity is high when   25 are at the normal and high levels but the sensitivity is low when is at the low level.This may result from the joint control of  and   25 on   estimates.
When   25 is low, the photosynthetic rate is limited by   25 , thus limiting the modeled stomata aperture.Figure 5a and Figure 5c show that GPP is sensitive to   25 (350~500 g m -2 year -1 to +50%   25 ) but the sensitivity decreases when  is low.This is likely because low  limits the change of modeled   to

25
, which oppositely affects the modeled photosynthetic rate.Figure 5b and Figure 5d show that ET is sensitive to   25 (20~70 mm year -1 to +50%   25 ). Figure 4 and Figure 5 show the sensitivities of GPP and ET to  and   25 .GPP is more sensitive to   25 than to , and ET is more sensitive to  than to   25 .
However, both GPP and ET are still influenced by both  and   25 , suggesting that regarding The comparison of the three optimization scenarios demonstrates that using GPP and ET as the constraints simultaneously is the best method to estimate  and   25 .By incorporating variable  and   25 into BEPS, the temporal patterns of simulated GPP and ET are greatly improved and the biases are considerably reduced.Figure 6 shows the differences of modeled  Figure 7   represents the photosynthetic capacity normalized to 25 °C and has been found to be linearly correlated with leaf chlorophyll content (Croft et al., 2017;Luo et al., 2018b).The nitrogen content in leaves determine the enzyme amount for the carboxylation processes and the enzymatic activities are affected by air temperature (Chen et al., 2022b;Smith et al., 2019).The correlations between   at this forest site.However, the four sample plant species that were reachable from the tower did not account for all species in the footprint of the EC tower, so there could be some uncertainties of upscaled measurements of  and   25 .
Second, the monthly  and   25 were optimized using GPP and ET with measured site-level EC fluxes, measured LAI, measured site-level meteorological forcing data.The Borden site lies within in an ecotone that contains both temperate and boreal species, which can to some extent represent the forest ecosystems in the temperate and boreal regions.Our method achieved good performance for boreal mixed forests.However, further research is needed to extend the framework in this study to other ecosystems and validate the accuracy in those settings.
Third, as the two pivotal parameters in photosynthesis and stomatal conductance models,  and   25 can be physiologically affected by the environment.In this study,  was found to be well correlated with the product of shortwave radiation and soil water content, while   25 was found to be well correlated with air temperature.However, there could be more environmental factors that should be included in the correlation and it is possible that the number of optimized  and   25 datapoints may not be sufficient to generate clear correlations.Further research will emphasize on exploration of additional environmental factors to improve the applicability of the method in this study.
Understanding how  and   25 change over space and time is crucial to guide their dynamic parameterization within TBMs.We appeal for retrieval of seasonally variable  and   25 , which can lead to accurate simulations of GPP and ET time series by TBMs, to better understand the GPP-ET decoupling trends (Cheng et al., 2017;Keenan et al., 2013) and to quantify GPP and ET under climate change (Chen et al., 2022a;Keenan et al., 2021).

Conclusion
Coupling photosynthesis and stomatal conductance models is a widely used method for simulating carbon and water cycles.
TBMs,  or   25 are often prescribed as constant values.In this study, we explore the seasonal variations of  and   25 and assess how seasonally varied validated  and   25 improve the accuracy of GPP and ET simulations.A TBM is utilized to optimize  and   25 at a mixed forest eddy-covariance flux tower site, and to conduct GPP and ET modeling with optimized  and   25 .The aims of this study are to: 1) propose a parameter optimization framework to estimate seasonally varied  and   25 at the ecosystem level and validate the optimized  and   25 photosynthesis, thus not likely to change much in  and   25 optimizations.

Figure 1 .
Figure 1.Scatter plots of optimized  and   25 compared with measured  and   in the growing seasons of 2014-2018 at the Borden site.(a) and (b) refer to the carbon scenario (Scenario 1), considering only the GPP constraint.(c) and (d) refer to the water scenario

Figure 1
Figure 1 shows the comparisons of measured  and   25 to optimized  and   25 in three scenarios: (1) carbon flux only, (2) water flux only, and (3) coupled carbon and water fluxes.
Figure 2. Scatter plots of the predicted and measured GPP of the growing seasons in 2014-2018 at the Borden site.(a) and (b) refer to the hourly GPP validation with fixed  and   25 and optimized  and   25 with seasonal variation, respectively; (c) and (d) refer to the daily GPP validation with fixed  and   25

Figure 3 .
Figure 3. Scatter plots of the predicted and measured ET of the growing seasons in 2014-2018 at the Borden site.(a) and (b) refer to the hourly ET validation with fixed  and   25 and optimized  and   25 with seasonal variation, respectively; (c) and (d) refer to the daily ET validation with fixed  and   25 , and optimized  and   25 with seasonal variation,

3. 3
Impacts of  and   25 on estimated GPP and ET According to the modeling results, optimized and seasonally variable  and   25 markedly improved the estimation of GPP and ET during the growing season in 2014-2018 at the Borden site.To understand how optimized  and   25 affects GPP and ET estimates, we explored the sensitivities of GPP and ET to changing  and   25 .To access the sensitivities of GPP and ET, each of  and   was changed from -50% to 50% by 2% steps to drive BEPS, while other variables were unchanged.To further reveal the interaction between  and   on simulations of GPP and ET, in the sensitivity test of ,   25 was prescribed as low (  25 = 30 µmol m -2 s -1 ), normal (  25 = 60 µmol m -2 s -1 ), and high values (  25 = 90 µmol m -2 s -1 ), while in the sensitivity test of   25 ,  was also prescribed as low ( = 4), normal ( = 8), and high values ( = 12).

Figure 4 .
Figure 4. Sensitivities of ET and GPP to  at the Borden site.(a) GPP by percentage (%); (b) ET by percentage (%); (c) annual GPP (g m -2 year -1 ); (d) annual ET (mm year -1 ).Red, green, and blue lines refer to the sensitivities of ET and GPP to  when   25 values are low (  25 =

Figure 5 .
Figure 5. Sensitivities of ET and GPP to   25 at the Borden site.(a) GPP by percentage (%); (b) ET by percentage (%); (c) annual GPP (g m -2 year -1 ); (d) annual ET (mm year -1 ).Red, green, and blue lines refer to the sensitivities of ET and GPP to   25 when  values are low ( = Figure 6.The 5-year average differences of modeled stomatal conductance (  ) simulated by BEPS using optimized  and   25 and prescribed  and   25 over the growing seasons of 2014-2018.The blue semitransparent dots refer to the modeled hourly   differences.The red lines refer to the 3-day average of the modeled hourly   differences.(a) sunlit leaves; (b) shaded leaves.

Figure 7 .
Figure 7. Sensitivities of sunlit and shaded   to  and   at the Borden site.(a) sunlit   by percentage (%) to ; (b) shaded   by percentage (%) to ; (c) sunlit   by percentage (%) to   25 ; (d) sunlit   by percentage (%) to   25 .Red, green, and blue lines in (a) and (b) refer to the sensitivities of sunlit and shaded   to  when   25 values are low (  25 = 30 may be physiologically determined by environmental conditions which can influence photosynthetic activities.To determine how  and   25 are affected by the environment,  and   25 were correlated with five environmental variables: air temperature, shortwave radiation, soil water content, relative humidity, and CO2 concentration.All environmental variables were normalized.Then  and   25 were also correlated with the product of each two environmental variables among the five, as shown in Figure 8 and Figure 9.The product of monthly mean shortwave radiation (SW) and soil water content (SWC) shows the best correlation with optimized  (R = 0.53, P < 0.01).This is likely because  is influenced by water availability and environmental conditions that affects the photosynthesis.Air temperature (TA) shows the best correlation with optimized   25 (R = 0.73, P < 0.01).This is likely because   25

𝑚
and SW n * SWC n , and between   25 and TA n imply the possibilities for determining  and   25 directly from environmental variables thus to better quantify GPP and ET over large scales.
Among the model parameters,  and   25 are the most important two variables that can significantly affect the GPP and ET estimates through their joint influences on modelled stomatal behavior.However, state-of-the-art TBMs consider  as a constant even though it has been found to vary across a growing season.To address the lack of research focusing on the retrieval of  at the canopy level, in this study, we designed a parameter optimization framework for monthly  and   25 retrievals and validated their accuracy with measurements at a boreal forest site.The following conclusions are drawn: 1.  and   25 optimized with GPP and ET as joint constraints show good accuracy with measured  and   25 at a boreal forest site (R 2 =0.70 for both  and   25) .Both  and   25 show great seasonal variations.2. The two-leaf TBM with optimized and seasonally variable  and   25 is capable of more accurately simulating GPP and ET than using traditional fixed  and   25 values.The biases of GPP and ET are considerably reduced for both daily and hourly simulations and the temporal correlations between modeled and measured fluxes are considerably improved.By incorporating variable  and   25 values, the RMSE for daily GPP simulation is reduced from 2.579 to 2.088 g m -2 d -1 and RMSE for daily ET simulation is reduced from 1.151 to 1.137 mm d -1 , with R 2 increasing from 0.76 to 0.85 for GPP and from 0.65 to 0.71 for ET.This study, for the first time, estimates seasonally variable  and   25 at the ecosystem level across the growing season, and presents a unique comparison with a large temporal dataset of measured  and   25 values.The findings demonstrate that with variable  and   25 parameterization, a TBM can generate carbon and water flux estimates for forest ecosystems with improved accuracy, particularly for low and high GPP and ET values.This work points to the need to incorporate dynamic, seasonally variable  and   25 to more accurately quantify the dynamic nature of the global carbon and water cycles under climate change with TBMs.
-2 d -1 for hourly and daily results, respectively.With optimized variable  and   25 , ET modeled by BEPS show the better correlation with measured ET (R 2 =0.81 and R 2 =0.71 for hourly and daily results, respectively), smaller derivation from the identical line of measured ET (y=0.92x+0.021mmh -1 and y=0.95x+0.407mmd - for hourly and daily results, respectively), and lower RMSE of 0.093 mm h -1 and 1.137mm d -1 for hourly and daily results, respectively.Although BEPS with prescribed constant  and   25 generates satisfactory GPP and ET simulations for intermediate values of GPP and ET, small GPP and ET values tend to be overestimated while large GPP and ET values tend to be underestimated.This is likely because constant  and   25 of sunlit (   ) and shaded leaves (  ℎ ) by BEPS using optimized  and   25 and prescribed  and   25 .The discrepancies between the two modeled   mainly result from the   estimates of sunlit leaves.   is more affected by  and   25 and there are only slight differences in modeled   ℎ with variable  and   25 .At the beginning and ending of the growing season,    modeled with constant  and   25 are higher than    modeled with variable  and   25 , which results in the underestimation of GPP and ET during those periods.At the peak of the growing season,    modeled with constant  and   25 are lower than    modeled with variable  and   25 , which results in the overestimation of GPP and ET.This implies that the enhancement in modeling   of sunlit leaves contributes to the improvements in GPP and ET simulations over the full growing season.To explore how  and   25 affect the modeling of    and   ℎ , the same sensitivity test as Section 3.3 was conducted for    and   ℎ .
shows that    and   ℎ are sensitive to the changes in .   is the most sensitive to  but the sensitivity low (+15%    to +20%  ) when      to +20% ).ℎ is also sensitive to changes in  and the sensitivity is smaller at lower   25 .This likely results from the model linkages between   and , indicating that the sensitivity of   to  is also dependent on the influence of   on  .Although the sensitivities of    and   ℎ to   25 are lower than those to , the patterns of    and   ℎ variations to changing   25 are similar to the pattern of GPP to changing   25 .This is also because of the model linkages between   and , in which   25 determines the photosynthetic capacity and indirectly affects   via its direct impact on  .The lower sensitivities of    and   ℎ to   Since  is multiplied with  in the BWB model, the influence of  on   depends on , and  depends on available energy for photosynthesis and the CO2 carboxylation rate, represented by   25.When  is low, the multiplicative impact of  on   is minimal while when  is high, the multiplicative impact of  on   will be also high, and both  and   jointly influence .Plants regulate transpiration via their stomata, so   is influenced by the source, the availability of water to plants, and the demand, the humidity and heat stresses in the ambient environment of plants.Therefore,  may be physiologically determined by environmental variables which can depict water availability and photosynthetic activities, and