dataesr/scientific-paragraphs-categorization huggingface dataset:
Model setup To simulate atmospheric CO 2 concentrations that can be compared to observations, the transport model has to be forced not only by terrestrial biospheric fluxes (calculated by ORCHIDEE), but also by other natural (e. g. , ocean) and anthropogenic CO 2 fluxes. We imposed a net emission due to land use change (i. e. , deforestation) of 1. 1 GtC yr -1 , although we also accounted for a larger flux from biomass burning but which was compensated for partly by forest regrowth (see  Peylin et al. , 2016, for more details). The global maps of biomass burning emissions were taken from the Global Fire Emission Database version 3 dataset (  Van der Werf et al. , 2006; Randersen et al. , 2013)  over the period 1997-2010 at a monthly time step and gridded at 0. 5 • × 0. 5 • resolution. The global fossil fuel CO 2 emission products used here were developed by the University of Stuttgart/IER based on EDGAR v4.
In situ flux measurements (F) The NEE and LE measurements come from the FLUXNET global network. We used harmonized, quality-checked and gap-filled data (Level 4) at 68 sites from the La Thuile global synthesis dataset  (Papale, 2006). The site locations are presented in Fig. 1. These ecosystem measurements cover very different time spans, ranging from 1 single year at some sites up to 9 years. They constrain seven PFTs among the 12 natural vegetation types represented in ORCHIDEE: tropical evergreen broadleaf forest -TrEBF (3 sites corresponding to 6 site years), temperate evergreen needleleaf forest -TeENF (16 sites, 45 site years), temperate evergreen broadleaf forest -TeEBF (2 sites, 4 site years), temperate deciduous broadleaf forest -TeDBF (11 sites, 37 site years), boreal evergreen needleleaf forest -BoENF (12 sites, 44 site years), boreal deciduous broadleaf forest -BoDBF (3 sites, 6 site years), and C 3 grassland -C3GRA (21 sites, 56 site years).
Atmospheric CO 2 measurements (CO2) The surface atmospheric CO 2 concentration data come from the NOAA Earth System Laboratory (ESRL) GLOBALVIEW-CO2 collaborative product (GLOBALVIEW-CO2, 2013). The data include in situ measurements, made by automated quasi-continuous analyzers, and air samples collected in flasks and later analyzed at central facilities. In this study, we used monthly-mean values of these measurements  (Peylin et al. , 2016). A period of 10 years of observations over the 2000-2009 period was used from a total of 53 stations located around the world (Fig. 1 ).
Data assimilation framework The data assimilation system associated with the OR-CHIDEE model (ORCHIDAS) has been described in previous studies regarding the assimilation of these data streams alone  (Kuppel et al. , 2012; Santaren et al. , 2014; MacBean et al. , 2015; Bastrikov et al. , 2018)  or their combinations  (Bacour et al. , 2015; Peylin et al. , 2016). The assimilation system relies on a variational Bayesian framework that optimizes ORCHIDEE parameters gathered in a vector x, by finding the minimum of a global misfit function J (x) iteratively. J (x) is a linear combination of the misfit functions associated with each data stream. It is assumed that the errors of observations and on the model parameters are Gaussian and that the data stream errors are independent from each other: J (x) = 1 2 H LMDz • H ORCH (x) -y CO2 T • R -1 CO 2 • H LMDz • H ORCH (x) -y CO2 + H ORCH (x) -y F T • R -1 F H ORCH (x) -y F + H ORCH x -y VI T • R -1 VI H ORCH (x) -y VI + x -x b T B -1 • x -x b , (1) where y o denotes the observation vectors (with o = F (flux), VI (satellite NDVI) or CO 2 (CO 2 concentration); H ORCH and H LMDz are the observational operators of the ORCHIDEE and LMDz models, respectively.
, 1994; Todd-Brown et al. , 2013). Given that we do not have access to that information, neither at the site scale (for assimilation of NEE measurements) nor at the global scale (for assimilation of atmospheric CO 2 concentrations), we use a steadystate assumption where ORCHIDEE has been brought to near equilibrium with a long spin-up of the soil carbon pools. To correct for this bias, the initial state of the soil carbon reservoirs is optimized using a multiplicative parameter of both the slow and passive pools as in  Peylin et al. (2016). The use of these correction factors is a handy way to correct any issues related to the use of our soil organic C model and the soil carbon disequilibrium. Two multiplicative parameters are used depending on the type of data considered (and their associated spatial scale): for in situ flux measurements, we considered site-specific parameters K soilC,site ; for atmospheric CO 2 concentration data, instead of resolving the initial conditions for all LMDz grid cells we scaled the carbon pools for 30 large-scale regions K soilC,reg.
Note that having correct soil carbon pools is less important when assimilating satellite NDVI data because these are more closely related to carbon uptake rather than net carbon flux. In total, up to 182 parameters are optimized depending on the data streams considered. The prior values x b of the parameters are set to the standard values of ORCHIDEE (Table  S3 ). Not all parameters are constrained by all three data streams. In particular, satellite FAPAR/NDVI products inform the timing of phenology of plant vegetation (start and end of the growing season) rather than photosynthesis or respiration with our DA system  (Bacour et al. , 2015; MacBean et al. , 2015). The dependency of each parameter with respect to the assimilated data streams is indicated in Table  1.
Data assimilation experiments Different data assimilation experiments were tested in order to understand the respective constraint brought by each data stream and evaluate their compatibility with each other and with the model. First, each data stream was assimilated separately, and then its combinations with the other two were considered. Second, the three data streams were assimilated all together. The various experiments are described in Table  2 , with the number of data points assimilated and the number of parameters optimized. Indeed, the number of optimized parameters differs with the type of data assimilated as described in Sect. 2. 3. 2 and in Table  1. The assimilations have a high computational cost, with an average value for joint assimilations using all three data streams of about 50 000 h central processing unit time on AMD Rome compute nodes, at 2.
6 GHz with 256 GB memory per node. Two assimilation experiments combining the three data streams were tested: one experiment (F + VI + CO2) with all parameters optimized in a single step; and an additional experiment following a two-step optimization (F + VI + CO2-2steps), as described hereafter. In the first step, the global soil carbon reservoirs were constrained by assimilating atmospheric CO 2 data only and optimizing the two main parameters controlling soil respiration, K soilC,reg and Q 10. In the second step, all parameters but K soilC,reg were optimized from the three data streams: K soilC,reg was retained from the first step, and Q 10 was optimized, but the prior uncertainty Table  1. List of the ORCHIDEE parameters to be optimized and data streams that constrain them (F for in situ flux measurements, VI for normalized satellite NDVI data and CO 2 for atmospheric CO 2 concentration data).
Observation error statistics Like in previous studies with ORCHIDAS, we defined R o as diagonal and computed the variances from the root mean square difference (RMSD) between the data and the a priori ORCHIDEE simulations (i. e. , performed with the model default parameter values) for fluxes and satellite observations. However, it is worth noting that this approach overestimates the variances in order to compensate for any neglected correlations. For atmospheric CO 2 measurements, we followed a different methodology given the large discrepancy in the modeled a priori concentrations with respect to the observed data (i. e. , large bias that increases over time due to biases in the land net carbon sink (too small)). The errors were determined at each site as the standard deviation of the observed temporal concentrations  (Peylin et al.
(7) In principle, the tuning of B and R needs to be performed iteratively for successive values of x a and of the corresponding residuals, until convergence, which is prohibitive in terms of computing time. The estimation of the covariance matrices depends on the mathematical expectation (E), which would require several realizations of the residuals to diagnose the error statistics  (Desroziers et al. , 2005; Cressot et al. , 2014). In this study, only one optimization was performed using one set of a priori parameters for each dataset. We therefore calculated these metrics by averaging the diagonals of the matrices described by both sides of the equations for all available observations  (Kuppel et al. , 2013). This way, both sides are scalar values  (Cressot et al. , 2014). The standard deviation of the errors was determined after a few trials considering the three single data stream assimilation experiments independently: for each DA experiment we started from an initial parameter error set at 40 % of the variation interval for each parameter (as in  Peylin et al.
Optimization performance We measured the efficiency of any assimilation by quantifying the reduction of the cost function as the ratio of the prior to posterior values. It should be noted that the minimum value of the cost function is not expected to be zero given the uncertainty in both the data and model and the limited number of degrees of freedom (number of optimized parameters) allowed. We also looked at the ratio of the norm of the gradient between the prior and posterior misfit functions, as it illustrates the progression towards the expected optimum, for which the gradient is null. The decrease of the norm of the gradient depends on the estimation problem (nonlinearities, number of observations versus number of optimized parameters, constraints of the data on the model processes, etc.
Overall fit to the observations The impact of assimilating one type of observation on all the data streams (including those that are not assimilated) was evaluated for the various assimilation experiments. The reduction of the model-data mismatch (i. e. , reduction in prior RMSD) after assimilation of each data stream (or any combination of them) is illustrated in Fig. 2. The length of the boxes (first and third quartiles) of the whisker plots highlights the spread in misfit reduction across sites/vegetation types. For fluxes, only the impact on NEE is shown, given the choice of the parameters optimized parameters is mostly related to the carbon cycle. Using the parameter values optimized in either the F or VI assimilations has a strong detrimental impact on the simulated atmospheric CO 2 data because the soil carbon pools were not adjusted in these DA experiments.
Therefore, we also analyzed the changes induced on the detrended seasonal cycles of atmospheric CO 2 concentrations (Fig. 2c ) (hence removing the trend using the time series decomposition based on the CCGCRV routine; Thoning et al. , 1989see Sect. S2 and Fig. S1  for representative comparisons of observed vs modeled time series of atmospheric CO 2 concentrations and their associated trend estimation). For a given data stream, the improvement is usually better for the experiment where that data stream is assimilated alone. One noteworthy exception is the assimilation of NDVI alone (VI experiment where only the phenology parameters are optimized) that results in a lower model improvement with respect to NDVI than when it is assimilated in combination with other data streams (where a higher number of parameters are optimized in these joint assimilations, hence improving the timing of phenology and the amplitude of the annual cycle when flux or atmospheric CO 2 data are also assimilated).
2a ). While using optimized phenological parameters retrieved from satellite data alone (experiment VI) degrades the modeled seasonality of NEE as compared to the measurements (median RMSD reduction of -3 %), the optimization with respect to in situ flux data (F ), with additional control parameters, leads to a general improved consistency between modeled FAPAR and satellite NDVI time series (median RMSD reduction of 8 %). The impact on LE is much lower for all DA experiments (median values close to 0 % in all cases, result not shown). One can also note the positive impact of the F and VI assimilations on the atmospheric CO 2 data with median RMSD reductions of 15. 8 % and 11. 2 % respectively for the detrended time series. Such an improvement after assimilation of in situ flux data corroborates the findings of  Kuppel et al.
The parameters retrieved from the CO2 experiment also have a small but positive impact at the site level with respect to NEE (median value of 3 %) and FAPAR (0. 8 %). For the joint assimilation experiment (F + VI, F + CO2, VI + CO2, or F + VI + CO2; Fig. 2 ), the model-data agreement is improved for all assimilated data streams, as expected, while the model degradation relative to the data not assimilated is generally not as severe as compared to the assimilation of individual data stream experiments described above, with the exception of the F + VI experiment. The latter experiment leads to enhanced model improvement compared to when flux and satellite NDVI data are assimilated alone (see Sect. S5). In the simultaneous assimilations involving atmospheric CO 2 data, most of the model improvement concerns CO 2 (Fig.
, 2014). Surprisingly, the use of satellite vegetation indices (VI) leads to a slightly lower residual amplitude bias than when atmospheric CO 2 data are assimilated, albeit a lower number of optimized parameters. For the length of the CUP, the relative model correction appears small for almost all experiments and is lower than what is achieved for the trend and amplitude. Some degradation (increased model-data bias) is even obtained for the cases F and F + CO2. This may be attributed to some inconsistency in the phasing of the CUP derived from the FLUXNET stations and from the atmospheric stations (given differences in the spatial-and temporal-scale constraints brought each data stream). Among the single data stream assimilations, the highest improvement is obtained for VI, where the optimization of the phenological parameters was the only improvement allowed for tuning the model.
Depending on the data assimilated in this study, the posterior GPP ranges from 147 GtC yr -1 (F + VI) to 170 GtC yr -1 (VI + CO2) at the global scale. The largest differences with the prior are obtained for the experiments involving flux and satellite data (alone or the two combined). This is directly linked to large corrections in photosynthesis and phenology parameters for these experiments (see  Sect. 3. 3). In comparison, the assimilations involving atmospheric CO 2 concentration data are more conservative with respect to GPP. Assimilating atmospheric CO 2 data alone lessens the GPP reduction by a factor of about 3 compared to assimilations with F and VI data, and the correction for the joint assimilations using CO 2 data is even lower (see Sect. S6 for details). By propagating the error on the parameters in the observation space (see Eq.
For GPP, the lowest posterior errors are found for the experiments combining F and CO2 data, while experiments F, CO2 and VI + CO2 lead to larger posterior errors. This is due to the fact that (i) F and CO 2 data provide a stronger constraint on the annual mean photosynthesis than VI data and that (ii) F and CO 2 data provide cross constraints on photosynthesis. Experiment VI, in which about 10 times fewer parameters are optimized and targeting primarily the timing of phenology, results in the highest posterior GPP errors (although still a reduction from the prior). Finally, one can observe that the posterior errors are higher in the tropics for both NEE and GPP (and the reduction compared to the prior error is lower), which is even more prominent in the experiments using in situ flux data alone or with satellite data, a direct consequence of the lower data availability (eddy-covariance measurements) to constrain the model parameters for tropical PFTs.
Parameter estimates and associated uncertainties Figure  6  shows the impacts of the different assimilation experiments on a subset of the retrieved parameter values and their associated uncertainties (the remaining parameters are shown in Fig. S2 ). While the stepwise study showed only few changes in the parameter estimates between the sequential steps (and hence as a function of the data stream from which the parameters were constrained)  (Peylin et al. , 2016) , our results show a large variability between the assimilation experiments. For most parameters, the highest departures from the prior values are obtained for the single data stream assimilations. Higher changes are obtained for flux or satellite data as compared to the estimates retrieved with atmospheric CO 2 data alone which remain closer to the prior values.
6 ), with below 0. 5 % difference, and consistent with the estimates of the other joint assimilation experiments. For some parameters/PFTs, the direction of the departure with respect to the prior value (increase or decrease) may differ depending on the data stream assimilated. At the first order, the estimated parameter uncertainties decrease with the number of observations assimilated, as expected from Eq. (  3 ), and given that the observations are treated as independent data. However, given that the estimated parameter errors strongly depend on the setup of B and R matrices and that we did not use error correlations in these matrices, we should only focus on the relative error reduction between experiments. The uncertainty reduction achieved through the assimilation of atmospheric CO 2 data is usually lower than when flux and satellite data are assimilated alone and typically varies between 10 % and 60 % for most photosynthetic and phenological parameters.
Relative constraints brought by the different datasets We now quantify the impact of each of the three data streams on the analysis using the global average observation influence (quantified by OI) and information content (DFS) metrics defined in Sect. 2. 4. 4. We recall that OI (i. e. , trace of S normalized by the number of observations) gauges the aver-  age influence that each single observation has on the analysis, while the relative DFS measures the overall weight of one data stream in the optimization (the difference between OI and DFS is due to the number of observations assimilated;  Cardinali et al. , 2014). OI and DFS are determined for the joint assimilation experiments combining the three data streams. Because of the very large number of observations (above 300 000) involved in the assimilation, only the diagonal elements of the influence matrix (Eq.
10) can be calculated. The trace of S measures the equivalent number of parameters and is equal to 132. Such a value, lower than the number of parameters (182), indicates that the optimized parameters may not be fully independent (although parameter error correlations have been ignored in our B matrix) as already reported in  Kuppel et al. (2012)  or that some are not constrained during the optimization process (as for instance LAI MAX , whose estimate remains at its a priori value for some PFTs; Fig. S2 ). The values of OI are provided in Table  4  for flux, NDVI and atmospheric CO 2 data. With about the same number of observations considered (Table  2 , last column), one in situ flux measurement has about 10 times more weight than one NDVI observation. This is a consequence of the larger number of parameters constrained by flux measurements than by NDVI data in our setup.
The highest influence is found for at-mospheric CO 2 data, the relative weight of one atmospheric CO 2 measurement being 4 times larger than that of one flux observation, despite the much lower number of data assimilated. Again, this is a consequence of the large weight of the mismatch between the a priori-simulated trend and the observed trend in the atmospheric CO 2 data, which is drastically reduced through the optimization. However, the smaller number of atmospheric CO 2 data assimilated, compared to flux and NDVI datasets, reduces the overall constraint on the analysis provided by atmospheric CO 2 data, as gauged by its relative DFS. Hence, our optimization is mainly controlled by flux data, which have an overall contribution of about 75 %, which is about 5 times larger than the constraint brought by atmospheric CO 2 data and 7 times larger than that of satellite NDVI.
Realism of the regional-to global-scale C fluxes The overarching objective of the study was more about assessing how to make the best of a synergistic exploitation of different data streams within a consistent assimilation framework rather than achieving an up-to-date reanalysis of the global carbon fluxes, especially since we focused on a limited dataset both in terms of temporal coverage (no atmospheric CO 2 data nor satellite data after 2010, no in situ flux data beyond 2007) and of informational constraint. Indeed, we did not assess the potential of other data that can bring relevant (and possibly more direct) additional constraints on the dynamics of terrestrial carbon stocks and fluxes, such as aboveground biomass  (Thum et al. , 2017)  or solar-induced fluorescence  (Bacour et al. , 2019b) , which have already been investigated with ORCHIDAS and with an updated version of the ORCHIDEE model.
The expansion of the assimilated datasets to provide the most up-to-date constraint on modeled carbon fluxes will be the subject of future work. In spite of these limitations, we saw that the regional-/global-estimated NEE and GPP budgets are realistic and in agreement with independent estimates. There are still important differences in the model predictions for the different assimilation experiments (and we have not attempted to identify what was the most reliable optimized model, which would require the use of an ensemble of independent data, an effort beyond the scope of this paper). Still, our optimized simulations allow for a more in-depth exploration of the partitioning of the land carbon budget between the northern extratropics and the tropics. From the global carbon budget, a discrepancy exists between the partition estimated by the atmospheric CO 2 inversions and by the terrestrial biosphere models  (Kondo et al.
Exploiting TBM databases more directly related to regional soil carbon contents, such as the Harmonized World Soil Database (HWSD) (FAO/IIASA/ISRIC/ISSCAS/JRC, 2012), the International Soil Carbon Network  (Nave et al. , 2016)  or the global soil respiration database  (Jian et al. , 2021) , is not straightforward because of the errors associated with these datasets  (Todd-Brown et al. , 2013)  and inconsistencies between the estimated quantities and the model state variables and underlying processes (as for instance the depth of the soil carbon). In any case, what is sorely needed is data that track changes in C stocks over long time periods. Still, it is of primary importance for the science community to endeavor to bridge the gap between state-of-the art estimates of soil carbon stocks and the quantities that TBMs simulate over the historical period.
We investigated how the different combinations of data streams constrain the parameters of the ORCHIDEE land surface model and by consequence the simulated historical spatial and temporal distribution of the net and gross carbon fluxes (NEE and GPP), as well as FAPAR and atmospheric CO 2 concentrations. We quantified how the combination of these data streams (two by two or all together) impacts the reliability of the model predictions. Although it leads to lower fitting performances with respect to the assimilation of any individual dataset (because the optimization seeks a tradeoff solution between all data streams), the simultaneous assimilation of the three data streams is found to be the most consistent approach. In particular, it avoids model overfitting, which can degrade the model predictions with respect to data streams not assimilated.