Exploring alternative SMAP Level-4 carbon model formulations for the North American Arctic–Subarctic growing season

The Soil Moisture Active Passive Level-4 Terrestrial Carbon Flux model (hereafter referred to as the L4C model) provides daily estimates of net ecosystem CO₂ exchange (NEE), gross primary production (GPP), and ecosystem respiration (ER) at a global scale. The model is based on direct mechanistic forcing–response relationships between CO₂ fluxes and energy proxies (absorbed photosynthetically active radiation and temperature) and moisture proxies (soil moisture and vapor pressure deficit). Although the L4C model aims to provide a representative estimation of the CO₂ budget of Arctic and Subarctic (AS) environments, a deeper understanding of carbon cycle processes and targeted refinements are needed to improve its accuracy. In this study, alternative model formulations are proposed for the North American AS regions during the growing season. These formulations are calibrated and evaluated using NEE-derived GPP and ER from 20 eddy covariance towers across western Canada and Alaska, covering the period from 2015 to 2022. Refinements in the representation of energy proxies resulted in greater improvements in model performance than adjustments to moisture proxies. Specifically, implementing a light-response curve in GPP estimation reduced unbiased root mean squared error and bias, while incorporating growing degree days improved correlation. Adjustments to rootzone and surface soil moisture in GPP and ER estimation, respectively, did not yield conclusive performance improvements. Vapor pressure deficit showed limited importance as a driver of GPP in upland tundra and wetlands, whereas it had a stronger impact in taiga forests. Finally, the litterfall scheme used to represent SOC dynamics in the L4C ER model formulation in version 8 demonstrated improved performance relative to version 7. Although some adjustments in ER and GPP formulations yielded strong performance gains, improvements in NEE were more modest than for the individual components. Overall, the results highlight opportunities to enhance the accuracy of the L4C model for the North American AS growing season and underscore the need for further research on CO₂ flux modeling.

1 Introduction

Arctic and Subarctic (AS) environments store nearly half of the global soil organic carbon (SOC) pool (Tarnocai et al., 2009; Hugelius et al., 2014; Mishra et al., 2021) and are experiencing accelerated warming (Rantanen et al., 2022). Rising temperatures increase photosynthetic activity and extend the growing season, leading to higher CO₂ uptake by vegetation (Myneni et al., 1997; Jia et al., 2003; Euskirchen, E. S. et al., 2009; Natali et al., 2012; Forkel et al., 2016; Fisher et al., 2018). In addition, rising temperatures enhance autotrophic respiration (AR) as well as heterotrophic respiration (HR) in two pathways: directly, by stimulating microbial activity, and indirectly, by thawing permafrost and exposing frozen SOC to decomposition. The combined increase in AR and HR intensifies CO₂ release to the atmosphere (Natali et al., 2019; Turetsky et al., 2020; Virkkala et al., 2025). Consequently, estimating the net CO₂ budget of AS regions is essential for understanding their role in global climate system feedbacks (Oechel et al., 1993; Hayes et al., 2011; Turetsky et al., 2011; Bell et al., 2013; Schuur et al., 2013; Schaefer et al., 2014; Zona et al., 2016). Nevertheless, our understanding of CO₂ fluxes in the AS environments remains limited. This is due to the inherent complexity and high cost of measuring CO₂ fluxes, the scarcity of such measurements, and the seasonal variability in the dominant processes controlling CO₂ fluxes (Baldocchi et al., 2007; Fisher et al., 2018; Pallandt et al., 2022; Mavrovic et al., 2023b).

Net ecosystem CO₂ exchange (NEE) represents the overall balance between CO₂ uptake by photosynthesis, called gross primary production (GPP), and CO₂ release through ecosystem respiration (ER), as follows (Chapin et al., 2006)

$$\begin{array}{}\text{(1)}& \text{NEE}=\text{HR}+\text{AR}-\text{GPP}=\text{ER}-\text{GPP}\end{array}$$

GPP is a light-driven process whose efficiency is modulated by air temperature, soil moisture availability within the plant root zone and vapor pressure deficit, which can induce stomatal closure and thereby reduce CO₂ uptake (Davis et al., 2014; Bao et al., 2022). Comparatively, HR is governed by SOC availability, soil temperature, and surface soil moisture, whereas AR primarily depends on air temperature, plant metabolic activity and GPP rate (Reichstein et al., 2005; Davis et al., 2014; Zona et al., 2023). When soil temperature drops near 0 °C, the soil starts freezing and GPP and AR progressively ceases, following a soil freezing characteristic curve (Salmabadi et al., 2025). Under fully frozen conditions, NEE is equal to HR, which is controlled by soil temperature and SOC availability (Natali et al., 2019; Mavrovic et al., 2023b).

Although global terrestrial carbon flux (TCF) models, atmospheric inversions (which infer surface CO₂ fluxes from atmospheric CO₂ concentrations), and data-driven flux-upscaling approaches are available to estimate the CO₂ budget of AS regions, they often disagree on whether these regions are CO₂ sources or sinks (McGuire et al., 2012; Fisher et al., 2018; López-Blanco et al., 2019; Virkkala et al., 2021; Ramage et al., 2024; Virkkala et al., 2025; Foster et al., 2024). In recent decades, satellite-based microwave remote sensing (300 MHz–100 GHz) has provided a valuable approach for monitoring land–atmosphere interactions and carbon cycle dynamics through the retrieval of key surface variables (Fisher et al., 2018; Lees et al., 2018; Mavrovic et al., 2023a; Pulliainen et al., 2024) such as soil moisture (Kerr et al., 2012; Colliander et al., 2017), snow properties (Lievens et al., 2019), aboveground biomass (Mialon et al., 2020), and freeze-thaw state (Rautiainen et al., 2016; Derksen et al., 2017; Prince et al., 2019). In 2015, the Soil Moisture Active Passive (SMAP) satellite was launched to monitor surface soil moisture and freeze-thaw dynamics using L-band brightness temperature observations (Entekhabi et al., 2010). One of the science objectives of the SMAP mission is to improve our understanding of interconnected water, energy, and carbon cycles, as well as to quantify the boreal landscape CO₂ budget (Entekhabi et al., 2014). In this context, the SMAP Level-4 Global Daily 9 km EASE-Grid Carbon Net Ecosystem Exchange (SPL4CMDL) product currently provides global, daily estimates of NEE and GPP, as well as ER (indirectly derived from HR, GPP and NEE). These estimates are derived from a TCF model (hereafter referred to as the L4C model), which is notably informed by the SMAP Level-4 Global 9 km EASE-Grid Surface and Root Zone Soil Moisture Geophysical Data (SPL4SMGP) product (Jones et al., 2017; Endsley et al., 2022; Kimball et al., 2025; Reichle et al., 2025a). Although the L4C model achieves an unbiased root-mean-square error of NEE within the targeted accuracy of 1.6 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in AS environments, recent studies have reported that it fails to capture the amplitude of GPP and ER during the green-up phase (Endsley et al., 2022; Madelon et al., 2025). The authors also reported discrepancies in annual CO₂ budgets when compared with eddy covariance (EC) measurements, leading to uncertainties in classifying AS environments as net CO₂ sources or sinks (Madelon et al., 2025). From April to July 2025, the L4C model transitioned from version 7 to version 8 (Kimball et al., 2025), featuring a major update partly due to (i) the upgrade of the SMAP SPL4SMGP product, which transitioned from its own version 7 to version 8 (Reichle et al., 2025a), and (ii) changes to the litterfall estimation scheme used for modeling SOC dynamics and ER (Sect. 3).

The goal of this study is to better characterize how key environmental drivers influence GPP and ER, and to refine their modeling for the North American AS growing season. To achieve this, we:

• explore alternative formulations of the L4C model (hereafter referred to as the AS-adapted models) that adjust GPP and ER responses to absorbed photosynthetically active radiation, air and soil temperature, rootzone and surface soil moisture, and vapor pressure deficit;

• calibrate and evaluate these formulations using GPP and ER data from 20 EC towers across western Canada and Alaska from 2015 to 2022;

• identify and interpret the specific model adjustments that yield the greatest performance improvements in terms of Pearson correlation, unbiased root-mean-square error, and bias;

• provide recommendations for more accurate satellite-derived estimates of GPP and ER, with indirect benefits for NEE and CO₂ budget estimation.

2 Eddy covariance measurements

