Partitioning net carbon dioxide fluxes into photosynthesis and respiration using neural networks

Abstract The eddy covariance (EC) technique is used to measure the net ecosystem exchange (NEE) of CO2 between ecosystems and the atmosphere, offering a unique opportunity to study ecosystem responses to climate change. NEE is the difference between the total CO2 release due to all respiration processes (RECO), and the gross carbon uptake by photosynthesis (GPP). These two gross CO2 fluxes are derived from EC measurements by applying partitioning methods that rely on physiologically based functional relationships with a limited number of environmental drivers. However, the partitioning methods applied in the global FLUXNET network of EC observations do not account for the multiple co‐acting factors that modulate GPP and RECO flux dynamics. To overcome this limitation, we developed a hybrid data‐driven approach based on combined neural networks (NNC‐part). NNC‐part incorporates process knowledge by introducing a photosynthetic response based on the light‐use efficiency (LUE) concept, and uses a comprehensive dataset of soil and micrometeorological variables as fluxes drivers. We applied the method to 36 sites from the FLUXNET2015 dataset and found a high consistency in the results with those derived from other standard partitioning methods for both GPP (R 2 > .94) and RECO (R 2 > .8). High consistency was also found for (a) the diurnal and seasonal patterns of fluxes and (b) the ecosystem functional responses. NNC‐part performed more realistic than the traditional methods for predicting additional patterns of gross CO2 fluxes, such as: (a) the GPP response to VPD, (b) direct effects of air temperature on GPP dynamics, (c) hysteresis in the diel cycle of gross CO2 fluxes, (d) the sensitivity of LUE to the diffuse to direct radiation ratio, and (e) the post rain respiration pulse after a long dry period. In conclusion, NNC‐part is a valid data‐driven approach to provide GPP and RECO estimates and complementary to the existing partitioning methods.


