Improved estimation of net ecosystem CO₂ exchange over North America using LSTM-based flux upscaling (2001–2021)

Accurate estimation of regional-scale terrestrial carbon budgets is of great importance but remains challenging. With particular advantages, the Long Short-Term Memory (LSTM) network method shows potential for improving regional carbon budget upscaling estimations. Here, using LSTM, we upscaled regional net ecosystem carbon exchange (NEE) with available flux tower measurements and satellite land surface observations in North America. With well-established ecosystem-specific LSTMs, we produced monthly NEE at a spatial resolution of 0.1°×0.1° over 2001–2021 (labelled as MemoryFlux). Unlike existing upscaling estimates, our dataset properly identified the Midwest Corn Belt as a region of large seasonal carbon uptake during the peak growing season, a feature revealed by previous top-down studies and recognized as a model benchmark. Moreover, the seasonal variations in NEE estimated by MemoryFlux were strongly correlated with those from independent atmospheric inversions, including the ensemble mean of the Orbiting Carbon Observatory-2 Model Intercomparison Project (OCO-2 v10 MIP; r=0.96, p<0.001) and CarbonTracker2022 (CT2022) (r=0.97, p<0.001). The mean annual NEE was estimated at −1.27 ± 0.12 Pg C yr⁻¹, which was closer in magnitude to inversions (−0.83 to −0.70 Pg C yr⁻¹) than existing upscaling estimates (−3.30 to −1.68 Pg C yr⁻¹). In addition, MemoryFlux captured spatial NEE anomaly patterns associated with six selected severe drought and flood events. We further found that explicitly incorporating historical predictor information improved the representation of NEE interannual variability and spatial anomalies associated with climate extremes. MemoryFlux provides an improved bottom-up estimation of North American NEE and shows greater consistency in magnitude with independent atmospheric inversion estimates than several existing EC-based upscaling products. The MemoryFlux dataset is available at https://doi.org/10.5281/zenodo.20482274 (Huang and He, 2026).

1 Introduction

Terrestrial ecosystems, as important carbon sinks on Earth, absorb approximately one fourth of the carbon dioxide (CO₂) released by human activities, which significantly slows down the rising CO₂ in the atmosphere and the pace of global warming (Friedlingstein et al., 2022). Net ecosystem CO₂ exchange (NEE), defined as the difference between the CO₂ released by total ecosystem respiration and the CO₂ sequestered by photosynthesis (Ballantyne et al., 2021), indicates the current net carbon balance of terrestrial ecosystems. Climate extremes, e.g., droughts, floods, heatwaves, and wildfires, strongly impact terrestrial CO₂ exchange processes, leading to non-negligible interannual variations (IAVs) in ecosystem carbon budgets (Le Quéré et al., 2018). Accurate quantification of NEE at regional and global scales is essential for advancing our understanding of carbon-climate feedbacks (Yao et al., 2018).

NEE is only measurable at the plot scale with the eddy covariance (EC) technique (Wofsy et al., 1993), and cannot be observed directly over large scales. Feasible ways to estimate large-scale NEE include top-down inversion of atmospheric CO₂concentration measurements and bottom-up estimation using process-based terrestrial biosphere models (TBMs) and data-driven machine learning (ML) upscaling. In recent years, the availability of atmospheric CO₂ observations has been greatly enhanced by space-based observations (Byrne et al., 2019), making it possible to infer the large-scale space-time distribution of carbon fluxes by using transport models and atmospheric CO₂ data. Compared to bottom-up estimates, top-down inversions constrained by atmospheric CO₂ data provide an independent perspective on regional carbon exchange and a seasonal cycle of NEE that is more consistent with tower observations by assimilating the loss of CO₂ from various carbon release and transport pathways (Byrne et al., 2018). However, confined by sampling density, current atmospheric CO₂ inversions only provide reliable estimates at relatively coarse spatial resolutions (Byrne et al., 2023). Process-based TBMs, which explicitly take the physical processes of energy, carbon, and water cycles into account, have become important tools for studying ecosystem carbon budgets (White et al., 2000). Nevertheless, TBM-based flux estimates are subject to substantial uncertainties arising from model structure, parameterization, and forcing data (Foster et al., 2024; Huntzinger et al., 2012), and lack effective constraints from multi-source observations (Ciais et al., 2022). Recently, the rapid development of ML approaches has paved new avenues for modelling the terrestrial carbon cycle. Data-driven ML approaches are straightforward and efficient for assessing NEE since they allow for easy integration of data and full use of information from various sources (Xiao et al., 2008, 2011, 2014). These methods derive the NEE estimates through spatial extrapolation without making any assumptions about ecosystem patterns (Jung et al., 2019; Xiao et al., 2014).

However, compared with top-down inversions, which assimilate various carbon flux transport pathways, large-scale bottom-up ML-based NEE estimates contain considerable uncertainty, primarily due to the difficulty of capturing the complex and highly nonlinear ecosystem responses to disturbances and climate extremes (Vicca et al., 2014). Characterizing the potential effects of climate variability on terrestrial carbon cycling remains challenging because the spatiotemporal dynamics and underlying response mechanisms are not yet fully understood. Observations suggest that ecosystem CO₂ flux responses are influenced not only by contemporary environmental conditions but also by memory effects of previous climate variability, disturbances, and their interactions (Monger et al., 2015). Such memory effects can operate across multiple time scales through dynamic interactions between biomes and their environment (Kraft et al., 2019), because terrestrial ecosystems may require months to years to recover from disturbances associated with climate extremes (Ogle et al., 2015). These temporal dependencies pose challenges to the advancement of research on the carbon cycle and climate change. Incorporating information on antecedent vegetation states and environmental conditions can therefore improve the representation of current ecosystem response to climate variability at large spatial scales (Seddon et al., 2016). Moreover, the dependence of carbon cycle dynamics on climate variability is more complex, as the IAVs in temperature and precipitation can substantially influence carbon flux estimates in terrestrial ecosystems (Piao et al., 2020).

To understand the memory effects of climate extremes on ecosystem carbon fluxes, a considerable amount of research has been conducted (Aubinet et al., 2018; Liu et al., 2018; Shen et al., 2016; Wu et al., 2015). Generally, the link between NEE IAVs and the memory effects of these climate extremes has not been sufficiently characterized by ML-based models. ML methods typically do not directly consider carbon fluxes and stocks at previous time steps, while reflecting the effects of past disturbances to some extent through the use of certain independent datasets (e.g., vegetation indices) (Tao et al., 2020). Thus, these approaches cannot fully represent the inherent temporal interactions of ecosystem carbon cycle processes (Papagiannopoulou et al., 2017). Poor understanding and inadequate representation of memory effects result in substantial uncertainty in estimating regional carbon budgets and considerable differences between top-down inversions and bottom-up estimates (Deng et al., 2013). Consistency between bottom-up estimates and independent atmospheric inversions across fundamentally different methodologies is nevertheless a key way to increase the confidence in quantifying terrestrial ecosystem carbon fluxes at large scales. Therefore, developing the capacity of bottom-up models to explicitly incorporate historical information contained in climate and vegetation time series is urgently needed to better represent temporal variability in ecosystem carbon fluxes.