NEE, GPP, and ER data were collected from 20 EC towers located in AS environments (Fig. 1), all within the NASA Arctic-Boreal Vulnerability Experiment (ABoVE) study domain (https://above.nasa.gov/sites.html, last access: 15 January 2026). The dataset spans from April 2015 through December 2022 (Table 1) and includes

• half-hourly fluxes from 13 EC towers in Alaska, downloaded from the AmeriFlux Network website (https://ameriflux.lbl.gov, last access: 21 August 2025).

• half-hourly fluxes from 7 EC towers in western Canada, provided directly by the principal investigators to ensure the use of the most up-to-date records; some of these sites are not yet available on the AmeriFlux Network website.

EC towers measure NEE by quantifying the turbulent vertical exchange of CO₂ in the surface layer of the atmospheric boundary layer (typically within the lowest tens of meters), where turbulence dominates the airflow (Aubinet et al., 2012; Burba, 2013, 2022). The spatial footprint can extend up to 1 km or more, but remains complex to estimate and varies with wind direction, wind speed, and tower height (Leclerc and Thurtell, 1990; Schuepp et al., 1990; Aubinet et al., 2012; Webb et al., 2016). GPP and ER are derived from NEE using flux-partitioning methods selected by the investigators at each EC tower. In this study, we placed confidence in these flux-partitioning methods and considered the resulting partitioned GPP and ER values to be credible representations of the underlying processes and suitable for use as reference data for model calibration and evaluation. Additional information on NEE, and the derived GPP and ER is provided in Appendix A. EC NEE, GPP, and ER were averaged from 30 min intervals to daily time steps, using at least 24 out of the theoretical 48 data points available per day. Hereafter, NEE_EC, GPP_EC, and ER_EC refer to the daily means of EC NEE, GPP, and ER.

Figure 1 Locations of the 20 eddy covariance (EC) towers providing measurements of net ecosystem CO₂ exchange (NEE_EC), NEE-derived gross primary production (GPP_EC) and ecosystem respiration (ER_EC) from April 2015 to December 2022. The High Arctic, Low Arctic and Subarctic zones were delineated following the Conservation of Arctic Flora and Fauna (CAFF) working group of the Arctic Council, using Moderate Resolution Imaging Spectroradiometer (MODIS) and Landsat imagery (Potapov et al., 2008). The permafrost extent is estimated in percent areal coverage (Brown et al., 2002): continuous (>90 %–100 % areal extent), discontinuous (>50 %–90 %), sporadic (10 %–50 %) and isolated patches (<10 %). Due to overlapping, a single color dot may represent up to four EC towers on the map. Information for each EC tower is listed in Table 1. The figure was inspired by Madelon et al. (2025) and Mavrovic et al. (2023a).

Madelon et al. (2025)Madelon et al. (2025)Sonnentag and Marsh (2021a)Sonnentag and Quinton (2021)Sonnentag and Quinton (2018)Sonnentag (2021)Sonnentag and Marsh (2021b)Euskirchen (2022d)Euskirchen (2022c)Euskirchen (2022b)Euskirchen (2022a)Bracho et al. (2021)Euskirchen et al. (2022a)Euskirchen et al. (2022b)Euskirchen et al. (2022c)Iwahana et al. (2023)Ueyama et al. (2023b)Ueyama et al. (2023a)Natali (2024)Natali (2025)

Table 1 List of the 20 eddy covariance (EC) tower sites used in this study. Ecosystem types were assigned based on site descriptions while plant functional type (PFT) classes were assigned using the Moderate Resolution Imaging Spectroradiometer (MODIS) MCD12Q1 type 5 product (Friedl and Sulla-Menashe, 2019). GRA, SHR, and ENF stand for grassland, shrubland, and evergreen needleleaf forest, respectively. Column 7 shows the number of data points (NDP) within the EC net ecosystem CO₂ exchange (NEE_EC), gross primary production (GPP_EC), and ecosystem respiration (ER_EC) time series after daily averaging and data filtering (see Sects. 2 and 4.2).

All EC towers considered in this study are located in the tundra and taiga biomes, in areas underlain by sporadic, discontinuous, or continuous permafrost (Fig. 1). Permafrost (ground that has remained frozen for more than 2 years) lies beneath an active layer that thaws during the growing season, enabling plant growth, as roots can only establish in thawed soil (Blume-Werry et al., 2019). Permafrost also restricts surface drainage, promoting water saturation and slowing decomposition rates (Wania et al., 2009; Robinson and Moore, 2000; Rouse et al., 1997; Maltby and Immirzi, 1993). These conditions can favor the formation of wetlands, including peatlands as well as other seasonally or permanently waterlogged ecosystems, across both tundra and taiga biomes (Treat et al., 2022). Based on site descriptions, EC towers were grouped into three distinct ecosystem types (Table 1):

• Taiga forests: 7 EC towers are located in taiga forests, characterized by a vertically stratified vegetation structure, with an open canopy of coniferous trees and an understory of shrubs, mosses, and lichens (Crawford, 2013; Juday, 2025).

• Upland tundra: 5 EC towers are located in upland tundra, which may exhibit lower vegetation density and diversity, as well as reduced soil biological activity, compared with taiga forests (Crawford, 2013; Hagedorn et al., 2025). The landscape is treeless and dominated by dwarf shrubs, grasses, sedges, mosses, and lichens, as plant growth is constrained by cold temperatures, short growing seasons, and the shallow depth of the permafrost active layer (Crawford, 2013; Hu and Bliss, 2025; Juday, 2025; Péwé, 2025).

• Wetlands: 8 EC towers are located in wetlands, where the term “wetland” refers to a wide range of types, including peatlands, bogs, fens, marshes, wet meadows, and shrub swamps, present in both tundra and taiga biomes. Compared with taiga forests and upland tundra, wetlands may exhibit higher species richness (McPartland et al., 2019) and localized microtopography, such as hummocks and hollows, whose characteristics depend on water table depth (Rouse et al., 1997; Zhang et al., 2024).

3 L4C model

The L4C model provides global, daily estimates of NEE, GPP, and ER at 9 km resolution from 31 March 2015, to the present (Jones et al., 2017; Kimball et al., 2025). It takes as inputs

• A static global plant functional type (PFT) classification at 500 m resolution, retrieved from the Moderate Resolution Imaging Spectroradiometer (MODIS) MCD12Q1 Type 5 product (Friedl and Sulla-Menashe, 2019).

• Eight-day fraction of photosynthetically active radiation (canopy-intercepted FPAR) and leaf area index (LAI) data at 500 m resolution, retrieved from the Visible Infrared Imaging Radiometer Suite (VIIRS) VNP15A2H product (Myneni and Knyazikhin, 2018).

• Daily means of three-hourly data at 9 km resolution, retrieved from the SMAP SPL4SMGP product version 8 (Reichle et al., 2025a), including 10 cm deep soil temperature (ST₁₀), surface skin temperature, incident shortwave solar radiation (SW_in), surface soil moisture (SSM), and rootzone soil moisture (RZSM). SSM and RZSM estimates are obtained by assimilating SMAP L-band brightness temperature observations into the Goddard Earth Observing System Version 5 Catchment Land Surface Model (GEOS-5 CLSM) (Reichle et al., 2019).

• Daily vapor pressure deficit (VPD) and minimum air temperature (MNT) at 0.25° (approximately 25 km resolution), retrieved from the GEOS-5 Forward Processing (FP) product (Lucchesi, 2018).

SW_in is combined with FPAR to compute canopy-absorbed photosynthetically active radiation (APAR), assuming that PAR constitutes 45 % of SW_in, as follows:

$$\begin{array}{}\text{(2)}& \text{APAR}=\mathrm{0.45}\cdot {\text{SW}}_{\text{in}}\cdot \text{FPAR}=\text{PAR}\cdot \text{FPAR}\end{array}$$

RZSM is rescaled using a normalized logarithmic transformation (Jones et al., 2017), and SSM is converted from volumetric units to relative wetness units. With the exception of APAR, the variables MNT, VPD, RZSM, ST, and SSM affect GPP and ER estimation only after being converted into stress scalars, denoted as S_MNT, S_VPD, S_RZSM, S_ST, and S_SSM. The derivation of these stress scalars is described later in Eq. (7a–e).

The L4C model runs at a daily time step and is defined as follows:

$$\begin{array}{ll}\text{(3a)}& {\displaystyle }\text{NEE}\left(t\right)& {\displaystyle }=\text{ER}\left(t\right)-\text{GPP}\left(t\right){\displaystyle }\text{GPP}\left(t\right)& {\displaystyle }={\mathit{ϵ}}_{\text{max}}\cdot \text{APAR}\left(t\right)\cdot S_{\text{MNT}}\left(t\right)\cdot S_{\text{VPD}}\left(t\right)\\ \text{(3b)}& {\displaystyle }& {\displaystyle }\cdot S_{\text{RZSM}}\left(t\right){\displaystyle }\text{ER}\left(t\right)& {\displaystyle }=\text{AR}\left(t\right)+\text{HR}\left(t\right)=\mathit{\alpha }\cdot \text{GPP}\left(t\right)+[k_{\mathrm{1}}\cdot {\text{SOC}}_{\mathrm{1}}(t)\\ {\displaystyle }& {\displaystyle }+(\mathrm{1}-\mathit{\eta })\cdot k_{\mathrm{2}}\cdot {\text{SOC}}_{\mathrm{2}}\left(t\right)+k_{\mathrm{3}}\cdot {\text{SOC}}_{\mathrm{3}}\left(t\right)]\\ \text{(3c)}& {\displaystyle }& {\displaystyle }\cdot S_{\text{ST}}\left(t\right)\cdot S_{\text{SSM}}\left(t\right)\end{array}$$

GPP is modeled using a light-use efficiency approach (Jones et al., 2017; Xiao et al., 2013), where ϵ_max represents the bulk environmental reduction in PAR conversion efficiency. AR is modeled as a fixed proportion of GPP, determined by the coefficient α. HR is estimated using a cascading three-pool SOC decomposition model (Ise and Moorcroft, 2006; Kimball et al., 2008; Jones et al., 2017), assuming that carbon fixed from atmospheric CO₂ through GPP enters the SOC pools as litterfall (L_fall). The daily SOC change for each of the three SOC pools is specified as:

$$\begin{array}{ll}{\displaystyle }{\text{SOC}}_{\mathrm{1}}\left(t\right)& {\displaystyle }={\text{SOC}}_{\mathrm{1}}(t-\mathrm{1})+[\mathit{\lambda }\cdot L_{\text{fall}}(t)-k_{\text{l}}\\ \text{(4a)}& {\displaystyle }& {\displaystyle }\cdot {\text{SOC}}_{\mathrm{1}}(t-\mathrm{1})\cdot S_{\text{ST}}\left(t\right)\cdot S_{\text{SSM}}\left(t\right)]\cdot \text{dt}{\displaystyle }{\text{SOC}}_{\mathrm{2}}\left(t\right)& {\displaystyle }={\text{SOC}}_{\mathrm{2}}(t-\mathrm{1})+\left[\right(\mathrm{1}-\mathit{\lambda })\cdot L_{\text{fall}}(t)-k_{\text{s}}\\ \text{(4b)}& {\displaystyle }& {\displaystyle }\cdot {\text{SOC}}_{\mathrm{2}}(t-\mathrm{1})\cdot S_{\text{ST}}\left(t\right)\cdot S_{\text{SSM}}\left(t\right)]\cdot \text{dt}{\displaystyle }{\text{SOC}}_{\mathrm{3}}\left(t\right)& {\displaystyle }={\text{SOC}}_{\mathrm{3}}(t-\mathrm{1})+[\mathit{\eta }\cdot k_{\text{s}}\cdot {\text{SOC}}_{\mathrm{2}}(t)\cdot S_{\text{ST}}(t)\\ {\displaystyle }& {\displaystyle }\cdot S_{\text{SSM}}\left(t\right)-k_{\text{r}}\cdot {\text{SOC}}_{\mathrm{3}}(t-\mathrm{1})\cdot S_{\text{ST}}\left(t\right)\\ \text{(4c)}& {\displaystyle }& {\displaystyle }\cdot S_{\text{SSM}}\left(t\right)]\cdot \text{dt}\end{array}$$

SOC₁, SOC₂, and SOC₃ represent the labile, structural, and recalcitrant SOC pools, respectively, with corresponding decay rates k₁, $k_{\mathrm{2}}=\mathrm{0.4}\cdot k_{\mathrm{1}}$, and $k_{\mathrm{3}}=\mathrm{0.01}\cdot k_{\mathrm{1}}$. The model parameters λ and η account for the fraction of L_fall allocated to the SOC₁ and SOC₂ pools, and the fraction of material transferred from the SOC₂ pool to the SOC₃ pool, respectively. The parameters ϵ_max, k_l, λ, and η are treated as free parameters estimated during the optimization process. The model integration step dt is set to one day.

In the L4C model version 7, L_fall was derived as a constant daily fraction of the mean annual estimated net primary productivity (NPP), as follows:

$$\begin{array}{}\text{(5)}& L_{\text{fall}}\left(t\right)=\stackrel{\mathrm{‾}}{{\text{NPP}}_{\text{annual}}}\cdot \mathrm{1}/\mathrm{365}=(\mathrm{1}-\mathit{\alpha })\cdot \stackrel{\mathrm{‾}}{{\text{GPP}}_{\text{annual}}}\cdot \mathrm{1}/\mathrm{365}\end{array}$$

$\stackrel{\mathrm{‾}}{{\text{NPP}}_{\text{annual}}}$ and $\stackrel{\mathrm{‾}}{{\text{GPP}}_{\text{annual}}}$ denote the mean annual NPP and GPP, respectively, and the daily fraction is set to $\mathrm{1}/\mathrm{366}$ for leap years. Instantaneous NPP is initially derived as $\text{NPP}\left(t\right)=\text{GPP}\left(t\right)-\text{AR}\left(t\right)$. In the L4C model version 8, the allocation timing was changed from constant to dynamic, and is now determined using a leaf-loss function (L_loss), derived from climatological LAI, as follows:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{L}_{\text{fall}}\left(t\right)& =\stackrel{\mathrm{‾}}{{\text{NPP}}_{\text{annual}}}\cdot \left[f_{\text{E}}\cdot \text{dt}+(\mathrm{1}-f_{\text{E}})\cdot {\displaystyle \frac{L_{\text{loss}}\left(t\right)}{\sum L_{\text{loss}}\left(t\right)}}\right]\\ & \text{with}{f}_{\text{E}}={\displaystyle \frac{min\left(\text{LAI}\right)}{max\left(\text{LAI}\right)}}\end{array}\end{array}$$

Here, f_E represents the proportion of the canopy that is evergreen. L_loss is computed using a triangular moving average centered on the current time step, where weights increase linearly toward the center. It represents the difference between lagged and leading climatological LAI values (Endsley et al., 2022).

In Eqs. (3a–b) and (4a–c), the stress scalars S_MNT, S_VPD, S_RZSM, S_ST, and S_SSM represent the ecosystem responses to their respective environmental variables and are defined as follows:

$$\begin{array}{}\text{(7a)}& {\displaystyle }{S}_{\text{MNT}}\left(t\right)=min\left(\mathrm{1},max\left(\mathrm{0},{\displaystyle \frac{\text{MNT}\left(t\right)-{\text{MNT}}_{min}}{{\text{MNT}}_{max}-{\text{MNT}}_{min}}}\right)\right)\text{(7b)}& {\displaystyle }{S}_{\text{VPD}}\left(t\right)=min\left(\mathrm{1},max\left(\mathrm{0},\mathrm{1}-{\displaystyle \frac{\text{VPD}\left(t\right)-{\text{VPD}}_{min}}{{\text{VPD}}_{max}-{\text{VPD}}_{min}}}\right)\right)\text{(7c)}& {\displaystyle }{S}_{\text{RZSM}}\left(t\right)=min\left(\mathrm{1},max\left(\mathrm{0},{\displaystyle \frac{\text{RZSM}\left(t\right)-{\text{RZSM}}_{min}}{{\text{RZSM}}_{max}-{\text{RZSM}}_{min}}}\right)\right)\text{(7d)}& {\displaystyle }{S}_{\text{SSM}}\left(t\right)=min\left(\mathrm{1},max\left(\mathrm{0},{\displaystyle \frac{\text{SSM}\left(t\right)-{\text{SSM}}_{min}}{{\text{SSM}}_{max}-{\text{SSM}}_{min}}}\right)\right)\text{(7e)}& {\displaystyle }{S}_{\text{ST}}\left(t\right)=min\left(\mathrm{1},\exp \left({\mathit{\beta }}_{\mathrm{0}}\cdot \left({\displaystyle \frac{\mathrm{1}}{{\mathit{\beta }}_{\mathrm{1}}}}-{\displaystyle \frac{\mathrm{1}}{\text{ST}\left(t\right)-{\mathit{\beta }}_{\mathrm{2}}}}\right)\right)\right)\end{array}$$

Each stress scalar ranges from 0 to 1, where a value of 0 indicates that the environmental variable fully constrains model estimates, while a value of 1 indicates no constraint. The thresholds MNT_min, MNT_max, VPD_min, VPD_max, RZSM_min, RZSM_max, SSM_min, and SSM_max are free parameters estimated during the optimization process. These are used as model thresholds and do not correspond to the actual minimum or maximum values within the time series. Similarly, β₀ is a free parameter, while β₁ and β₂ are fixed at 66.02 and 227.13 K, respectively (Kimball et al., 2025). The behavior of the ecosystem response functions is shown in Appendix B, Fig. B1(A1–C1). The GPP formulation originally includes a stress scalar based on the freeze–thaw state, computed using surface skin temperature from the SMAP SPL4SMGP product (Jones et al., 2017; Kimball et al., 2025). However, it is not shown here, as this study focuses on the growing season.

L4C model estimates are initially derived at a 1 km sub-grid resolution for up to 8 MODIS MCD12Q1 PFT classes, and then averaged to 9 km resolution. In this study, only the L4C model estimates corresponding to the PFT class in which the EC towers are located were considered (Table 1). A total of two towers are located in the grassland (GRA) class, 8 in the shrubland (SHR) class, and 10 in the evergreen needleleaf forest (ENF) class. Hereafter, NEE_L4C, GPP_L4C, and ER_L4C refer to daily estimates derived from the L4C model, which were retrieved from the SMAP SPL4CMDL product version 8. Frequently used abbreviations throughout this study are listed in Table 2.

Table 2 Summary of frequently used abbreviations.

4 Method

4.1 Arctic-Subarctic adapted model formulations

In this study, we explored alternative GPP and ER model formulations that retain the core GPP and ER equations of the original L4C model version 8. We aimed to preserve the structure and variable set of the original L4C model while enabling the incorporation of constraints or additional flexibility guided by literature-based findings on ecosystem responses and flux-partitioning methods. Five different formulations are presented for both GPP and ER, with each formulation building incrementally on the previous one by incorporating earlier modifications along with additional adjustments. Testing modifications incrementally, rather than independently, allowed us to determine whether their interactions improved or degraded model performance.

The GPP formulations, labeled GPP₁ through GPP₅, mainly adjust ecosystem responses to environmental variables and are defined as follows (Table 3):

• GPP₁: Under sub-freezing air temperatures, photosynthetic activity is expected to be severely reduced, approaching cessation (Schaefer et al., 2012; Ensminger et al., 2004; Bowling et al., 2018; Parazoo et al., 2018). This behavior is not well represented in the original L4C model, where S_MNT still remains near 0.5 at 270 K (Figs. 2, 3, 4(A2)), indicating that GPP capacity is reduced by only half at this temperature. In the proposed formulation, GPP is ensured to cease when MNT is equal to or below 273.15 K by fixing the minimum threshold (MNT_min) of S_MNT to 273.15 K (Eq. 7a).

• GPP₂ (defined as GPP₁ with additional adjustments): Some flux-partitioning methods use a nonlinear light-response curve to partition NEE_EC into GPP_EC and ER_EC (Lasslop et al., 2010b; Runkle et al., 2013), capturing the saturation of leaf-level photosynthesis at high solar irradiance. In the original L4C model, APAR directly scales the dynamic range of GPP and is not transformed through a transfer function into a stress scalar, as is the case for the other environmental variables (Eq. 3a). In GPP₂, this linear dependence is replaced by a nonlinear stress scalar, S_APAR, inspired by the light-response curve, which is defined as follows:

  $$\begin{array}{ll}{\displaystyle }\text{GPP}\left(t\right)& {\displaystyle }={\text{GPP}}_{\text{max}}\cdot S_{\text{APAR}}\left(t\right)\cdot S_{\text{MNT}}\left(t\right)\\ \text{(8a)}& {\displaystyle }& {\displaystyle }\cdot S_{\text{VPD}}\left(t\right)\cdot S_{\text{RZSM}}\left(t\right)\text{(8b)}& {\displaystyle }{S}_{\text{APAR}}\left(t\right)& {\displaystyle }={\displaystyle \frac{\text{APAR}\left(t\right)}{\text{APAR}\left(t\right)+{\text{APAR}}_{\text{crit}}}}\end{array}$$

  At low APAR, GPP increases rapidly, but as APAR increases, the rate of increase slows down, and GPP asymptotically approaches a maximum value, GPP_max. This behavior is shown in Appendix B, Fig. B1(D1). GPP_max and APAR_crit are free parameters estimated during the optimization process.

• GPP₃ (defined as GPP₂ with additional adjustments): In the original L4C model, GPP responds linearly to MNT, VPD, and RZSM (Eq. 7a, b, c). In GPP₃, GPP responses to these variables are modeled with increased flexibility, with no assumed shape other than monotonicity, to allow varying rates of change across different ranges. To achieve this, they are redefined as logistic ramp functions:

  $$\begin{array}{}\text{(9a)}& {\displaystyle }\begin{array}{rl}& S_{\text{MNT}}\left(t\right)={\displaystyle \frac{g\left(\text{MNT}\right(t\left)\right)-g\left({\text{MNT}}_{min}\right)}{g\left({\text{MNT}}_{max}\right)-g\left({\text{MNT}}_{min}\right)}}\text{with }g\left(\text{MNT}\right(t\left)\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp (-{\mathit{\gamma }}_{\text{MNT}}\cdot (\text{MNT}\left(t\right)-{\text{MNT}}_{\text{crit}}\left)\right)}}\end{array}\text{(9b)}& {\displaystyle }\begin{array}{rl}& S_{\text{VPD}}\left(t\right)={\displaystyle \frac{g\left(\text{VPD}\right(t\left)\right)-g\left({\text{VPD}}_{max}\right)}{g\left({\text{VPD}}_{min}\right)-g\left({\text{VPD}}_{max}\right)}}\text{with }g\left(\text{VPD}\right(t\left)\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp ({\mathit{\gamma }}_{\text{VPD}}\cdot (\text{VPD}\left(t\right)-{\text{VPD}}_{\text{crit}}\left)\right)}}\end{array}\text{(9c)}& {\displaystyle }\begin{array}{rl}& S_{\text{RZSM}}\left(t\right)={\displaystyle \frac{g\left(\text{RZSM}\right(t\left)\right)-g\left({\text{RZSM}}_{min}\right)}{g\left({\text{RZSM}}_{max}\right)-g\left({\text{RZSM}}_{min}\right)}}\text{with }g\left(\text{RZSM}\right(t\left)\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp (-{\mathit{\gamma }}_{\text{RZSM}}\cdot (\text{RZSM}\left(t\right)-{\text{RZSM}}_{\text{crit}}\left)\right)}}\end{array}\end{array}$$

  In this formulation, the thresholds MNT_min, MNT_max, VPD_min, VPD_max, RZSM_min, and RZSM_max are used to scale the stress scalars between 0 and 1, and are fixed to 273.15, 293.15 K, 0, 2.5 kPa, 0, and 1 m³ m⁻³, respectively. The parameters MNT_crit, γ_MNT, VPD_crit, γ_VPD, RZSM_crit, and γ_RZSM are treated as free parameters and are estimated during the optimization process. The behavior of the logistic ramp functions is illustrated in Appendix B, Fig. B1(A2, B2).

