dataesr/scientific-paragraphs-categorization huggingface dataset:
Precipitation data were used wherever available, but considered optional in order to not further reduce the dataset. Long-term mean annual temperature (MAT), precipitation (MAP, Supplementary Fig. 2 ), coordinates (Supplementary Fig. 1 ), average number of days with snow cover, and International Geosphere-Biosphere Programme (IGBP) ecosystem type were provided by site operators with the respective dataset and via the European Fluxes Database Cluster (http://www. europe-fluxdata. eu/home/ sites-list). Where MAT or MAP were missing, the mean annual temperature and precipitation of the time series were used. The dataset comprises 176 sites with 1 to 25 valid years each, totalling to 1167 site-years. The dataset covers 12 out of 17 IGBP land-cover classes, missing only the classes urban, mosaic, barren, water body and deciduous needleleaf, and covers the global biome variability in temperature-precipitation space (Supplementary Material 1.
1). Imputation values for data gaps in SW in and NEP, and its partitioning into gross primary production (GPP) and ecosystem respiration (R eco ), were provided along with the above-mentioned data product, with methodology as described by Pastorello, et al. 73. The dataset does not include imputation values for data gaps in SW out. A bias in the resulting albedo α s was avoided by applying an adaptive window imputation method screening the joint time series of SW in and SW out for time periods of approximately constant α s (Supplementary Methods 2,  76  ). Albedo and (half-hourly) minimum temperature of these periods were also used to identify the possible presence of snow or hoar frost (Supplementary Fig. 3 ). Only full calendar years of successfully gap-filled data were retained in order to ensure that resulting time averages of the above-mentioned variables for each station were unbiased by seasonality.
For any time interval of aggregation (e. g. year), α s was computed as the ratio of the interval's average SW out and SW in , which is equivalent to a flux-weighted average of α s and differs from a physically misleading arithmetic time average of instantaneous α s. The fraction of absorbed photosynthetically active radiation (fAPAR) was taken at 500 m spatial and 4 day temporal resolution from the MODIS fAPAR/LAI (MCD15A3H v6. 1) product. Data were downloaded through the AppEEARS web portal  77  for the tower pixels and, after quality control (MODLAND_QC bits equal to zero) averaged to a monthly time scale (gaps were linearly interpolated if necessary). For five sites (DE-Akm, DE-Rus, FR-Tou, US-Snf, US-Wpt) a neighbour pixel had to be downloaded as the tower pixel was not processed correctly.
Out of the potential partner sites identified for each single site according to these criteria, we selected for: Scenario 1 (NEP optimization): The site with the highest (multi-) annual average NEP. Scenario 2 (albedo optimization regardless of NEP): The site with the highest (multi-) annual average α s. Scenario 3 (balanced): The site with the largest relative improvement parallel to the axis connecting the global minima and maxima of α s and NEP, i. e. x ¼ 1 ffiffi ffi 2 p α new À α site α max À α min þ 1 ffiffi ffi 2 p NEP new À NEP site NEP max À NEP min ; ð1Þ where subscripts max and min indicate the maximum and minimum value in the dataset shown in Fig. 1 , such that the denominators normalize albedo and NEP to their total range. Table  1  Albedo α and net ecosystem productivity NEP of land cover types and species during the snow-free growing season (NEP > 0, index act) with standard errors for those groups including more than one site, and deviations from it (Δ) during snowfree dormancy (NEP ≥ 0, index dorm) and snow or hoar-frost -affected periods (index snow, see methods for details).
Surface albedo effect at the top of the atmosphere in scenarios. The modification of net shortwave radiation in response to an albedo change at the same level (e. g. the land surface or top-ofatmosphere TOA) under unchanged incoming radiation is given by SW in Δα at this level  14, 24. In our case, SW in and Δα are measured at the land surface, while effects at the TOA are of interest for comparison to CO 2 effects. This would require an additional correction either for cloud masking of the surface (if using SW in at TOA), or for atmospheric absorption (if using SW in Δα in at the surface). We therefore replace SW in by model-derived, publicly available surface albedo kernels  87. We extracted monthly kernels for all sites from four available kernel datasets 31,87-89 and compared them, their ensemble mean, and monthly measured SW in to each other (Supplementary Methods 3).
All kernels were highly correlated (R² ≥ 0. 92) to our measured SW in and related to it by factors between 1. 03 and 0. 81 (ensemble mean 0. 89, Supplementary Fig. 4 ). We used the largest-kernel  88  and lowest-kernel 87 dataset in the following, adding a strong albedo effect and a weak albedo effect ensemble member to the ensemble of assumptions on CO 2 -based effects described in the previous subsection. While an annual resolution is sufficient to describe the slow evolution of net radiative effect shown in Fig. 2  and the cumulative effect of CO 2 , annual albedo changes multiplied with annual kernels can lead to systematic errors if Δα s and k co-vary over the course of a year  27, 30. We therefore use a monthly sub-loop, corresponding to the highest temporal resolution at which globally distributed albedo kernel datasets are available, to compute the TOA albedo effect for each site and scenario from k and α s for each month of the year.
Where a site observation or kernel dataset covered multiple years, monthly climatological k and α s were determined by averaging across all available years. Monthly effects were integrated over the year before combining with CO 2 -based radiative forcings from the previous subsection to yield the net TOA radiation change, ΔR, for the respective site, year and ensemble member. The final net effect of surface albedo further differs from kernel-based instantaneous shortwave effects at the TOA, because any land use change modifies surface temperature, dust emission  70  , sensible and latent heat fluxes and resulting atmospheric temperature, humidity, particle and cloud profiles, which in turn affect SW in and longwave radiation. Quantification of the resulting effective radiative forcing (ERF  69  ) would require a global circulation model run with the exact spatiotemporal pattern and magnitude of all surface property changes accompanying the albedo change for each scenario.
To provide a rough estimate of the importance of the adjustments distinguishing ERF from IRF, we refer to a study comparing among others the net effect of the same land use change across 14 Coupled Model Intercomparison Project (CMIP6) models  31. Methods for Table  1. For each site, the time windows from the gap-filling process (defined by a temporally stable albedo and of variable length, see Methods: Dataset) used to determine the presence of snow were further distinguished into active (NEP > 0) or dormant (NEP ≤ 0) when snow-free. Sites were grouped according to IGBP land cover type, dominating ( ≥ 75% coverage) species where available, and the two contrasting pairs needle-vs. broadleaf and deciduous vs. evergreen for forests. Species information was compiled from Flechard, et al. 90  and Musavi, et al.
Data availability This study is based on publicly available data from the following sources: FLUXNET2015 product  73  : https://fluxnet. org/data/fluxnet2015-dataset/. ICOS Warm Winter 2020 product (in FLUXNET2015 format)  74  : https://doi. org/10. 18160/2G60-ZHAK. AMERIFLUX FLUXNET2015 compatibility product  75  : https://ameriflux. lbl. gov/ introducing-the-ameriflux-fluxnet-data-product. MODIS fAPAR  77  : https://appeears. earthdatacloud. nasa. gov. CAM5 radiative kernels  89  : https://zenodo. org/record/997902  92. HadGEM2 radiative kernels  88  : https://doi. org/10. 5518/406  93. CACK 1. 0 radiative kernels  87  : https://doi. org/10. 6073/pasta/d77b84b11be99ed4d5376d77fe0043d8  94. HadGEM3-GA7. 1 radiative kernels  31  : https://doi. org/10. 5281/zenodo. 3594673  95. These data were quality filtered and aggregated as described in the methods section.
Acknowledgements We thank  Talie Musavi ,  Max Planck Institute for Biogeochemistry, Jena, Germany , for providing species information originally collected for Musavi, et al. 91 ,  Manuel Acosta ,  CzechGlobe, Czech republic , for provision of the data of site CZ-Krp, and the  ICOS research infrastructure  for data provision. L. Š. and N. K. acknowledge support by the  Ministry of Education, Youth and Sports of the Czech Republic  within the  CzeCOS program  (grant number  LM2018123 ) and SustES-Adaptation strategies for sustainable ecosystem services and food security under adverse environmental conditions ( CZ. 02. 1. 01/0. 0/0. 0/16019/0000797 ). In 2009 funding for the Australian sites was provided to the  Australia Terrestrial Ecosystem Research Network (TERN)  (http://www. tern. org.
au)through the  Australian government's National Collaborative Research Infrastructure Strategy (NCRIS) , which provides support for many OzFlux sites along with other capabilities; an overview of all Australian sites is given in Beringer, et al. 97. Data acquisition for FR-Lam were mainly funded by the  Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS-INSU)  through the  ICOS  and OSR SW observatories. Facilities and staff were also funded and supported by the  University Toulouse III -Paul Sabatier , the  CNES  and 732  IRD (Institut de Recherche pour le Développement). We are grateful to CESBIO technical team for their technical support at the field and to  Aurore BRUT  for data processing. We extend special thanks to  Ecole d'Ingénieur de Purpan  for accommodating the measurement devices in the plot at FR-Lam.
Results The reconstruction of the glacier-wide MBs of these glaciers from our simple temperature-index model shows good agreement with data (Fig. 1 ). Using these data, we calculated the MB sensitivities to temperature and winter precipitation at 2,750 metres and 3,100 metres on the Argentière and Sarennes glaciers respectively. These altitudes were selected because they correspond to the approximate center of the glaciers. For each day of each series, we calculated an annual MB anomaly by adding a temperature anomaly or a precipitation anomaly. We report the results in Figure  2  to mirror Figure  3  of  Bolibar et al. (2022)  and make the comparison easier. We also ran these numerical experiments at different altitudes and over the entire glacier surface (Fig. 3 ). From our experiments, we found first that the response of MB to temperature, using a temperature-index model, is not linear, contrary to the conclusions of  Bolibar et al.
Acknowledgements This study was funded by  Observatoire des Sciences de l'Univers de Grenoble (OSUG)  and  Institut des Sciences de l'Univers (INSU-CNRS)  in the framework of the  French GLACIOCLIM (Les GLACIers , un Observatoire du CLIMat) program. We thank all those who conducted the field measurements. We are grateful to  Cary Bartsch  for reviewing the English.
Conclusion The Lower Sandstones of the Ennedi is suggested to be essentially Ordovician in age. Here, the Palaeozoic succession illustrates the depositional conditions that occurred across one the most proximal segment of the north-facing Gondwana platform. It shows in particular that the Ordovician transgressions, well-characterized in North Africa  [Ghienne et al. , 2007] , must have reached northern Chad. This archive also expands the area of glacially-controlled deposition to the south and east with respect to the pre-existing dataset. It confirms the eastward prolongation of the North African palaeoglacial framework in the eastern Sahara and, most likely, up to eastern Arabia and Ethiopia  [Le Heron and Howard, 2012]. As far south as the Ennedi, glaciomarine sedimentation and northward-flowing ice flows are recognized.
Acknowledgments The authors wish to thank the  Governor of the Ennedi-Ouest  and the  Prefect of Fada , as well as the staff of the  Ennedi African Park  for their hospitality and for facilitating the fieldwork. They also gratefully acknowledge the  University of N'Djaména , the  CNRD (Centre National de Recherche pour le Développement) , and the  French Embassy in N'Djaména  for their financial and technical support. The authors are grateful to  G. Meinhold  and  J. N. Proust  for their careful reviews, which helped improve an earlier version of the manuscript.
, 2014; Duggen et al. , 2009) , the Algerian Hoggar  (Aït-Hamou et al. , 2000; Allègre et al. , 1981; Maza et al. , 1998) , the Northeastern African volcanism (Sudan and Egypt;  Lucassen et al. , 2008 , the Canary, Lustrino & Wilson, 2007) , and the Cape Verde islands (GEOROC database; "https:// georoc. eu"). According to magma source modeling (Figure  8 ), the generation of subalkaline and transitional magmatism for the Middle Atlas  (Duggen et al. , 2009; El Azzouzi et al. , 2014 , 2010)  would require higher melting degrees (2%-5%) of a mantle source similar to that inferred for Saghro and Middle Atlas alkaline lavas. The anorogenic Rif magmas are compatible with a genesis using similar source and melting conditions of the Sahgro and Middle Atlas alkaline lavas (Figure  7 ; Guilliz, Oujida, Gourougou;  Duggen et al.
e. , thickness and density) provided in Table  S1  and Figure  S1a  in Supporting Information S1. To solve the typical convection equations for the conservation of mass and momentum in the infinite Prandtl number regime, we use the ASPECT finite element code (ver. 2. 2. 0;  Bangerth et al. , 2020; Clevenger & Heister, 2021; Fraters et al. , 2019; He et al. , 2017; Heister et al. , 2017; Kronbichler et al. , 2012; Rose et al. , 2017) , considering visco-plastic rheology and incompressible fluids. Details of the numerical calculations including the rheological formulations are provided in the Appendix A. Any numerical model has to make simplifications, sometimes because of technical limitations, and more importantly and generally, to reduce complexity, and to be able to extract the basic underlying physics in hopefully general ways.
Acknowledgments The authors warmly thank  Paolo Ballato ,  Frank Pazzaglia  and  Maria Laura Balestrieri  for the inspiring discussions concerning the evolution of the Atlas system over the years. We acknowledge the reviewer  Ioannis Baziotis  and an anonymous reviewer for their comments and criticism that significantly improved the manuscript and the Editor  Paul Asimow  for constructive suggestions. TWB was partially supported by  NSF EAR  Grants  2045292  and  1925939. We thank the  Computational Infrastructure for Geodynamics  (https://geodynamics. org/), which is funded by the  National Science Foundation  under award  EAR-0949446  and  EAR-1550901  for supporting the development of ASPECT.
We have compiled chronological and glaciological information including novel orbital age markers from new data on EDC ice core as well as accurate firn modeling estimates in a Bayesian dating tool to construct the new AICC2023 chronology. The average uncertainty of the ice chronology is reduced from 2,500 years to 1,800 years in AICC2023 over the last 800 kyr. The new timescale diverges from AICC2012 and suggests age shifts reaching 3,800 years towards older ages over Marine Isotopes Stages (MIS) 5, 11 and 19. But, the coherency between the new AICC2023 timescale and independent chronologies of other archives (Italian Lacustrine succession from Sulmona Basin, Dome Fuji ice core and northern Alpine speleothems) is improved by 1,000 to 2,000 years over these time intervals.
To determine EDC age scale, they used various orbital dating constraints including: 39 tie points attached to a 6,000 years uncertainty derived from δ 18 O atm tuning to 5,000 years delayed precession between 800 and 300 ka BP, 20 tie points associated with a 4,000 years uncertainty from δO 2 /N 2 alignment to local summer solstice insolation between 800 and 300 ka BP, and 14 tie points linked to an uncertainty between 3,000 and 7,000 years using TAC synchronised to integrated summer insolation between 430 and 0 ka BP. However, due to the lack of data for the orbital dating approach, AICC2012 uncertainty is of 2,500 years on average, reaching 8,000 years at the bottom of the core. The origins of AICC2012 uncertainty can be divided in the following points: (i) discrepancy between δ 18 O atm ,δO 2 /N 2 and TAC series and their orbital target; (ii) discontinuity and poor quality of the δO 2 /N 2 and TAC records; (iii) uncertainty on the phasing between δ 18 O atm and precession; (iv) poor constraint on the LID scenario due to a disagreement between δ 15 N data and firn modeling estimates  (Bréant et al.
(2017)  can be employed to estimate LID evolution in the past when δ 15 N data are still missing. In this work, we implement new absolute age constraints spanning the last 800 kyr derived from 81 Kr measured in air trapped in EDC ice core as well as new orbital age constraints obtained by synchronizing up-to-date EDC records with their orbital target. We combine these data to recent volcanic matching and methane records synchronization which provide additional stratigraphic links, relating EDC to other ice cores over the past 122 kyr  (Svensson et al. , 2020; Baumgartner et al. , 2014). Finally, we propose the new chronology AICC2023 with reduced chronological uncertainties.
Dating strategy The Paleochrono Python software is a probabilistic dating tool similar to Datice and Icechrono1 with improved mathematical, numerical and programming capacities  (Parrenin et al. , 2021). The dating strategy of Paleochrono relies on the Bayesian inference of three glaciological functions forming the input background scenario: accumulation rate (𝐴), thinning of annual ice layers (𝜏) and Lock-In-Depth (𝐿𝐼𝐷). The three variables evolve along the ice core depth 𝑧 and are used to estimate the ice (𝜓) and gas (𝜒) age profiles as follows: 𝜓(𝑧) = ∫ 𝐷(𝑧 ′ ) 𝜏(𝑧 ′ )𝐴(𝑧 ′ ) 𝑑𝑧 ′ 𝑧 0 (1) 𝜒(𝑧) = 𝜓(𝑧 -Δ𝑑𝑒𝑝𝑡ℎ(𝑧)) (2) ∫ 𝐷(𝑧 ′ ) 𝜏(𝑧 ′ ) 𝑑𝑧 ′ = 𝑧 𝑧-∆𝑑𝑒𝑝𝑡ℎ(𝑧) 𝐿𝐼𝐷(𝑧) × 𝐷 𝜏 | 𝑓𝑖𝑟𝑛 0 ( 3 ) where 𝐷 is the relative density of the snow/ice and in the firn when the air particle was at the lock-in-depth (this parameter is usually ~0.
As a consequence, and in addition to our new extensive δ 15 N dataset, we have chosen to use here the firn model approach of  Bréant et al. (2017). In order to make a correct calculation of uncertainties linked to firn modeling at EDC, we ran two tests of the model with and without including the impurity concentration parameter (Table  3 ). Q depending on temperature  (Bréant et al. , 2017)  Q depending on temperature  (Bréant et al. , 2017 ).
81 Kr age constraints Three ice samples from the bottom part of EDC have been analyzed and provide three age estimates displayed in  V and VI (Fig. 1 ). As these novel δO 2 /N 2 measurements have been performed on ice samples stored at -50°C, there is no storage effect and they can directly be merged with the 800 kyr long record of  Extier et al. (2018c)  (Table  1 ). The new dataset improves the resolution of the long EDC record, reaching sub-millennial scale accuracy over MIS 5, 7, 9 and in particular over MIS 11 and MIS 13, periods of sparsity in the ancient record  (Extier et al. , 2018c). Although the two datasets agree well over recent periods (last 350 kyr), they show some discrepancies during older periods (between 550 and 375 ka BP, see Fig. 1 ). Such dissimilarities are observed over MIS 11 (between 424 and 374 ka BP) where the sampling resolution of the previous dataset is particularly low (2,500 years).
In addition, the MIS 11 is a period characterised by a low eccentricity in the Earth orbit, inducing subdued variations of insolation, causing δO 2 /N 2 changes of smaller magnitude and leading to lower signal to noise ratio. Following a data processing treatment consistent with the method described in  Kawamura et al. (2007) , the compiled dataset is linearly interpolated every 100 years, and then smoothed using a finite-duration impulse response (FIR) filter with a KaiserBessel20 window (cut-off from 16. 7 to 10. 0 kyr period, number of coefficients of 559 for the 800 kyr long record) designed with the software Igor Pro, in order to reject periods inferior to 10,000 years and erase the noise present in the data. Note that using a low-pass (rejecting periods below 15 kyr) or a bandpass filter (keeping periods between 100 and 15 kyr periods, used by  Bazin et al.
Total Air Content The new TAC record is continuous over the last 800 kyr with a mean sampling resolution of 2,000 years (Fig. 3 ). The raw data between 800 and 440 ka BP are not shown here and will be published in a separate study  (Capron et al. , in prep). The TAC series shows a good resemblance with the integrated summer insolation (ISI, obtained by a summation over a year of all daily insolation at 75° South above a chosen threshold). After comparison in the frequency domain of the TAC record of EDC ice core with ISI obtained using different thresholds, the ISI curve calculated for a threshold of 375 Wm -2 (ISI375) exhibits the finest spectral agreement with the TAC record of EDC ice core over the past 800 kyr. The coherency between the TAC record and ISI is deficient over MIS 11 (between 430 and 370 ka BP) and in the deepest part of the core (prior to 700 ka BP) where the signal to noise ratio is low.
Following a data processing treatment consistent with the method described by  Lipenkov et al. (2011) , the 800 kyr long TAC dataset is interpolated every 100 years, and then filtered with a band-pass filter rejecting periods below 15,000 and above 46,000 years (IgorPro FIR filter with a KaiserBessel20 window: cut-off from 15 to 14 kyr period and from 46 to 47 kyr period, number of coefficients of 559). Outliers in the raw TAC dataset are discarded if they show an anomaly superior to 1. 0 mL/kg (standard deviation of TAC record) when compared to the band-pass filtered signal. 45 outliers are rejected out of 399 datapoints (among which 16 outliers are identified between 100 and 0 ka BP). The TAC is interpolated and filtered again after removal of the outliers. Tie points are mostly determined by matching variations extrema of TAC and integrated summer insolation at 75°S (see Fig.
3 ). Indeed, in case of a non-linear relationship between TAC and insolation, extrema are better indicators of TAC response to insolation forcing. Moreover, filtering the dataset induces a bias in the mid-slope position. The method employed to determine extrema position is the same as for δO 2 /N 2 insolation tie points. Only one of the tie points is identified by matching mid-slopes (i. e. derivative extremum) at 362 ka BP rather than minima at 375 ka BP due to the flatness of the insolation minimum which precludes to identify an accurate tie point. All extrema are not tuned to the target due to the poor resemblance between the signal and insolation and 42 unambiguous tie points were kept out of 64 detected by the automated method. The tie point uncertainty finds its origin in the age errors associated with the filtering (~700 years), tie point identification and outlier rejection (~900 years).
𝛅 𝟏𝟖 𝐎 𝐚𝐭𝐦 In this work, new highly resolved δ 18 O atm data on EDC ice core are presented over Terminations II, III, IV, V and VI (Fig. 4 ). The available δ 18 O atm data can be sorted out in two groups: new δ 18 O atm data  (Grisart, 2023)  at high temporal resolution (between 333 and 1,375 years, see Table  1 ) and old measurements compiled by  Extier et al. (2018b) , characterised by a lower sampling resolution (between 1,000 and 1,500 years, see Table  1 ), except between 374 and 153 ka BP (resolution between 160 and 700 years, see Table  1 ). The new dataset allows us to improve the resolution of the long EDC record over MIS 5, 7, 9 and in particular over MIS 11 and 13, periods of  Between 810 and 590 ka BP, we updated the following approach of  Bazin et al. (2013) : EDC δ 18 O atm and 5 kyr delayed climatic precession are synchronized using mid-slopes of their variations.