| INTRODUC TI ON
The eddy covariance (EC) technique offers a unique opportunity for monitoring carbon and energy exchanges between land ecosystems and the atmosphere (Baldocchi, 2003) allowing near-continuous measurements integrated at the ecosystem scale. The number of study sites equipped with EC systems has increased over the years (Baldocchi, 2014;Chu, Baldocchi, John, Wolf, & Reichstein, 2017;Pastorello et al., 2017), and now we estimate that more than 1,400 globally distributed sites  are monitoring the most representative land ecosystems in different climate conditions. Most of the EC study sites are organized in regional networks such ICOS, AmeriFlux, and AsiaFlux, and contribute to the global FLUXNET network (Baldocchi, 2014).
The EC method allows for the measurement of the net ecosystem exchange (NEE) which is the difference between two larger flux components: gross primary production (GPP) and ecosystem respiration (RECO). GPP is the gross amount of carbon uptake by photosynthesis from vegetation while RECO is the total carbon efflux by the respiration processes of all organisms. Estimating GPP and RECO is a key step to better understand the underlying mechanisms constraining ecosystem function. Moreover, GPP and RECO estimates from EC are useful for modeling, supporting process-based model parameterization and validation, data assimilation, plant trait retrieval by model inversion (Dutta, Schimel, Sun, van der Tol, & Frankenberg, 2019;Pacheco-Labrador et al., 2019), upscaling Tramontana et al., 2016), as well as photosynthesis estimates based on remote sensing (e.g., Arnone et al., 2008;Verrelst et al., 2016;Zhang et al., 2016).
The EC technique does not directly measure GPP and RECO, and for this reason, numerical methods (termed partitioning methods) have been proposed for estimating GPP and RECO from NEE measurements (e.g., Desai et al., 2008;Keenan et al., 2019;Lasslop et al., 2010;Reichstein et al., 2005;Sulman, Tyler Roman, Scanlon, Wang, & Novick, 2016). The most widely used approaches are based on the use of NEE measurements for fitting simple physiologically based nonlinear relationships to estimate GPP and RECO using few meteorological drivers. These functional relationships are in general either light response functions linking global incoming radiation and GPP (Gilmanov et al., 2003), or respiration functions based on temperature (Reichstein et al., 2005), or also combinations of the two approaches (Keenan et al., 2019;Lasslop et al., 2010). The simple implementation and the robustness of the results (Lasslop et al., 2010) have led to their adoption as standard processing tools in the FLUXNET community (Pastorello et al., 2020, accepted).
However, these partitioning methods rely on important assumptions about the flux dynamics and their relationship with environmental drivers. Importantly, the assumed functional relations used can be similar to those applied in ecosystem models that make use of these data for their validation, creating a risk of circularity. In addition, although the functional relationships used in standard partitioning approaches are known to be valid at the organ level (where they can be measured), their direct application at the ecosystem spatial scale should be evaluated carefully (Medlyn, 1998;Medlyn et al., 2017;Musavi et al., 2016). Furthermore, to guarantee their wide applicability, partitioning methods use a limited number of drivers (i.e., air temperature, vapor pressure deficit, and global radiation).
However, fluxes dynamics are also potentially affected by many important environmental factors, such as soil moisture, soil temperature, or the ratio of diffuse to direct radiation, among others, that are generally not considered in the partitioning methods (Lasslop et al., 2012;Wohlfahrt & Galvagno, 2017). These limitations are partially compensated by the use of short temporal moving windows for parameters estimation, which takes more slowly changing factors indirectly into consideration (such as phenology, water, or substrate availability), but does not handle fast ecosystem responses well, like for example respiration pulses following rain events (Williams, Hanan, Scholes, & Kutsch, 2009).
Some studies have attempted to solve the limitation of FLUXNET standard partitioning methods by developing more comprehensive approaches. For instance, Scanlon and Kustas (2010) coupled CO 2 and H 2 O fluxes dynamics in the flux variance similarity approach for simultaneously partitioning carbon and water fluxes. The method, although interesting, requires canopy scale estimates of water-use efficiency (WUE), which introduces assumptions and uncertainty.
More recently, the estimation of gross CO 2 fluxes from NEE has been approached using the EC method in combination with parallel measurements of a trace gas such as carbonyl sulfide (COS; Commane et al., 2015) or 13 C isotopes (Ogée, Peylin, et al., 2003;Oikawa et al., 2017;Wehr et al., 2016;Wehr & Saleska, 2015) with the aim to disentangle the photosynthesis signals from respiration in daytime NEE measurements. Both methods are promising and are starting to be applied in the field; however, they currently require extensive and expensive instrumentation and the uncertainty in the results is still large (Oikawa et al., 2017;Whelan et al., 2018), limiting, for now, their application to a restricted number of study sites.
These methods, however, are all subject to assumptions that could affect the resulting estimates of GPP and RECO. A possible alternative approach, based on existing measurements and not subject to the limitation of the standard methods, could be provided by machine learning methods.
Machine learning methods are generally used to model underlying complex relationships linking a set of predictors with one or more outputs. The core characteristics of machine learning methods are that predictors are not prescribed and the relationships between input (drivers) and output (fluxes) are inferred from the data. Several studies have reported the capability of machine learning to reproduce complex patterns in ecological studies, and in particular in relation to EC measurements (e.g., see Moffat, Beckstein, Churkina, Mund, & Heimann, 2010;Moffat et al., 2007;Papale & Valentini, 2003;Reichstein et al., 2019;Tramontana et al., 2016).
In this study, we develop and test an empirical machine learning-based approach for retrieving GPP and RECO from NEE. The main motivation was to evaluate if a purely empirical approach could provide estimates of the two components without any predefined relationships and with the flexibility to use multiple sets of meteorological and biological forcings.
In this experiment, we used artificial neural networks (ANNs), which have already been tested as a partitioning tool in a crosssite intercomparison among several partitioning methods (Desai et al., 2008;Oikawa et al., 2017), without appreciable differences with respect to the other methods. In those studies, the ANNs were trained using only nighttime measurements and then applied to extrapolate RECO during daytime. Here, we revisited the ANN implementation by proposing a new scheme in which the ANN is constrained using the measured net CO 2 fluxes to directly estimate the two components GPP and RECO.
Our aim was to design a machine learning-based method flexible enough to: (a) use a large set of drivers, accounting for the complexity of the GPP and RECO responses, (b) ensure general applicability of the method to allow ANN-based partitioning in the context of FLUXNET thanks to the use of generally measured variables. In order to evaluate the results, we compared the GPP and RECO produced by the ANN with the estimates obtained by two partitioning methods used in FLUXNET, by analyzing: (a) the consistency of the estimates at yearly, seasonal, daily, and (half) hourly time steps; (b) the dynamics of the seasonal and diurnal cycles of the estimated gross CO 2 fluxes; (c) the functional relationships between micrometeorological inputs and fluxes in output as reproduced by the ANN; and (d) the effectiveness of the proposed method to predict additional patterns of gross CO 2 fluxes not (or less) accounted for by the standard partitioning methods.