• GPP₄ (defined as GPP₃ with additional adjustments): Growing degree days (GDD) are widely used in agricultural and ecological studies as a proxy to plant development (Fotouo Makouate and Zude-Sasse, 2025), and have recently been used to develop a phenology scheme that improved GPP modeling in a temperate bog (He et al., 2025). In the present formulation, GDD is incorporated into GPP modeling to capture the vegetation green-up and senescence phases through an additional stress scalar (S_GDD). GDD is first derived from mean air temperature (AT) using a base temperature of 273.15 K. It is then normalized for each site and each year using the annual minimum and maximum values, resulting in a normalized range from 0 to 1. This normalization ensures that GDD acts as a seasonal shape or trend driver, rather than a magnitude driver, which is instead represented by the instantaneous variables (APAR, MNT, VPD, and RZSM). Hereafter, GDD refers to normalized GDD and is used to derive S_GDD, as follows:

  $$\begin{array}{ll}{\displaystyle }{\text{GPP}}_{\mathrm{4}}\left(t\right)& {\displaystyle }={\text{GPP}}_{\text{max}}\cdot S_{\text{APAR}}\left(t\right)\cdot S_{\text{MNT}}\left(t\right)\\ \text{(10a)}& {\displaystyle }& {\displaystyle }\cdot S_{\text{VPD}}\left(t\right)\cdot S_{\text{RZSM}}\left(t\right)\cdot S_{\text{GDD}}\left(t\right){\displaystyle }{S}_{\text{GDD}}\left(t\right)& {\displaystyle }={\displaystyle \frac{g\left(\text{GDD}\right(t\left)\right)}{\text{g}\left(\frac{a}{a+b}\right)}}\text{with }g\left(\text{GDD}\right(t\left)\right)\\ \text{(10b)}& {\displaystyle }& {\displaystyle }=\text{GDD}(t{)}^a\cdot (\mathrm{1}-\text{GDD}\left(t\right){)}^b\end{array}$$

  S_GDD is defined as a beta-like, bell-shaped function, normalized between 0 and 1 (Appendix B, Fig. B1(C2)). The parameters a and b are treated as free parameters and are estimated during the optimization process.

• GPP₅ (defined as GPP₄ with additional adjustments): Some studies have reported that water-saturated soil conditions limit oxygen and nutrient availability to plant roots, restrict cellular respiration, and consequently hinder photosynthetic activity (Kreuzwieser et al., 2004; Nawaz et al., 2025). Additionally, Peng et al. (2024) showed that the relationship between GPP and soil moisture may follow a bell-shaped curve at EC tower sites with PFTs similar to those in the present study (e.g., ENF, SHR, GRA; Table 1). A similar response has also been reported for peatlands (Valkenborg et al., 2023). In GPP₅, it is similarly assumed that GPP peaks at an optimal RZSM level, beyond which excessive moisture reduces efficiency. This assumption is tested by redefining S_RZSM as a beta-like, bell-shaped function, analogous to S_GDD (Eq. 10b and Fig. B1(C2) in Appendix B):

  $$\begin{array}{}\text{(11)}& \begin{array}{rl}{S}_{\text{RZSM}}\left(t\right)& ={\displaystyle \frac{g\left(\text{RZSM}\right(t\left)\right)}{\text{g}\left(\frac{a}{a+b}\right)}}\text{with }g\left(\text{RZSM}\right(t\left)\right)\\ & =\text{RZSM}(t{)}^a\cdot (\mathrm{1}-\text{RZSM}\left(t\right){)}^b\end{array}\end{array}$$

  The parameters a and b are treated as free parameters and are estimated during the optimization process.

Figure 2 Ecosystem responses used to compute modeled gross primary production (GPP) from the L4C model (GPP_L4C, column A) and from the Arctic-Subarctic (AS) adapted formulations (GPP₁ through GPP₅, columns B–F) in upland tundra during the growing season. The ecosystem responses are expressed as stress scalars, S_x, for each environmental variable x, where x represents absorbed photosynthetically active radiation (APAR), minimum air temperature (MNT), vapor pressure deficit (VPD), root zone soil moisture (RZSM), and normalized growing degree days (GDD). In the background of each subplot, the histogram of the corresponding environmental variable is shown in light grey. All ecosystem responses of a given model formulation are calibrated jointly using eddy covariance GPP (GPP_EC) as reference, and not based on the histogram values. In subplot A4, RZSM^* refers to the normalized RZSM used as input in the original L4C model. Subplots A1 and B1 are omitted because APAR is used as a direct input, rather than through a stress scalar, in GPP_L4C and GPP₁ formulations. Similarly, subplots A5-D5 are omitted because GDD is not used as input in the GPP_L4C, GPP₁, GPP₂, and GPP₃ formulations. Refer to Sects. 3 and 4.1 for a detailed description of each model formulation. The last row displays scatter plots between modeled GPP against GPP_EC. It is a spatiotemporal comparison, accounting for both temporal and spatial variability. It includes all available GPP_EC data points for upland tundra after filtering, including those used for calibration (1650 data points; Sect. 4.3).

Figure 3 Same as Fig. 2, but for taiga forests. The spatiotemporal comparison includes all available GPP_EC data points for taiga forests after filtering, including those used for calibration (4632 data points; Sect. 4.3).

Figure 4 Same as Fig. 2, but for wetlands. The spatiotemporal comparison includes all available GPP_EC data points for wetlands after filtering, including those used for calibration (3653 data points; Sect. 4.3).

Table 3 Summary of the specificity of each Arctic-Subarctic (AS) adapted model formulation of the L4C model. Column 1 lists the gross primary production (GPP) and ecosystem respiration (ER) formulations, labeled GPP₁ through GPP₅ and ER₁ through ER₅, respectively. The free parameters estimated during the optimization process are provided in column 2. Each GPP and ER formulation is incrementally built on the previous one by incorporating its modifications along with an additional change summarized in column 3. APAR, MNT, VPD, RZSM, SSM, GDD, NPP, L_fall, SOC, R_base, denote absorbed photosynthetically active radiation, minimum air temperature, vapor pressure deficit, rootzone soil moisture, surface soil moisture, normalized growing degree days, net primary productivity, litterfall, soil organic carbon, and baseline heterotrophic respiration, respectively. Refer to Sects. 3 and 4.1 for a description of each model formulation and parameter.