The Long Short-Term Memory (LSTM) network is a temporal dynamic statistical approach (Besnard et al., 2019; Hochreiter and Schmidhuber, 1997), which effectively addresses long-term dependence issues and thus can incorporate the memory effects of climate change on the terrestrial carbon cycle. LSTM has exhibited an outstanding power in characterizing the memory effects of climate and vegetation on NEE and improving the overall prediction accuracy at the site level (Besnard et al., 2019; Huang et al., 2024; Liu et al., 2023; Ma et al., 2026). As a pioneer such study over large scales, Liu et al. (2023) proved the importance of the memory effects for assessing IAVs in NEE by conducting an upscaling estimation of global gridded NEE using LSTM. In our previous work, we established plant functional type (PFT)-specific LSTM models for NEE prediction using flux tower measurements from more than 80 sites and multiple satellite land surface datasets over North America, and found that the LSTM framework provided better site-level predictive performance and representation of NEE IAVs than a Random Forest model using contemporaneous predictors (Huang et al., 2024). Yet, the ability of LSTM models to characterize regional spatiotemporal NEE patterns, particularly with the aid of multi-source satellite land surface data, has not been sufficiently examined. For example, Liu et al. (2023) used only one satellite vegetation variable, the normalized difference vegetation index (NDVI), and did not investigate the performance of the upscaling result in capturing the impacts of large-scale climate extremes.

In this study, we focus on North America, a region with the most extensive eddy flux observations globally and a relatively well-studied carbon budget. In particular, the substantial carbon uptake in the Midwest Corn Belt has been widely reported as a continental-scale phenomenon since the start of the CarbonTracker inversion (Peters et al., 2007) and during the Mid-Continent Intensive campaign (Schuh et al., 2013), and it has been revisited as a model benchmark feature. However, existing global flux upscaling datasets generally fail to capture this phenomenon. To upscale NEE, we employ an LSTM DL algorithm that fully utilizes in situ CO₂ flux observations in combination with diverse remote sensing and climate variables. From this, we generate a long-term gridded NEE dataset at a spatial resolution of 0.1° and monthly intervals from 2001 to 2021. The resulting dataset, named MemoryFlux, is publicly available for research purposes. We aim to answer a key scientific question: Does explicitly incorporating historical predictor information and integrating multi-source satellite data improve the representation of North American NEE? An ensemble of flux upscaling datasets and atmospheric inversions is used for comparison. We evaluate the consistency and differences among MemoryFlux, previously reported empirical knowledge, existing bottom-up products, and independent top-down inversion estimates, with particular attention to seasonal, interannual, and spatial variability. We also assess whether MemoryFlux shows greater consistency with inversions than existing bottom-up upscaling products and whether it provides a coherent representation of NEE anomalies associated with selected climate extremes. We anticipate that MemoryFlux, which explicitly incorporates historical ecosystem information within an LSTM framework, provides a robust representation of ecosystem carbon fluxes and serves as a valuable resource for investigating terrestrial carbon dynamics and ecosystem–climate interactions.

2 Data and methods

2.1 Datasets

2.1.1 CO₂ eddy covariance flux data

The study area encompassed the major terrestrial land-cover types across North America, extending from 5 to 80° N. We used observations from 84 flux sites (Table S1 in the Supplement), totalling 7711 monthly NEE records (Table S2), drawn from FLUXNET2015 (Pastorello et al., 2020) and the AmeriFlux network (Novick et al., 2018). These observations were combined with various remote sensing and climate reanalysis datasets. The target variable was monthly average NEE (NEE_VUT_REF; $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), with negative NEE values indicating net carbon uptake by terrestrial ecosystems and positive values indicating net carbon release to the atmosphere. The database covered 10 major plant functional types (PFTs) in North America (Fig. 1), following the International Geosphere-Biosphere Programme (IGBP) classification: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), mixed forest (MF), cropland (CRO), grassland (GRA), woody savanna (WSA), open shrubland (OSH), closed shrubland (CSH), savanna (SAV), and permanent wetland (WET).

Figure 1 Distribution of the eddy covariance flux towers used in this study. The background information reveals the land cover types across North America, derived from the MCD12C1 dataset. MODIS Land Cover Type Product (MCD12C1) © NASA and U.S. Geological Survey.

2.1.2 Climate reanalysis data

Information about the predictive features used in this study for upscaling NEE is listed in Table 1. These variables were selected to correspond as closely as possible to the climate and environmental drivers used at the site level (Huang et al., 2024). We extracted seven variables from ERA5-Land, including downward shortwave radiation (DSR), air temperature (T_a), dew point temperature (T_d), the 10 m u-component of wind (U₁₀), the 10 m v-component of wind (V₁₀), soil water content in the first soil layer (0–7 cm; SWC), and precipitation (P). ERA5-Land is the land component of the fifth generation of European Reanalysis (ERA5) (Muñoz-Sabater et al., 2021), developed by the Copernicus Climate Change Service at the European Centre for Medium-Range Weather Forecasts. It integrates various data sources, including both remote sensing observations and ground-based measurements, ultimately providing monthly data at 0.1° resolution.

Table 1 Basic information about the predictive datasets used in this study. All gridded predictors were converted to the same units as their corresponding variables used for site-level model training.

The U₁₀ and V₁₀ components from ERA5-Land were used to calculate wind speed (WS), which served as a predictor of NEE:

$$\begin{array}{}\text{(1)}& \text{WS}=\sqrt{(U_{\mathrm{10}}{)}^{\mathrm{2}}+(V_{\mathrm{10}}{)}^{\mathrm{2}}}\end{array}$$

According to Eqs. (2)–(6), vapor pressure deficit (VPD; hPa) was calculated using the T_a and T_d retrieved from ERA5-Land:

$$\begin{array}{}\text{(2)}& \text{VPD}=\text{SVP}-\text{AVP}\end{array}$$

where SVP and AVP represent saturated and actual vapor pressure (hPa), respectively. These values can be computed using the following equations:

$$\begin{array}{}\text{(3)}& {\displaystyle }\text{SVP}=\mathrm{6.112}\times f_w\times \exp \left({\displaystyle \frac{\mathrm{17.67}{T}_{\mathrm{a}}}{T_{\mathrm{a}}+\mathrm{243.5}}}\right)\text{(4)}& {\displaystyle }\text{AVP}=\mathrm{6.112}\times f_w\times \exp \left({\displaystyle \frac{\mathrm{17.67}{T}_{\mathrm{d}}}{T_{\mathrm{d}}+\mathrm{243.5}}}\right)\text{(5)}& {\displaystyle }{f}_w=\mathrm{1}+\mathrm{7}\times {\mathrm{10}}^{-\mathrm{4}}+\mathrm{3.46}\times {\mathrm{10}}^{-\mathrm{6}}{P}_{\text{mst}}\text{(6)}& {\displaystyle }{P}_{\text{mst}}=P_{\text{msl}}\times {\left({\displaystyle \frac{T_{\mathrm{a}}+\mathrm{273.16}}{T_{\mathrm{a}}+\mathrm{273.16}+\mathrm{0.0065}Z}}\right)}^{\mathrm{5.625}}\end{array}$$

where T_a is the air temperature (°C), while T_d denotes the dew point temperature (°C). P_mst and P_msl represent the air pressure (hPa) and the air pressure at the mean sea level (1013.25 hPa), respectively. Z represents the altitude (m) and was provided by the Advanced Land Observing Satellite (ALOS) (Tadono et al., 2016). The 30 m data were resampled to 0.1° to align with the resolution of the ERA5-Land reanalysis dataset. Thereby, monthly VPD values with consistent spatial and temporal coverage were generated.