| ANN algorithm
ANNs are nonlinear and nonparametric methods for regression and classification that artificially emulate the functioning of a biological brain (Haykin, 1999). The base unit of an ANN is the "neuron" where the numeric information in the input is weighted, condensed, and transformed (by a linear/nonlinear transfer/activation function) to be transferred to other neurons. Neurons are organized in layers that are interconnected: the outputs of m neurons in one layer are the inputs for n neurons of the next layer and the signals are transferred through connections associated with multiplicative weights.
The learning (or training) process of an ANN consists of adjusting the weights of the network on the basis of specific examples provided as input (supervised training).
In this experiment, we developed a new ANN architecture designed to emulate 1 the ecosystems processes driving NEE. This customized neural network (hereafter NN C-part ) is based on the concept that NEE, measured by the EC system, is the difference between RECO and GPP (Equation 1): and that these two fluxes have individual drivers and dependencies.
The uniqueness of our approach consists of imposing physical constraints in the proposed neural network, on the basis of the known properties of RECO and GPP. The overall structure is based on two subnetworks ( Figure 1): one for retrieving RECO (sub neural network, SNN RECO ) and the other for retrieving GPP (SNN GPP ). The output of these two subnetworks (S RECO and S GPP for SNN RECO and SNN GPP, respectively) is combined in the last neuron of the overall structure, where NEE is calculated and compared with the measurements in order to optimize the network's weights (network training). The output signals of both subnetworks are constrained to be always positive, except for SNN GPP , which is constrained to be 0 during nighttime, as photosynthesis requires light. However, RECO and GPP have opposite signs; therefore, the connection weight of RECO is fixed to be positive (w RECO = 1) and that of GPP to be negative (w GPP = −1), mirroring the sign convention of the measured NEE. Since the last transfer function is unbounded and linear, and the bias is fixed to 0, the last node effectively reproduces the equation of NEE as: where and The two subnetworks have different input drivers but a similar structure including two hidden layers: the first hidden layer has n hidden neurons (different for the two fluxes) and a hyperbolic tangent sigmoid activation function ("Tansig"), which improved the performance of network optimization; the second has only one neuron and a log-sigmoidal activation function ("Logsig") that constrains the output to positive values, allowing the implementation of Equation (2) in the last node. The SNN GPP , used to simulate GPP, has also a third layer where the measured incoming shortwave radiation enters as input in a node where a product is applied. This last node has a positive linear transfer function ("Poslin"); thus, the output from the SNN GPP is 0 at nighttime and positive during the day. With this structure, the last node of the SNN GPP mimics the light-use efficiency (LUE) approach: the output coming from the two previous hidden layers is an LUE that is then multiplied by the available shortwave radiation to estimate GPP. In this way, the drivers are used to estimate an instantaneous proxy of the LUE which takes into consideration seasonal and diurnal variability of carbon uptake, including photosynthesis saturation at high light, while shortwave radiation mainly defines the magnitude of 1 The use of this term is not incidental: surrogate modeling, also known as emulation, is about developing statistical models that learn to mimic costly physical models, such as radiative transfer or climate models, using a representative dataset of simulation/ (1) NEE = RECO − GPP, the flux. This decreases the weight given to the incoming shortwave radiation in relation to other driving variables. Up to this step, the output of the two subnetworks (S RECO and S GPP ) have positive signs; the sign convention of the NEE equation is finally restored in the last neuron of the overall structure by associating a positively signed weight to S RECO and a negatively signed weight to S GPP .

| Dataset used
For the purpose of this experiment, we used data from the FLUXNET2015 includes meteorological and EC measurements that were quality checked and processed with standard tools (Papale et al., 2006;Pastorello et al., 2020) and provided with per-variable quality flags. For more details, see the webpage http://fluxn et.fluxd ata.org/data/fluxn et201 5-datas et/.
Data used for training and validation of the neural network were taken from the "FULLSET" "TIER 1" collection. Among the sites available in this collection, we selected a subset of 36 study sites (listed in Table 1) on the basis of the data quality and data availability in order to ensure the best conditions for the partitioning methods comparison. Site-years were selected if the percentage of meteorological gap-filled data was less than 20% and the measured NEE covered at least 10% of both daytime and nighttime periods.
Among the variables stored in the FLUXNET2015 dataset, we used the GPP and RECO obtained with the two standard methods: the nighttime method (NT) from Reichstein et al., 2005, and the daytime method (DT) from Lasslop et al., 2010. In the NT method, the Arrhenius-type temperature-response curve of respiration as modeled by the Lloyd and Taylor equation (Lloyd & Taylor, 1994), driven by air temperature, is used for estimating RECO. This method makes use of the nighttime (when global incoming radiation <20 W/ m 2 ) NEE observations (assumed as representative of RECO given the absence of photosynthesis at night) to fit the Lloyd and Taylor equation. There are two parameters to fit in the Lloyd and Taylor equation: the activation energy (E0) and the respiration rate at the reference temperature (R ref ). In the NT method, E0 is estimated at an annual scale by calculating E0 values every 15 days and then averaging the three with smaller uncertainty. Once E0 is fixed, the R ref parameter is estimated using short-term moving windows (4 days).
On the basis of Equation (1), GPP is finally calculated by subtracting NEE from the daytime RECO extrapolated by the fitted model.
The DT method combines a rectangular hyperbola light response curve (Falge et al., 2001) for estimating GPP and the Lloyd and Taylor equation for estimating RECO (as in the NT method). In the F I G U R E 1 The scheme of the customized neural network applied for retrieving gross primary production (GPP) and ecosystem respiration (RECO) from net ecosystem exchange (NEE) measurements. The two subnetworks (SNN GPP left, in green, SNN RECO right, in brown) are connected in the last step to estimate NEE. Inputs for RECO are air and soil temperature; soil water content; wind speed and wind direction; the day of the year (sine and cosine values of angular day of year); the average value of nighttime NEE. Input GPP are, in addition to the shortwave incoming radiation used in the product, air temperature; vapor pressure deficit; soil water content; potential and actual shortwave incoming radiation; wind speed and wind direction; a proxy of the mean seasonal GPP dynamic derived from the nighttime and daytime averaged NEE. See Section 2.3 of the main text for details DT method, nighttime NEE is used only to estimate the E0 parameter of the Lloyd and Taylor equation, while the remaining unknown parameters, including R ref , are fitted using daytime measurements.
In DT, the light response curve is driven by the incoming shortwave radiation (SW_IN), but it is adjusted for the effect of stomata closure due to the atmospheric evaporative demand using VPD as an additional driver of GPP (Lasslop et al., 2010). Similar to NT: (a) the Lloyd and Taylor model is driven by air temperature, (b) the model TA B L E 1 List of the subset of FLUXNET2015 study sites used in this experiment for net ecosystem exchange partitioning. The study sites used also for NN C-part validation (see section S1) are marked with the symbol (*) in the "Validation" column.  Archibald et al. (2009) parameters are fitted on short-term moving time windows to account for slow changing factors.