Regarding the ER modeling, we tested different formulations for the HR component while leaving the AR component unchanged. The ER formulations, labeled ER₁ through ER₅, are described below (Table 3):

• ER₁: Rather than using the the L_fall estimation scheme from the baseline L4C model version 8, ER₁ instead adopts the one from version 7 (Eq. 5). The L4C model transitioned from version 7 to version 8 over the course of the present study was conducted, during which the L_fall estimation scheme was modified. The version 7 scheme was retained to enable comparison with the one introduced in version 8. Additionally, HR response to SSM (S_SSM) is redefined as a logistic ramp, analogous to S_RZSM in GPP₃ (Eq. 9c and Fig. B1(A2)). The original linear response (Eq. 7d) was directly replaced because the logistic ramp can reproduce a linear behavior if the relationship between SSM and HR is actually linear. The ER response function to ST is unchanged but was recalibrated (Eq. 7e).

• ER₂ (defined as ER₁ with additional adjustments): The L_fall estimation scheme is reverted to that of the baseline L4C model version 8. Hence, ER₁ and ER₂ differ only in the L_fall estimation scheme (version 7 vs. 8), isolating its impact.

• ER₃ (defined as ER₂ with additional adjustments): None of the established flux-partitioning methods requires SOC data to derive GPP and ER from NEE (Reichstein et al., 2005; Lasslop et al., 2010b; Runkle et al., 2013; Helbig et al., 2017a). Consequently, in ER₃, we aimed to capture the added value of incorporating SOC dynamics in ER modeling. To this end, SOC dynamics are replaced by a single constant representing a baseline heterotrophic respiration rate, (R_base), as follows:

  $$\begin{array}{}\text{(12)}& \text{ER(t)}=\mathit{\alpha }\cdot \text{GPP(t)}+R_{\text{base}}\cdot S_{\text{ST}}\left(t\right)\cdot S_{\text{SSM}}\left(t\right)\end{array}$$

  The parameter R_base is treated as a free parameter and is estimated during the optimization process.

• ER₄ (defined as ER₃ with additional adjustments): In the present formulation, we aimed to mimic the flux-partitioning methods that derive R_base every few days (Reichstein et al., 2005; Lasslop et al., 2010b; Runkle et al., 2013; Helbig et al., 2017a). R_base is then redefined as:

  $$\begin{array}{}\text{(13)}& \begin{array}{rl}{R}_{\text{base}}\left(t\right)& =(\mathrm{1}-\mathit{\omega })\cdot R_{\text{base}}(t-\mathrm{1})+\mathit{\omega }\cdot R_{\mathrm{0}}\\ & \cdot \stackrel{\mathrm{‾}}{S_{\text{ST}}\left(t\right)\cdot S_{\text{SSM}}\left(t\right)}\end{array}\end{array}$$

  This follows a first-order auto-regressive (AR1) approach, in which the current value depends on the previous one and the 7 d backward mean of the product of temperature and moisture stress scalars ($\stackrel{\mathrm{‾}}{S_{\text{ST}}\left(t\right)\cdot S_{\text{SSM}}\left(t\right)}$), weighted by the parameter ω (ranging from 0 to 1). This approach introduces acclimation behavior by smoothing short-term variability in environmental conditions. The parameters R₀ and ω are treated as free parameters and are estimated during the optimization process.

• ER₅ (defined as ER₄ with additional adjustments): Endsley et al. (2022) showed that including an O₂ diffusion limitation in the HR response to SSM, which penalizes HR rates under high SSM conditions, improved seasonal ER performance. In ER₅, we adopted the beta-like, bell-shaped function, analogous to S_RZSM in GPP₅ (Eq. 11), to model HR response to SSM. This avoids to collect or estimate O₂ concentration data while still representing diminishing returns on HR under high SSM conditions.

In this study, daily VPD and MNT were retrieved from the Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2), M2T1NXSLV Version 5.12.4 product (Gelaro et al., 2017), instead of the GEOS-5 FP product, because it is sparsely documented. MERRA-2 re-analysis dataset is better constrained by observations, exhibits a climatology comparable to GEOS-5 FP and is used in the L4C model calibration due to its longer period of record. AT, required to derive GDD, was also retrieved from the MERRA-2 M2T1NXSLV product. Finally, RZSM and SSM from the SMAP SPL4SMGP product were retained in volumetric units (m³ m⁻³), unlike in the original L4C model.

4.2 Growing season timing and data filtering

GPP_EC was used as an indicator to identify the growing season (Gonsamo et al., 2013). For each EC tower, GPP_EC values below the noise threshold of 0.05 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ were first attributed to the winter season and removed. Among the remaining values, those below the 10th percentile were considered part of the shoulder seasons (i.e., transitional periods between fully frozen and fully thawed states) and were excluded. The growing season is more commonly defined using fixed fractions of annual maximum GPP_EC rather than percentile-based thresholds (Panwar et al., 2023; Luo et al., 2025). However, such approaches rely on the assumption that the annual maximum GPP_EC is robustly captured each year, which may not hold due to data gaps for instance. In addition, a previous study has shown that the growing season identification is sensitive to the chosen fraction of annual maximum GPP_EC, indicating a lack of consensus on an optimal threshold (Panwar et al., 2023). In contrast, the percentile-based approach used here provides a distribution-based criteria that is less sensitive to extreme values and inter-annual variability, ensuring consistency across years. For instance, two years with the same growing season timing (i.e., identical start and end dates) but different photosynthetic peaks are treated identically, which would not be the case when using the annual maximum GPP_EC. Although the 10th percentile threshold may not be optimal for each GPP_EC time series, it is consistent with the structure of the mean seasonal cycle of GPP_EC, where it separates winter and shoulder seasons from the growing season (Fig. B2). In addition, similar separations are obtained when using thresholds within the 5th–15th percentile range, indicating low sensitivity to the chosen percentile. Finally, GPP_EC and ER_EC values above the 99th percentile were treated as outliers and removed.

Complementary filtering flags were applied to ensure biophysical plausibility of root-level soil activity and photosynthesis during the growing season from a modeling perspective. The specific criteria for these flags are as follows:

• ST_10 cm ≥ 275.15 K (i.e., 2 °C above freezing)

• ST_20 cm ≥ 275.15 K

• ST_39 cm ≥ 275.15 K (applied to ENF sites only, see Table 1)

• MNT ≥ 275.15 K

ST_20 cm and ST_39 cm refer to soil temperature at 20 and 39 cm depths, respectively, and were retrieved from the SMAP SPL4SMGP product. For the remainder of this study, ST refers to ST_10 cm as ST_20 cm and ST_39 cm are not used further.

After filtering, a total of 1650 data points (23 %) remained for the upland tundra ecosystem, 4632 data points (33 %) for the taiga forest ecosystem, and 3653 data points (30 %) for the wetland ecosystem (Table 1).

4.3 Model formulation calibration

The AS-adapted GPP and ER formulations were calibrated separately for each ecosystem type (upland tundra, taiga forests, and wetlands), using GPP_EC and ER_EC as reference targets. The optimization framework for calibration used least-squares minimization via the MATLAB lsqcurvefit function (MathWorks, Inc., 2023), which minimizes the sum of squared differences between the model outputs and target values. This is an unconstrained optimization, with no additional penalty terms applied. To mitigate overfitting, the optimization process was repeated 100 times using different random training (70 %) and testing (30 %) subsets of the data for each ecosystem type: 1155 (495) data points for upland tundra, 3243 (1389) for taiga forests and 2557 (1096) for wetlands. In each optimization run, model parameters were estimated using the corresponding training subset, and final parameter values were defined as the median of the parameter estimates across all 100 runs. This approach aimed to capture diverse data combinations and promote a more stable and representative parameterization by smoothing out the influence of outliers or any individual biased subset. A 70 % subset size was arbitrarily chosen to balance between providing sufficient data for robust model calibration and retaining enough data variability across the 100 iterations (Martinez Molera, 2025). Contrary to the global calibration of the original L4C model, no reference SOC data were used to constrain the recalibration over the AS environments. As a result, in ER₁ and ER₂, only the SOC₁ pool was modeled to avoid potential parameter compensation and to prevent unrealistic SOC distribution across the original three pools. This reduction in the number of SOC pools constitutes an additional adjustment for ER₁ and ER₂ that is confounded with the change in L_fall estimation scheme (Sect. 4.1). An initial guess was necessary for SOC₁ on 31 March 2015, to explicitly solve the SOC dynamics, since the system is formulated recursively and requires a starting value to iterate forward in time (Eq. 4). 31 March 2015, was chosen as the start of the simulation because it precedes the first date of the period of study. The initial guess was set to 0 gC m⁻² for the first spin-up iteration, providing a neutral starting point to avoid biasing the early simulation. It was subsequently updated using the SOC value on 31 March 2022, corresponding to the last 31 March within the study period. A total of 20 spin-up iterations were performed.

4.4 Model formulation evaluation

For each ecosystem type, the AS-adapted GPP and ER formulations were evaluated spatiotemporally, capturing the combined effects of spatial and temporal variability, with GPP_EC and ER_EC used as reference targets. For each of the 100 optimization runs (Sect. 4.3), the Pearson correlation coefficient (r), the unbiased root mean square error (ubRMSE), and the bias (B) were computed separately for the training and testing subsets, enabling the assessment of model calibration and generalization, respectively. Median values across the 100 runs were reported. The trade-off between goodness of fit and model complexity was quantified using the Akaike Information Criteria (AIC) and Bayesian Information Criteria (BIC), which were applied on the training subsets. The best-performing ER and GPP formulations were identified based on the lowest ΔAIC and ΔBIC values, defined as AIC − AIC_L4C and BIC − BIC_L4C, respectively. Median values across the 100 runs were reported. As a complementary analysis, site-level (temporal-only) performance was evaluated for each EC tower by computing r, ubRMSE, and B using model outputs obtained with the median parameter values across the 100 runs. For this analysis, training and testing data were pooled, and median metrics across EC towers were reported.

As the ER formulations require GPP as an input, all 25 possible ER configurations were systematically evaluated (ER₁ through ER₅ combined with GPP₁ through GPP₅). This approach also yields 25 corresponding NEE configurations, as NEE is defined as the difference between ER and GPP (Eq. 3). Evaluating all combinations enables assessment of whether the best-performing GPP formulation is associated with the best-performing ER formulation, and whether the combination of the best-performing ER and GPP formulations also results in the most accurate NEE configuration. It additionally enables identification of cases where GPP and ER formulations that perform less well individually (relative to ER_EC and GPP_EC) nevertheless yield a more accurate NEE configuration (relative to NEE_EC). Such behavior may arise because errors in modeled ER and GPP can either compensate or accumulate when computing NEE.

The 25 NEE configurations were evaluated following the same approach as for the GPP and ER formulations, using NEE_EC as the reference target. Although NEE configurations were not directly calibrated, the same training and testing subsets were retained for consistency with the GPP and ER formulation evaluations, ensuring that each NEE_EC value was paired with its corresponding GPP_EC and ER_EC values.

Because the evaluation of 25 ER formulations produces many metrics, only ER₁ through ER₅ using the GPP input that yields the best performance are presented in the main text. Results for all configurations are provided in Appendix D. This avoids bias arising from the choice of GPP input when comparing ER formulations. For NEE, summarized performance for all configurations is presented in the main text, while detailed results are provided in Appendix E. Hereafter, NEE_i, j denotes modeled NEE computed as ER_i − GPP_j.

5 Results

5.1 Gross primary production

5.1.1 Upland tundra

Based on the spatiotemporal evaluation of the testing splits for upland tundra (Table 4), GPP₁ performs better than GPP_L4C in terms of B and ubRMSE (−0.33 vs. 0.41 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and 1.19 vs. 1.46 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively). However, GPP_L4C achieves a higher r value (0.56 vs. 0.64). Introducing a nonlinear light-response in GPP₂ (Eq. 8) leads to better performance compared with GPP₁, reducing B from −0.33 to −0.02 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. In addition, ubRMSE decreases, and r increases, approaching the r observed for GPP_L4C. Replacing linear ramps with logistic ramps to simulate ecosystem responses in GPP₃ (Eq. 9) increases model complexity but provides limited improvement over GPP₂. GPP₄ incorporates GDD through an additional stress scalar (Eq. 10). This results in improved r and ubRMSE compared with GPP₃ (0.75 vs. 0.64, and 0.84 vs. 0.96 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively). In GPP₅, a bell-shaped function is used to simulate the influence of RZSM (Eq. 11). This further improves r and ubRMSE, though B slightly increases (−0.03 vs. −0.01 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$).

