Modelling decadal trends and the impact of extreme events on carbon fluxes in a temperate deciduous forest using a terrestrial biosphere model

Changing climatic conditions pose a challenge to accurately estimating the carbon sequestration potential of terrestrial vegetation, which is often mediated by nitrogen availability. The close coupling between the nitrogen and carbon cycles controls plant productivity and shapes the structure and functional dynamics of ecosystems. However, how carbon and nitrogen interactions affect both carbon fluxes and plant functional traits in dynamic ecotones, which are experiencing biotic and abiotic changes, remains unclear. In this work, we use in situ measurements of leaf chlorophyll content (Chl_Leaf, 2013–2016) and the leaf area index (LAI, 1998–2018) to parameterize the seasonal dynamics of the QUINCY (QUantifying Interactions between terrestrial Nutrient CYcles and the climate system) terrestrial biosphere model (TBM) to simulate the carbon fluxes at the Borden Forest Research Station flux tower site, Ontario, Canada, over 22 years from 1996 to 2018. Our goals are to assess the additional value of using Chl_Leaf in the model parameterization, to study how well QUINCY can capture observed trends related to the carbon cycle at the site, and to investigate how well the processes associated with a drought year and its legacy effects are captured by the model.

QUINCY was able to simulate leaf-level maximum carboxylation capacity (V_c(max),25), Chl_Leaf and leaf nitrogen in a manner quite consistent with observations. The model with the improved parameterization captured observed daily gross primary production (GPP) well (r²=0.80, root mean square error (RMSE) = 2.2 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). Nevertheless, we found that although observed GPP increased significantly during the study period (22.4 $\mathrm{gC}{\mathrm{m}}^{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) and net ecosystem exchange (NEE) shifted towards a stronger sink, these trends were not captured in the model. Instead, QUINCY showed a significant increasing trend for total ecosystem respiration (TER) that was not present in the observations. The severe drought in 2007 strongly affected observed carbon fluxes, lowering both GPP and TER in the following year as well. QUINCY was able to capture some of the decrease in GPP and TER in 2007. However, the legacy effect of the drought in 2008 was not captured by the model. These results call for further work on representing legacy effects in TBMs, as these can have long-lasting impacts on ecosystem functioning.

1 Introduction

Climate change impacts the exchange of carbon (C), water and energy between vegetation and the atmosphere, as well as the biogeochemical cycles and carbon storage potential of ecosystems (Canadell et al., 2022). The future land carbon balance can be predicted by terrestrial biosphere models (TBMs), which are versatile tools for studying the effects of climate on biogeochemical cycles (Blyth et al., 2021). However, TBMs exhibit significant inconsistencies in their simulated results over space and time due to divergent representations of important biogeochemical processes (O'Sullivan et al., 2022). One way to address these inconsistencies and understand the most important development needs is to compare these models to data to gain a better understanding of their performance. The variables related to the nitrogen (N) cycle and long-term observations are particularly relevant (Zaehle and Dalmonech, 2011), as the role of the N cycle can constrain the C storage capacity of terrestrial ecosystems (Zaehle et al., 2015) and long-term observations can provide information on how a changing climate affects the ecosystem. These long-term observations may also include anomalous years, such as years with droughts, one of the major stressors of extreme events that can have profound impacts on the carbon cycle (Piao et al., 2019). Droughts can also have legacy effects for years to come, and these legacy effects can vary between forests according to species and structure (Yu et al., 2022).

The C and N cycles are closely interconnected, with N being a significant component of plants and a vital macronutrient. N is required for growth, development and metabolic processes and is a fundamental constituent of DNA and various plant structural and photosynthetic components, such as light-harvesting complexes and the electron transport chain (ETC). Additionally, it is an integral component of many enzymes involved in the Calvin cycle. The N deficit therefore limits photosynthesis, which ultimately decreases plant productivity (LeBauer and Treseder, 2008). The future carbon pools and budgets in the coming decades depend, in part, on the N cycle and availability of N to vegetation (Arora et al., 2020; Huntingford et al., 2022; Zaehle, 2013).

In recent decades, many TBMs have incorporated elements of the N cycle to varying degrees (Thornton et al., 2009; Sokolov et al., 2008; Zaehle and Friend, 2010). However, this inclusion presents many challenges, e.g. in terms of N uptake (Davies-Barnard et al., 2022) and losses (Meyerholt and Zaehle, 2018) and particularly in regards to how models represent the N limitation of photosynthesis (Medlyn et al., 2015; Thomas et al., 2015; Walker et al., 2021). N limitation may directly affect photosynthesis rates or its effects may be buffered via different stoichiometry-related implementations (Thomas et al., 2015), and the different hypotheses and parameter values related to N cycle processes lead to differences between models (Medlyn et al., 2015). Various modelling approaches to the N cycle have resulted in different ecosystem responses e.g. to carbon dioxide (CO₂) fertilization (Arora et al., 2020; Meyerholt et al., 2020a; Thomas et al., 2013). A model intercomparison study of five CMIP6 (Coupled Model Intercomparison Project) models showed a wide range of response in net primary production for increased atmospheric CO₂ and atmospheric N deposition (Davies-Barnard et al., 2020).

Given the importance of N in physiological, biochemical and structural processes, novel data sources linked to the N cycle and its connections to the C cycle are necessary for model development (Kou-Giesbrecht et al., 2023; Meyerholt et al., 2020b). Remotely sensed observations of ecosystem traits are an important means of parameterizing the temporal dynamics of these models in a spatially explicit manner (Rogers et al., 2017). Satellite observations have shown that the leaf area index (LAI) has increased globally and has contributed to an increase in the land carbon sink (Chen et al., 2019) and a change in the Bowen ratio of energy fluxes (Forzieri et al., 2020). Simultaneously, lowering of leaf N content has been observed through intensive long-term monitoring plots in European forests (Jonard et al., 2015). Climate-change-induced changes led to increases in LAI (Chen et al., 2019) and tree productivity (Jonard et al., 2015), and the changes in N availability and demand have led to declining N availability relative to demand in terrestrial ecosystems (Mason et al., 2022).

One candidate that brings N observations and remote sensing together is leaf chlorophyll content (Chl_Leaf), which can be used as a proxy for the photosynthetic N component (Croft et al., 2017) and can be accurately retrieved at ecologically relevant temporal and spatial resolutions from remote sensing data, due to the presence of large absorption features in spectral bands typically sampled by optical sensors (Croft and Chen, 2018). There is a large body of literature on leaf chlorophyll retrieval from remote sensing data (Sims and Gamon, 2002; Dash and Curran, 2004), leading to the creation of large-scale national and global products spanning several years (e.g. Croft et al., 2020). Most TBMs simulate photosynthesis through the Farquhar–von Caemmerer–Berry (FCB) kinetic enzyme model (Farquhar et al., 1980), where photosynthetic capacity is represented by the maximum carboxylation capacity (normalized to 25 °C; V_c(max),25) and the maximum rate of electron transport (normalized to 25 °C; J_(max),25). Integrating physiological information through Chl_Leaf data has led to developments in improving modelling gross primary production (GPP). For example, Houborg et al. (2013) developed a semi-empirical relationship between Chl_Leaf and V_c(max),25 and used remotely sensed Chl_Leaf to replace V_c(max),25 in the Community Land Model (CLM) to improve simulations of GPP for a maize field. Lu et al. (2022) used observations of Chl_Leaf and V_c(max),25 to create linear relationships dependent on the plant functional type (PFT) and successfully retrieved the V_c(max),25 parameter at several ecosystems. At the site level, at the Borden Forest Research Station (hereafter referred to as Borden Forest), Luo et al. (2018) improved GPP modelling by directly linking the seasonal cycle of the V_c(max),25 parameter to Chl_Leaf. Our present study differs from previous literature in that we model Chl_Leaf explicitly and it is based on a predicted nitrogen cycle.

In this research, we use the long time series of in situ LAI and Chl_Leaf observations at Borden Forest together with eddy covariance observations of carbon fluxes over a decadal timescale. The long-term nature of the continuous flux record at Borden Forest with 22 years of near-continuous data provides an almost unparalleled opportunity to examine longer-term trends in ecosystem processes at a deciduous forest site due to warming temperature against a background of temperature, drought variability and extremes, and their legacy effects. The seasonal cycle of LAI has several implications for the climate–vegetation exchange, including influencing canopy conductance and fluxes of water, energy and carbon dioxide (Richardson et al., 2013). The development of Chl_Leaf and canopy structural parameters, such as LAI, decouple during the shoulder seasons in deciduous forests, necessitating the separate parameterization of leaf-level physiological processes and LAI in TBMs (Croft et al., 2015). Continuous long-term ground-based observations of LAI at the site scale (Rogers et al., 2021) provide a means to assess these phenomena in one forest.

In our research, the observational data are combined with a terrestrial biosphere model, QUINCY (QUantifying Interactions between terrestrial Nutrient CYcles and the climate system) (Thum et al., 2019). QUINCY simulates fully coupled carbon–nitrogen–phosphorus cycles of terrestrial ecosystems along with water and energy budgets. QUINCY is one of the few TBMs that explicitly represents relationships between time-varying foliar nitrogen content, which vary given the ecosystem's nutrient availability and carbon uptake capacity, and Chl_Leaf and photosynthetic activity (such as V_c(max),25). The model treats the impact of leaf chlorophyll and its vertical distribution on leaf- and canopy-level photosynthesis using an extension of the FCB model (Kull and Kruijt, 1998). QUINCY includes a representation of plant growth separating sink and source processes, acclimation of many ecophysiological processes to meteorological and/or nutrient availability, and the explicit representation of vertical soil processes. Whilst the impacts of climatic events and longer-term climatic shifts are complex to model, we hypothesize that QUINCY can capture changes in ecosystem function if they are related to meteorological conditions and atmospheric CO₂. Additionally QUINCY has the potential to capture legacy effects of extreme events via its carbohydrate pool structure. This improved modelling capacity will enable us to better understand the nitrogen and carbon cycles under both episodic events and over interannual and decadal timescales in a temperate deciduous–boreal ecotone. The diverse observations available at the site allow us to evaluate the model's performance in several aspects. The following four research questions are addressed in this work:

• How does the decoupling of LAI and Chl_Leaf seasonal development in the model affect the estimation of annual carbon fluxes?

• Is the QUINCY model able to simulate any long-term changes in seasonal shifts in carbon fluxes and LAI values?

• Is there a nitrogen constraint on carbon fluxes at Borden Forest, and does it change over the 23-year period?

• Can the QUINCY model simulate the effects of drought events on the carbon cycle?

2 Materials and methods

2.1 Study site

Borden Forest (44°19^′ N, 79°56^′ W) is a mixed forest situated in southern Ontario, Canada. This forest is located in the Great Lakes–St. Lawrence forest ecotone, which is a transition zone that includes both southern temperate forest species and northern boreal species (Froelich et al., 2015). Based on the 2006 vegetation survey (Teklemariam et al., 2009), the forest species composition was primarily red maple (Acer rubrum) and eastern white pine (Pinus strobus) with 52 % and 14 % respectively. Other species included large-tooth aspen (Populus grandidentata, 8 %), white ash (Fraxinus americana, 7 %) and trembling aspen (Populus tremuloides, 3 %). The forest has 15 %–20 % evergreen coniferous vegetation. The understory consists of short ferns, small shrubs and saplings (Halliday, 2010). The forest has been naturally regrown from farmland that was abandoned in the early 20th century, with a canopy height of approximately 22 m (Froelich et al., 2015). The soil consists mainly of sand with a thin layer of organic matter. The mean annual temperature at the site is 7.4 °C, and the mean annual precipitation is 784 mm over the 2000–2014 period (Froelich et al., 2015). The site is a member of the AmeriFlux network (site ID: CA-Cbo).

2.2 Site-level observations

Site-level measurements of carbon fluxes, meteorology and soil moisture, LAI, leaf nitrogen, and Chl_Leaf and biochemical model parameters (J_(max),25 and V_c(max),25) were used in the study. In this study, we define LAI as half of the total (all-sided) leaf area per unit of ground area (Chen and Black, 1992). Three datasets were used, two of which were long-term, (i) carbon fluxes and meteorology (1996–2018) and (ii) LAI (1998–2018), along with a leaf-level biochemical dataset from 2013 to 2016. Meteorological data were used to force our model simulations, while soil moisture and temperature observations (from 2005–2015) were used to evaluate the model's performance.

2.2.1 Net ecosystem exchange (NEE) of CO₂, meteorological and soil moisture observations

CO₂ flux data from half-hourly eddy covariance measurements sampled at the Borden Forest tower at 33 m height between 1996 and 2018 were used for model evaluation. The instrumentation is described in detail in Froelich et al. (2015). No observations were made in 2004 due to instrument and tower replacement. The fluxes were determined on a half-hourly timescale using a programme developed at the State University of New York at Albany (SUNY Albany) (Froelich et al., 2015) up to 2013 and using EddyPro (Fratini and Mauder, 2014) thereafter. The vegetation remains uninterrupted from 1.5–4 km towards the southeast and southwest and 1 km towards the northeast. However, there is cropland less than 400 m away in the northwest direction. Data for wind directions between 285 and 20° were excluded from analysis due to insufficient fetch and gap-filled (Luo et al., 2018). Also observations recorded when the friction velocity was less than 0.3 m s⁻¹ were removed according to Froelich et al. (2015). The data were gap-filled, and the measured carbon dioxide net ecosystem exchange (NEE) flux was partitioned into half-hourly GPP and total ecosystem respiration (TER) according to Barr et al. (2004). This is the standard method used for the FLUXNET Canada sites (Pierrat et al., 2021) and has been used in the previous publications from this site. This procedure first derives the component fluxes from the NEE and then uses simple empirical models constrained by the measured data (Barr et al., 2004; Rogers, 2022). The other empirical relationship is between TER and soil temperature at a shallow depth, and the other is between GPP and photosynthetically active radiation (PAR) above the canopy, also including a time-varying parameter (Barr et al., 2004; Rogers, 2022). These time-varying parameters are determined using a flexible moving-window approach, including 100 points (Barr et al., 2004; Rogers, 2022). The exact formulations used are as in Rogers (2022).

Air temperature, relative humidity, air pressure, longwave and shortwave radiation, and wind speed and direction were also measured by instruments on the flux tower (Froelich et al., 2015). For the air temperature, relative humidity, radiation and wind speed we used observations made at 42 m; for air pressure observations were at 2 m. Soil temperature (from 5 to 100 cm depth) and profiles (depths from 2 to 100 cm) were measured at two locations (Froelich et al., 2015), one 40 m southwest of the flux tower and the other 50 m west of it. Precipitation data were obtained from the nearby Egbert weather station (44°23^′ N, 79°78^′ W), which has provided hourly observations since 2014. Prior to 2014, the hourly precipitation was obtained from the ERA5-Land product (Muñoz Sabater et al., 2021) and scaled to match the average precipitation values as estimated from the hourly Egbert observations (more information in Sect. S1). In addition, the ERA5-Land product was used in the gap filling of observed meteorological data.

2.2.2 LAI observations

We used a daily LAI time series of 1999–2018 from Rogers et al. (2021), estimated from photosynthetically active radiation (PAR) observations collected above and below the canopy. PAR was measured above the canopy with a LI-COR LI-190SA (LI-COR, Lincoln, Nebraska) sensor and below the canopy with a LI-COR LI-191 sensor. We estimated daily LAI values from half-hourly observations of above-canopy and transmitted PAR using the Miller integral (Miller, 1967), as recommended by Rogers et al. (2021). To improve spatial representativeness of the daily LAI estimates at the site, the values were then calibrated to match effective LAI (L_e) measured along a 100 m transect using a handheld LI-COR LAI-2000 plant canopy analyser (LI-COR, Lincoln, Nebraska) using a linear relationship. True LAI was estimated from all the observations as

$$\begin{array}{}\text{(1)}& \mathrm{LAI}={\displaystyle \frac{\left[\right(\mathrm{1}-\mathit{\alpha }\left)L_{\mathrm{e}}{\mathit{\gamma }}_E\right]}{{\mathrm{\Omega }}_E}},\end{array}$$

where α is the ratio of woody area to total area and γ_E is the ratio of needle area to shoot area (taken as unity in deciduous forests). The value of α (0.17) was taken from the literature (Gower et al., 1999), similar to Croft et al. (2015). The clumping index (Ω_E) is 0.95 and was measured using a TRAC (Tracing Radiation and Architecture of Canopies) instrument (Huiming Instrumentation Limited, Nanjing, China) (Rogers et al., 2021).

2.2.3 Leaf-level trait measurements

In situ leaf-level data included measurements of Chl_Leaf and nitrogen content, maximum electron transport and carboxylation rates (J_max,25, V_c(max),25), and specific leaf area (SLA). The measurements were collected at an average interval of 9 d during the growing seasons (day of year 130–290) from 2013 to 2016. In 2013, only leaf chlorophyll and SLA data were collected. In 2014, all six variables were measured. In the following years, only Chl_Leaf, J_max,25 and V_c(max),25 were measured. For these biochemical measurements, leaves were sampled from the top of the canopy from the flux tower. Leaf-level gas exchange measurements were made using an LI-6400 portable infrared gas analyser (LI-COR, Lincoln, Nebraska, USA) (Croft et al., 2017). Chl_Leaf and leaf nitrogen content were destructively analysed in the laboratory from five sampled leaves per measurement date. The methodology and the measurements are described in detail in Croft et al. (2013, 2014) and Wellburn (1994). The SLA was calculated as the ratio of leaf area to leaf dry mass. In this study, we used a species-weighted canopy average of the leaf-level parameters, based on the species composition of the forest (i.e. red maple, 60.4 %; large-tooth aspen, 12.9 %; trembling aspen, 12.4 %; and ash, 14.2 %) (Croft et al., 2015). The values for the different tree species are given in Table S2, and the calculation of the species-weighted average is explained in Sect. S2.

2.3 The QUINCY model

2.3.1 General description

The QUINCY model (Thum et al., 2019; Zaehle et al., 2019) was used to simulate ecosystem functioning in the study area. QUINCY simulates fully coupled carbon, nitrogen and phosphorus cycles as well as water and energy balances in vegetation and soil on a half-hourly timescale. Here we give a brief description of the model, focusing on the parts relevant to this study. A more detailed description can be found in Thum et al. (2019).

Vegetation is grouped by plant functional type (PFT) and represented as an average individual composed of structural pools (leaves, sapwood, heartwood, coarse roots, fine roots and fruits), a labile pool (a fast overturning and respiring non-structural pool) and a reserve pool (seasonal, non-respiring and non-structural storage pool). The non-structural pools represent storage pools for non-structural carbohydrates and associated nutrients. Trees are also characterized by height (m), diameter (m) and stand density (m⁻²). The tree canopy is composed of 10 canopy layers, which increase in depth (of LAI) exponentially with layer depth (LAI_cl). Photosynthesis and stomatal conductance are calculated separately for the sunlit and shaded leaves in each layer, as estimated by a radiative transfer model (see below) (Zaehle and Friend, 2010) using a Farquhar-model-based scheme from Kull and Kruijt (1998). According to this scheme, the role of leaf chlorophyll is explicitly taken into account in the photosynthesis calculation and determines the proportion of the leaf area at each canopy layer that is light-saturated. Photosynthesis for the non-light-saturated part is calculated using the light-limited rate of photosynthesis (relying on the J_max,25 parameter). Photosynthesis in the light-saturated part is calculated as a co-limited rate of the electron transport capacity (determined using J_max,25) and carboxylation capacity (via Rubisco; determined using V_c(max,25)). Therefore both biochemical model parameters and leaf chlorophyll influence photosynthesis. The determination of these variables from leaf nitrogen is described in Sect. 2.3.3. The leaf stoichiometry can be set as fixed or it can be dynamic, when it is varied in response to nutrient demand and supply.

Stomatal conductance is modelled after Medlyn et al. (2011), and in addition to stomatal conductance, soil moisture can limit photosynthesis directly through a modifier in the calculation of the biochemical parameters (Egea et al., 2011). Photosynthesis can also be down-regulated by the sink limitation (example in Fig. S1 in Thum et al., 2019).

Maintenance respiration is a linear function of the N content of each pool. Temperature acclimation for photosynthesis is as in Friend (2010) and for maintenance respiration as in Atkin et al. (2014). Tissue growth is defined by allometric equations, and the allometric relationship between leaves and fine roots responds to N and water limitation by increasing the uptake capacity under nutrient limitation.

Soil biogeochemistry is largely based on a CENTURY-style (Parton et al., 1993) approach, except that the vertical soil profile of biogeochemical pools, including metabolic, structural and woody litter as well as fast- and slow-overturning soil organic matter (SOM), is explicitly represented. Each soil layer also includes N pools of ammonium (NH₄) and nitrate (NO₃). The soil profile consists of 15 layers, reaching a depth of 9.5 m. The depth of each layer increases exponentially as it goes deeper. The stoichiometry of the litter is determined by the stoichiometry of the plant pool it comes from. The fast pool's stoichiometry is dependent on the availability of inorganic nutrients, while the slow pool has a fixed stoichiometry. Plants and microbes compete for nutrients based on their respective demand and uptake capacity. NH₄ is oxidized to NO₃ through nitrification in the aerobic part of the soil, and NO₃ is reduced to diatomic nitrogen (N₂) through denitrification in the anaerobic part of the soil (Zaehle et al., 2011). Both processes also produce nitrogen oxide (NO_y) and nitrous oxide (N₂O). Biological nitrogen fixation (BNF) is considered to be an asymbiotic and symbiotic process (Meyerholt et al., 2016).

Soil temperature and moisture are calculated for each layer based on soil physical characteristics, as well as the transport and atmospheric exchange of energy and water. The radiative transfer scheme has been developed following the two-stream approach by Spitters (1986), and the original implementation in the terrestrial biosphere model O-CN (Zaehle and Friend, 2010) has been extended to include diagnostic canopy albedo, clumping and attenuation of the shortwave backscatter from the soil. The radiative transfer is calculated separately for the visible and near-infrared radiation bands. It estimates the sunlit and shaded leaves for each canopy layer and separates the incoming radiation into direct and diffuse components. Snow dynamics is a later addition to the model (Lacroix et al., 2022), and the model can be run with or without the snow. The snow dynamics are represented in a five-layer scheme that accounts for flows of heat and water within and between the snow layers.

2.3.2 Phenology of the deciduous trees

The seasonal development of leaf biomass is affected by the plant's ability to grow new tissues and fractional allocation to plant organs. The start and end of the growing season are determined by meteorological and soil moisture values, which are averaged over 7 d to mitigate the impact of daily climate variability. The start of the growing season is determined by the accumulated growing degree days (GDD_acc), representing the current number of growing degree days above the temperature threshold ($t_{\mathrm{air}}^{\mathrm{GDD}}$) since the beginning of last dormancy period as

$$\begin{array}{}\text{(2a)}& {\displaystyle }{\mathrm{GDD}}_{\mathrm{acc}}& {\displaystyle }>{\mathrm{GDD}}_{\mathrm{req}}^{\max }\times {\exp }^{-k_{\mathrm{dormancy}}^{\mathrm{GDD}}\times \mathrm{NDD}},\mathrm{where}\text{(2b)}& {\displaystyle }& {\displaystyle \frac{{\mathrm{GDD}}_{\mathrm{acc}}}{\mathrm{d}t}}={\mathrm{GDD}}_{\mathrm{acc}}+\mathrm{MAX}(t_{\mathrm{air}}-t_{\mathrm{air}}^{\mathrm{GDD}},\mathrm{0.0}),\end{array}$$

where NDD is the number of dormancy days, taken as days since the last growing season; $k_{\mathrm{dormancy}}^{\mathrm{GDD}}$ (value of 0.007 d⁻¹) relates dormancy to the maximum growing degree-day requirement (GDD${}_{\mathrm{req}}^{\max }$, 800 degree days) to account for the chilling requirements of the buds (Krinner et al., 2005); and dt denotes the time step in days. The growing season ends when the decreasing average air temperature falls below the temperature threshold ($t_{\mathrm{air}}^{\mathrm{sen}}$, 8.5 °C) and then senescence occurs.

2.3.3 Leaf N partitioning

While the overall amount of leaf nitrogen is largely driven by phenological development, the leaf N concentration per leaf area (N_leaf) responds to soil nutrient availability, as plants take up mineral N from the soil pools. This uptake is determined by the amount of N in each soil pool and fine root biomass and is further modulated by plant N demand. Leaf nitrogen has a vertical gradient that decreases exponentially towards the bottom of the canopy, in accordance with observations (Niinemets et al., 1998). Leaf N in each layer (N_leaf,cl) is divided into structural and photosynthetic parts (Friend et al., 1997). The fraction of structural N (fN_struc,cl) is calculated for each canopy layer as a function of the total leaf N in the respective layer (Zaehle and Friend, 2010):

$$\begin{array}{}\text{(3)}& {\mathrm{fN}}_{\mathrm{struc},\mathrm{cl}}=k_{\mathrm{0}}^{\mathrm{struc}}-k_{\mathrm{1}}^{\mathrm{struc}}{N}_{\mathrm{leaf},\mathrm{cl}},\end{array}$$

where $k_{\mathrm{0}}^{\mathrm{struc}}$ is the maximum fraction of structural leaf N (0.63 for deciduous forest; Friend et al., 1997; Kattge et al., 2011) and $k_{\mathrm{1}}^{\mathrm{struc}}$ is the slope of structural leaf N with total N (7.14×10³ gN⁻¹) (Friend et al., 1997).

The photosynthetic N pools have three compartments: Rubisco associated (fN_rub), electron transport associated (fN_et) and chlorophyll associated (fN_chl). The photosynthetic fractions all play a role in the calculation of photosynthesis (Kull and Kruijt, 1998). The fractions are used directly in the calculation of the photosynthetic parameters V_c(max),25 and J_max,25, where the leaf N in each content is multiplied with these fractions and some other modifiers (Eqs. S7 and S10 in Thum et al., 2019). According to Zaehle and Friend (2010), the fraction of leaf N in chlorophyll is calculated to increase with canopy depth:

$$\begin{array}{}\text{(4)}& {\mathrm{fN}}_{\mathrm{chl}}={\displaystyle \frac{k_{\mathrm{0}}^{\mathrm{chl}}-k_{\mathrm{1}}^{\mathrm{chl}}{e}^{-k_{\mathrm{fn}}^{\mathrm{chl}}{\mathrm{LAI}}_{\mathrm{c}}}}{a_{\mathrm{chl}}^n}},\end{array}$$

where $k_{\mathrm{0}}^{\mathrm{chl}}$ (value of 6.0; Zaehle and Friend, 2010), $k_{\mathrm{1}}^{\mathrm{chl}}$ (value of 3.6; Zaehle and Friend, 2010) and $k_{\mathrm{fn}}^{\mathrm{chl}}$ (value of 0.7; Friend, 2001) are empirical parameters. $a_{\mathrm{chl}}^n$ is the molecular N content of chlorophyll (25.12 mol mmol⁻¹; Evans, 1989). LAI_c is the cumulative leaf area.

The photosynthetic parameters V_c(max),25 and J_max,25 are assumed to have a fixed ratio of 1.97 (Wullschleger, 1993). Based on this ratio, fN_rub and fN_et are calculated, using the calculated values of the structural and chlorophyll fractions.

2.3.4 Modelling protocol

The QUINCY model requires half-hourly meteorological forcing, including short- and longwave radiation, air temperature, precipitation, air pressure, humidity, and wind speed as well as atmospheric CO₂ concentration (hereafter referred to as [CO₂]) and N and P deposition rates. Meteorological forcing was measured at the site, and the [CO₂] was obtained from Friedlingstein et al. (2019) and N deposition data from Lamarque et al. (2010, 2011). Borden Forest is described as a broadleaf deciduous forest PFT, and we performed point scale simulations.

To obtain a near-equilibrium state of soil and vegetation, the model was spun up for 500 years, using atmospheric CO₂ concentration from a randomly selected year from the 1901–1930 period and a random year of observed meteorological data. This was followed by a transient simulation starting in 1901, which was using meteorological data randomly picked from the observed meteorology. The transient simulation used the atmospheric [CO₂] and N deposition values derived from data sources mentioned above and, from 1996 onwards, the measured site-level meteorology for the respective years. For the purposes of this study, we ran simulations where only the C or both the C and N cycles were active. In the simulation where only C was active, the plants had access to all the N that they needed. The P concentration was kept at a level where it did not limit plant uptake or SOM decomposition. The temperature response of the BNF (Bytnerowicz et al., 2022) was set to have an optimum temperature of 18 °C, replacing the default value of 32 °C. The default value is based on observations in the tropics, and the shape of the curve predicts very low BNF for more northern regions with the default optimum temperature. Lowering this value to typical air temperature at the Borden site provides more realistic BNF for the site and assumes local temperature acclimation. The SLA has a constant value of 320 cm g⁻¹ for the broadleaf deciduous forest PFT in the model simulations. To compare to the leaf-level observations (Chl_Leaf, V_c(max),25, J_max,25) made at the top of the canopy, we used only the top canopy layer values from the model.

The model was run with several different parameterizations to study the influence of using LAI and Chl_Leaf on the parameterization and how their decoupling influences the results. For the carbon-cycle-only (C) simulation, we show results from the original model formulation (orig), then from the simulation using LAI to tune phenology (LAI tuning, C:LAI), and finally from a simulation using both LAI and Chl_Leaf for tuning (C:LAI&chl). These simulations have been done with dynamic leaf stoichiometry, which is basically showing forest at N saturation. This simulation is shown in the first sections of this paper, since the nitrogen-cycle-enabled version showed too low a GPP at the site and we also wanted to show the influence of model tuning with LAI and Chl_Leaf with GPP levels comparable to the observations. For comprehensiveness, we also report the values for the carbon-cycle-only simulations with fixed stoichiometry, which reflects the C cycle with average N availability but with no N limitation on growth and soil processes. We made an additional simulation of adding snow (CN:LAI&chl&snow) to assess its influence on the soil temperature and carbon fluxes.

For the C and CN simulation, we use the results after both LAI and Chl_Leaf have been used in the tuning in the main analysis. The abbreviations of the model simulations are found in Table 1. A schematic figure (Fig. S1) showing the workflow of the study, with parameter tuning and a comparison to the observations, is found in the Supplement.

Table 1 Abbreviations used for the model runs.

The original parameter values and the tuned parameter values are shown in Table 2. We adjusted the parameters by comparing the modelled LAI and Chl_Leaf to the observations (Fig. 1b, e) and tried to match those. The autumn decline in the simulated LAI was adjusted to match the observations by modifying the parameter controlling leaf senescence ($t_{\mathrm{air}}^{\mathrm{sen}}$) from the default value of 8.5 to 15.0 °C. It is likely that leaf senescence at the site is partially controlled by light availability, a process which is not yet present in the model; therefore the higher temperature threshold accounts for this missing factor. To adjust the summertime magnitude of Chl_Leaf, the parameters $k_{\mathrm{0}}^{\mathrm{struc}}$ (Eq. 3) and $k_{\mathrm{0}}^{\mathrm{chl}}$ (Eq. 4) were adjusted. In the model the whole leaf nitrogen is initially allocated to the structural nitrogen during the initial stages of the growing season. Subsequently, the amount allocated to the photosynthetic compartments (including leaf chlorophyll, V_c(max),25 and J_max,25) begins to increase as the season progresses. The model was modified to incorporate a delay in the transition from structural nitrogen to photosynthetic nitrogen. We added in a delay of 20 d by introducing a leaf age factor to the simulations (equations for this change: Eqs. S2–S4) in the C simulations and delay of 15 d in the Cfix:LAI&chl and CN:LAI&chl simulations.

Table 2 Parameter values in different simulations. The unit of the parameter is in parentheses after the parameter name.

2.3.5 Estimation of seasonal metrics and trends

We estimated the growing-season metrics separately using GPP or LAI. The start and end of the season (SOS and EOS) estimated from GPP were calculated from the first and last pass of the threshold, which was defined as the 30 % of the year's 90th-percentile value (an example year of 2013 in Fig. S2a). For LAI the threshold was calculated as being 20 % of the difference between the summer and winter values, starting from the winter value (Fig. S2b). Length of the growing season (LOS) is the time between SOS and EOS. These calculations were made on smoothed data using an averaging weekly window to minimize anomalies. The trend assessment was carried out with a particular focus on statistically significant trends, which were identified through the application of Student's t test on slope values obtained from the linear regression (p<0.05).

3 Results

3.1 Dynamic parameterization of LAI and leaf chlorophyll content improves modelled GPP estimates

The average modelled GPP, LAI and Chl_Leaf values for the three C-cycle-only model simulations (C:orig, C:LAI and C:LAI&chl) across a growing season from 1996 to 2018 are shown alongside the measured data in Fig. 1a, c and e. The observed GPP starts to increase already after DOY (day of year) 100, whereas in all the simulations, the increase begins later and at a faster rate. The increase to maximum summer time values in the simulations happens rapidly, and the maximum summer values occur early in the season, around DOY 160 (early June). The observations show a more gradual increase in GPP and a sharp increase near the peak, with maximum summertime values occurring around DOY 200 (mid-July).

At the end of the season, the inclusion of LAI (C:LAI) and Chl_Leaf (C:LAI&chl) data in the tuning improved the representation of senescence at the end of the season and the consequent decline in GPP, compared to the C:orig simulation (Figs. 1a, S3). The summertime average was accurately simulated by the model, with an average overestimation of 1.1 % for June, July and August with the C:orig and C:LAI simulations and 5.0 % overestimation with the C:LAI&chl simulation (Fig. 1a). The annual carbon flux values together with root mean square error (RMSE) and R squared (r²) values are shown in Table 3. The different C model simulations (i.e. C:orig, C:LAI and C:LAI&chl) using dynamic stoichiometry (i.e. not including simulations using fixed stoichiometry, Cfix) did not largely impact RMSE and r² values (Table 3).

Table 3 Annual average carbon fluxes for different model parameterizations, together with their standard deviations and root mean square error (RMSE) and r² values. RMSE and r² have been calculated from the daily averages for the time period 1996–2018. Also shown is the discrepancy in percent between the simulation result and the observation, with positive values denoting an overestimation and negative denoting an underestimation.

Figure 1 Average yearly cycles of (a) GPP, (c) LAI and (e) leaf chlorophyll from the C simulations and for the CN simulations GPP in (b), LAI in (d) and leaf chlorophyll in (f). The shaded regions show the standard deviations between the years. The data represent the mean values for (a–b) 1996–2018, (c–d) 1999–2018 and (e–f) 2013–2016. The lines have been smoothed with a 7 d averaging window, except for the observed leaf chlorophyll which has been smoothed with a 3 d window.

The simulations with the N constraint on carbon fluxes (CN simulation) demonstrate a reduction in summertime GPP values, with the average GPP during July–August underestimated by 14 % compared to observations (Fig. 1b). The same modification to the phenology parameter was made as in the C simulation, with the objective of improving the fit of the simulated LAI (CN:LAI simulation). This resulted in a more accurate representation with better timing of senescence of the observed seasonal cycles of GPP, LAI and Chl_Leaf in the simulations (Figs. 1b, d, f; S3).

The use of LAI in the model tuning had a more pronounced impact on the GPP fluxes in both the C and CN simulations because of the shortening of the growing season, even though the changes in r² (on the order of 0.01) or RMSE remain relatively minor. After LAI and Chl_Leaf tuning, the simulation had a 3.6 % overestimation in the annual GPP in the C simulations and an underestimation by 17 % in the CN simulations (Table 3). The impact of the N constraint on the annual GPP was 302 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Table 3), representing a 20 % decrease in the annual GPP value relative to the C simulation. When considering the Cfix:LAI&chl simulation, the estimated GPP was very similar to the CN simulations, with only a 1.4 % larger value (Table 3).

From this point onward in the paper, when we refer to the C simulations, we refer to the results from the simulations with the tuning of LAI and Chl_Leaf (C:LAI&chl), and the same is done for CN simulations. CN:LAI&chl was the most successful of the simulations in terms of r² and RMSE when compared against observed GPP (Table 3).

3.2 Other carbon fluxes: TER and NEE

TER was decreased compared to the original simulation after altering the parameters (Table 2) for both C and CN simulations (Table 3). This is connected to the declining GPP, as the litter input influences the amount of soil carbon. The seasonal cycle was not influenced by these changes (Fig. S4). The CN simulations had better r² values for TER than the C simulations (Table 3). This occurred because the magnitude of TER was better captured by these simulations (Fig. S5b). The most pronounced underestimation of TER in the CN simulations by the model compared to the observations occurs during the summer months of July and August (Figs. 2b, S5b).

The observed annual TER was overestimated by the C simulation by 18 % (Table 3). The CN simulation yielded a more accurate representation of the annual estimate (Table 3). The Cfix:LAI&chl simulation gave similar values to the CN simulation (Table 3). The nitrogen cycle was found to constrain annual TER by 300 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, representing a 20 % decrease from the C simulation.

The observations indicated that the forest acted as a sink, with a net carbon uptake of −205 ± 140 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ over the measurement period. The simulations indicated that it was a weak sink, for both the C and the CN simulations, with a very small difference (Table 3). However, it should be noted that the interannual variation in the simulations was considerable. A steep increase in GPP in the C simulations led to a short-term overestimation of the ecosystem carbon sink in early summer (Figs. S5c, S6), but generally during the summer the ecosystem carbon sink was underestimated in the simulations. Despite the GPP summertime magnitude being underestimated in the CN simulations, the performance of the model as estimated by r² and RMSE was better for NEE and the component fluxes with the CN than C simulations (r² better by 0.08 for GPP, 0.12 for TER, 0.15 for NEE) (Table 3, Fig. S5).

Figure 2 Average seasonal cycles of net ecosystem exchange (NEE) (a), total ecosystem respiration (TER) (b), soil temperature at 5 cm depth (c) and soil moisture at 5 cm depth (d) averaged over the period 2005–2015. The observations are represented in black, the CN:LAI&chl model results are in orange and the standard deviation is shown by the shaded regions. Both the observations and simulations are smoothed with a 7 d window.

The early-season pattern observed in the simulated NEE is attributed to too late an onset of GPP and too early an onset of heterotrophic respiration (Figs. 2a, b; S7), which is regulated by the soil moisture and soil temperature in the model. QUINCY simulates these soil conditions based on the meteorological conditions and soil texture. The majority of the simulated heterotrophic respiration originates from the uppermost soil layers, as most of the organic material is located in these layers in the model. The soil temperature of the upmost layer is significantly underestimated during the winter months, and an increase to the summertime levels occurs earlier in the simulations than in the observations (Fig. 2c). This phenomenon occurs during the period between DOY 70 and 150. It is during this time that TER is being overestimated in the spring (Fig. 2b, c). The maintenance respiration is also activated during this time period, although its increase is less pronounced due to the temperature acclimation implemented in the model (Fig. S7). In autumn, the simulated drawdown of TER occurs simultaneously with the observations, despite the simulated soil temperature decreasing at a faster rate than that observed (Fig. 2).

The soil conditions are relevant to the estimation of TER, as heterotrophic respiration is an important part of it (Fig. S7). Too early an increase in simulated soil temperatures in spring occurs also at deeper depths (Fig. S8b, c). Also at deeper levels at 10 and 20 cm the wintertime soil moisture is underestimated (Fig. S8e, f). The observed summertime variability in soil moisture is better captured by the model in deeper layers than in the 5 cm depth, even though the summertime magnitude is overestimated at 10 cm depth (Fig. S8).

The simulation with snow (CN:LAI&chl&snow) had an influence on the spring soil conditions and thus on TER (Figs. S8, S9). Overall the simulations successfully captured the observed snow depth, with some underestimation in some years, which also resulted in snow clearance dates that were too early (Fig. S10). Snow influenced both soil temperature and soil moisture (Fig. S8). In winter, the soil temperatures were colder in the CN:LAI&chl&snow simulation than in the CN:LAI&chl simulation, and despite a fairly synchronous increase from the winter levels, the simulations with snow were slower to reach above-zero levels and therefore better matched the observations (Fig. S8a–c). Soil moisture was also affected, with snow delaying the increase to the summer levels (Fig. S8d–f). At 10 cm depth, the simulation without snow appears to be more consistent with spring observations of soil moisture (Fig. S8e). Snow in the simulations caused a slight delay in the increase in TER during spring (Fig. S9b), and the RMSE of TER during mid-March to mid-April (DOY 74–105) improved by 0.1 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for the 2005–2015 period. Overall, GPP and TER do not begin to increase to summer levels before snow clearance (Fig. S9).

3.3 Simulated structural and biochemical parameters

The continuous observations of LAI provide values of 3.78 ± 0.43 m² m⁻² in summer (averaged over the June–July period, along with the standard deviation). The simulated values were closer to the LAI-2000 observations than to the values obtained from the continuous observations. The LAI-2000 observations have a summertime average of 4.63 ± 0.71 m² m⁻² for the years 2013–2018. The C simulation overestimated continuous observations by 39 %, while the overestimation in CN simulations was 26 % (Table 4). The LAI from the CN simulations was closer to the value estimated from the LAI-2000 observations, with only 3 % overestimation. In the model the specific leaf area (SLA) is an important factor in converting the leaf biomass to LAI and is held constant. In the observations the SLA showed a dynamic change, with higher values (303 cm g⁻¹) observed in the early season and a subsequent decrease to a summer value of 162 cm g⁻¹ within about 1 month.

The C simulation overestimated observed leaf N by approximately 3-fold (Table 4). The CN simulation overestimated the observed value by 59 %. The C simulation estimated a 2.5-fold overestimation of observed J_max,25 (Fig. 3a, Table 4). Both the CN and C simulations overestimated the values of J_max,25 and V_c(max),25, with a considerably higher overestimation taking place in the C simulations (Fig. 3, Table 4). These high predicted values are not unexpected, given the direct link between plant N and photosynthetic parameters in the model and the implicit unlimited N availability in the C model.

Table 4 Observed and simulated (C:LAI&chl and CN:LAI&chl simulations) leaf traits. Values are estimated for the June–July means with standard deviations. The unit of the parameter is in parentheses after the parameter name. For LAI the values are for 1998–2018, for leaf N and the biochemical parameters the values are for 2014, and for Chl_Leaf the values are for 2013–2016.

The tuning of the model did not have a pronounced effect on the summertime magnitude of the photosynthetic parameters (Fig. 3) because the changes in the nitrogen allocated to leaf chlorophyll were derived from the structural nitrogen. The springtime delay imposed on the leaf chlorophyll also influenced the photosynthesis parameters, resulting in an improved seasonal cycle compared to observations (Fig. 3).

Figure 3 The seasonal cycle of J_max,25 (a, c) and V_c(max),25 (b, c) in 2014 with C simulations (a, c) and CN simulations (b, d). The observations are represented in black, the original QUINCY results from the C simulations are in light brown, LAI tuning is in blue, and both LAI and leaf chlorophyll tuning are in dark red. The CN:orig simulation results are in dark yellow, the CN:LAI simulation results are in magenta and CN:LAI&chl simulation results are in orange. The modelling results have been averaged with a 7 d smoothing window.

3.4 The influence of drought on carbon fluxes

A severe drought occurred at the site in 2007, when the precipitation was approximately 20 % lower than in a regular year, at 608 mm yr⁻¹ (Fig. S11). Figure 4 depicts the averages over the time period 2005–2015 (with the exception of 2007 and 2008) for GPP, TER and soil moisture at 5 cm, as the soil moisture observations were available for this time period. The following year, 2008, was characterized by precipitation 21 % above average (923 mm yr⁻¹) (Fig. S11).

The annual observed TER in 2007 was clearly below the average annual TER by 37 % (i.e. 794 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), and the level remained low throughout the summer (Fig. 4a). Furthermore, the summertime maximum TER values in 2008 were below average (Fig. 4b). Additionally, the values of TER exhibited a slower rate of increase and a more rapid decrease to winter levels after midsummer compared to regular years (Fig. 4b). This resulted in the observed TER for 2008 being 36 % below the annual average TER (i.e. 804 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). In contrast to the measurements, the simulations did not predict low TER for the beginning of the season in 2007 or 2008. Instead, the behaviour in the CN simulations was similar to that observed in other years (Fig. 4a–b). Only when there was a pronounced decrease in soil moisture around DOY 160 in 2007 did a decrease in the simulated TER result. Following the precipitation event around DOY 200, the TER values exhibited a recovery to typical summertime levels for several days. After this a decline to lower levels occurred. The precipitation events that occurred around DOY 240 resulted in TER returning to a regular level. In 2008 the simulated TER followed a pattern similar to regular years. The annual simulated TER was 9 % lower in 2007 compared to regular years and 10 % higher in 2008.

Observed GPP exhibited a decline in 2007, with the annual value being 24 % lower (i.e. 1112 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) compared to the annual average GPP. In 2008, the observed GPP was found to be lower than in 2007, with a value of 1086 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, representing a 26 % decline below the average. A later increase in GPP to summer values in spring and an earlier decrease to winter values were responsible for the reduction in the annual value, with peak season GPP being higher than in 2007 (Fig. 4c–d). The simulated GPP exhibited a regular pattern in both years until DOY 160 in 2007 (Fig. 4c). The annual simulated GPP was 14 % lower in 2007 compared to regular years and 15 % higher in 2008.

The measured soil moisture was consistently below the typical values throughout 2007, and that continued until early summer 2008 (DOY 150, Fig. 4e, f). In 2007 the simulated soil moisture was generally at a higher level compared to observations. However, the responses of the simulated soil moisture to precipitation events were similar to those in the observations (Fig. 4e). Overall, the simulated soil moisture exhibits a narrower range than the observed soil moisture data in regular years and declines less in the drought year than observations.

Figure 4 The average seasonal cycles for observations and CN simulations for total ecosystem respiration (TER) (a in 2007 and b in 2008), gross primary production (GPP) (c in 2007 and d in 2008) and soil moisture at 5 cm depth (e in 2007 and f in 2008). For these years, the observations are a thick black line and simulations are a thick orange line. For 2005–2015 without 2007 and 2008, the average seasonal cycle is a thin grey line and simulations are a thin red line. Standard deviation of the average seasonal cycle is shown as shaded regions. Both the observations and the model results have been smoothed using a 7 d averaging window.

The model demonstrated a less pronounced effect of drought on the carbon fluxes than was observed in 2007. The water potential of soil lowered the simulated GPP (Fig. S12), but the drawdown was not large enough (Fig. 4). Furthermore, the likely legacy effects of drought that were observed in 2008 were not replicated by the model. The non-structural carbohydrate pools were affected by the drought. The labile pool was 25 % lower in 2007 than in normal years but had fully recovered by 2008. The reserve pool was in 2007 at 18 % lower in 2007 than in normal years and was still 10 % lower in 2008. However, this did not significantly affect the LAI values or the annual GPP levels. Because the reduction in GPP was not pronounced enough, the legacy effect was not seen in the TER values. Overall, the simulated soil moisture at 5 cm depth did not reach the lower values in the observed range (Fig. S13a), so although the simulated carbon fluxes do respond to soil moisture (Fig. S13), the legacy effects were not captured.

In regular years a hysteresis effect of TER versus the soil temperature relationship was observed, with the values in the later half of the year having lower values (Fig. S14a). In drought year 2007 this effect was not visible, with values staying at a low level. In 2008 the hysteresis effect was more pronounced than in 2006. The model simulations were not able to replicate this kind of behaviour. The observed soil moisture at 5 cm depth did not explain the observed hysteresis effect (data not shown). As a possible explanation to the hysteresis effect, one could think about a strong connection between photosynthesis and respiration. The observations indeed show a higher GPP in early-summer months compared to later months in the year, especially in 2008 (Fig. S15). This phenomenon does not take place in the simulations (Fig. S15b, d, f).

3.5 Interannual variability and longer-term trends in annual carbon fluxes

The growing-season metrics were estimated based on both LAI and GPP from the observations and CN:LAI&chl simulation. The start of the season (SOS) takes place almost at the same time according to the GPP and LAI observation-based estimates. The end of the season (EOS) is estimated to be some days later in LAI-based estimates (Table 5). The simulations agree on both observation-based estimates of SOS by 1 to 5 d, and there is a larger difference between the observed and simulated EOS estimates. The simulated EOS based on LAI is earlier than observed (Table 5). In the simulations GPP and LAI are tightly coupled. However, in the observations the EOS estimated from LAI takes place on average 13 d later than the EOS according to GPP. The observed LAI remains high, despite decreasing GPP. LAI is therefore not so tightly coupled to the seasonality of GPP in autumn as it is in the model. C simulations provide similar estimates for these growing-season metrics.

The time series of the metrics indicates that the EOS from simulations has a larger range of variability (38 d from GPP-estimated EOS) than the observations (27 d from GPP-estimated EOS, 18 d from LAI-estimated EOS) (Fig. 5b). Additionally, the simulations demonstrated a greater interannual variation in the SOS estimation (range of 34 d for LAI-based estimation) than seen in the observations (range of 20 d for GPP-based estimation, 32 d for LAI-based estimation) (Fig. 5a). Stronger interannual variability in model estimates of LAI is driven by the phenological parameters governed by air temperature. The real forest with several species may be more resilient to different environmental conditions, and therefore different species may be able to benefit from different environmental conditions during the shoulder seasons. No discernible trends are evident in any of the time series under consideration. The cold spring of 2017 resulted in the simulation estimates of SOS occurring at a later point in time, although the impact was not as pronounced in the observations of GPP (Fig. 5).

Figure 5 The start of the season (SOS) (a) and end of the season (EOS) (b) as estimated from observed and simulated (CN:LAI&chl) time series of GPP (observation in black, simulation in orange) and LAI (observation in grey, simulation in light brown).

Table 5 The average start of the season (SOS), end of the season (EOS) and length of the season (LOS) as determined from the observations and the QUINCY model simulation (CN:LAI&chl) of GPP and LAI.

The next step involved the calculation of trends for the annual values of GPP, TER and NEE (Fig. S16, Table S3) derived from the CN:LAI&chl simulations. In order to assess the ability of QUINCY to capture observed interannual fluctuations in the annual carbon balance, the changes were investigated over time. This was done with a particular focus on statistically significant trends, which were identified through the application of Student's t test (p<0.05). During the time period 1996–2018 the observed GPP exhibited a significant increasing trend (Fig. S16, Table S3), with a magnitude of 22.4 $\mathrm{gC}{\mathrm{m}}^{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. A significant trend was also present in the summertime and autumn observations (Table S3). QUINCY showed a minor and non-significant upward trend for GPP (Table S3). When the final 5 years of the data series were arbitrarily excluded, the observed GPP trend was no longer statistically significant and had decreased to 11.6 $\mathrm{gC}{\mathrm{m}}^{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. This was comparable to the QUINCY estimation for the same period, which was 8.6 $\mathrm{gC}{\mathrm{m}}^{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$.

The observed TER in the period from 1996 to 2018 showed a non-significant trend of 9.0 $\mathrm{gC}{\mathrm{m}}^{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, which was a bit larger than the significant trend seen in TER from QUINCY, 6.9 $\mathrm{gC}{\mathrm{m}}^{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. The observations indicated a significant trend for NEE towards a larger sink, with a rate of −13.4 $\mathrm{gC}{\mathrm{m}}^{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ over the same time period. The simulations instead proposed an increasing source, with a magnitude of 5.2 $\mathrm{gC}{\mathrm{m}}^{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. The continuous observations of summertime LAI (averaged over June, July and August) showed a significant trend of 0.033 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Fig. S16d). The QUINCY model estimated a minor non-significant negative trend (−0.001 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). Tuning done for LAI and Chl_Leaf did not influence the simulated trends.

The GPP estimates derived from the CN simulation were found to be in close agreement with the observations until 2010 (Fig. S16a). However, a divergence is observed between the simulations and observations from that point onwards. The observations indicate an increase, whereas the simulations remain at a consistent level before declining in the last 3 years (Fig. S16a). In the observations the interannual variation (in terms of standard deviation) is found to be in similar magnitude for both TER (220 gC m² yr⁻¹) and GPP (224 gC m² yr⁻¹) (Table 3). In contrast, the simulations indicate that the interannual variation is greater in GPP (163 gC m² yr⁻¹) than TER (100 gC m² yr⁻¹) (Table 3). The nitrogen deposition that was input to the model showed a continuous decline in the study period (Fig. S17).

4 Discussion

In this study we employed a variety of observational data sources in conjunction with a terrestrial biosphere model. Our objective was to assess the utility of these data in enhancing and evaluating the model's performance, as well as to ascertain the model's capacity to simulate the biogeochemical cycles within Borden Forest. Our focus extends beyond the carbon cycle, as our modelling approach incorporates the nitrogen cycle, for which the available observational data provide valuable insights.

4.1 Using the continuous LAI observations

The site-level continuous LAI observations provide a valuable data source for model evaluation and development. Our simulations showed a discrepancy between the simulated absolute value of LAI and these observations. The LAI estimated from litterfall at the site was 5.1 m² m⁻² (Neumann et al., 1989) and in close agreement with the model estimate of 5.3 m² m⁻² derived from the C simulations and 4.8 m² m⁻² derived from the CN simulations. The changes that the biochemical model parameters cause to GPP are more influential than LAI values alone. The winter LAI will remain overestimated in our current modelling setup because only deciduous forest is being simulated. Changing the allocation pattern in QUINCY would have allowed for a reduction in simulated LAI, but this would have resulted in lower GPP values.

Tuning of the senescence parameter (Table 2) resulted in a substantially higher temperature threshold for the start of the senescence period at the site compared to the standard parameterization. As described in Sect. 2.3.2, the phenology model of QUINCY is a simple growing-degree-based formulation and does not take into account other environmental conditions such as light availability or day length, which might be contributing to early leaf senescence. Day length in particular has been shown to be another important variable controlling senescence in these regions (Bowling et al., 2024). Future work will address whether implementing such a dependency in QUINCY improves the seasonality of LAI estimation without adjusting the temperature threshold for senescence.

The continuous measurements are of particular value for assessing model seasonality, as they provide continuous data, unlike the point values obtained from LAI-2000 or litterfall. Therefore, using these data to improve the seasonality of the modelled LAI is a logical approach. By adjusting the senescence in accordance with the LAI observations, it is possible to improve the seasonal cycle of GPP, as the senescence was occurring too late in the default model at this site (Fig. 1). In this context, the continuous measurement of LAI is of particular significance, as it represents an independent measure from the carbon fluxes that can also be observed from space.

The parameterization of QUINCY often relies on multiple sites in order to perform successful large-scale simulations. An ongoing study evaluating the QUINCY seasonal cycles with flux tower and remote sensing data (Tuuli Miinalainen, personal communication, 2024) has shown that a similar bias in autumn phenology for temperate broadleaf deciduous trees occurs at several other sites, so this PFT could benefit from parameter tuning. However, the satellite observations have some issues in predicting the autumn phenology (Wang et al., 2024), and the long-term in situ observations at the Borden site can serve as a valuable verification resource.

4.2 Other leaf-level observations

In addition to LAI, there were several other leaf-level observations at the site, which we used in our model evaluation. In the original QUINCY formulation the development of leaf chlorophyll is fully coupled to the LAI development (Thum et al., 2019). However, observations reveal a clear decoupling taking place during the early season (Croft et al., 2017). The approach we have adopted here is to delay the development of leaf chlorophyll from the structural, i.e. non-photosynthetic, nitrogen. Previous studies at the site indicated that leaf chlorophyll lags 30 d behind LAI in reaching its maximum summertime value (Croft et al., 2015). However, the difference observed in the modelling was less pronounced, as the objective was to accurately represent the seasonal development of leaf chlorophyll correctly, while the observed maximum value of Chl_Leaf was reached later than the simulated summertime maximum level (Fig. 1c). The observed Chl_Leaf values exhibit elevated levels during the midsummer period, spanning approximately 1 month. In the modelling conducted, the objective was not to replicate this specific behaviour but rather to capture the average summertime level. The observed midsummer peak in Chl_Leaf was not simultaneously accompanied by increases in the biochemical model parameters V_c(max),25 and J_max,25 (Fig. 3), which would have a stronger effect on the simulated photosynthesis than Chl_Leaf alone. We could test the effect of increased midsummer Chl_Leaf on photosynthesis using the canopy module of QUINCY, which calculates only the canopy part of the model. Given the structure of the model, this effect is likely to be small. Luo et al. (2018) used the fact that Chl_Leaf and V_c(max),25 have a linear relationship at the Borden site to improve their model results for GPP and evapotranspiration. In QUINCY we currently do not have this linear relationship, although it has been observed in some studies (e.g. Qian et al., 2021). A further study with remote sensing data and more sites will explore the implementation of such a description of leaf N partitioning in QUINCY. At the Borden site the linear relationship between Chl_Leaf and V_c(max),25 was found to be off due to their different rates of light acclimation (Yu et al., 2024). Once we have a linear relationship between Chl_Leaf and V_c(max),25 implemented in QUINCY, we can use the canopy module to estimate what effect this would have on photosynthesis in our model. Currently, the way we have described the decoupling of LAI and Chl_Leaf in QUINCY does not lead to pronounced differences in the annual C balance (Table 3).

Another variable that undergoes changes during the spring season in the observations is the specific leaf area (SLA). Rather than delaying the onset of chlorophyll development from the structural nitrogen, it would be possible to introduce a dynamically changing SLA to QUINCY, as the leaves are thinner immediately following budburst. The SLA in QUINCY is set to a PFT-specific, time-invariant constant representative of the average observed value in Kattge et al. (2011), which is almost double the observed summertime value at Borden Forest. This overestimation may contribute to the large LAI predicted for the site by QUINCY, as the SLA is used in the leaf area calculation (Thum et al., 2019). Testing a dynamically changing SLA in the model falls outside the scope of our current study but is an important future step for improving predictions of leaf process seasonality.

4.3 Limitations of the model

It should be noted that, like with any other modelling study, our modelling approach is subject to certain limitations, due to the necessity of making certain simplifications. Borden Forest is a mixed forest, and the different deciduous species differ in their leaf chlorophyll contents and SLA values (Croft et al., 2017) and likely in their responses to climate variability. Since we do not have the ability to model different tree species, we model a deciduous forest composed of trees with identical traits. The rationale behind our approach can be attributed to the large scale that we are aiming to model and the focus on the utilization of parameters that can be derived from remote sensing observations which are unavailable at a scale where individual trees (and their species) can be resolved.

Another challenge in characterizing the forest as a deciduous PFT is the exclusion of the coniferous trees at the site. This leads to discrepancies in our wintertime estimates of LAI, as the simulated deciduous trees lack leaves during this season. Another effect can be seen in the delay of the simulated GPP increase in the spring, which occurs in the observations after DOY 100 (Fig. 1a, b). This is partly due to the understorey vegetation, which we do not simulate, and partly due to the coniferous trees at the site, which can start photosynthesis as soon as meteorological conditions and the release of possible winter acclimation allow. The shape of the seasonal cycle of GPP is also different in the simulations compared to the observations (Fig. 1b). The increase in the simulations is more abrupt to the summer levels and decline from early-summer values occurs quite early, probably due to dry periods occurring during summer. This could be due to the fact that the observed transitions are more smooth in time as several tree species contribute to the trend, which is represented only by one functional type in the model. Issues with species mixtures are common in TBMs, and while this is certainly an area that needs further improvement, it is not an issue unique to QUINCY.

Testing the performance of a TBM designed for large-scale simulation at the site level is challenging. The model must necessarily make generalizations in the process representation due to the limited knowledge and data required for large-scale parameterization. Of particular importance to this study is the parameterization of the partitioning of nitrogen into different compartments, which has a theoretical basis (Evans, 1989). However, the approach applied in QUINCY is simple and relies on PFT-specific parameterization, with very limited consideration of leaf-level data. The dataset available for the Borden site is very valuable in this respect, as it allows for evaluating the division of nitrogen to different compartments. The parameterization used in QUINCY does not take all factors into account. For example, the phosphorus content has been shown to influence the relationship between V_c(max),25 and leaf nitrogen (Walker et al., 2014). Furthermore, the description of the canopy nitrogen gradient in QUINCY appears to be sound, but it may not account for all possible variation (Niinemets et al., 2015). Borden Forest is located in the temperate boreal forest ecotone, and many species are close to the limits of their temperature and moisture ranges (Froelich et al., 2015). The tree species composition has undergone changes at the site during our study period; e.g. the red maple was reported to have coverage of 36 % in 1995 (Lee et al., 1999) and 52 % in 2006 (Teklemariam et al., 2009). The impacts that these changes in the tree composition have on the carbon fluxes could be studied by a demographic model with sufficient granularity in the description of tree functional diversity (see Fisher et al., 2018, for a review).

4.4 Legacy effects of drought

Carbon fluxes at the site were strongly influenced by the 2007 drought, which also led to legacy effects visible in 2008 (Fig. 4). Often the effect of drought on GPP is stronger than on TER (Schwalm et al., 2010; Piao et al., 2019), but here a stronger effect on TER was observed. GPP can be reduced by drought through both physiological and structural effects (van der Molen et al., 2011). The inclusion of non-structural carbohydrate pools in models has shown improvements in the ability of models to capture drought responses (Jones et al., 2020).

No decrease in LAI was observed at the site in 2007 or 2008. One of the recognized mechanisms for legacy effects is that the drought-induced decrease in GPP can lead to a decrease in the carbohydrate pool and therefore influence the LAI development in the following year (Yu et al., 2022). QUINCY has an explicit reserve pool and could theoretically simulate this type of behaviour but did not suggest that such legacy effects affected GPP in 2008. Although the simulated soil moisture shows a dynamic behaviour similar to that of the observations (Fig. 4), simulated topsoil moisture did not show the same dynamic range of the seasonal cycle of soil moisture in QUINCY when compared to the in situ observations, possibly explaining the underestimated response to soil moisture stress by the model. However, we note that the modelled range was more similar to the observations at deeper depths than to those at the top layer (Fig. S8e, f). At present, it is unclear whether this reflects shortcomings in the representation of soil physical processes or if it is a result of the lack of in situ precipitation observations during the study period (see “Materials and methods”).

There are many possible explanations for the strong legacy effect observed in TER in 2008. Soil microbial activity is dependent on soil moisture (Gaumont-Guay et al., 2006; Liu et al., 2009; Orchard and Cook, 1983), and drought can thus strongly reduce soil respiration directly and indirectly through several different mechanisms (von Buttlar et al., 2018). Direct effects include the dependence on the presence of water films for substrate diffusion and exo-enzyme activity (Davidson and Janssens, 2006) as well as microbial dormancy and death (Orchard and Cook, 1983). Indirect effects affect microbial activity through, for example, changes in soil nutrient retention and availability (Bloor and Bardgett, 2012) or changes in microbial community structure (Frank et al., 2015). In addition, GPP and TER fluxes are tightly coupled, as heterotrophic respiration is driven also by the recently assimilated carbon and not only by environmental conditions (Ruehr et al., 2012). QUINCY was not able to capture the drought-induced decrease in TER in 2007, which could either be due to too low an impact of the drought on soil moisture or the fact that the version of QUINCY applied here does not simulate microbial activity and root exudation (Yu et al., 2020). This also contributes to the failure to simulate potential legacy effects in the observed TER (Fig. 4).

The drought response at the site could potentially be improved by calibrating the soil moisture response functions in the model, but probably some structural changes in the description of soil physics would also be required. These changes may include modifications to the water-retention curve, pedotransfer functions (Weber et al., 2024) or infiltration properties (Vereecken et al., 2019). Soil moisture is a challenge for many models and often in need of improvement (De Pue et al., 2023). Future research will investigate whether a more sophisticated soil biogeochemical model can better represent the effects on microbial communities and through them the legacy effects on respiration (Yu et al., 2020). Overall, terrestrial biosphere models are not yet well equipped to capture the legacy effects (Bastos et al., 2021), and more work is needed to better understand the processes governing ecosystem recovery in order to improve models in this respect.

4.5 Seasonality in total ecosystem respiration

QUINCY is generally capable of modelling the observed magnitude and seasonal amplitude of observed TER, based on empirical responses of soil organic turnover to soil temperature and moisture (Table 3). The premature increase in the simulated soil respiration occurs during average years due to too early an increase in the soil temperature in the simulations (Fig. 2). This behaviour suggests that the coupling between the atmosphere and the soil in the model is too strong, which may be associated with parameters controlling heat diffusion in the soil.

One interesting feature observed in TER is the strong seasonality with respect to soil temperature in normal years, which cannot be explained by the soil moisture (Fig. S14). This behaviour was first discussed by Lee et al. (1999) at this site, and they called it the hysteresis effect. It is most pronounced in the year following the drought, 2008 (Fig. S14e). Based on the data available to Lee et al. (1999), they found that the early season had lower respiration values than the late season and speculated that this difference might be due to warmer soil temperatures in deeper soil layers as the season progressed, as well as greater litter accumulation. The data available to date for a longer period showed an “inverse” hysteresis effect, in which the later season has lower TER than the early season (Fig. S14a and e). Soil moisture does not provide an explanation for such a shift. Rather, this behaviour could be driven by the seasonality of photosynthesis, as the early-season GPP coincides with lower TER values compared to the late season (Fig. S15). A fast coupling between GPP and soil respiration, e.g. through photosynthesis supplying carbohydrates to rhizosphere respiration (Zhang et al., 2018), observed already early in a girdling experiment (Högberg et al., 2001), could explain the observed hysteresis effect. Hysteresis effects on soil respiration versus soil temperature are quite common (Zhang et al., 2018), and further explanations for dynamics taking place at the Borden site could be caused by substrate depletion late in summer (Kirschbaum, 2006) or by greater root productivity in the early season (Oe et al., 2011). Models generally describe soil respiration as a function of soil temperature responses and would not capture the hysteresis effects (Zhang et al., 2018). This is also the case for QUINCY, which is not able to capture this hysteresis effect (Fig. S14). Future work should evaluate whether including for instance increased vegetation–soil coupling via root exudates would improve the representation of the interannual variability in TER.

4.6 Trends and growing-season length

There were no significant changes in the growing-season metrics SOS, EOS and LOS over the period studied, similar to Gonsamo et al. (2015). QUINCY was generally successful in simulating these metrics, but the end of the season as estimated from GPP and LAI differed in the observations, whereas it was coupled in the model. The previous studies at the site that have assessed the growing-season metrics (Froelich et al., 2015; Gonsamo et al., 2015) also used the carbon uptake period (CUP). We did not assess CUP because its onset would have been biased in the simulations due to the premature increase in heterotrophic respiration caused by the premature increase in soil temperature (Fig. 2).

Froelich et al. (2015) found a significant increase in summertime GPP, and Gonsamo et al. (2015) found a significant increase in carbon uptake between 1996 and 2012 at the Borden site. These are consistent with our observational results, which additionally also showed a small but significant increase in the summertime LAI. The increase in net carbon uptake is attributed to increased PAR (photosynthetically active radiation; 400–700 nm), which leads to increased photosynthetic activity at this site (Gonsamo et al., 2015). When we looked at these PAR data in more detail, we noticed a significant increasing trend that had been reported. However, its magnitude was less than 1 % yr⁻¹, which would not lead to large increases in annual GPP according to our model. Also, when the time period was extended to 2018, the significant trend disappeared. The input to QUINCY was shortwave radiation, and it did not show an increasing trend for the time period, instead showing a very small decrease in annual values.

There have been reductions in atmospheric sulfur (S), nitrogen oxides and total nitrates and ozone deposition since 1992 at the Borden site (as measured at the nearby Egbert station), and the brightening seen in the PAR observations has been attributed to reductions in gaseous and particulate emissions, while declines in ozone emissions reduce the damages to the leaves (Gonsamo et al., 2015). The N and S deposition was found to reduce net primary production in southern Ontario (Aherne and Posch, 2013), where Borden Forest is also located. Decreasing S deposition was found to increase forest productivity in the northeastern US region (Dalton et al., 2024; Phelan et al., 2024), which is close to Borden Forest. Decreasing S and N deposition will also affect the soil pH (Dalton et al., 2024), which can potentially affect the N dynamics. The N deposition may be useful for forests under N deficit (Horn et al., 2018); however according to our model results, Borden Forest would not suffer from N deficiency. Recovery from S deposition may be one of the reasons for the increasing trend in GPP. Long-term S exposure could also reduce the drought resilience of sites (Dalton et al., 2024), which may be another explanation for why our model failed to capture the drought effects at this site.

Furthermore, the QUINCY model does not take into account potential damage to the leaves caused by ozone. Ozone influences both photosynthesis and stomatal conductance and can cause them to become decoupled (Novak et al., 2005; Lombardozzi et al., 2015). Decreases in photosynthesis by 21 % and in stomatal conductance by 11 % after chronic ozone exposure have been estimated (Lombardozzi et al., 2013), but also lower estimates have been presented (6 %–10 % for Europe) (Franz et al., 2017). Ozone exposure can also have an impact on the N cycle (Simpson et al., 2014). The impact of ozone has been modelled by direct influence on V_c(max),25 and stomatal conductance (Lombardozzi et al., 2012) or on photosynthesis, which then has feedbacks on stomatal conductance (Franz et al., 2017; Lombardozzi et al., 2015). Decreasing ozone deposition at the site could further affect the increase in photosynthesis.

QUINCY approximately reproduces the leaf-level photosynthetic parameters in the CN version of the model but at the same time overestimates LAI (compared to the continuous observations) and underestimates GPP. The uncertainty in the annual flux estimates by the eddy covariance method is usually around 10 %–20 % (Loescher et al., 2006), so QUINCY's estimates are within the uncertainty in the observations. Another source of uncertainty in the observations is the use of gap-filling (Mahabbati et al., 2021) and partitioning methods (Desai et al., 2008), which introduce uncertainty in the annual carbon balance estimates.

Possible reasons for the discrepancy between QUINCY and the observations are the lack of understorey representation in the model; the simplified representation of the mixed forest with a single deciduous plant functional type; and possible biases introduced by the assumed within-canopy gradient of leaf nitrogen, which might not hold for this diverse forest. It is interesting to note that the differences between the measured and simulated annual GPP are not apparent in the early years of the record but emerge in the later years (Fig. S16). The underestimation of annual GPP is 13 % by the CN:LAI&chl simulations for the whole period but only 8 % for the years 1996–2010. Therefore, the increasing trend in observed GPP that the model fails to reproduce contributes strongly to the model–data discrepancy.

4.7 Impact of the nitrogen cycle on carbon fluxes

Including the nitrogen cycle in the simulations did not cause a change in the net carbon balance of the ecosystem, as both GPP and TER were attenuated by approximately the same amount, about 20 % when comparing to the C version with N saturation. With the fixed stoichiometry the C model gave values similar to the CN simulations (Table 3). This means that there was no nitrogen constraint on the carbon fluxes according to the QUINCY model. The leaf C:N was only 9.7 in the C:LAI&chl-saturated case simulation, which is an unrealistically low value, but we chose to show these results here to assess the influence of parameterization with simulations with magnitudes comparable to the observed GPP. The foliar C:N ratio was 20.4 in the CN:LAI&chl simulation and 22.4 in Cfix:LAI&chl, showing that they are similar.

Nitrogen availability limits the carbon cycle (Du et al., 2020), especially in the boreal region (Högberg et al., 2017). The estimated effect of the N cycle on the carbon fluxes is not as large as some previous estimates (Thornton et al., 2007). The summertime variation in GPP values is more pronounced in the C simulations with dynamic stoichiometry than in the CN simulations, highlighting the more stable behaviour of the model when including the N constraint and compared to the N-saturated case. Although the annual GPP was underestimated with the inclusion of the N cycle, after tuning by both LAI and Chl_Leaf, the CN simulation gave the best r² and RMSE metrics for both GPP and TER, in line with the C cycle simulations with fixed stoichiometry (Table 3).

4.8 Using site-level observations in model development

The different site-level observations available at the site provided means to evaluate the model performance from different aspects. The QUINCY model is a large-scale model and cannot capture all the small-scale variations. Furthermore, the different tree species at the Borden site complicate simulating the forest with QUINCY, as the model needs to rely on a general PFT description. However, to better understand the processes occurring at the site some further observations would be useful. To capture the forest structure and to facilitate estimation of the radiative transfer inside the canopy, lidar observations (Balestra et al., 2024) would be beneficial, and if done on a temporally continuous scale (such as in StrucNet; see Calders et al., 2023), valuable information on the allocation of annual net primary production could also be obtained. Soil chamber observations of respiration would help to separate the role of soil in total ecosystem respiration. The use of isotopes would enable revealing the processes behind the observed hysteresis behaviour of the soil respiration. A rain gauge at the site would help to study potential biases of using precipitation data at a nearby site. The ozone profile concentration observation at the site could help in estimating the potential ozone damage to vegetation that could be addressed by a model. Observations of leaf and soil C:N ratios would help to better understand the nitrogen status of the forests.

4.9 Revisiting the research questions and future responses

Our first research question was whether the decoupling of LAI and Chl_Leaf affects the estimation of annual carbon fluxes. Delaying the development of Chl_Leaf did affect the seasonal development of GPP, but the effect on annual fluxes was very small. If the Chl_Leaf was more linearly related to the biochemical model parameters in the model, the influence might be somewhat greater.

Our second research question was whether QUINCY is able to simulate any long-term changes in seasonal shifts in carbon fluxes and LAI values. The observed trends in GPP found at the site were quite large and were not captured by our model, even though it responds to climate and increasing atmospheric [CO₂]. Possible reasons for this could be processes that are not included in our model, such as the role of the understory, changes in tree composition, and recovery from S and ozone deposition. The observations showed a small increasing trend in LAI that was not captured by the model. If the reason for the increase is a change in tree composition, our model would not include the process behind it.

Our third research question was whether nitrogen limitation exists in Borden Forest and whether it changes over the 22-year period. To answer this question, we compared the nitrogen-cycle-enabled version to the carbon-cycle-only version with fixed stoichiometry. The differences in the annual carbon fluxes between these two simulations were very small, and there was no change in time between these results. Therefore Borden Forest would not be nitrogen limited. This is consistent with the finding that observed GPP continues to increase, even though the N deposition at the site has decreased since 1995.

Our fourth research question was whether QUINCY could simulate the effects of drought events on the carbon cycle. In general, the model captures many environmental responses that occur at the site, as the model performance is generally good. However, the effects of the 2007 drought on carbon fluxes in that year and the following year were not captured by the model. The top layer of the soil contains most of the organic matter, and most of the heterotrophic respiration in the model comes from this layer. The soil model does not have a wide enough range of soil moisture values in the top layer compared to observations. This may be one reason for the discrepancy. QUINCY has non-structural carbohydrate pools, and they get lower values in 2007 and also in 2008 in the reserve pool compared to regular years. These changes are not reflected in the summer LAI values, which would be one way to carry out the legacy effect of the drought into other years, but this was also not detected in the observations. The legacy effect of respiration is likely due to a process not currently represented in QUINCY.

QUINCY was generally able to capture the environmental responses of the forest. However, the long-term trends and extreme events pose challenges to the model that would require a more detailed representation of some processes. These challenges highlight the value of long-term data series. The combination of data and simulations shows a decoupling of LAI development and Chl_Leaf as well as too strong a coupling of LAI and GPP in the model in autumn compared to the observations.

In the future Borden Forest is likely to have warmer temperatures and more precipitation in winter (Bush et al., 2022). A warmer climate is expected to increase the severity of heat waves, which in turn may contribute to increased droughts (Bush et al., 2022). The 2007 drought likely had a legacy effect on the 2008 component fluxes. However, because TER was reduced more than GPP in 2007, the forest was a net sink (−318 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) in that year and also in the following year (−282 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). Thus, the drought did not endanger the sink capacity of the forest but rather increased it. If droughts become more frequent, they could threaten the sink capacity of the forest. Estimating the effect of warmer temperatures is challenging because the observed trend in GPP is not captured by the model. QUINCY would predict both increased photosynthesis and respiration at warmer temperatures, but whether this would lead to a higher sink or source would depend on the overall conditions.

5 Conclusions

In this work we used several data streams measured at the Borden Forest Research Station, some extending over 2 decades, and aimed to improve and evaluate the terrestrial biosphere model QUINCY. This work demonstrated the usefulness of using different data sources and the importance of observational long time series. The use of leaf chlorophyll content and LAI in parameterizing the model improved simulated GPP in the CN simulations. These changes also decreased the RMSE for TER. Generally the model did capture the average seasonal cycle of GPP (daily r²=0.80) and TER (daily r²=0.75). QUINCY was also successful in its estimation of the growing-season metrics, even though the ending of the season was more coupled between LAI and GPP than in the observations.

Two important data sources used in this work, leaf chlorophyll content and the leaf area index (LAI), can also be measured from space. Therefore our work paves the way toward combining terrestrial biosphere models with remote sensing data for parameterization, as has been proposed by Rogers et al. (2017). Work on this front has been done by combining leaf chlorophyll with the photosynthesis parameters of models (Lu et al., 2022). In this work we explicitly model the leaf chlorophyll, which links this variable directly to the nitrogen cycle. The unique dataset accessible from the Borden site permitted the assessment and enhancement of the parameterization employed to divide leaf nitrogen to different compartments. In addition to utilizing LAI and leaf chlorophyll, sun-induced chlorophyll fluorescence (SIF) represents a pivotal variable observed from space that is linked to the carbon cycle (Sun et al., 2023). SIF is currently being implemented in QUINCY and will, in the future, provide a means of conducting global-scale assessments of the carbon cycle together with nitrogen-cycle-related metrics.