| Input variables and data preparation
The variable used as target values in the NN C-part training was the half hourly NEE (µmol CO 2 m −2 s −1 ) measured with the EC technique. In particular, we used the NEE_CUT_USTAR50 variable (Pastorello et al., 2020). We used a comprehensive subset of micrometeorological variables measured at the flux towers as candidate drivers for NN C-part , and additional variables derived from NEE, with the aim to explain variability of carbon fluxes at seasonal, daily, and half hourly resolutions. In particular, in terms of meteorological drivers, we used the used in order to provide information about seasonality and length of the day. The SW_IN_POT was also aggregated at daily resolution to better track the seasonality of light conditions; the first derivatives of half hourly and daily SW_IN_POT were also calculated to add specific information about the seasonal and diurnal dynamics of light (Bodesheim, Jung, Gans, Mahecha, & Reichstein, 2018;Papale, Black, et al., 2015;Papale, Migliavacca, et al., 2015).
Finally, since plant photosynthesis and RECO change seasonally also due to substrate availability and management (in case of managed sites), daily average of nighttime NEE (NEE NIGHT ) and a proxy of GPP (GPP prox ) were calculated and used in input. In this case, gapfilled (Reichstein et al., 2005) half hourly NEE was used in order to have a continuous time series. Proxy of GPP was calculated by using NEE NIGHT , and the average of daytime NEE (NEE DAY ) as follow: where k is the fraction of daytime hours for each day.
Some of the drivers are periodic, which means that extreme values have similar meanings (e.g., the DOY 0 and 365 or the WD 0 and 360). To take this into consideration, we used their sine and cosine transformations instead of the original variable's value.
The two subnetworks that estimate RECO and GPP use a different set of drivers, selected on the basis of the expected role and tests on performances and that are listed in GPP prox , and sine and cosine values of WD. Note that SW_IN is then used both as driver and as single input in the GPP subnetwork for the product with the LUE (see Figure 1).

| ANN training
Because our goal is to partition NEE into its two fluxes components, the NN C-part was trained at the site level and year-by-year; we used only high quality measured half hourly NEE (see Section 2.2) as training target. The training of the NN C-part was carried out using both nighttime and daytime measurements of NEE.
As commonly used in ANN training and application, inputs and outputs were normalized in the range +1/−1 by Equation (6): where X min and X max are estimated as the maximum of the absolute value of X then negative signed in the case of X min and positive signed in the case of X max ; this preserved the zero value in both normalized

TA B L E 2
List of the variables used as drivers to estimate gross primary production (GPP) and ecosystem respiration (RECO) pacity of the ANN to simulate NEE was also analyzed (see Data S1, section S1).

| Statistical analysis and evaluation
Results of this experiment were evaluated by comparing GPP and RECO from NN C-part with the NT-and DT-based partitioning methods in the FLUXNET2015 collection. As complementary information, we repeated the same analysis for the gross CO 2 fluxes derived from the synthetic NEE dataset provided by the MUSICA model importance of the different drivers (see Data S1).

| Cross consistency between NN C-part and standard partitioning methods (DT and NT)
The half hourly GPP and RECO retrieved by NN C-part were consistent with the output obtained by the NT and DT methods, with a slightly higher correlation (R 2 and lower RMSE) with the DT method for GPP estimation and with the NT method for RECO, in both cases higher than the correlation between DT and NT methods ( Table 3).
The mean squared error of partitioned gross CO 2 fluxes between NN C-part and the NT/DT methods was <1.53 µmol CO 2 m −2 s −1 on average and comparable to the average value of the estimated random uncertainty of NEE ("NEE_CUT_USTAR50_RANDUNC," variable provided by FLUXNET 2015). Looking at the bias, the magnitude of both GPP and RECO retrieved by NN C-part was slightly lower than the DT and NT methods, with the predicted values by NN C-part closer to the DT method (though the spread among methods was very small).
The average values of bias among methods were low on average (<27 g C m −2 year −1 ) and comparatively lower with respect to the NEE random uncertainty reported in the FLUXNET2015 database (<10% of the reported "NEE_CUT_USTAR50_RANDUNC"). The agreement between NN C-part and the two FLUXNET standard methods was significantly lower in the case of RECO suggesting that the differences occurred in the diurnal and seasonal cycles of RECO between NN C-part and both DT and NT methods were comparatively larger than GPP.
High correlation and consistency among NN C-part and the DT and NT methods in both GPP and RECO were found for daily, seasonal, and yearly values and also for the daily average anomalies with respect to the seasonal values ( Figure 2). Daily, growing season, and yearly average values from NN C-part were highly correlated with the two standard methods (R 2 > .83), in particular for GPP (R 2 > .97).
In terms of bias, results by NN C-part were more consistent with the DT results confirming the findings from the subdaily comparison.
High consistency was also found in the case of GPP daily anomalies (R 2 > .88 and RMSE < 0.52 µmol CO 2 m −2 s −1 ), but it decreased in the case of RECO (Figure 2).
This general consistency between NN C-part and standard meth-  Figures S5 and S6).

