dataesr/scientific-paragraphs-categorization huggingface dataset:
8 wt% and 2 wt%, 234]ectively. Note that a water content in the lithospheric mantle of 2 wt% corresponds to approximately 15% serpentinization, which is consistent with serpentinization degrees inferred from geophysical observations (see ref. 33). The phase relations are computed with the thermodynamic code Perple_X[24]using the thermodynamic database of ref23. More details on the thermodynamic calculations are given in Supp. Info. For the complex peridotite model, we combine the phase assemblages predicted by thermodynamic calculations to experimental data on phase relations at high PT obtained from natural serpentinites containing magnetite (see below).
Results Magmas are differentiated, with whole rock compositions in the andesitic and dacitic eld (Fig. 1c-d ; Supplementary Note and Supplementary Table  1 ). Opx are present in the dacite and andesite. 70% are unzoned (Supplementary Fig. 2 ). Zoned crystals were classi ed according to the type of zonation observed: either a single one (normal or reverse zonation) or a multiple one (Fig. 1g-j ; Supplementary Fig. 2 ). In the dacite, the majority of opx are normal-zoned while, in the andesite, reverse-zoned opx are more abundant (Supplementary Fig. 2 ). Few opx show reactional borders or remelting textures (Supplementary Fig. 2 ). Out of the 1832 opx crystals mounted, 112 zoned opx compositions and 108 unzoned opx compositions were investigated in the dacite and andesite (Supplementary Data 1).
In the experimental study of D we used to model our timescales, they hypothesized an fO 2 dependence similar to olivine with an exponent of 1/6 60. This fO 2 dependence has been incorporated in the formulation of D in other works  57, 61  but this dependence has never been veri ed experimentally for opx with En contents of our dataset  62, 63. For En-rich opx of the mantle (En 98 and En 91 ), between 870-1100°C and oxygen fugacity controlled (fO 2 = 10 - 11 -10 - 7 Pa), an fO 2 dependence of an exponent between 0 and 1/20 was found to be appropriate from their data 63. The estimate for D Fe-Mg used here 60 , when extrapolated to higher temperatures of the study on En-rich opx without considering any fO 2 dependence appear to be consistent with their dataset  63. The oxygen fugacity dependence on Fe-Mg interdiffusion coe cient is still under discussion as incorporating it can give longer timescales  47  but no study has yet examined the appropriate En compositions of the opx and appropriate oxygen fugacity ranges.
Therefore, the parametrization of the interdiffusion coe cient D used is the one determined on opx compositions close from our dataset without an fO 2 correction 60 , as done in previous studies  25, 55. D is de ned as follows in Eq. (  1 ) 60 : logD = -5. 54 + 2. 6X Fe -12530 T #(1) The main assumption of the model is that the initial pro le between two zones follows a step function, which is modi ed by diffusion to form sigmoidal concentration gradients 64. This steep pro le becomes broader with time, due to interdiffusion of the different elements. With the eruption, the diffusion is stopped. Maximum timescales are then estimated 64. Pro les with smooth compositional variations of sigmoid shape are studied and those displaying growth dominant pro les are isolated. These pro les display linear trends or curves and bumps on the pro les and are not consistent with a diffusion sigmoid calculated for the given initial compositional contrast  25  (Supplementary 5).
The Al and Ca gradients are compared to the Mg#, as Al and Ca diffuse slower than Fe-Mg 47,57,65,66 (Supplementary Data 1; Supplementary Fig. 5 ). The model that we used consists in the intercalibration of chemical composition pro les of zoned opx, measured with an electron microprobe, with grayscale pro les on high-resolution north-south oriented images taken with SEM (spatial resolution of 755 nm)  25, 55  , following the methodology developed in two other studies  57, 59. The Mg number is the parameter that seems to best control the grayscale of the images 50,59. The ImageJ software (https://imagej. nih. gov/ij/; version 1. 52a) is used to produce the grayscale pro les. They are averaged over a certain width to reduce the noise (between 20 to 50 µm, depending on the crystal rim size)  25, 55.
After this intercalibration, the pro les are modelled to obtain the timescales of interdiffusion. As the diffusion coe cient used is composition-dependent, the curve shape of the pro le is determined by the degree of compositional contrast but the curve width is dependent on time  25, 59  Based on a library of pro les of particular values of compositional contrasts, the shape of the curve that ts the best to the one obtained by intercalibration will be searched in this library according to the compositional data  25. The curve will then be scaled in width to solve for time and the timescale corresponding to the zonation is obtained, as the scaling factor is an expression of crystal composition, the diffusivity used for the library of the pro les, the diffusivity in the crystal and the ratio between the timescales of the database and the timescale associated with the zoning in the crystal  25  (with the timescale of the crystal being solved as the only unknown).
Acknowledgments We would like to thank  S. Hidalgo  for her help with sample preparation,  M. Fialin  and  N. Rividi  for assistance during EPMA analyses, as well as  O. Boudouma  and  S. Borensztajn  for SEM imaging. We also thank  D. M. for providing the les to model the timescales of the opx. We also acknowledge that the image from Figure 1a is from a  Digital Elevation Model  provided by the  Polar Geospatial Center  under NSF-OPP awards 1043681, 1559691, and  1542736  and Figure 1b from  NASA Earth Observatory  images by  Robert Simmon , using Landsat 8 data from the  USGS Earth Explorer. Thank you also to  R. Bajou  for the discussion on the cumulative number of earthquakes. This project was funded by the  Institut de physique du globe de Paris  ( Doctoral grant  from the  French Ministry of Higher Education and Research and Innovation)  and  ANR V-Care-18-CE03-0010.
, 2015) , have also been used to constrain the model parameters at various spatial scales. Finally, atmospheric CO2 mole fraction measurements have been assimilated to provide valuable information on large-scale net ecosystem exchange (NEE)  (Rayner et al. , 2005; Koffi et al. , 2012). In the early days of DA studies, most focused on the assimilation of a single data stream (e. g. , only NEE). Then, assimilations with multiple different C cycle related datasets have soon been considered  (Moore et al. , 2008; Richardson et al. , 2010; Ricciuto et al. , 2011; Keenan et al. , 2013; Thum et al. , 2017; Knorr et al. , 2010; Kaminski et al. , 2012; Kato et al. , 2013; Bacour et al. , 2015; Peylin et al. , 2016). The underlying motivation behind assimilating multiple data streams is that using a greater number and diversity of observations should provide stronger constraints on model parameters, including a wider range of processes, hence resulting in a greater reduction in model uncertainty.
However, many previous studies that assimilated multiple datasets hardly considered potential incompatibilities between the model and the observations  (Bacour et al. , 2015; Thum et al. , 2017). Besides, a few have quantified the actual benefit of assimilating multiple data-sets compared to the single data stream assimilations, in particular in the context of global scale C cycle DA experiments. The assimilation of multiple data streams can be done either sequentially, in which one observation type is assimilated at a time, or simultaneously (joint assimilation approach or "batch" strategy as defined in  Raupach et al. , 2005) , where the model is calibrated with all data included in the same https://doi. org/10. 5194/bg-2022-109 Preprint. Discussion started: 31 May 2022 c Author(s) 2022.
CC BY 4. 0 License. optimization (e. g. Richardson et al. , 2010; Kaminski et al. , 2013; Schürmann et al. , 2016). While the latter joint assimilations are more optimal as it maximizes the consistency of the model with the whole of the datasets considered  (Richardson et al. , 2010) , sequential approaches remain appealing for modelers: They require less initial technical investment and enable easier assessment of the impact of the data stream assimilated successively onto the optimized variables. Both approaches however face similar challenges, like defining the model-data uncertainty and hence the weight of that each dataset has on the optimization outcome (see, e. g. ,  Richardson et al. , 2010; Keenan et al. , 2013; Kaminski et al. , 2012; Bacour et al. , 2015; Thum et al. , 2017; Peylin et al.
, 2006)  based on pre-calculated transport fields for assimilating atmospheric CO2 concentration data. ORCHIDEE and LMDz are the terrestrial and atmospheric components of the IPSL Earth System Model  (Dufresne et al. , 2013). As a preliminary step, we tuned prior and observation error statistics: we enhanced their realism by making them consistent with the differences between prior model simulations and observations. We then conducted different assimilation experiments in which each data stream was assimilated alone on in combination (for all combinations of datasets) to assess what the impact of each assimilation experiment was on: a) the fit to each dataset (including datasets included and excluded from the assimilations); b) on the resultant regional to global scale gross and net CO2 flux budget (NEE and https://doi.
org/10. 5194/bg-2022-109 Preprint. Discussion started: 31 May 2022 c Author(s) 2022. CC BY 4. 0 License. GPP); and c) on the optimized parameter values and uncertainties. We complemented our analysis by comparing our regional to global scale C budgets with independent process-based model and atmospheric inversion estimates from the Global Carbon Project's 2020 Global Carbon Budget  (Friedlingstein et al. , 2020). We completed our assessment of the impact of each assimilation scenario by quantifying the observation influence of each of the three data streams on the joint assimilation in which all three datasets were included in the optimization. By determining the overall constraint brought by each data set on the inversion, this final analysis allowed us to highlight the importance of atmospheric CO2 data in the optimisation of soil C pools despite the smaller number of observations assimilated.
Model Set-up The set-up of the simulations performed with ORCHIDEE depends on the data assimilated. The model is run at site scale for the assimilation of eddy-covariance measurements, at spatial resolution 0. 72°173 for the assimilation of the satellite NDVI data, and at the resolution of the atmospheric transport model LMDz (3. 75°x2. 5°) for the assimilation of atmospheric CO2 measurements. The Olson land cover classification at 5 km is used to derive the PFT fractions at each spatial resolution, but for the flux tower simulations where the proportion of each PFT is set based on expert knowledge. For satellite pixels and global simulations, ORCHIDEE is forced using the 3-hourly ERA-Interim gridded meteorological forcing fields  (Dee et al. , 2011 ) (aggregated at 3. 75°x2. 5°when assimilating atmospheric CO2 concentrations).
For the flux tower simulations, the model is forced by local measurements of the meteorological variables at a half-hourly time step. For each spatial resolution, a prior spin-up simulation was performed by recycling available forcing data. The objective was to bring the different soil carbon reservoirs to "realistic" values, albeit the spin-up runs result in neutral net carbon flux by construction. Each spin-up simulation was then followed by a transient simulation (starting from the first year of measurement for each data stream) and accounting for the secular increase of atmospheric CO2 concentrations.
Model set-up To simulate atmospheric CO2 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 CO2 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 compensated 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)   (Peylin et al. , 2016).
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 location is presented in Figure  1. These ecosystem measurements cover very different time spans, ranging from one single year at some sites up to nine years. They constrain seven PFTs among the twelve 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 sites-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 C3 grassland -C3GRA (21 sites, 56 site-years).
Atmospheric CO2 measurements (CO2) The surface atmospheric CO2 concentration data come from three databases: The NOAA Earth System Laboratory (ESRL) archive (ftp://ftp. cmdl. noaa. gov/ccg/co2/), the CarboEurope IP project (http://ceatmosphere. lsce. ipsl. fr/database/index_database. html), and the World Data Centre for Greenhouse Gases of the World Meteorological Organization Global Atmospheric Watch Programme (http://gaw. kishou. go. jp). The data include in situ measurements, made by automated quasicontinuous analysers, 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). Ten years of observations over the 2000-2009 period were used from a total of 53 stations located around the world (Figure  1 ).
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). 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 KsoilC,site; for atmospheric CO2 concentration data, instead of resolving the initial conditions for all grid cells we scaled the carbon pools for 30 large scale regions KsoilC,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.
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 are assimilated altogether. 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 §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 hr 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 2-step optimization (F+VI+CO2-2steps), as described hereafter. In the first step, the global soil carbon reservoirs are constrained by assimilating atmospheric CO2 data only, and optimizing the two main The results of these assimilations were compared to the companion study of  Peylin et al. (2016)  in which the same data streams were assimilated in a sequential/stepwise approach: NDVI data were assimilated first, then in situ flux measurements, and finally atmospheric CO2 concentration measurements. In addition, there are a few differences in the set-up compared to the present study, in particular the use of only three years of atmospheric CO2 data (cf.
Observation error statistics Like in previous studies with ORCHIDAS, we defined R 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 CO2 measurements, we followed a different methodology given the large discrepancy in the modelled a priori concentrations: the errors were determined as a function of the observed and modelled temporal concentration variability at each site  (Peylin et al. , 2016) , thus neglecting the structural errors of the terrestrial model.
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 were 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.
Optimisation 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 (non-linearities, number of observations versus number of optimized parameters, constraints of the data on the model processes, etc.
Assessment of the information content of each data stream For the joint assimilations using the three different data streams, we further analyzed the influence matrix S that quantifies their leverage on the model-data fit  (Cardinali et al. , 2004 ): = -. .  (10) A diagonal element Sii is the rate of change of the simulated observable i with respect to variations in the corresponding assimilated observation i. Sii is referred to as "self-sensitivity" of "self-influence". A zero self-sensitivity indicates that this i th observation does not contribute to improving its simulation We analyzed the trace (i. e. the sum of all diagonal elements) of S that quantifies a measure of the amount of information that can be extracted from all observations / all data streams. We used two derived quantities: the global average observation influence (OI) and the relative degrees of freedom for signal (DFS) associated with the data stream o, which measures its relative contribution to the fit.
For the joint assimilation experiment (F+VI, F+CO2, VI+CO2, or F+VI+CO2; Figure  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 simultaneous assimilation of flux measurements and satellite NDVI data (F+VI) leads to enhanced model improvement compared to when these data are assimilated alone (cf. Supplementary Text S5). In the simultaneous assimilations involving atmospheric CO2 data, most of the model improvement concerns CO2 (Figure  2c ) while the benefit for the fluxes and FAPAR/NDVI is weak (RMSD reduction below 3%). Noteworthy, the 2-step assimilation F+VI+CO2 (see Section 2.
Specific improvements at CO2 stations With the default parameter values, the simulated fluxes by ORCHIDEE lead with LMDZ to overestimates the (trend) by about 1 ppm. yr -1. When assimilating atmospheric CO2 data, most of the parameter correction aims at reducing this bias. This is mostly achieved by tuning the regional KsoilC_reg parameters: the net land carbon sink is increased globally in order to match the observed trend at most stations (reducing the bias from around 1 ppm. yr -1 to 0. 1 ppm. yr -1 ). Compared to the improvement in the bias in the trend, the improvements (reduction in bias) in the amplitude of the CO2 seasonal cycle and in the CUP length (Figures  3b  and c ) are marginal. Note that our joint DA experiments lead to significantly lower trend biases compared to the stepwise approach, which is probably due to the longer period of the atmospheric CO2 data considered (10 years vs 3 years for the stepwise).
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 deduced 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 optimisation of the phenological parameters was the only improvement allowed for tuning the model. For the joint assimilations, those combining the three data streams provide the best performance and perform better than the stepwise approach. Among the joint assimilations with three data streams, the 2-step approach results in the largest reduction in amplitude and CUP bias, but, on the other hand, the larger trend bias.
The correction of the CO2 trend in the first step of the 2-step approach, with the optimisation of the soil disequilibrium (KsoilC_reg parameters), tends to favor a tropical C sink, which has a direct impact on the atmospheric CO2 stations of both hemispheres. Although such an approach offers the advantage of correcting more substantially the seasonal cycle of atmospheric CO2 (see section 3. 1. 3), it may still be suboptimal for partition of the terrestrial C sink (between northern and tropical regions). With a global scale budget of 171 GtC. yr -1 for GPP, the prior ORCHIDEE model is on the high range of recent estimates of the global GPP, as synthesized in  Anav et al. (2015) , the mean value of which being around 140 GtC. yr -1. Depending on the data assimilated in this study (note that we include all assimilation experiments in estimating the posterior GPP), the posterior GPP ranges from 147 GtC.