2.1.3 Satellite data

We further selected remote sensing variables to represent the growth status of vegetation. NDVI was utilized in this study, which can indicate disturbance by displaying spectral browning signals linked to drought, harvesting, and fire (Myers-Smith et al., 2020). We employed the new global seamless 0.1°, monthly NDVI product, which was derived from Moderate Resolution Imaging Spectroradiometer (MODIS) surface reflectance data using an LSTM network method (Xiong et al., 2023). This product provides long-term information on global vegetation dynamics.

Two MODIS-based Global Land Surface Satellite (GLASS) product suites (Liang et al., 2021) were used: leaf area index (LAI) and fraction of absorbed photosynthetically active radiation (FAPAR). Both components have a 0.1° resolution and 8 d intervals, covering the period from 2000 to 2021. To match the monthly temporal resolution of NEE, we generated monthly composites by averaging the 8 d values. For 8 d windows spanning two months, the values were allocated to each month by weighting them according to the number of days within that month.

We also represented the space-time variations of NEE by incorporating novel satellite datasets. Solar-induced chlorophyll fluorescence (SIF) has brought considerable advancement in measuring terrestrial photosynthesis and is reported to have a strong correlation with NEE changes (Huang et al., 2024; Shiga et al., 2018). We employed the monthly Global Orbiting Carbon Observatory-2 (OCO-2) based SIF (GOSIF), a reconstructed SIF dataset with a high resolution (0.05°), developed using OCO-2 SIF, MODIS vegetation index, and MERRA-2 reanalysis data (Li and Xiao, 2019). Subsequently, we resampled the SIF data to a 0.1° resolution using the averaging approach.

For each 0.05° grid cell, the modal MCD12C1 land-cover class over 2001–2021 was designated as the representative class, such that interannual land-cover changes were not considered (Jung et al., 2019). This procedure produced a static PFT map and reduced potential noise associated with year-to-year classification error in the MODIS land-cover product. The resulting map was then resampled to a spatial resolution of 0.1° using the majority aggregation method.

2.1.4 NEE data products for comparison

Due to the sparse distribution of EC flux towers and the mismatch between product grids and tower footprints, large-scale evaluation of the upscaling NEE estimates required comparisons with published regional carbon flux products. Atmospheric inversions infer surface fluxes by combining atmospheric CO₂ observations with atmospheric transport models. This makes them useful for characterizing large-scale seasonal and interannual variations, and they therefore serve as independent references for evaluating bottom-up products. We selected two monthly atmospheric inversions for comparison: the OCO-2 Model Intercomparison Project (MIP) (Byrne et al., 2023) and CarbonTracker2022 (Peters et al., 2007). For OCO-2 v10 MIP, we used the ensemble-mean posterior land flux from the “LNLGIS” experiment (“land”), which assimilates OCO-2 Land Nadir and Land Glint XCO₂ retrievals together with in situ CO₂ measurements and provides 1° estimates for 2015–2020. To improve comparability with EC-based ecosystem NEE, we removed the prescribed wildfire-emission component using the National Oceanic and Atmospheric Administration (NOAA) CarbonTracker “fire_flux_imp” estimate. We also included CarbonTracker2022 (CT2022), the latest version of CarbonTracker developed by NOAA (Jacobson et al., 2023; Peters et al., 2007). For CT2022, we used the optimized terrestrial biosphere flux (“bio_flux_opt”), which provides 1° global land–atmosphere CO₂ exchange for 2001–2020, with wildfire emissions represented separately from the optimized biosphere flux. Importantly, closer agreement between MemoryFlux and the atmospheric inversions was not interpreted as direct evidence of higher accuracy. Instead, the inversions were treated as independent references based on fundamentally different methodologies, providing complementary information for assessing the consistency of regional NEE estimates in the absence of large-scale direct observations.

We also compared MemoryFlux with representative ML-based bottom-up NEE products at regional and global scales. At the regional scale, we used the EC-MOD2012 dataset for North America, which provided NEE estimates at 1 km spatial and 8 d temporal resolution for 2000–2012 and was derived from EC flux measurements and MODIS observations (Xiao et al., 2014). At the global scale, we included the FLUXCOM 2020 ensemble products (RS+CRUJRA1.1, 2001–2017; RS+ERA5, 2001–2018) developed by the Max Planck Institute for Biogeochemistry (MPI-BGC) (Jung et al., 2019). These products provide 0.5° monthly NEE estimates with broad spatial coverage and continuous temporal records. We also included FLUXCOM-X, the latest generation of FLUXCOM products, which provided global 0.05° NEE estimates (2001–2021) with improved methodology and data integration (Nelson et al., 2024). In addition, we used NIES2020 (1999–2019, 0.1° every 10 d) from the National Institute for Environmental Studies, Japan. Unlike FLUXCOM, which used static PFT information, NIES2020 derived vegetation-type information from LAI-related variables, thereby improving the representation of spatial heterogeneity and enhancing model generalization across ecosystems (Zeng et al., 2020).

Together, these complementary datasets provided a broad comparison framework, covering regional to global upscaling products and atmospheric inversion results. Comparisons with these datasets allowed us to assess the consistency and differences among MemoryFlux, previously reported regional NEE patterns, existing bottom-up products, and independent top-down inversion estimates, as well as to evaluate how effectively MemoryFlux represents regional-scale seasonal, interannual, and extreme-event-related variability in NEE.

2.2 The LSTM model for NEE upscaling

Conventional non-recurrent models generally treat observations at each time step independently and do not contain a dedicated architecture for explicitly representing temporal dependencies between current NEE and antecedent predictor conditions. The Long Short-Term Memory (LSTM) network (Hochreiter and Schmidhuber, 1997), a dynamic statistical model, can use previous variable data of a certain look-back window as input to predict the NEE output at the current time. NEE is a variable influenced by multiple factors, such as climate and environment. Considering that this influence changes over time, it is essential to predict NEE using historical data that capture its response to climate and disturbance conditions over specific periods.