TA B L E 3
Cross consistency among the retrieved gross primary production (GPP) and ecosystem respiration (RECO) by NN C-part , daytime method (DT), and nighttime method (NT) methods at half hourly time step, also after removing the mean daily value of the fluxes. Statistics reported in the

F I G U R E 2
Scatter density plot showing the cross-consistency between partitioned gross primary production and ecosystem respiration by NN C-part (x-axis) and the ones by the daytime method and nighttime method methods (y-axis) aggregated at daily time step (a, e, i, m), yearly (b, f, j, n),), and looking for the seasonal cycle (c, g, k, o) and daily anomalies (d, h, l, p) 3.2 | Consistency among the seasonal and diurnal dynamics of partitioned gross CO 2 fluxes

| Seasonal cycle
The seasonal pattern of the retrieved gross CO 2 fluxes by NN C-part closely matched that of the DT and NT methods, with the only exception being that predicted RECO by NN C-part was slightly lower in late spring/early summer (Figure 3, left). The predicted seasonal patterns were also similar among the represented International Geosphere-Biosphere Programme vegetation classes (Figure 3, right). Results showed a large degree of agreement among the three methods, in particular between NN C-part and NT.
The seasonal dynamic of gross CO 2 fluxes derived from NEE MUSICA by NN C-part , NT, and DT agreed consistently. However, all methods underestimated both GPP MUSICA and RECO MUSICA (in particular the DT method) during the growing season (between April and July, see section S4.2 and Figure S7 in Data S1).

| Mean diurnal cycle of GPP
We found a good agreement between the pattern of the mean diurnal cycle of GPP predicted by NN C-part and FLUXNET standard methods ( Figure 4). The GPP diurnal dynamic followed an asymmetrically bell-shaped curve with the maximum value reached before midday (Figure 4, top panels). Analyzing the differences more in depth, the mean diurnal cycle of GPP predicted by NN C-part systematically diverged from the one predicted by the DT standard method, following a pattern that was consistent across the seasons (Figure 4, bottom panels). In particular, the GPP estimated by NN C-part was lower during the first hours of the morning (just after the dawn) and in the late afternoon; instead, the GPP estimated by NN C-part was slightly higher than DT from the morning until mid-  Figure S8).

| Mean diurnal cycle of RECO
In comparison to GPP, the diurnal cycle of RECO modeled by NN C-part showed larger differences with respect to the DT and NT methods, in particular regarding the magnitude of the flux. All the methods showed an increased respiration when moving toward midday, in particular during the growing season ( Figure 5, top plots); nevertheless, the NN C-part estimated values were lower than predictions from the other methods. Conversely, the differences between the two

F I G U R E 3
The mean seasonal cycles of gross primary production (GPP; a) and ecosystem respiration (RECO; c) obtained by NN C-part , daytime method, and nighttime method. The consistency among methods was evaluated by the determination coefficients (R 2 ) among the seasonal patterns of GPP (b) and RECO (d) per International Geosphere-Biosphere Programme (IGBP) vegetation class (the number in brackets after each IGBP vegetation class stands for the number of sites in each category). Only FLUXNET2015 study sites at northern hemisphere (latitude > +15°N) were used for that analysis. The vegetation classes are the same as Table 1 standard methods were negligible with a slightly higher RECO by NT in comparison to DT ( Figure 5, lower plots, calculated as NT-DT).
This finding was not fully mirrored in the patterns of RECO MUSICA (see section S4.3 and Figure S9). In that case, NN C-part predicted higher RECO in comparison to both DT and NT, but more consistent with RECO MUSICA . Conversely, NN C-part predicted lower RECO than the NT method (and more consistent with DT) during the morning; the differences between DT and NT fitted on NEE MUSICA were also enhanced.

F I G U R E 4
The trends of the mean diurnal cycle of gross primary production (GPP; a-f) retrieved by NN C-part (black line), daytime method (DT; red line), and nighttime method (NT; green line). The differences between diurnal cycle of retrieved GPP are also shown (g-l): NN C-part -DT (red line), NN C-part -NT (green line), and NT-DT (dashed blue line). Only FLUXNET2015 study sites at northern hemisphere (latitude > +15°N), having at least 2 years of data, were used for the analysis F I G U R E 5 The trends of the mean diurnal cycle of ecosystem respiration (RECO; a-f) retrieved by NN C-part (black line), daytime method (DT; gray line), and nighttime method (NT; dashed gray line). The differences between diurnal cycle of retrieved RECO (g-l) are also shown: NN C-part -DT (continuous black line), NN C-part -NT (dashed black line), and NT-DT (gray line). Only FLUXNET2015 study sites at northern hemisphere (latitude > +15°N), having at least 2 years of data, were used for that analysis These differences increased during the nighttime hours, while the predicted RECO by the two standard methods became closer during the daytime hours (particularly when air temperature reached its maximum value). It is interesting to note that the peak of the maximum RECO by NN C-part was shifted compared to the DT and NT methods and closer to the one of RECO MUSICA , suggesting a better capacity to reproduce diurnal patterns when multiple drivers were involved ( Figure S9).