Table 4 Performance of modeled gross primary production (GPP) against daily averaged eddy covariance (EC) GPP (GPP_EC) for upland tundra, taiga forests, and wetlands during the growing season. GPP_L4C refers to outputs from the original L4C model, while GPP₁ through GPP₅ correspond to the five Arctic–Subarctic (AS) adapted formulations (Sects. 3 and 4.1). The Pearson correlation coefficient is denoted by r (dimensionless), and ubRMSE and B denote the unbiased root mean square error and bias, respectively (in $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). A positive (negative) B indicates overestimation (underestimation) of GPP_EC. (a) Spatiotemporal performance metrics derived from 100 random splits (70 % training, 30 % testing), accounting for both spatial and temporal variability. The total number of data points used for training (testing) is 1155 (495) for upland tundra, 3243 (1389) for taiga forests, and 2557 (1096) for wetlands (Sect. 4.3). Median values across the 100 splits are reported. (b) Site-level (temporal-only) performance metrics computed for each EC tower using model outputs obtained with the median parameter values across the 100 splits. Metric values are reported as the median across towers. For this analysis, training and testing data are pooled.

Considering the site-level evaluation (Table 4), GPP₄ and GPP₅ exhibit the highest r (0.77 vs. 0.76) and the lowest ubRMSE (0.65 vs. 0.71 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). In contrast, GPP₅, GPP₃, and GPP₂ show the lowest B with 0.04, 0.05, and 0.08 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively.

Across all formulations, S_VPD remains equal to 1 throughout the entire range of VPD variability (Fig. 2(B3–F3)). The use of a non linear light-response in GPP₂ introduces an early saturation (Fig. 2(B6–F6)), where modeled GPP peaks are lower than those of GPP_EC. This premature flattening is progressively reduced in GPP₄ and GPP₅, due to the incorporation of GDD and the use of a bell-shaped function for simulating RZSM influence.

The lowest ΔAIC and ΔBIC are obtained for GPP₅, followed by GPP₄, GPP₃, GPP₂ and GPP₁ in last place (Table C1). For GPP₂ and GPP₃, ΔAIC and ΔBIC are equal within each subset (training or testing) because both formulations have the same number of free parameters as GPP_L4C.

5.1.2 Taiga forests