In our previous work, we used an LSTM model to simulate NEE and its IAVs at the site level by linking monthly CO₂ fluxes with environmental and vegetation predictors (Huang et al., 2024). The model was trained with a six-month look-back window to predict current NEE, enabling it to capture temporal dependencies in vegetation dynamics and meteorological conditions. This window length corresponds to the empirical memory effects of most ecosystems reported by Liu et al. (2023) and was further validated to provide good predictive performance at the site scale (Huang et al., 2024; Ma et al., 2026). Key inputs for site-scale model training included meteorological and environmental variables measured at EC flux towers (TA_F, WS_F, P_F, DSR_F, SWC_F_MDS_#, and VPD_F), together with vegetation variables derived from remote sensing datasets (NDVI, LAI, FAPAR, and SIF). For SWC, the shallowest available SWC measurement at each site was used to maximize data availability. We used the quality-controlled and gap-filled variables provided by the flux networks, removed values flagged as −9999 as well as obvious outliers, and retained only sites with at least 18 consecutive valid monthly observations. For sites represented in both FLUXNET and AmeriFlux, the record with the longer and more complete time series was used. The LSTM models were trained using tower-observed meteorological variables rather than ERA5-Land data. This design was intended to make full use of co-located EC-site observations when establishing empirical relationships between NEE and its environmental drivers. We adopted an improved temporal cross-validation strategy to minimize overfitting and enhance model robustness (Fig. S1 in the Supplement). Specifically, the initial 70 % of data from each site were assigned to the training set, while the remaining 30 % were reserved as an independent testing set. During model training, we further employed nested ten-fold cross-validation to internally validate the training set and optimize model hyperparameters using a combination of early stopping and grid search (Bergstra and Bengio, 2012). This design effectively reduced the risk of data leakage and enhanced the model's generalization to unseen data while strictly maintaining the independence of the training, validation, and testing datasets. Loss curves (Fig. S2) confirmed the robustness of the LSTM models.

Building on this foundation, we applied the LSTM framework to regional carbon flux upscaling across North America. To reduce inconsistencies between site-training and regional prediction, ERA5-Land variables were harmonized with the corresponding EC-site variables in terms of physical definitions and units before being used as model inputs. Remote sensing and meteorological reanalysis gridded datasets were then processed into pixel-level predictor time series and fed into the corresponding PFT-specific LSTM models, yielding monthly NEE estimates at 0.1° spatial resolution. Predictions were generated for 10 PFTs (ENF, DBF, MF, CRO, GRA, WSA, OSH, CSH, SAV, and WET) represented in the site-level training dataset. For the SAV PFT, SWC was not included as a predictor because reliable site-level SWC observations were unavailable. To maintain consistency between site-level model training and regional upscaling, SWC was also excluded from the corresponding gridded SAV model. Despite the reduced predictor set, the SAV LSTM model still demonstrated strong predictive performance (R²=0.83; Fig. S3).

3 Results

3.1 Spatial patterns of NEE in North America

Figure 2 illustrated the spatial patterns of long-term mean NEE from MemoryFlux over 2001–2021. The long-term mean highlighted the dominant spatial distribution of NEE while minimizing the influence of interannual variations. Annually, the strongest carbon uptake was observed in the southeastern United States, followed by the northern Pacific Coast and tropical regions, whereas major carbon sources were concentrated in the tundra of northern and central Canada, as well as in the arid central U.S., where precipitation is limited. During the peak growing season (July–August), the Corn Belt in the Upper Midwest U.S. (Fig. 2b and c) emerged as the dominant carbon-uptake region. As the peak growing season progressed, the spatial extent of strong carbon uptake contracted in August, though it remained largely consistent with the North American Corn Belt. By the end of the peak growing season (September), the carbon uptake area further diminished (Fig. 2d). From July to September, CO₂ sources were mainly concentrated in the western U.S. and southwestern Canada, with both their spatial extent and magnitude increasing over time.

Figure 2 Spatial distribution of NEE from MemoryFlux dataset, (a) on an annual basis, (b, c) during the peak growing season, and (d) in September, averaged over 2001–2021 at 0.1° resolution across North America.

To further evaluate spatial patterns, we compared MemoryFlux with seven other NEE datasets (Fig. 3). At the annual scale, MemoryFlux and the other ML-based flux-upscaling products consistently identified the southeastern United States as the region of strongest carbon uptake, whereas the OCO-2 v10 MIP and CT2022 inversions showed stronger annual uptake centred over the Midwest Corn Belt. This discrepancy may reflect methodological differences and uncertainties inherent in both EC-based upscaling and atmospheric inversion approaches. During the peak growing season, MemoryFlux reproduced pronounced seasonal carbon uptake over the Midwest U.S. croplands, a feature also evident in OCO-2 v10 MIP, CT2022, and EC-MOD2012. In contrast, other global flux-upscaling products showed their strongest seasonal uptake in the southeastern United States. These results indicated that the seasonal and annual spatial patterns represented by MemoryFlux were more consistent with previous empirical evidence and with most published estimates. Latitudinal gradients were generally consistent across datasets, although the inversions exhibited weaker NEE variations with latitude (Fig. 3b and d). Furthermore, only the regional NEE estimates from MemoryFlux, EC-MOD2012, and the global FLUXCOM-X consistently identified annual CO₂ sources in northern and central Canada, as well as in the central U.S., and captured growing season CO₂ release in the western U.S. and southwestern Canada. Among these datasets, only the OCO-2 v10 MIP inversions indicated the northern Pacific Coast and tropical regions as annual carbon sources.

Figure 3 Spatial distribution of NEE in North America on an annual basis and during the peak growing season (July–August) derived from OCO-2 v10 MIP, CT2022, MemoryFlux, EC-MOD2012, NIES2020, FLUXCOM2020 RS+CRUJRA1.1, FLUXCOM2020 RS+ERA5, and FLUXCOM-X. The shades around the lines in (b) and (d) represent the range of standard deviation.

We further analysed the contributions of major ecoregions to annual total NEE using the vegetation-type classification of Olson et al. (2001), as implemented in the CarbonTracker project (Fig. S4). On average, MemoryFlux indicated that Forests/Wooded ecoregions served as the dominant carbon sink in North America, with a mean uptake of −0.76 Pg C yr⁻¹, whereas Grass/Shrubs contributed only marginally, with uptake close to zero. Although all products agreed on the leading role of Forests/Wooded ecosystems as carbon sinks, the magnitude of their contributions varied across datasets. MemoryFlux and FLUXCOM-X attributed more than 85 % of the total carbon sink to forests; EC-MOD2012 and the other global upscaling products estimated about 75 %, while OCO-2 v10 MIP and CT2022 reported substantially lower shares of 65 % and 38 %, respectively. Moreover, only a few products (EC-MOD2012 and FLUXCOM-X) suggested that Grass/Shrubs acted as a weak carbon source.

3.2 Seasonal variations of NEE in North America

To evaluate the seasonal variations represented by MemoryFlux, we compared it with other NEE datasets (Fig. 4 and Table S3). This comparison was limited to 2015–2017 when all products except EC-MOD2012 were available. The results revealed that MemoryFlux exhibited similar seasonal patterns to other NEE datasets, demonstrating carbon uptake during warm seasons and carbon release during cold seasons. The estimated NEE seasonal cycle matched well with the OCO-2 v10 MIP (RMSE=0.11 Pg C per month, MAE=0.09 Pg C per month, r=0.96, p<0.001) and CT2022 inversions (RMSE=0.13 Pg C per month, MAE=0.11 Pg C per month, r=0.97, p<0.001). However, it had weaker net carbon uptake in June, particularly relative to the global upscaling products. Notably, all upscaling NEE products exhibited earlier spring onset and later autumn senescence compared with atmospheric inversions. Furthermore, these global upscaling NEE products generally estimated stronger carbon uptake than both the atmospheric inversions and MemoryFlux, particularly during the growing seasons. NIES2020 exhibited the largest amplitude in the seasonal NEE cycle, particularly in July (Fig. 4a).

Figure 4 Comparison of the seasonal patterns of MemoryFlux estimated by LSTM models against other products. (a) NEE mean seasonal cycle variations (Pg C per month) in 2015–2017 derived from MemoryFlux (the shaded area which boundary is standard deviation, the grey line represents the NEE seasonal cycle), NIES2020, OCO-2 v10 MIP, CT2022, FLUXCOM2020 RS+CRUJRA1.1, FLUXCOM2020 RS+ERA5, and FLUXCOM-X. (b–d) Statistics of NEE over the season scale in whole, boreal, and temperate North America. Error bars indicate the standard deviation across years, representing interannual variability. Note that EC-MOD2012 was excluded from this comparison due to the data unavailability.

When we compared these NEE products across the four seasons, we observed consistent patterns but notable differences in magnitude (Fig. 4b and Table S4). In MAM (March–May) and JJA (June–August), most global upscaling products indicated stronger net carbon uptake than both the atmospheric CO₂ inversions and MemoryFlux. Even the recently updated FLUXCOM-X showed relatively strong JJA carbon uptake, second only to NIES2020. In contrast, MemoryFlux was closer in magnitude to the top-down inversions during MAM and JJA; for example, its JJA uptake was more similar to OCO-2 v10 MIP and CT2022 than to the estimates of any other upscaling product. During SON (September–November), FLUXCOM2020 RS+CRUJRA1.1 and NIES2020 indicated net carbon uptake over North America, whereas the atmospheric inversions, MemoryFlux, and the other globally upscaling estimates suggested net carbon release. In DJF (December–February), the FLUXCOM2020 ensembles and NIES2020 exhibited weaker carbon release than the other products. Consequently, MemoryFlux, which explicitly incorporated historical predictor information, showed greater consistency with atmospheric inversions than the other upscaling estimates.

We further compared seasonal NEE patterns between the boreal and temperate regions of North America based on the Transcom regional partitioning (Fig. 4c and d). Overall, the boreal region showed better agreement among products, whereas the larger discrepancies in carbon uptake between the global upscaling products and both MemoryFlux and the atmospheric inversions were mainly concentrated in the temperate region. During SON, the net carbon uptake reported by NIES2020 and FLUXCOM2020 RS+CRUJRA1.1 was mainly associated with the temperate region, while the boreal region remained a net carbon release. FLUXCOM-X exhibited net carbon uptake in MAM in both regions, whereas CT2022 and MemoryFlux showed this pattern only in the temperate region. Additionally, MemoryFlux estimated stronger carbon release in the boreal region during MAM and DJF, but weaker carbon release during SON, compared with the other products.

3.3 Annual total and interannual variations of NEE in North America

We first investigated annual total NEE and its IAVs across the eight datasets (Fig. 5 and Table S5). MemoryFlux estimated annual NEE ranging from −1.50 to −1.05 Pg C yr⁻¹ over 2001–2021, with a mean of −1.27 Pg C yr⁻¹ and an interannual standard deviation of 0.12 Pg C yr⁻¹. The atmospheric inversions yielded weaker mean carbon uptake (CT2022: −0.70 ± 0.20 Pg C yr⁻¹; OCO-2 v10 MIP: −0.83 ± 0.11 Pg C yr⁻¹), whereas the bottom-up upscaling products estimated stronger uptake (Table S5). We further compared annual NEE over common periods. For 2015–2020, mean annual NEE was −1.38 Pg C yr⁻¹ for MemoryFlux, −0.83 Pg C yr⁻¹ for OCO-2 v10 MIP, −0.69 Pg C yr⁻¹ for CT2022, and −1.75 Pg C yr⁻¹ for FLUXCOM-X. For 2001–2012, the corresponding estimates were −1.21 Pg C yr⁻¹ for MemoryFlux, −0.66 Pg C yr⁻¹ for CT2022, −1.81 Pg C yr⁻¹ for EC-MOD2012, −3.19 Pg C yr⁻¹ for NIES2020, −3.05 and −2.78 Pg C yr⁻¹ for the two FLUXCOM2020 ensemble products, and −1.65 Pg C yr⁻¹ for FLUXCOM-X. These common-period comparisons showed that MemoryFlux was closer in magnitude to the atmospheric inversions than several regional and global EC-based flux-upscaling products, including EC-MOD2012, NIES2020, the FLUXCOM2020 ensemble products, and FLUXCOM-X. Among all datasets, NIES2020 estimated the strongest carbon sink at −3.30 ± 0.22 Pg C yr⁻¹.

Figure 5 Annual total and interannual variations of North American NEE estimated by OCO-2 v10 MIP, CT2022, MemoryFlux, EC-MOD2012, NIES2020, FLUXCOM2020 RS+CRUJRA, FLUXCOM2020 RS+ERA5, and FLUXCOM-X. Dashed lines represent the long-term trend calculated using the Mann–Kendall test and Theil–Sen slope method. Due to the limited data, we did not estimate a slope for OCO-2 v10 MIP.

To further evaluate model performance, we calculated the correlation coefficients (r) for interannual variations in NEE between MemoryFlux and the other products. MemoryFlux showed positive correlations with OCO-2 v10 MIP (r=0.70, p=0.12), CT2022 (r=0.41, p=0.07), EC-MOD2012 (r=0.68, p<0.05), and FLUXCOM-X (r=0.69, p<0.05), whereas weaker correlations were found with the other datasets. Although the correlations with OCO-2 v10 MIP and CT2022 were not statistically significant, the overall comparison indicated that MemoryFlux, estimated using the LSTM DL approach, captured both the magnitude and IAVs of NEE reasonably well. Long-term NEE trends were assessed using the Mann–Kendall test and Theil–Sen slope method. Due to limited data, no slope was estimated for OCO-2 v10 MIP. MemoryFlux, CT2022, NIES2020, and FLUXCOM-X showed tendencies toward increasing carbon uptake (MemoryFlux: $k=-\mathrm{0.009}$ Pg C yr⁻², p=0.09; CT2022: $k=-\mathrm{0.009}$ Pg C yr⁻², p=0.35; NIES2020: $k=-\mathrm{0.032}$ Pg C yr⁻², p<0.05; FLUXCOM-X: $k=-\mathrm{0.006}$ Pg C yr⁻², p=0.19), whereas EC-MOD2012 and the FLUXCOM 2020 ensembles indicated a slight decrease, although most trends were not statistically significant. Regional analysis further found that the temperate region accounted for much of the inter-product divergence in annual carbon sink magnitude (Fig. S5). Compared with the atmospheric inversions and MemoryFlux, the upscaling products tended to indicate stronger carbon uptake in the temperate region, whereas the boreal region showed relatively good agreement across datasets.

To evaluate whether explicitly incorporating historical climate and environmental information improved the representation of NEE IAVs, we generated two LSTM-based NEE datasets, MemoryFlux and nonMemoryFlux. The two configurations were identical with respect to predictors, model architecture, data partitioning, and nested cross-validation, differing only in their temporal input structure. MemoryFlux updated its internal state using information from previous time steps, thereby modelling temporal dependencies and potential memory effects between antecedent environmental conditions and current NEE. In contrast, nonMemoryFlux was driven only by synchronous explanatory variables. Both configurations may retain some implicit legacy information through contemporaneous vegetation-state predictors (e.g., NDVI and LAI), whereas MemoryFlux additionally incorporated historical predictor sequences explicitly through its temporal look-back structure. Therefore, the comparison of the two datasets provided a controlled assessment of the added value of explicitly incorporating historical predictor information. We found that nonMemoryFlux reproduced the IAVs of NEE reasonably well (r=0.82, p<0.001), but clear differences from MemoryFlux remained (Fig. 6). In particular, nonMemoryFlux tended to exhibit a stronger increasing trend (−0.013 Pg C yr⁻²) and estimated stronger carbon uptake in recent years (since 2014), when the 2015/16 El Nino events and frequent droughts and heatwaves hit North America, suggesting that it failed to capture the legacy effects of climate extremes on terrestrial ecosystems. The maximum difference reached 186 Tg C yr⁻¹ in 2014 (14.6 % of the mean carbon uptake), and the mean absolute difference was 64 Tg C yr⁻¹ (5.0 %), highlighting the importance of accounting for memory effects in carbon flux upscaling.

Figure 6 Comparison of interannual variations in North American NEE estimated by MemoryFlux and nonMemoryFlux. MemoryFlux explicitly incorporates historical predictor information through a temporal look-back structure, whereas nonMemoryFlux uses only contemporaneous predictors. The dashed lines represent long-term trends calculated using the Mann–Kendall test and Theil–Sen slope method.

We further compared, at the PFT level, the predictive performance of the LSTM models used to generate MemoryFlux and nonMemoryFlux (Fig. S6). Overall, R² increased from 0.63 for nonMemoryFlux to 0.79 for MemoryFlux, while RMSE decreased from 1.12 to 0.83 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and MAE decreased from 0.69 to 0.52 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Improvements were observed across all PFTs, with particularly pronounced gains for ENF, for which R² increased by 0.39 (Tables S6 and S7). We also evaluated whether the temporal input structure improved the representation of seasonal and interannual variations in carbon fluxes at the site scale by comparing NEE estimates from MemoryFlux and nonMemoryFlux with long-term flux-tower observations across PFTs (Fig. S7). Clear differences were evident in their representation of temporal variability. MemoryFlux generally reproduced the observed seasonal cycle and IAVs more closely than nonMemoryFlux, particularly at lower NEE values, whereas nonMemoryFlux showed weaker IAVs. Together, these results suggested that explicitly incorporating historical predictor information improved both overall predictive performance and the representation of temporal variations in NEE.

3.4 Capturing climate extremes impacts on regional NEE

We further evaluated the ability of the LSTM-based model to represent NEE IAVs at the grid scale by examining its representation of spatial NEE anomalies associated with selected climate extremes. Given that the largest impacts typically occurred from June to August (summer, JJA), we focused our analysis on regional carbon flux anomalies during summer. Regional NEE anomalies were calculated as the difference between mean summer NEE during the event year and the corresponding multi-year mean summer NEE for each dataset. Soil moisture (SM) anomalies were used to characterize the severity and spatial extent of drought and flood conditions. The Global Land Evaporation Amsterdam Model (GLEAM) v3.7a dataset at 0.25° resolution was employed for this study (Martens et al., 2017). Furthermore, regional anomalies in SIF from the GOSIF dataset were utilized to provide information on concurrent changes in vegetation activity. The SIF and SM datasets were used as supporting ecological and hydroclimatic indicators and were not intended as direct proxies or independent observations of NEE. In the analysis, we chose six events, including the 2011 drought in Mexico and the southern U.S. (Dobler-Morales and Bocco, 2021; Weiss et al., 2012), the 2012 North American drought (Boyer et al., 2013), the 2013 California drought (Swain et al., 2014), the 2017 flash drought across the northern U.S. and Canadian Prairie (Hoell et al., 2020), the 2019 Midwest American floods (Yin et al., 2020), and the southwestern North American drought and associated wildfires during 2020–2021 (Mankin et al., 2021). Across these climate extremes, drought and flood conditions caused negative and positive SM anomalies, respectively, and negative SIF anomalies, leading to anomalous CO₂ release from land to the atmosphere, as reflected by positive NEE anomalies.

Figure 7 showed the areas affected by selected climate extremes across North America during JJA 2011–2013, JJA 2017, JJA 2019, and JJA 2021, based on multiple datasets. Black rectangles indicated the regions most strongly affected by the selected extreme climate events. For the 2020–2021 drought in southwestern North America, the NEE anomaly was evaluated for JJA 2021 relative to the corresponding climatological baseline, as this period represents a more mature and severe phase of the event. In arid regions, MemoryFlux exhibited positive NEE anomalies during JJA that spatially coincided with negative SM and SIF anomalies indicated by the environmental proxy datasets. During the 2019 Midwest American floods, SM showed positive anomalies (highlighted by the red rectangle and the reversed colour bar), which spatially coincided with negative SIF and positive NEE anomalies. The spatial extent of the NEE anomalies represented by MemoryFlux was also broadly consistent with the reported extent of the 2019 flood, with widespread flooding across the Midwest and more limited impacts farther south (Fig. S8). The spatial correspondence among MemoryFlux NEE anomalies, environmental proxy anomalies, and the documented extent of the selected extreme events suggested that MemoryFlux provided a coherent representation of regional NEE anomaly patterns associated with these events. We further examined the added value of explicitly incorporating historical predictor information by comparing MemoryFlux with nonMemoryFlux. The nonMemoryFlux model showed weaker spatial correspondence with the anomaly patterns associated with several extreme events, including the 2011–2012 droughts and the 2019 flood, whereas MemoryFlux represented these regional anomaly patterns more clearly (Fig. 7). These results suggested that explicitly incorporating historical predictor information can improve the representation of NEE anomalies associated with climate extremes.

Figure 7 Spatial anomalies of NEE derived from MemoryFlux and nonMemoryFlux, together with environmental proxy anomalies from GOSIF and GLEAM, during JJA for selected major drought and flood events in North America over the study period. The third row shows the difference between the MemoryFlux and nonMemoryFlux NEE anomalies. Black rectangles indicate the regions most strongly affected by the selected extreme climate events.

To further assess how MemoryFlux represents NEE anomalies associated with climate extremes, we compared its spatial anomaly patterns with those from several contemporary NEE products (Fig. 8). It can be observed that, in general, the regional flux upscaling estimates from EC-MOD2012 captured the climate extremes–induced regional flux anomalies well, yet they indicated a much larger spatial extent of drought impacts. Among the global flux upscaling estimates, NIES2020 and FLUXCOM-X exhibited anomaly patterns broadly similar to those of MemoryFlux, whereas the two FLUXCOM2020 datasets generally showed weaker anomaly magnitudes. MemoryFlux also showed broad spatial consistency with the OCO-2 v10 MIP in indicating the climate extreme events in 2017 and 2019. Furthermore, MemoryFlux generally agreed with the CT2022 inversion for several of the selected events but with spatial differences, possibly due to the limited capacity of in situ based atmospheric CO₂ inversions to capture flux anomalies at finer spatial scales. To improve comparability, we recalculated the anomalies using a common 2011–2017 reference period for products with overlapping temporal coverage, including MemoryFlux, FLUXCOM-X, CT2022, NIES2020, and the FLUXCOM2020 ensemble products (Fig. S9). The recalculated results showed that the spatial extent and overall patterns of the regional NEE anomalies remained broadly similar to those obtained using the original product-specific reference periods, although the anomaly magnitudes changed slightly. Overall, among the NEE datasets considered, MemoryFlux provided a coherent representation of regional carbon flux anomaly patterns across the selected extreme events, including relatively moderate events such as the 2013 California drought (Fig. 8).

Figure 8 Spatial anomalies of JJA NEE associated with selected major drought and flood events in North America during the study period, as represented by MemoryFlux, EC-MOD2012, OCO-2 v10 MIP, CT2022, NIES2020, FLUXCOM2020 RS+CRUJRA1.1, FLUXCOM2020 RS+ERA5, and FLUXCOM-X. Black rectangles indicate the regions most strongly affected by the selected extreme climate events.

4 Discussion

4.1 North American terrestrial carbon budget: spatial distribution, magnitude, and IAV

Our estimate reasonably identified the spatial features of NEE on an annual basis and during the peak growing season. Annually, the most spatially extensive carbon sinks were located in the southeastern U.S. and characterized by mild temperatures and ample precipitation, agreeing with previous reports (Guanter et al., 2014; Huntzinger et al., 2012). Additionally, large annual carbon sinks were observed along the North Pacific coast and in tropical regions of North America, where evergreen forests prevailed (Anthoni et al., 2002). During the peak growing season, the net carbon uptake of croplands was highest in the central and western regions, leading to a large regional carbon uptake in the Midwest Corn Belt (Schuh et al., 2013). Previous studies have demonstrated the peak global photosynthesis in this region (Guanter et al., 2014; Hilton et al., 2017; Mueller et al., 2016). However, the seasonal carbon uptake was transitory. As the peak growing season advanced, both the spatial extent and magnitude of the strong carbon uptake decreased, because most of the absorbed carbon ended up in crop products that were consumed by humans and livestock or converted into bioethanol. The detection of this seasonal carbon uptake further validated the ability of MemoryFlux to diagnose spatial patterns of NEE. The large carbon sources on an annual basis were observed in the tundra regions of northern and central Canada (Virkkala et al., 2021), as well as in the arid regions of the central U.S. Warming in these regions could lead to the release of soil carbon, thereby accelerating climate change (Natali et al., 2019; Schuur et al., 2021). Many areas in the western U.S. were almost carbon neutral owing to limited vegetation distribution and precipitation (Xiao et al., 2014). Furthermore, according to the carbon sink statistics over boreal and temperate North America, all products indicated that Forests/Wooded biomes constituted the largest ecosystem carbon sink. These biomes were characterized by relatively low disturbance and high productivity (Yadav et al., 2022). Several studies supported forests as the main carbon sink across North America (Kurz et al., 2013; Stinson et al., 2011).

The mean annual NEE estimated by MemoryFlux was −1.27 ± 0.12 Pg C yr⁻¹, which was closer in magnitude to top-down atmospheric inversions than to the bottom-up upscaling estimates considered in this study. Several previous studies provide context for our estimates. King et al. (2015) revealed a mean annual net land-atmosphere exchange for North America of −0.89 ± 0.40 Pg C yr⁻¹ based on 11 inverse models during 2000–2009. Foster et al. (2024) reported a mean ±1 standard deviation annual NEE for North America of −0.81 ± 0.48 Pg C yr⁻¹ across 8 inverse models from 2007 to 2010. The Second State of the Carbon Cycle Report (SOCCR2) reported estimates of −0.61 Pg C yr⁻¹ with an uncertainty of ±75 % using a bottom-up approach and of −0.70 Pg C yr⁻¹ with an uncertainty of ±12 % using a top-down method over the period 2004–2013 (Hayes et al., 2018). Our slightly higher NEE estimates compared to those from previous top-down inversions may partly reflect differences in system boundaries and flux components represented by EC-based upscaling and atmospheric inversions. For example, fluxes associated with harvested products, fires, inland waters, and other lateral or disturbance-related carbon transfers may not be fully represented by tower-based ecosystem NEE but can influence atmospheric inversion estimates (Ballantyne et al., 2021). Over the past two decades, MemoryFlux showed a non-significant tendency toward increasing carbon uptake in North America, similar in direction to CT2022, NIES2020, and FLUXCOM-X and broadly consistent with previous studies reporting a strengthening of terrestrial carbon uptake (Keenan and Williams, 2018; Walker et al., 2021).

4.2 Capturing climate extremes-induced carbon flux anomalies over North America

The IAVs in atmospheric CO₂ concentrations can primarily be ascribed to the response of terrestrial ecosystem carbon fluxes to climate anomalies (Battle et al., 2000). Therefore, in principle, it is reasonable to use top-down approaches (e.g., OCO-2 v10 MIP and CT2022) to represent carbon anomalies in these events. OCO-2 XCO₂-based inversions (OCO-2 v10 MIP) have been proven to effectively capture regional carbon flux anomalies associated with climate extremes and have shown significantly better performance than in situ CO₂-based inversions (Chen et al., 2024; He et al., 2023a, b). In this study, the regional NEE anomaly patterns represented by MemoryFlux showed broad consistency with those from atmospheric inversion products and were spatially coherent with concurrent anomalies in SM and SIF. This agreement suggested that MemoryFlux provided a reasonable representation of regional NEE anomaly patterns associated with the selected extreme events. However, this assessment was primarily qualitative. Because continuous regional-scale NEE observations that could serve as an independent reference were unavailable, rigorous quantitative evaluation of model performance during extreme events remains difficult. Moreover, SM and SIF were not direct observations of NEE and therefore provided supporting ecological and hydroclimatic information rather than direct validation of carbon-flux estimates. Accordingly, the present analysis should be interpreted as a diagnostic assessment of spatial consistency rather than as evidence of superior quantitative accuracy.

Here, we further discuss how the implementation of an LSTM architecture in MemoryFlux improved the model estimates of IAVs and climate extremes-induced carbon flux anomalies, i.e., what the roles of explicitly incorporating historical predictor information and employing critical predictive variables play in this issue. Our findings showed that excluding the historical climate and environmental information from the LSTM models resulted in changes to the predicted spatial and temporal variability of NEE, which showed lower consistency with top-down inversions and contemporary upscaling products. At the PFT level, MemoryFlux also achieved better predictive performance than nonMemoryFlux, with the overall R² increasing from 0.63 to 0.79, accompanied by corresponding reductions in RMSE and MAE (Fig. S6; Tables S6 and S7). Substantial differences between MemoryFlux and nonMemoryFlux were also observed in modelling the seasonality and IAVs of NEE at the pixel scale. In particular, nonMemoryFlux showed a weaker representation of NEE IAVs (Fig. S7). However, the annual total NEE estimated by nonMemoryFlux was relatively close to that estimated by MemoryFlux. This may be because the differences between MemoryFlux and nonMemoryFlux diminished as the carbon flux was aggregated from the pixel scale to the annual total at the regional scale of North America. This scale effect reduced the differences between the two models. Therefore, the improvement in the magnitude of total NEE in North America can be primarily attributed to the inclusion of additional predictive variables, such as SIF and SWC. In contrast, the effective characterization of NEE IAVs and climate extremes-induced carbon flux anomalies resulted from explicitly incorporating historical information of vegetation and the environment.

Overall, MemoryFlux showed a more coherent representation of subcontinental NEE anomaly patterns across the selected events than nonMemoryFlux. However, it should be noted that nonMemoryFlux also roughly captured the major spatial patterns of NEE anomalies associated with climate extremes in many cases, although differences in anomaly magnitude and spatial extent were evident relative to MemoryFlux. This could be explained by the fact that the impacts of climate extremes are typically most severe during the concurrent period, while their impact gradually weakens over time (Piao et al., 2019). Therefore, a model using only contemporaneous predictors may capture the dominant concurrent response while representing potential legacy effects less clearly. In addition, nonMemoryFlux included contemporaneous vegetation variables that may themselves retain information related to antecedent environmental conditions and therefore partially encode legacy signals. Accordingly, the comparison between MemoryFlux and nonMemoryFlux should not be interpreted as a strict contrast between models with and without memory, but rather as an assessment of the added value of explicitly incorporating historical predictor information through the LSTM look-back structure.

4.3 Limitations and Future Perspectives

Several limitations should be acknowledged. Firstly, the reliability of flux upscaling is tightly linked to the number of observations, meaning that regions and ecosystems with few observations may suffer from larger uncertainties (see Fig. S2g and h) (Zou et al., 2024). Although we included as many flux-tower sites as possible across North America, observational coverage remains highly uneven among regions and PFTs. For example, SAV, MF, WSA, and CSH were each represented by only two sites, with fewer than 300 monthly NEE observations (Table S2). For the SAV PFT, SWC was excluded from the predictor set. Although the SAV model still showed good predictive performance, the reduced predictor set may limit further improvements in model performance. Predictions for sparsely sampled PFTs and regions with limited flux-tower coverage should therefore be interpreted with greater caution. Furthermore, although we provide quantitative PFT-level model performance metrics (Tables S6 and S7), the current version of MemoryFlux does not include a spatially explicit uncertainty framework that fully characterizes uncertainties arising from model structure, predictor data, and spatial extrapolation. Developing such a framework will be an important direction for future improvements to the dataset.

Second, the site-level models were trained using tower-observed meteorological variables, whereas regional upscaling relied on ERA5-Land gridded predictors, potentially introducing a predictor-domain mismatch despite harmonization of variable definitions and units. Comparisons at flux-tower locations showed strong agreement for T_a, DSR, VPD, and P, but weaker agreement for WS and SWC (Fig. S10). These discrepancies may reflect differences in spatial representativeness between point-scale observations and 0.1° gridded reanalysis data, as well as differences in variable definitions, local topography, surface roughness, and soil-moisture measurement depth, and therefore represent an additional source of uncertainty in regional extrapolation.

Third, the accuracy of PFTs or land cover type classifications may substantially impact the flux upscaling results, as the modelled relationships depend on these classifications. Although the static PFT map used here helps reduce interannual classification noise, it cannot capture actual land-cover transitions, including cropland expansion or abandonment, forest disturbance and regrowth, and fire-induced vegetation changes. Such mismatches may introduce additional uncertainty into regional NEE estimates. How to balance the representativeness of observations across ecosystems is an important issue for upscaling flux from sites to regions. Developing cross-PFT methods could be a promising way, for example, by replacing them with parameter-environmental relationships, as proposed by Famiglietti et al. (2023).

Lastly, while this study suggests that explicitly incorporating historical information improved the temporal and spatial representation of upscaling NEE, differences between MemoryFlux and top-down atmospheric CO₂ inversions remained. In particular, the large discrepancies occurred in the Midwest U.S., which may be linked to the pronounced impact of human management on croplands (Marcolla et al., 2017). Incorporating disturbance history and management information into ML modelling of carbon fluxes is expected to further improve the representation of these processes and enhance comparability among NEE estimates derived from different approaches. Furthermore, although these inversions provide valuable independent references for evaluating seasonal and interannual variations and large-scale carbon budgets, they are not direct measurements of NEE and are subject to uncertainties arising from factors such as atmospheric transport characterization and prior assumptions. Given the limited availability of large-scale direct flux observations, comparisons with atmospheric inversions provide a practical complementary means of evaluating regional NEE estimates. Future improvements in EC networks and inversion methods are expected to further enhance the accuracy and reliability of regional carbon flux modelling. Therefore, our results are preliminary, and we anticipate that refining our methodology will yield improved estimates in the future, particularly in terms of seasonality during the non-growing season.

5 Code and data availability

The datasets used in this paper are available from open resources. Eddy covariance data for the 35 FLUXNET sites utilized in this study are available from the FLUXNET2015 dataset (https://fluxnet.org/data/fluxnet2015-dataset/, last access: 17 October 2023). Eddy covariance data for the remaining 49 sites are available from the AmeriFlux website (https://ameriflux.lbl.gov/data/download-data/, last access: 17 October 2023). GLASS LAI and FAPAR products are available from https://glass.hku.hk/download.html (last access: 17 October 2023). NDVI dataset is available from https://doi.org/10.6084/m9.figshare.22267048.v1 (Xiong, 2023). GOSIF is available from http://data.globalecology.unh.edu/data/GOSIF_v2/ (last access: 21 October 2023). ERA5-Land data are available from https://doi.org/10.24381/cds.68d2bb30 (Muñoz Sabater, 2019). DEM is available from: https://www.eorc.jaxa.jp/ALOS/en/dataset/aw3d30/aw3d30_e.htm (last access: 19 October 2023). EC-MOD2012 is available at the Global Ecology Data Repository (http://globalecology.unh.edu/data.html, last access: 29 February 2024). The data for the OCO-2 MIP project are available at https://gml.noaa.gov/ccgg/OCO2_v10mip/ (last access: 18 November 2024). The CarbonTracker2022 carbon flux data are publicly available at https://gml.noaa.gov/aftp/products/carbontracker/co2/ (last access: 30 November 2023). The upscaling FLUXCOM-X products are open access at https://doi.org/10.18160/5NZG-JMJE (Nelson et al., 2023). The MemoryFlux dataset, which explicitly incorporated historical predictor information, along with the main related code, is available at https://doi.org/10.5281/zenodo.20482274 (Huang and He, 2026).

6 Conclusion

This study presents a new regional-scale upscaled estimate of NEE from a bottom-up perspective. The LSTM-based upscaling, combined with selected remote sensing variables and environmental factors, produced results that were consistent with established knowledge and with most NEE products, whereas other datasets exhibited consistency only in specific aspects. MemoryFlux reproduced major growing-season and annual spatial patterns of NEE at a fine resolution of 0.1°×0.1° and showed strong seasonal consistency with the top-down atmospheric CO₂ inversions from OCO-2 v10 MIP (r=0.96, p<0.001) and CT2022 (r=0.97, p<0.001). The mean annual total and standard deviation of NEE from 2001 to 2021 were −1.27 ± 0.12 Pg C yr⁻¹, which was closer in magnitude to top-down inversions (OCO-2 v10 MIP and CT2022) than to the regional (EC-MOD2012) and global (NIES2020, FLUXCOM2020 RS+CRUJRA1.1, FLUXCOM2020 RS+ERA5, and FLUXCOM-X) flux upscaling estimates. Annual NEE ranged from −1.50 to −1.05 Pg C yr⁻¹ over the study period. MemoryFlux also showed positive interannual correlations with OCO-2 v10 MIP (r=0.70, p=0.12) during 2015–2020, CT2022 during 2001–2020 (r=0.41, p=0.07), EC-MOD2012 during 2001–2012 (r=0.68, p<0.05), and FLUXCOM-X during 2001–2021 (r=0.69, p<0.05). The annual total NEE exhibited a non-significant tendency toward increasing carbon uptake in North America ($k=-\mathrm{0.009}$ Pg C yr⁻², p=0.09). In addition, MemoryFlux represented pronounced NEE anomaly patterns associated with the selected droughts during 2011–2013, 2017, and 2020–2021, as well as the 2019 Midwest floods. When historical predictor information was not explicitly incorporated during model training, the resulting NEE estimates showed markedly different spatial and temporal variability, showing less agreement with top-down inversions and other upscaling products, particularly at smaller scales.

Overall, the results indicate the added value of explicitly incorporating historical predictor information for representing NEE IAVs, particularly under climate extremes. MemoryFlux improves the representation of the spatial and temporal characteristics of NEE across North America and therefore provides a useful regional dataset for investigating terrestrial carbon dynamics and ecosystem–climate interactions.