| Functional relationships between partitioned fluxes and meteorological drivers
The functional response obtained by artificially varying the drivers in NN C-part was consistent with the current knowledge. The response of GPP to SW_IN was as expected, with an increased photosynthesis due to the increased light that resulted as saturated at high radiation values (Figure 6a-d). The response to VPD clearly showed the effect of VPD limitation on photosynthesis (Körner, 1995) that was evident also at relatively low values of VPD (Figure 6e-h).
The response of photosynthesis to air temperature showed a strong and positive response at the start of the growing season that became flat in the summer period (Figure 6i-l). The functional relationships between RECO and TA followed the expected increase of RECO with TA (Figure 6m-p) except for the dry ecosystems (red points in Figure 6) where TA was not the main driving factor. These general patterns of relationships were consistent across the seasons and in the different vegetation types, despite the expected variability due to the ecosystem-specific properties.
When applied using the original measurements, the functional ecosystem response to the micrometeorological forces was largely consistent with the ones obtained using the DT and NT standard methods, although affected by multiple co-acting and confounding factors ( Figure 7). The shape of the relations was slightly different in the four periods analyzed, but it was important to recall that these were average responses from sites in different climatic conditions. It could be noted that an almost perfect match was present in the case of GPP at this broad scale; in the case of RECO, the NN C-part showed in general lower respiration fluxes at high temperatures (so mainly daytime) as opposed to the other methods, although the pattern of response was basically the same. This could be an indication of a water resources limitation during the summer period and at higher temperatures, that could be detected by the NN C-part since it also used SWC as a driver. In fact, the functional relationships between RECO and SWC retrieved by NN C-part highlighted a clear direct relationship between SWC and RECO in drought conditions, in particular for non-forested ecosystems ( Figure S16).

| NN C-part instantaneous LUE and the ratio between diffuse to direct radiation
The instantaneous LUE estimated by NN C-part showed an increasing trend with respect to the increased diffuse to direct radiation ratio F I G U R E 8 Sensitivity of the light use efficiency (LUE) of NN C-part to the diffuse to direct radiation ratio. The trend of LUE by NN C-part with respect to the diffuse radiation ratio is reported in the panels (a-d; here LUE is normalized for comparison purpose). The difference between LUE (µmol CO 2 /W) as estimated by NN C-part and that from the daytime method (DT) standard partitioning method is also reported (e-h); the statistics were aggregated per diffuse to direct radiation class (here estimated by the proxy 1-SW_IN/SW_POT and reported as percentage).
We reported the following statistics of x (x = LUE(NN C-part ) − LUE(DT)): the median values (bar), the mean value (*), the interquartile range (by brackets). Only LUE at midday is used in this figure When compared against the LUE estimated by the DT method, NN C-part showed a moderate but consistent higher sensitivity to the fraction of diffuse to direct light with a higher LUE when the diffuse to direct radiation ratio was higher than 40% (Figure 8e-h).

| Hysteresis of the diel cycle of gross CO 2 fluxes
The diurnal cycle of gross CO 2 fluxes exhibited hysteresis, which we analyzed through the diurnal cycle of GPP with respect to SW_IN (Figure 9a-c) and the diurnal cycle of RECO with respect to TA (Figure 9e-g). We generally observed a clockwise cycle for GPP, with higher values of GPP in the morning in comparison to the afternoon. All the methods showed this pattern although it was more evident in the NN C-part and NT methods. We provided numerical estimates by calculating the integral area included in the hysteresis pattern (larger area indicates larger hysteresis) that was reported in Figure 9d.
For RECO, only NN C-part detected an appreciable hysteresis in the diurnal cycle ( Figure 9h); this is expected because both the DT and NT methods use the same TA-dependent and invariant F I G U R E 9 Hysteresis of the gross primary production (GPP) diel cycle predicted by NN C-part with respect to SW_IN in three sampled study sites (a-c). Three different International Geosphere-Biosphere Programme classes are reported: DBF (a), ENF (b), and GRA (c), while the period of interest was July-August. (d) Distribution of the GPP hysteresis integral area (then normalized) estimated by NN C-part , nighttime method (NT), and daytime method (DT) partitioning methods. The hysteresis of the diel cycle of ecosystem respiration (RECO; here logtransformed) with respect to the TA dynamic, in three sampled study sites, is reported (e-g). The distribution of the RECO hysteresis integral area (then normalized) estimated by NN C-part , NT, and DT partitioning methods is reported (h) F I G U R E 1 0 The sensitivity of half hourly ecosystem respiration (RECO) predicted by NN C-part to the SWC "pulse". The investigated study site was US-SRG. Fluxes are reported in panel (a) and more specifically we reported: RECO predicted by NN C-part (blue line), nighttime method (green) and daytime method (red) and the half hourly net ecosystem exchange (gray).
Meteorological variables are reported in panel (b); more specifically we reported the SWC dynamics (blue line) and the rain events (P [mm]; red line) relationships for the whole diel cycle. In general, the hysteresis of RECO is realized as a counterclockwise cycle, with the value of RECO, for the same TA, higher during the afternoon and nighttime hours and lower during the morning. However, the counterclockwise cycle observed in the RECO patterns was not always observed, and in some cases, the opposite hysteresis cycle was found.

| RECO pulses due to SWC variations
The evaluation of the capacity to correctly reproduce rapid responses of RECO (respiration pulse) to changes in SWC (e.g., a post rain event after a prolonged drought period in arid environments) is in general difficult due to the sporadic nature of the events and the noise associated with the measurements.
We evaluated the performances of the three models in a specific event that was clearly visible in the NEE time series (US-SRG site, Figure 10). The analysis of the predicted RECO shows that NN C-part was much more sensitive to the fast response of respiration pulses. The NT and DT methods, that did not have SWC as a direct input, showed a slow changing pattern with the peak of respiration significantly delayed compared to the measurements.