The spatiotemporal evaluation of the testing splits for taiga forests (Table 4) indicates that GPP₁ performs better than GPP_L4C in terms of r (0.62 vs. 0.58) and ubRMSE (1.75 vs. 1.93 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), although GPP₁ shows higher B (−0.44 vs. −0.37 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). As in upland tundra (Sect. 5.1.1), introducing a nonlinear light-response in GPP₂ (Eq. 9) leads to overall better performance compared with GPP₁, notably reducing B from −0.44 to −0.05 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. GPP₃ does not offer improvement over GPP₂, aside from a slight reduction in B (−0.02 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ vs. −0.05 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Due to the inclusion of GDD through an additional stress scalar (Eq. 10), GPP₄ outperforms GPP₃ (0.73 vs. 0.67 for r, and 1.30 vs. 1.41 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ for ubRMSE). Switching to a bell-shaped function for simulating RZSM influence in GPP₅ (Eq. 11) does not result in improved performance compared with GPP₄.

Considering the site-level evaluation (Table 4), the highest r are obtained for GPP₄ and GPP₅ (0.79 vs. 0.76). These two formulations also achieve the lowest ubRMSE, with 1.12 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ for GPP₅ and 1.15 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ for GPP₄. GPP₄ and GPP₄ exhibits the lowest B (−0.13 vs. −0.14 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$).

RZSM appears to be a negligible input in GPP_L4C, as S_RZSM remains equal to 1 throughout the entire range of RZSM variability (Fig. 3(A4)). However, RZSM gains more effect in GPP₁ through GPP₅, although its effect remains weaker than that of MNT, VPD, and GDD (Fig. 3(B4–F4)). As in upland tundra (Sect. 5.1.1), the use of a nonlinear light-response in GPP₂ introduces an early saturation, underestimating GPP_EC peaks (Fig. 3(B6–F6)). However, this premature flattening persists in GPP₄ and GPP₅, despite the incorporation of GDD and the use of a bell-shaped function for simulating RZSM influence.

The lowest ΔAIC and ΔBIC are obtained for GPP₅, closely followed by GPP₄, then GPP₃, GPP₂ and GPP₁ in last place (Table C1). For GPP₂ and GPP₃, ΔAIC and ΔBIC are equal within each subset (training or testing) because both formulations have the same number of free parameters as GPP_L4C.

5.1.3 Wetlands

Based on the spatiotemporal evaluation of the testing splits for wetlands (Table 4), GPP₁ outperforms GPP_L4C, notably exhibiting reduced ubRMSE and B (1.33 vs. 2.18 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, and −0.39 vs. 1.32 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively). As seen in upland tundra and taiga forests (Sects. 5.1.1 and 5.1.2), introducing a nonlinear light-response in GPP₂ (Eq. 8) results in improved r (0.63 vs. 0.53), reduced ubRMSE (1.04 vs. 1.33 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and reduced B (−0.05 vs. −0.39 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), compared with GPP₁. GPP₃ shows only a minor improvement over GPP₃. Due to the inclusion of GDD through an additional stress scalar (Eq. 10), GPP₄ outperforms GPP₃ (0.72 vs. 0.65 for r, and 0.92 vs. 1.01 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ for ubRMSE). In contrast to upland tundra and taiga forests (Sect. 5.1.1 and 5.1.2), using a bell-shaped function for simulating RZSM influence in GPP₅ (Eq. 11), result in degraded performance compared with GPP₄.

Considering the site-level evaluation (Table 4), GPP₄ and GPP₅ exhibit the highest r (0.79 vs. 0.76). GPP₄ and GPP₃ exhibits the lowest ubRMSE (0.76 vs. 0.79 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively), closely followed by GPP₅ (0.82 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and GPP₃ (0.85 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). B is similar for GPP₂ through GPP₅, with −0.22, −0.22, −0.23, and −0.19 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively.

Similar to upland tundra (Sect. 5.1.1), S_VPD remains equal to 1 throughout the entire range of VPD variability across all formulations (Fig. 4(B3–F3)). Likewise, as observed in taiga forests (Sect. 5.1.2), RZSM appears to be a negligible input in GPP_L4C, with S_RZSM staying equal to 1 throughout the entire range of RZSM variability (Fig. 4(A4)). However, RZSM gains a considerable effect in GPP₁ through GPP₅, especially under dry conditions (Fig. 4(B4–F4)). The introduction of the nonlinear light-response in GPP₂ triggers an early saturation (Fig. 4(B6–F6)), that persists throughout GPP₅.

In contrast to upland tundra and taiga forests (Sects. 5.1.1 and 5.1.2), the lowest ΔAIC and ΔBIC are obtained for GPP₄, followed by GPP₅, GPP₃, GPP₂, and GPP₁ in last place (Table C1). For GPP₂ and GPP₃, ΔAIC and ΔBIC are equal within each subset (training or testing) because both formulations have the same number of free parameters as GPP_L4C.

5.2 Ecosystem respiration

For upland tundra and taiga forests, ER formulations using GPP₅ as input are compared, as GPP₅ yielded the lowest ΔAIC and ΔBIC (Sect. 5.1.1, 5.1.2). For wetlands, the set using GPP₄ as input is compared instead, because GPP₄ achieved the lowest ΔAIC and ΔBIC for this ecosystem type (Sect. 5.1.3). Results for all 25 ER formulations are provided in Tables D1–D6.

5.2.1 Upland tundra

Based on the spatiotemporal evaluation of the testing splits for upland tundra (Table 5), ER₁, which uses the approach where mean annual NPP is allocated uniformly across the year to L_fall (Eq. 5), performs better than ER_L4C. It exhibits higher r (0.52 vs. 0.44), reduced ubRMSE (0.72 vs. 0.99 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and reduced B (−0.12 vs. 0.37 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Switching to the dynamic allocation in ER₂ (Eq. 6) leads to enhanced r (0.61) and ubRMSE (0.64 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) compared with ER₁. Using a constant R_base term instead of a SOC model in ER₃ (Eq. 12) results in r and ubRMSE close to those of ER₂, and lower B (−0.01 vs. −0.12 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Using a dynamic R_base in ER₄ (Eq. 13) shows better performance than ER₃ with greater r (0.72 vs. 0.60) and reduced ubRMSE (0.53 vs. 0.62 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). In ER₅, the use of a bell-shaped function to simulate SSM influence does not provide a clear benefit.

Table 5 Same as Table 4, but for ecosystem respiration (ER). GPP₅ was used as the GPP input for ER₁ through ER₅ for upland tundra and taiga forests. For wetlands, GPP₄ was used instead.

Considering the site-level evaluation (Table 5), the highest r are obtained for ER₂ and ER₅ (0.66 vs. 0.60). The lowest ubRMSE is observed for ER₄ and ER5 (0.41 vs. 0.42 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), closely followed by the other formulations (up to 0.51 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ for ER₁). In terms of B, ER₅, ER₄, and ER₂ show the smallest values with 0.00, −0.01, and 0.03 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$.

ST appears to have a stronger effect than SSM in ER₁ through ER₄ (Fig. 5(B1–E1) vs. (B2–E2)). S_ST mainly oscillates around 0.5, while S_SSM rapidly reaches 1 under dry conditions, and wet conditions do not constrain model outputs. In ER₅, the use of a bell-shaped function to simulate SSM increases its effect (Fig. 5(F2)), but without any improvement in performance, as previously observed.

Figure 5 Ecosystem responses used to compute modeled ecosystem respiration (ER) from the L4C model (ER_L4C, column A) and from the Arctic-Subarctic (AS) adapted formulations (ER₁ through ER₅, columns B–F) in upland tundra during the growing season. The ecosystem responses are expressed as stress scalars, S_x, for each environmental variable x, where x represents soil temperature (ST) and surface soil moisture (SSM). In the background of each subplot, the histogram of the corresponding environmental variable is shown in light grey. All ecosystem responses of a given model formulation are calibrated jointly using daily-averaged eddy covariance ER (ER_EC) as reference, and not based on the histogram values. In subplot A2, SSM^* refers to the SSM in relative wetness unit used as input in the original L4C model. Refer to Sects. 3 and 4.1 for a detailed description of each model formulation. The last row displays scatter plots between modeled ER against ER_EC. It is a spatiotemporal comparison, accounting for both temporal and spatial variability. It includes all available ER_EC data points for upland tundra after filtering, including those used for calibration (1650 data points; Sect. 4.3). GPP₅ was used as input for ER₁ through ER₅ (Sect. 5.2.1).

The lowest ΔAIC and ΔBIC are obtained for ER₄, closely followed by ER₅, then ER₃ and ER₂, and ER₁ in last place (Table D6). Overall, a similar pattern is observed regardless of the GPP input, and the performance of each ER formulation remains comparable or improves as the GPP input increases in complexity (lowest with GPP₁, highest with GPP₅; Tables D1–D6).

5.2.2 Taiga forests

Based on the spatiotemporal evaluation of the testing splits for taiga forests (Table 5), ER₁ outperforms ER_L4C, showing enhanced r (0.52 vs. 0.34), reduced ubRMSE (1.28 vs. 1.71 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), but slightly higher B (−0.17 vs. −0.12 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Switching from a constant to dynamic allocation of mean annual NPP to L_fall in ER₂ (Eq. 6) results in a minor improvement in performance over ER₁. Using a constant R_base term instead of a SOC model in ER₃ (Eq. 12) exhibits similar r (0.54 vs. 0.54) and ubRMSE (1.23 vs. 1.25 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), but lower B (0.02 vs. −0.12 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), compared with ER₂. Using a dynamic R_base in ER₄ (Eq. 13) shows better performance than ER₃ with greater r (0.59 vs. 0.54) and reduced ubRMSE (1.17 vs. 1.23 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). In ER₅, the use of a bell-shaped function to simulate SSM influence does not provide a clear benefit.

Considering the site-level evaluation (Table 5), r is similar across formulations, ranging from 0.61 to 0.65. The same pattern is observed for ubRMSE, which ranges from 0.95 to 1.01 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. The lowest B are obtained for ER₃, ER₄, and ER₂, with 0.00, −0.05, and −0.09 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively.

SSM appears to be a negligible input in ER_L4C, as S_SSM remains equal to 1 throughout the entire range of SSM variability (Fig. 6(A2)). However, SSM gains more effect in ER₁ through ER₅, although its effect remains weaker than that of ST (Fig. 6(B2–F2)).

Figure 6 Same as Fig. 5, but for taiga forests. The spatiotemporal comparison includes all available ER_EC data points for taiga forests after filtering, including those used for calibration (4632 data points; Sect. 4.3). GPP₅ was used as input for ER₁ through ER₅ (Sect. 5.2.2).

The lowest ΔAIC and ΔBIC are obtained for ER₅, nearly tied with ER₄, then ER₃, ER₂, and ER₁ in last place. As observed in upland tundra (Sect. 5.2.1), this ranking remains consistent regardless of the GPP input, and the performance of each ER formulation is comparable or improves as the GPP input increases in complexity (lowest with GPP₁, highest with GPP₅; Tables D1–D6).

5.2.3 Wetlands

Based on the spatiotemporal evaluation of the testing splits for taiga forests (Table 5), ER₁ outperforms ER_L4C, exhibiting enhanced r (0.49 vs. 0.24), reduced ubRMSE (0.81 vs. 1.80 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and reduced B (−0.05 vs. 1.77 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Switching from a constant to dynamic allocation of mean annual NPP to L_fall in ER₂ (Eq. 6) slightly increases r (0.53 vs. 0.49), and reduces ubRMSE and B (0.78 vs. 0.81 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, −0.03 vs. −0.05 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively). Using a constant or dynamic R_base term instead of a SOC model in ER₃ and in ER₄ (Eqs. 12 and 13) does not result in enhanced performance. In ER₅, the use of a bell-shaped function to simulate SSM influence does not provides any benefits either.

Considering the site-level evaluation (Table 5), r is similar across formulations, ranging from 0.63 to 0.67. The same pattern is observed for ubRMSE, with the highest value assigned to ER₁ (0.56 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and the lowest to ER₄ (0.48 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). ER₁ and ER₂ exhibiting the lowest B, with −0.12 and −0.09 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively (Fig. 8(C2)).

As observed in taiga forests (Sect. 5.2.2), SSM appears to be a negligible input in ER_L4C, with S_SSM staying equal to 1 throughout the entire range of SSM variability (Fig. 7(A2)). However, SSM gains a considerable effect in ER₁ through ER₅, especially under dry conditions (Fig. 7(B2–F2)). ER_L4C may overestimate ER_EC by up to a factor of two or three (Fig. 7(A3)). This magnitude discrepancy is largely removed in ER₁ through ER₅, but a pattern persists, in which ER_EC are systematically overestimated at low values (approximately 0 to 1 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) across all formulations (Fig. 7(B3–F3)).

Figure 7 Same as Fig. 5, but for wetlands. The spatiotemporal comparison includes all available ER_EC data points for wetlands after filtering, including those used for calibration (3653 data points; Sect. 4.3). GPP₄ was used as input for ER₁ through ER₅ (Sect. 5.2.3).

Figure 8 Spatiotemporal performance of modeled net ecosystem CO₂ exchange (NEE) against daily-averaged eddy covariance NEE (NEE_EC) for upland tundra, taiga forests, and wetlands during the growing season. NEE is computed as the difference between ecosystem respiration (ER) and gross primary production (GPP). GPP_L4C and ER_L4C refer to outputs from the original L4C model, while GPP₁ through GPP₅ and ER₁ through ER₅ correspond to the Arctic–Subarctic (AS) adapted formulations (Sects. 3 and 4.1). The (ER_i, GPP_j) grid cell corresponds to the NEE_i,j configuration. The Pearson correlation coefficient is denoted by r (dimensionless), and ubRMSE and $\left|B\right|$ denote the unbiased root mean square error and absolute bias, respectively (in $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Spatiotemporal metrics were derived from 100 random splits (70 % training, 30 % testing), accounting for both spatial and temporal variability. The total number of data points used for training (testing) is 1155 (495) for upland tundra, 3243 (1389) for taiga forests, and 2557 (1096) for wetlands (Sect. 4.3). Reported values correspond to median metrics computed over the testing splits (Tables E1–E5).

The lowest ΔAIC and ΔBIC are obtained for ER₂, then ER₃, ER₄, ER₅, and ER₁ in last place. Nevertheless, all formulations are nearly tied. A consistent tied ranking is observed regardless of which GPP is used as input (Tables D1–D6). The performance of each ER formulation (ER₁ through ER₅) remains comparable across formulations, although using GPP₁ as input lead to the lowest performance.

5.3 Net ecosystem CO₂ exchange

This subsection presents the performance of the AS-adapted NEE formulations (ER₁ through ER₅ combined with GPP₁ through GPP₅) and NEE_L4C relative to NEE_EC.

5.3.1 Upland tundra

Based on the spatiotemporal evaluation for upland tundra, NEE_{1–5, 4–5} exhibit higher r than the other NEE configurations, with values above 0.6 (Fig. 8(A1)). In contrast, NEE_L4C shows a r of 0.57. The lowest ubRMSE are obtained for NEE_{2–5, 5} and NEE_{2–3, 4} (Fig. 8(B1)), ranging from 0.65 to 0.70 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Fig. 8(B1)). By comparison, NEE_L4C yields a value of 0.84 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. NEE_{3–5, 2–5} and NEE_1,1 show the lowest absolute B, with values below 0.05 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, comparable to that of NEE_L4C (0.03 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; Fig. 8(C1)).

Overall, NEE_{2–5, 4–5} appears to perform best in terms of higher r, and lowest ubRMSE and B. This set of configurations includes the ER and GPP formulations associated with the lowest ΔAIC and ΔBIC (ER₄, ER₄, GPP₄, and GPP₅; Sect. 5.1.1 and 5.2.1). Although their performance is broadly comparable, some configurations may perform better depending on the metric considered (r, ubRMSE or B). A similar pattern is observed considering the site-level evaluation (Tables E1–E5).

5.3.2 Taiga forests

For taiga forests, the spatiotemporal evaluation reveals that NEE_{4–5, 4–5} achieve the highest r with values between 0.60 and 0.65, whereas NEE_L4C shows a r of 0.55 (Fig. 8(A2)). NEE_{4–5, 4} and NEE_{5, 5} exhibit the lowest ubRMSE, ranging from 1.05 to 1.10 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Fig. 8(B2)). NEE_{4, 5}, NEE_{3–5, 3}, and NEE_{3, 4} closely follow, with values between 1.10 and 1.15 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. In contrast, NEE_L4C shows a ubRMSE of 1.17 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. NEE_{3–4, 2–5} and NEE_{5, 3–5} exhibit absolute B below 0.05 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, lower than both the other NEE configurations and NEE_L4C (0.26 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; Fig. 8(C2)).

Overall, NEE_{4–5, 4–5} appears to perform best in terms of higher r, and lowest ubRMSE and B. This set of configurations includes the ER and GPP formulations associated with the lowest ΔAIC and ΔBIC (ER₄, ER₄, GPP₄, and GPP₅; Sects. 5.1.2 and 5.2.2). A similar pattern is observed considering the site-level evaluation (Tables E1–E5).

5.3.3 Wetlands

Based on the spatiotemporal evaluation for wetlands, NEE_{3–5, 4} exhibit the highest r, ranging from 0.45 and 0.5, comparable to that of NEE_L4C (0.47; (Fig. 8(A3))). NEE_{2–5, 4}, NEE_{5, 5}, NEE_{2, 5}, and NEE_{2, 3} exhibit the lowest ubRMSE with values between 0.90 and 0.95 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, slightly lower than that of NEE_L4C (0.99 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; (Fig. 8(B3))). In terms of absolute B, NEE_{1–5, 2–5} show values below 0.05 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, outperforming both NEE_{1–5, 1} and NEE_L4C (0.45 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; Fig. 8(C₃)).

Overall, NEE_{1–5, 4} appears to perform best in terms of higher r, and lowest ubRMSE and B. This set of configurations includes the ER and GPP formulations associated with the lowest ΔAIC and ΔBIC (ER₂, ER₃, and GPP₄; Sects. 5.1.3 and 5.2.3). A similar pattern is observed considering the site-level evaluation (Tables E1–E5).

6 Discussion

This section first discusses how the modifications implemented in the AS-Adapted GPP and ER model formulations affected their performance relative to GPP_EC and ER_EC. We specifically identify candidate GPP model adjustments for operational implementation, evaluate the contribution of incorporating SOC dynamics into ER modeling, and examine the influence of input variables. We then outline the implications for NEE estimation and highlight the limitations of our study.

6.1 Candidate GPP model adjustments

Implementing a nonlinear light-response function to represent the influence of APAR on GPP (GPP₂, Eq. 8) appears to be the most effective adjustment tested for reducing both ubRMSE and B across the three ecosystem types (Sect. 5.1, and Table 4, Fig. 8). However, our results indicate that this adjustment requires careful parametrization, as it can lead to underestimation of GPP peaks (Figs. 2–4).

Adding GDD into the GPP modeling (Eq. 10) complements the nonlinear light-response adjustment by further reducing ubRMSE and predominantly improving r across the three ecosystem types (Sect. 5.1, Table 4, Fig. 8). These results suggest that the current L4C model may lack a phenological proxy that accounts for the progressive functional adjustment of vegetation to environmental conditions over time (Maire et al., 2012), thereby complementing the instantaneous proxies currently used as inputs. Vegetation indices are assumed to capture vegetation phenology by tracking seasonal changes in canopy structure and greenness. Because GDD is derived from mean air temperature and minimum air temperature is used as a proxy for the instantaneous temperature response of GPP, replacing GDD with a vegetation index may help reduce redundancy (Huang et al., 2019; Pulliainen et al., 2024). In this context, LAI or normalized difference vegetation index (NDVI) could potentially serve this role. However, LAI is likely to introduce additional redundancy, as VIIRS LAI retrievals are already used to derive FPAR, which represents canopy phenology in the GPP formulation (Eq. 2). In contrast, NDVI may provide a more independent phenological proxy without duplicating existing model inputs. Nonetheless, MODIS and VIIRS vegetation indices exhibit large uncertainties at high northern latitudes, particularly during shoulder seasons, due to extensive cloud cover and snow contamination (Xu et al., 2018; Pu et al., 2024).

Overall, the nonlinear light-response adjustment appears to be the strongest candidate for correcting GPP magnitude discrepancies, while incorporating GDD emerges as the most effective adjustment to improve GPP seasonal dynamics. These two findings are consistent with McCallum et al. (2013), where the authors reported that the inclusion of temperature acclimation and nonlinear light-response in GPP modeling in Russian boreal forests improved model performance. Future studies could explore more complex light-response functions, such as a rectangular hyperbolic function, where the parameters may vary temporally with temperature, as suggested by Wang et al. (2014a). However, testing this approach would imply departing from the current multiplicative structure of the L4C model, in which direct mechanistic forcing–response behaviors are represented. Prior work also suggests that vegetation green-up onset is influenced by winter chilling accumulation and precipitation (Fu et al., 2014; Elmendorf and Hollister, 2023). Greater accumulation of chilling days may lead to earlier green-up, as vegetation exposed to colder winter temperatures requires less thermal accumulation (less GDDs) to initiate spring growth. In contrast, higher winter precipitation may contribute to delayed green-up through thicker snowpacks, cooler soil temperatures, and increased cloud cover that reduces incoming radiation. Future improvements could therefore integrate winter chilling days and precipitation into the normalized GDD-based phenological proxy to better represent early-season GPP dynamics.

Finally, it is noteworthy that GDD was normalized using annual minimum and maximum values, implying that early-season GDD values implicitly depend on values from later in the same year. This approach therefore requires the full annual air temperature cycle to be known, introducing a one-year lag in the computation of model outputs that is not appropriate for real-time forecasting. However, the primary objective of this study is not prediction in an operational setting, but rather the retrospective reconstruction of CO₂ fluxes and the estimation of CO₂ budgets. For potential forecasting applications, GDD can instead be normalized using long-term minimum and maximum values across all years, resulting in comparable or slightly lower performance (see GPP₄ and GPP₅ in Table 4 vs. Table C2). This alternative approach avoids any within-year information leakage and suggests that GDD normalization has a limited influence on the overall performance improvements associated with adding GDD as an input.

6.2 Comparison of SOC-based and empirical approaches for ER modeling

Updating the allocation of mean annual NPP to L_fall from a constant to a LAI-based formulation to represent SOC dynamics (ER₁ vs. ER₂; Eqs. 5 and 6) improves ER model performance. Both the spatiotemporal evaluation and the median metrics across EC towers indicate higher r and lower ubRMSE and B, with stronger improvements for upland tundra and weaker improvements for taiga forests and wetlands (Table 5). The benefits are limited relative to the added model complexity, especially in taiga forests and wetlands, compared with the simpler approach that replaces SOC dynamics with a single constant R_base (ER₃, Eq. 12). Based on the spatiotemporal evaluation, introducing temporal variability in R_base (ER₄, Eq. 13) leads to improved performance in upland tundra and taiga forests (Table 5). However, the median metrics across EC towers do not indicate a clear improvement across the three ecosystems (Table 5). Compared to upland tundra and taiga forests, all ER formulations are highly similar for wetlands, regardless of the metrics considered (r, ubRMSE, B, ΔAIC, and ΔBIC) and the type of evaluation (spatiotemporal vs. site-level; Tables 5 and D6).

Overall, using SOC dynamics with the L_fall estimation scheme from the L4C model version 8 to model ER appears to be the most suitable approach, as it performs better than version 7 and is physically grounded and mechanistically interpretable compared with the two empirical approaches. Unfortunately, the improved performance of ER₁ and ER₂ compared to ER_L4C is difficult to interpret, as it reflects the combined effects of changes in the L_fall estimation scheme, SOC pool structure, and GPP input (Sect. 4.1, 4.3). Therefore, this comparison does not constitute a clean test of the added-value of the L_fall allocation scheme version 8 compared to version 7 for the study regions.

Continuing to explore alternative ways to estimate L_fall may be a promising direction for future research. However, the assumption that mean annual NPP can serve as a proxy for the magnitude of L_fall may not be realistic (Sierra et al., 2022). In addition, the timing of NPP allocation to L_fall may not accurately represent litter production dynamics, particularly given the large uncertainties in LAI and FPAR retrievals at high northern latitudes (Xu et al., 2018; Pu et al., 2024). Furthermore, because NPP is derived from modeled GPP, any inaccuracies in GPP propagate directly into modeled NPP, L_fall, SOC, and ultimately ER. Finally, recent work in Alaska has shown that implementing vertical SOC transport to simulate depth-dependent L_fall, SOC distribution, and corresponding HR rates may further improve ER estimates (Yi et al., 2020).

6.3 Tested but unretained GPP and ER model adjustments

Implementing logistic ramps to represent GPP responses to MNT, VPD, and RZSM stress (GPP₃, Eq. 9) provides limited benefits based on both the spatiotemporal and site-level evaluation (GPP₂ vs. GPP₃ in Table 4). This suggests that MNT, VPD, and RZSM may exhibit nonlinear interactions with GPP, but this adjustment appears to be of secondary importance compared to the implementation of a light-response curve and the incorporation of GDD (Sect. 6.1).

The use of a bell-shaped function to represent RZSM influence on GPP (GPP₅, Eq. 11) provides only a limited performance improvement in upland tundra and no improvement in taiga forests (Table 4). One possible explanation is that, in upland tundra, RZSM exhibits both dry and wet conditions across years and EC tower sites (grey histogram in Fig. 2(F4)), whereas in taiga forests, conditions remain mostly dry with less seasonal variation (Fig. 3(F4)). This pattern is supported by the bimodal distribution of RZSM in upland tundra, in contrast to the unimodal distribution in taiga forests. Nevertheless, the RZSM distribution in wetlands is bimodal (Fig. 4(F4)), with both dry and wet conditions, but the bell-shaped function worsens the performance of modeled GPP (Table 4). The diminishing returns under wet conditions appear to penalize model calibration, indicating a different ecosystem response to RZSM in wetlands compared with upland tundra and taiga forests. This may reflect the adaptation of wetland vegetation to anaerobic conditions, where excessive RZSM does not hinder photosynthesis activity. Although several studies show that wetlands and peatlands are more sensitive to drought than to flooding (Churchill et al., 2015; Olefeldt et al., 2017; Heinzelmann et al., 2025), there is no clear evidence in the literature on whether GPP in wetlands continues to scale or levels off under high RZSM conditions. It is also noteworthy that the bell-shaped function does not improve model performance for any ecosystem when normalized RZSM is used (not shown), as is the case in the original L4C model (Sect. 3). Finally, the L4C model methodology focuses on direct mechanistic forcing–response behavior, where instantaneous RZSM data are used as input. However, a temporal lag in GPP response to RZSM saturation may be expected, as it can take several days to weeks for soil oxygen levels to become depleted to the point of restricting aerobic processes under saturation. A larger number of EC towers should also be included to increase RZSM variability during calibration before drawing conclusions about the value of this adjustment for North American AS regions.

As in the GPP modeling, the use of a bell-shaped function to represent SSM influence on ER provides no clear improvement, regardless of ecosystem type or whether dry and wet SSM conditions are included during calibration (Table 5 and Figs. 5–7). These results indicate that an unidirectional function ramp is more appropriate, with dry conditions limiting ER rates, and no diminishing returns under wet conditions. The same conclusion is drawn when SSM expressed in relative wetness units is used (not shown), as in the original L4C model (Sect. 3). This finding partly contrasts with Endsley et al. (2022), who reported improved seasonal ER performance after adding an O₂ diffusion limitation (also based on SSM) to the original monotonic linear response, thereby penalizing ER rates under high SSM. The differing behavior between studies may be attributed first to differences in SSM response functions and, second, to the fact that in Endsley et al. (2022), the SPL4SMGP product did not yet account for peatland hydrology (Reichle et al., 2023).

6.4 Key drivers in shaping model performance

As VPD rises, indicating atmospheric dryness, plants typically show stomatal closure to minimize water loss, which in turn reduces their photosynthetic activity (López et al., 2021). However, in upland tundra and wetlands, GPP appears insensitive to VPD as the corresponding stress scalar S_VPD remains equal to 1 across all five AS-adapted formulations (Figs. 2(B3–F3) and 4(B3–F3)). This suggests that either the VPD response is inadequately represented in the formulations, or that vegetation in these areas is inherently less responsive to stomatal closure than in taiga forests, where S_VPD strongly constrains the modeling (Fig. 3(B3–F3)). Indeed, VPD distributions are similar across the three ecosystem types, which supports the idea that the observed differences in S_VPD are not due to differing environmental conditions, but rather to ecosystem-specific sensitivity. In other words, for the same VPD values, the model applies a stronger constraint to vegetation photosynthetic activity in taiga forests, while in upland tundra and wetlands it remains unconstrained. These findings are consistent with those of Chen et al. (2023), where the authors observed that increasing VPD did not hinder vegetation growth in northern peatlands. Additionally, Zona et al. (2023) reported that VPD was not correlated to GPP at monthly scales in Arctic tundra, while Mirabel et al. (2023) found that tree growth in the Canadian boreal forest responded negatively to rising VPD.

Across all three ecosystem types, the most notable model improvement arises from revising the influence of APAR and AT (through GDD) on GPP (Table 4, Fig. 8, Sect. 6.1). Interestingly, APAR and AT are also the two drivers primarily used to partition NEE_EC into GPP_EC and ER_EC (Sect. 2). This indicates that model performance is inherently entangled with these two drivers, rather than to VPD, RZSM, SSM, SOC, or ST. It is important to note that drivers used in the L4C model formulations are provided at 9 and 25 km resolution (except FPAR), which is coarse relative to EC tower footprints (Sects. 2 and 3). Some of the discrepancies between EC measurements and model estimates may therefore be attributed to representativeness errors, as the coarse model resolution is expected to smooth spatial variability that is captured by the EC measurements. Coupled with the candidate GPP model adjustments (Sect. 6.1), using higher spatial resolution PAR and AT inputs could represent a promising avenue for improving GPP estimates and, consequently, ER estimates in future studies.

6.5 Implications for NEE estimation

Overall, the individually best-performing ER and GPP formulations, when combined, yield the most, or among the most, accurate NEE configurations (Sect. 5.3, Fig. 8, and Tables E1–E5). Among the 25 tested NEE configurations, those using GPP₄ and GPP₅ have the strongest influence on NEE performance, particularly for r, and to a lesser extent for ubRMSE. This is consistent with the inclusion of GDD as a seasonal driver into the GPP modeling (Eq. 10). In contrast, no clear and consistent pattern emerges regarding the impact of ER formulations on NEE performance, even though they exhibit distinct performance when evaluated against ER_EC (Sect. 5.2). The main exception is ER₁, which systematically leads to higher ubRMSE and absolute B in NEE. Finally, although the best-performing GPP and ER formulations showed strong performance gains relative to the original L4C model, particularly for GPP, the resulting improvement in NEE performance is more modest (Tables 4, 5 vs. Tables E1-E5 and Fig. 8). These results suggest that improvements in the representation of ER and GPP do not necessarily lead to comparable improvements in NEE (Fig. B3). Rather than compensating for each other, errors in modeled ER and GPP may accumulate when computing NEE, at least for the formulations tested.

6.6 Influence of temporal and spatial autocorrelation on model performance

The current model validation approach used random 70 %–30 % training–testing subsets for model calibration and evaluation (Sect. 4.3 and 4.4). The 70 %–30 % splits were applied to individual data points. Consequently, model performance may be partly influenced by temporal autocorrelation, since GPP_EC and ER_EC data points from the same EC tower site may be separated by only a few days while belonging to the training and testing subsets, respectively. Model performance may also be influenced by spatial autocorrelation, as data from the same site are included in both the training and testing subsets.

Consequently, two complementary sensitivity analyses were conducted using alternative validation approaches that provide stronger temporal and spatial separation, respectively, between the data subsets used for calibration and evaluation.

The first alternative validation approach consisted of withholding complete years from model calibration instead of using random 70 %–30 % training–testing subsets. For each ecosystem, the available data were assigned to 6-year training and 2-year testing periods. Because the study period spans 8 years, this yielded 28 possible unique train–test runs. The sizes of the training and testing subsets varied among runs because data availability differs among years and sites. This validation approach reduces the potential influence of temporal autocorrelation but introduces variability in subset size among runs, which may itself affect performance metrics. Nevertheless, results obtained with this complementary approach were comparable to those obtained with the random 70 %–30 % splits applied to the individual data points, indicating that the main findings of this study are robust to temporal autocorrelation (Tables C3 and D7 vs. Tables 4 and 5).

In the second alternative validation approach, the 70 %–30 % training–testing subsets were applied to the set of sites within each ecosystem rather than to the individual data points, such that all data from a given site were assigned exclusively to either the training or the testing subset. Because the number of sites differs among ecosystems, the number of possible unique train–test runs was 5, 21, and 28 for upland tundra, taiga forests, and wetlands, respectively. This validation approach reduces the potential influence of spatial autocorrelation while providing a stricter assessment of model ability to generalize to unseen sites. However, as in the first alternative approach, the sizes of the training and testing subsets varied among runs. Overall, results obtained with this complementary approach were comparable to those obtained with the other two validation approaches (Tables C4 and D8 vs. Tables 4 and 5 and vs. Tables C3 and D7). The most notable difference was that absolute B was higher for the testing splits in some cases. However, this does not alter the main findings of the study. This difference may reflect variations in the spatial representativeness of the training subsets, in subset size, in site composition among runs, or a combination of these factors. Nevertheless, these effects cannot be disentangled with the present dataset, particularly for upland tundra, where only five sites were available (Table 1).

6.7 Limitations

The reference GPP_EC and ER_EC used for calibration and evaluation are derived from NEE_EC partitioning, meaning they are not direct measurements but modeled outputs based on NEE_EC and structural assumptions (Appendix A). This creates a potential circularity in model evaluation, as the AS-adapted formulations may share the same structural assumptions as the flux-partitioning algorithms. Consequently, improvements in modeled ER and GPP may partly reflect the model formulations reproducing the behavior of the flux-partitioning algorithms, rather than independently improving the representation of carbon dynamics. In some cases, GPP_EC is constrained to follow a light-response curve, which is why a similar adjustment was tested in GPP₂ (Eq. 8). At first glance, this structural similarity likely explains why the adjustment enhanced performance (Sect. 6.1). However, in other cases, GPP_EC is not directly modeled, but derived as the residual between NEE_EC and ER_EC. Therefore, it is difficult to determine whether improvements from GPP₂ reflect better reproduction of specific flux-partitioning algorithms, a more accurate representation of the true flux dynamics, or a combination of both. In contrast, GDD is not used at all in flux-partitioning algorithms (Appendix A). Therefore, the improvements resulting from the inclusion of GDD in GPP₄ may capture true ecosystem state changes that are also well represented in GPP_EC.

Regarding ER, the situation is more complex. ER_EC is typically derived from a fitted power-based or exponential-based function (Appendix A). These functions depend solely on temperature and estimate the combined contribution of AR and HR as a single inseparable flux. In contrast, the L4C model and the tested AS-adapted formulations explicitly represent ER as the sum of AR and HR, with each component estimated separately using multiple drivers, including APAR, GDD, MNT, VPD, RZSM, ST, and SSM. This approach relies on assumed linkages between GPP and AR, and between GPP, L_fall, SOC, and HR (Kimball et al., 2008), resulting in a more mechanistic, interaction-rich framework than the flux-partitioning algorithms. Consequently, calibrating ER formulations is challenging, because the reference ER_EC is obtained using a simpler empirical approach, which may limit model performance. If the ultimate goal is to estimate the CO₂ budget accurately, rather than to predict the underlying GPP and ER components, it may be advantageous to calibrate the L4C model using NEE_EC as the reference, rather than relying on GPP_EC and ER_EC as intermediate targets. However, this approach prevents validating whether the modeled GPP and ER truly reflect the underlying processes and strongly limits the number of free parameters that can be estimated, since only a single reference is available instead of two. For future research, it could also be valuable to partition NEE_EC into GPP_EC and ER_EC using a more mechanistic approach similar to the L4C model, explicitly distinguishing between AR and HR.

In summary, when calibrating TCF models using GPP_EC and ER_EC as references, one attempts to explain variability in fluxes that originate from a reference framework with a relatively simple structure, a limited number of drivers, and parameters that may vary in space and time (Appendix A). In contrast, TCF models, such as the L4C model, rely on a more complex process representation, a larger set of environmental drivers, and parameters that are assumed to be constant in space and time within a given ecosystem. These fundamental differences inherently complicate model calibration, hinder the interpretation of model performance, and limit our ability to determine whether the underlying processes of ER and GPP are realistically represented when extrapolated to larger spatial and temporal scales.

Several studies have also shown that GPP responds to the ratio of leaf-internal to ambient CO₂ concentration (Wang et al., 2014b, 2017). Although this ratio is regulated by environmental conditions such as temperature and VPD, neither the original L4C model nor the tested AS-adapted formulations and flux-partitioning algorithms explicitly accounts for the response of GPP to changes in ambient CO₂ concentration. Because ambient CO₂ varies over time and may continue to increase in the future, this omission may limit the ability of the L4C model to accurately predict GPP over long temporal scales.

Finally, it is important to note that this study focuses exclusively on the growing season, whereas the ultimate objective of improving the L4C model for the North American AS is to better estimate the full annual CO₂ budget over recent years by integrating modeled NEE since 2015. Although CO₂ flux magnitudes are highest during the growing season, GPP and AR become minimal or absent during the shoulder seasons and winter, while HR persists, even under snow-covered and frozen soil conditions. Several studies have shown that the winter and shoulder seasons play a critical role in shaping the annual CO₂ budget (Kim et al., 2013; Natali et al., 2019), as they primarily constitute a CO₂ source due to the dominance of HR. Therefore, future work will focus on improving the L4C model for these periods to provide more reliable year-round estimates of NEE, GPP, and ER, as well as more representative annual CO₂ budgets.

7 Conclusions

The goal of this study was to refine the integration of energy and moisture proxies into the SMAP L4C GPP and ER modeling for the North American AS growing season. To this end, alternative GPP and ER model formulations were calibrated and evaluated against GPP_EC and ER_EC across upland tundra, taiga forests and wetlands, covering the period from 2015 to 2022.

Ultimately, we recommend two key adjustments related to energy proxies to enhance the L4C model ability to monitor the GPP process:

• Implementing a nonlinear light-response, particularly to reduce ubRMSE and B;

• Incorporating GDD to reflect vegetation green-up and senescence phases, thereby improving seasonal dynamics.

In contrast, model adjustments related to moisture proxies (VPD, RZSM, SSM) for both ER and GPP modeling do not currently emerge as essential for future operational implementation. Moreover, evaluating the benefits of integrating SOC dynamics into ER modeling remains challenging, even though the L4C version 8 approach represents an improvement over that of version 7, and therefore further research into ER modeling is recommended.

Overall, combining the best-performing ER and GPP formulations yields the most, or among the most, accurate NEE configurations. However, improvements in NEE are generally more modest than for the individual components.

We further encourage the scientific community to harmonize strategies between flux-partitioning algorithms and mechanistic modeling frameworks (such as the L4C model), particularly for estimating ER and its underlying HR and AR components. The alignment between partitioning and modeling frameworks is essential to enhance the reliability of spatial and temporal extrapolation of GPP_EC and ER_EC using satellite-based TCF models.

Finally, future work will aim to improve the representation of shoulder and winter seasons in the L4C model, as these periods are critical for accurately capturing year-round CO₂ fluxes and annual CO₂ budgets in North American AS regions.

Appendix A: EC tower measurements of NEE and derived GPP and ER

EC tower measurements of NEE are subject to systematic errors, which mostly arise from unmet assumptions, instrument design and calibration, physical phenomena (e.g. storage terms), and terrain-specific conditions. These errors are generally well characterized and are typically corrected using software, such as EddyPro, as part of the standard flux processing workflow (Aubinet et al., 2012; Burba, 2022). NEE measurements are also affected by random errors, notably turbulence sampling error, which arises when large eddies are not adequately captured within a 30 min window (Finkelstein and Sims, 2001). The standard deviation of this error tends to follow a consistent pattern across ecosystem types and increases linearly with the flux magnitude (Aubinet et al., 2012). Overall, random errors in NEE are difficult to quantify, but using simultaneous measurements from two collocated EC towers, they have been estimated at 15 % for a 30 min interval (Eugster et al., 1997; Dragoni et al., 2007).

Gaps in NEE are usually filled using the Marginal Distribution Sampling (MDS) method (Falge et al., 2001; Reichstein et al., 2005; Aubinet et al., 2012). In this approach, a missing measurement is replaced by the mean of valid values observed under similar meteorological conditions within a ±7 d time window. Meteorological conditions are considered similar when variables such as solar radiation, temperature, and vapor pressure deficit do not deviate beyond predefined thresholds within the time window. When no suitable meteorological analogues are available, missing measurements are replaced by the mean of valid values from adjacent days at the same time of day (±1 h).

There is no universal flux-partitioning algorithm to separate NEE into its ER and GPP components. The most established method assumes that nighttime NEE consists solely of the ER component, since photosynthesis, and therefore GPP, is considered negligible in the absence of light (Reichstein et al., 2005; Aubinet et al., 2012). Nighttime ER is modeled using an exponential function, with air or soil temperature as the primary driver (Lloyd and Taylor, 1994). Air temperature is generally preferred because it better represents the landscape surrounding the EC tower, whereas soil temperature varies spatially and with depth across heterogeneous terrain (Helbig et al., 2017a, b). Daytime ER is then extrapolated to isolate the GPP contribution from the NEE measurements. Alternative approaches fit a light-response curve combined with a Q₁₀ equation to NEE measurements, accounting for the effects of photosynthetically active radiation (PAR) on GPP and air temperature on ER (Falge et al., 2001; Gilmanov et al., 2003; Lasslop et al., 2010b; Runkle et al., 2013; Helbig et al., 2017a). ER and GPP respond to more than temperature and light, respectively, and depend on other drivers such as SOC availability or moisture limitations. However, establishing a model formulation with more variables is complex and such additional measurements are not systematically available (Aubinet et al., 2012). To overcome this problem, the model regression is generally performed over short time intervals of 4 through 15 d (Reichstein et al., 2005; Lasslop et al., 2010b; Runkle et al., 2013). This approach enables the effects of additional drivers to be implicitly incorporated into the estimation of free parameters and accounts for seasonal parameter variability, reflecting changes in ecosystem state that are not represented in the model (Lasslop et al., 2010b; Aubinet et al., 2012). As GPP and ER are modeled using additional data and rely on various assumptions, they have greater uncertainties than tower measurements of NEE (Lasslop et al., 2010a).

Appendix B: Additional figures

Figure B1 Behavior of the ecosystem response functions used in the L4C model and the Arctic–Subarctic adapted formulations. Each panel shows the same response function evaluated across a range of parameter values to illustrate how parameterization affects function shape. Refer to Sects. 3 and 4.1 for a detailed description of each model formulation. Panels (A1), (B1), and (C1) correspond to Eq. (7a, c, d), (7b), and (7e), respectively; panel (D1) corresponds to Eq. (8b). Panels (A2) and (B2) correspond to Eq. (9a, c) and (9b), respectively; panel (C2) corresponds to Eqs. (10b) and (11).

Figure B2 Mean seasonal cycle of daily averaged eddy covariance (EC) gross primary production (GPP_EC) for upland tundra, taiga forests and wetlands. The horizontal red line indicates the estimated separation between dormant and transitional periods with active photosynthetic periods (i.e., the growing season) using a percentile-based GPP_EC threshold (Sect. 4.2). This threshold corresponds to the mean 10th percentile computed across EC towers, with the 5th and 15th percentiles shown as variability bounds.

Figure B3 Mean seasonal cycle of daily averaged gross primary production (GPP), ecosystem respiration (ER), and net ecosystem CO₂ exchange (NEE) across eddy covariance tower sites for upland tundra, taiga forests, and wetlands. GPP_EC, ER_EC, and NEE_EC denote eddy covariance data. GPP_L4C, ER_L4C, and NEE_L4C represent outputs from the original L4C model, whereas GPP₄, ER₄, and NEE_4,4 (i.e., ER₄ − GPP₄) correspond to outputs from the Arctic–Subarctic adapted formulations using test case 4 for both GPP and ER (Sects. 3 and 4.1). The data shown represent the growing season after applying the filtering criteria described in Sect. 4.2.

Appendix C: Additional tables for GPP evaluation

Table C1 Median differences in Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) between GPP₁ through GPP₅ and GPP_L4C. GPP_L4C refers to GPP outputs from the original L4C model, while GPP₁ through GPP₅ represent the five Arctic–Subarctic (AS) adapted formulations (Sects. 3 and 4.1).

Table C2 Same as Table 4, but here GPP₄ and GPP₅ use growing degree days (GDD) as inputs, which were normalized by the long-term minimum and maximum across all years rather than by annual values (Sect. 6.1).

Table C3 Same as Table 4, but instead of 100 runs with 70 %–30 % training–testing splits, the data were split into 6-year training and 2-year testing periods. The period of study spans a total of 8 years, yielding 28 possible runs instead of 100.

Table C4 Same as Table 4, but the 70 %–30 % training-testing splits were applied to the set of eddy covariance (EC) tower sites rather than to the individual data points, such that all data from a given site were assigned to either the training or the testing split. Because the number of EC tower sites differs among ecosystems, this yielded 5, 21, and 28 possible runs for upland tundra, taiga forests, and wetlands, respectively, instead of 100.

Appendix D: Additional tables for ER evaluation

Table D1 Same as Table 4, but for ecosystem respiration (ER). GPP₁ was used as the GPP input for ER₁ through ER₅.

Table D2 Same as Table 4, but for ecosystem respiration (ER). GPP₂ was used as the GPP input for ER₁ through ER₅.

Table D3 Same as Table 4, but for ecosystem respiration (ER). GPP₃ was used as the GPP input for ER₁ through ER₅.

Table D4 Same as Table 4, but for ecosystem respiration (ER). GPP₄ was used as the GPP input for ER₁ through ER₅.

Table D5 Same as Table 4, but for ecosystem respiration (ER). GPP₅ was used as the GPP input for ER₁ through ER₅.

Table D6 Same as Table C1, but for ecosystem respiration (ER).

Table D7 Same as Table 5, but instead of 100 runs with 70 %–30 % training–testing splits, the data were split into 6-year training and 2-year testing periods. The period of study spans a total of 8 years, yielding 28 possible runs instead of 100.

Table D8 Same as Table 5, but the 70 %–30 % training-testing splits were applied to the set of eddy covariance (EC) tower sites rather than to the individual data points, such that all data from a given site were assigned to either the training or the testing split. Because the number of EC tower sites differs among ecosystems, this yielded 5, 21, and 28 possible runs for upland tundra, taiga forests, and wetlands, respectively, instead of 100.

Appendix E: Additional tables for NEE evaluation

Table E1 Same as Table 4, but for net ecosystem CO₂ exchange (NEE). Modeled NEE is computed as the difference between modeled ecosystem respiration (ER) and gross primary production (GPP). NEE_L4C denotes NEE from the original L4C model, while NEE_i, 1 refers to outputs from the Arctic–Subarctic (AS) adapted ER_i and GPP₁ formulations.

Table E2 Same as Table 4, but for net ecosystem CO₂ exchange (NEE). Modeled NEE is computed as the difference between modeled ecosystem respiration (ER) and gross primary production (GPP). NEE_L4C denotes NEE from the original L4C model, while NEE_i, 2 refers to outputs from the Arctic–Subarctic (AS) adapted ER_i and GPP₂ formulations.

Table E3 Same as Table 4, but for net ecosystem CO₂ exchange (NEE). Modeled NEE is computed as the difference between modeled ecosystem respiration (ER) and gross primary production (GPP). NEE_L4C denotes NEE from the original L4C model, while NEE_i, 3 refers to outputs from the Arctic–Subarctic (AS) adapted ER_i and GPP₃ formulations.

Table E4 Same as Table 4, but for net ecosystem CO₂ exchange (NEE). Modeled NEE is computed as the difference between modeled ecosystem respiration (ER) and gross primary production (GPP). NEE_L4C denotes NEE from the original L4C model, while NEE_i, 4 refers to outputs from the Arctic–Subarctic (AS) adapted ER_i and GPP₄ formulations.

Table E5 Same as Table 4, but for net ecosystem CO₂ exchange (NEE). Modeled NEE is computed as the difference between modeled ecosystem respiration (ER) and gross primary production (GPP). NEE_L4C denotes NEE from the original L4C model, while NEE_i, 5 refers to outputs from the Arctic–Subarctic (AS) adapted ER_i and GPP₅ formulations.