| Consistency of NN C-part functional response with theory on plant physiology and with the DT and NT standard methods
The retrieved GPP and RECO as obtained from NN C-part were consistent with the estimates by the DT and NT methods both when trained directly on EC measurements and using synthetic data with NEE MUSICA (Figure 2; Figures S3 and S4). This was valid also for the relationship between micrometeorological drivers and gross CO 2 fluxes calculated by NN C-part and those implemented in the DT and NT partitioning methods (Figures 6 and 7).
The GPP response to light used by NN C-part is curvilinear and consistent with the pattern of that used in DT method (despite some small but systematic divergences). The pattern of the GPP-VPD relationship captured by NN C-part is consistent with current knowledge about the protective (physiological) mechanisms of stomata closure carried out by plants in response to the increase in atmospheric evaporative demands (Körner, 1995). There is also a general consistency of the The largest mismatch we found between NN C-part and the two standard methods of partitioning (NT and DT) is in late spring/early summer, with lower fluxes of daytime RECO predicted by NN C-part compared to both standard methods. When evaluated on the synthetic dataset, the mismatch between NN C-part and the standard methods was significantly lower except for the anomalies, which were better reproduced by NN C-part . In addition, the estimates from NN Cpart better matched the diurnal dynamics of RECO MUSICA than the estimates from the standard methods ( Figure S9); conversely the method that showed the highest mismatch with respect to RECO MUSICA was the DT method. All the methods underestimated both the reference GPP MUSICA and RECO MUSICA (despite the differences were very tiny) without any relevant effect on NEE. This could be related to a certain limitation of these methods to infer the complex relationships of the MUSICA model with the reduced set of drivers used in this experiment, but we have also to consider that both NN C-part and standard methods of partitioning were trained on noisy NEE signals (see section S4 in Data S1 for details) while both GPP MUSICA and RECO MUSICA used for reference are noise-free.
Recent papers highlighted possible limitations of the standard methods for partitioning stemming from the fact that they do not consider the inhibition of daytime leaf respiration (the Kok effect; Keenan et al., 2019;Wehr et al., 2016) that produces an overestimation of daytime RECO and GPP in the NT method (Wehr et al., 2016) and an underestimation of nighttime respiration (see Keenan et al., 2019) in DT. Although some of the differences between the RECO estimated by NN C-part and the other methods (see Figures 3 and 5 and Figure S9) are consistent with experimental data focused on this topic (e.g., Keenan et al., 2019;Wehr et al., 2016), it is not possible to demonstrate that the NN C-part method, as implemented in this experiment, is able to reproduce the light inhibition of leaf respiration. Moreover, the consistency between the DT and NT partitioned fluxes used as reference in this experiment, and the pattern emerged by other modeling experiences (see e.g., Jung et al., 2020, where DT and NT partitioned fluxes were globally upscaled) suggests that the Kok effect could have a minor relevance on biases of RECO and GPP estimates compared to other sources of NEE measurement uncertainties.
In terms of ecological responses, the NN C-part outputs showed a clear and interesting response to the diffuse/direct light with an increased LUE in diffuse conditions and with hysteresis of the diel cycle of GPP and RECO that are only partially shown by the standard methods. The latter showed also strong limitation in reproducing the respiration pulse while NN C-part showed this capacity, although tested only in one case. In summary, these findings highlight the strength of this approach that is able to reproduce these patterns even without being trained specifically for this.

| The advantages of the proposed NN C-part approach
One relevant difference between NN C-part and the standard partitioning methods is the absence of prescribed relationships between drivers and fluxes. In fact, only a few constraints have been set in the NN C-part structure (RECO and GPP signals must be positive, weights in the last node have to be consistent with the sign convention of NEE) and the relationship between inputs and outputs is set on the basis of the data without any assumption (as common in all the machine learning methods). This property could be a great advantage when used to study a phenomenon that does not respond with the same invariant pattern. The synthetic datasets used in this experiment are generated using functional relationships which are more complex than the ones implemented in the FLUXNET standard methods. The better matching by NN C-part highlights the potential of this algorithm to capture additional patterns of carbon fluxes dynamics. As a direct consequence, this led to some systematic differences observed in the GPP and RECO patterns compared to standard methods. For instance, the shape of the GPP diurnal cycle predicted by the DT method is almost symmetrically bell shaped: GPP reached a maximum around midday and only VPD, used as a downregulation function, when higher than 10 hPa, can modify this pattern that is otherwise dictated by incoming radiation. This is also visible in the modeled NEE, by comparing NN C-part and DT outputs with measured data ( Figure S10).
Conversely, the pattern of GPP obtained by the NT method is less affected by prescribed relationships because GPP is calculated as the difference between RECO and the measured NEE. The diurnal dynamic of GPP estimated by NN C-part is more similar to the NT estimates, which suggests that the formulation of the DT approach could be improved, for example, by using the non-rectangular light response function (instead of the rectangular used in DT, see Gilmanov et al., 2003)  In the case of RECO, both DT and NT used the same modeling approach, based on the Lloyd and Taylor model driven by the air temperature, resulting in high consistency between their results.
However, it is known that soil temperature is also an important driver of RECO in particular to define its temporal (diurnal) dynamic (Lasslop et al., 2012;Wohlfahrt & Galvagno, 2017) because it is directly linked to soil respiration, which is a large contribution to the total RECO (Misson et al., 2007). The NN C-part method uses soil temperature (in addition to air temperature) as a direct driver of RECO and the analysis of the functional responses highlighted an important and direct effect ( Figure S3). By looking at the analysis of RECO driver importance (see section S8 and Figure S20 in Data S1), it seems that TS and TA had, in general, a similar importance, Another important driver of RECO that is not explicitly accounted for in the standard methods was SWC that is particularly important in dry sites (see section S8). The DT and NT methods use a dynamic parameter (R ref ) estimated using moving windows to indirectly consider the slow changes in water availability (Reichstein et al., 2005).
Soil water content, however, also has instantaneous effects in some ecosystems (e.g., the respiration pulses after short rain events, see Jarvis et al., 2007) that need the direct use of SWC as input, like in the NN C-part method, to be effectively modeled.
Other slow changing factors affecting the fluxes, such as phenological state, substrate availability, and management or other disturbance events, that are considered indirectly through R ref in the standard methods, are represented in the NN C-part approach through the use of averaged NEE-derived quantities as drivers.
Finally, the standard methods do not account for EC footprint variability, which can significantly affect the magnitude of the measured fluxes if the fetch is small and not fully homogeneous (see section S8). This could be considered, for example, by fitting the models per wind sectors and wind speed classes (or using footprint models), but the amount of available data could become critically low in less represented conditions. In the NN C-part , we used wind variables (direction, speed) as input in order to indirectly take into consideration footprint variability. Wind-related variables appeared to be important in a few study sites (see Figures S15, S16, and S17). By comparing the ANN trained by including/excluding wind variables, we detected an effect on performance (increase of RMSE), which highlights an important effect on the instantaneous estimated flux values. However, the role of wind variables is reduced in the yearly budget calculation of gross CO 2 fluxes (~6 g C m −2 year −1 , average value across study site) even though differences were detected site by site (site-specific differences ranged between −169 and +122 g C m −2 year −1 ). In summary, the implemented machine learning approach gives the possibility to use variables as drivers which would be difficult to consider in process-based approaches, which require that the role and effects of each single variable being known and correctly encoded (at the ecosystem scale).

| Uncertainty and limitation of the proposed method
Despite the encouraging results, there are also limitations and uncertainty sources to consider when a pure machine learning approach is used. The absence of prescribed relationships and the few constraints fixed made the method robust against incorrect or incomplete formulations used in the standard methods but also led to higher sensitivity to the uncertainty in the data used in the training phase of machine learning, particularly when both NEE and the meteorological variables used as drivers are affected by long gaps.
In the NN C-part , the weights of the ANN are optimized against NEE; thus, the uncertainty of NEE could have a significant impact on their estimation. This could be particularly relevant for nighttime NEE measurements that can be affected by higher uncertainty (e.g., advection, large footprint) that, even if filtered, lead to long gaps. In these situations, a model based on ecological responses is in general more robust, and needs less data for a proper calibration than a datadriven approach.
Another limitation is related to the availability of the large set of input data that this approach requires in order to fully exploit its advantages, which are often not fully available at all FLUXNET sites.
In this situation, the model can be used with a reduced set of drivers, but this could also lead to higher uncertainty and lower performance in specific sites, while models based on ecological processes like the NT and DT may be more robust because the response pattern is already prescribed.

| CON CLUS IONS
In this study, we propose a machine learning method for NEE partitioning (NN C-part ) using an ANN approach with a tailored structure to simultaneously retrieve GPP and RECO. From a methodological point of view, this approach is an example of simple hybrid modeling , combining a neural net with a "theoretical" ecosystem equation based on LUE.
The NN C-part is designed to use as input a comprehensive set of drivers selected among the most commonly measured variables at EC sites. This property of NN C-part (and of machine learning in general) can be exploited in future studies, for example, using as input other variables such as diffuse radiation, COS, 13 C isotopes, or SIF.
Additionally, NN C-part does not make use of prescribed relationships between inputs and outputs and only a few constraints are introduced. Results in terms of the consistency of the gross CO 2 estimated fluxes by NN C-part with the estimates from FLUXNET standard methods were good and the proposed model is able to reproduce magnitudes, patterns, and functional responses. However, the systematic differences emerged by the cross-consistency analysis and the ecological patterns captured by NN C-part suggest that the implemented ANN reproduces fluxes dynamics more realistic than standard methods of partitioning. Being that said, further investigation are required in order to clarify the impact of the additional drivers used in NN C-part , or possible missing functions and relationships in the standard methods.
Altogether, this new method provides GPP and RECO estimates that are based on purely data-driven empirical relationships without any assumption of the driver-output relations. This product would be of high interest for process model parameterization and validation as it avoids unwanted circularity, and for this reason, we propose it as a complement to the standard processing in large networks such ICOS, AmeriFlux, or FLUXNET.

DATA AVA I L A B I L I T Y S TAT E M E N T
All the data used in the paper are included in the FLUXNET2015 col-