dataesr/scientific-paragraphs-categorization huggingface dataset:
PlioMIP2 models PlioMIP  (Haywood et al. , 2011)  is a project aimed at studying climate and environments of the late Pliocene. PlioMIP is part of the Palaeoclimate Modelling Intercomparison Project (PMIP). In consideration of the development of models and proxy data in the past years, this project has progressed to the next phase (PlioMIP2)  (Haywood et al. , 2016b (Haywood et al. , , 2020)). In the second phase, there are more models (17 models in comparison to 11 models that participated in PlioMIP phase 1). Most of the newly participating models are descendants of existing models. For example, CCSM4 participated in PlioMIP1. The models derived from it, CCSM4-Utr, CCSM4-UoT, CESM1. 2 and CESM2, participated in the second phase; IPSLCM5A from the Institute Pierre-Simon Laplace (IPSL) participated in PlioMIP, and then joined PlioMIP2 with IPSLCM5A2 and IPSL-CM6A-LR; NorESM-L, which has been developed from CCSM4, participated in the first phase, and then participated in the second phase with NorESM1-F.
From this map it is clear that in the mid-Pliocene, several modern ocean gateways were absent, such as the Bering Strait, the connection between Canada and Greenland (Canadian Arctic Archipelago), and some straits in the MC. Details of the experimental design are described by  Haywood et al. (2016a). Note that HadGEM3 did not change the land sea mask in the Eoi400 experiment, and therefore for this model the various gateways are the same as in E280. Data from the experiments E280 and Eoi400 used in this study are from the PlioMIP2 dataset, part of which is available at https://esgf-node. llnl. gov/search/cmip6/. In this study, the majority of the results are calculated across the final 100 model years of each simulation. Data has been regridded to a 1 ⇥ 1 grid in order to calculate ensemble means; results for individual models are illustrated based on the original grids.
During the calculation of ocean currents, we find that there are some disagreements between results derived from regridded data and analyses performed on individual models' native ocean model grids. Therefore, to reduce errors and improve the accuracy, we use native grid data to calculate the oceanic flow intensity. Fig. 1  shows the gateways (red rectangle), and differences in the land sea distribution between the E280 and the Eoi400 simulations. In the modelled mid-Pliocene, some regions, which are islands in the preindustrial, are still a peninsula that is connected to the Asian mainland in some models (CESM2, EC-Earth3. 3, CESM1. 2, MIROC4m, HadCM3, GISS-E2-1-G, CCSM4, CCSM4-Utr and CCSM-UoT). In these models except for HadGEM3 which adopted the same land-sea mask of the MC in both experiments, the Malacca Strait, and some other small straits in the northern part of the MC, are closed as shown in Fig.
Acknowledgements. This research is funded by the  European Union 's  Horizon 2020 research and innovation programme  under the  Marie Skłodowska-Curie  grant agreement No  813360 4DREEF  and the  NERC SWEET  grant (grant no. NE/P01903X/1 ). Wing-Le Chan  and  Ayako Abe-Ouchi  acknowledge funding from  JSPS  ( KAKENHI grant  no. 17H06104  and  MEXT KAKENHI  grant no. 17H06323 ) and computational resources from the Earth Simulator at  JAMSTEC, Yokohama, Japan. Development of GISS-E2. 1 was supported by the  NASA Modeling, Analysis , and  Prediction (MAP) Program. CMIP6 simulations with GISS-E2. 1 were made possible by the  NASA   High-End Computing (HEC) Program  through the  NASA Center for Climate Simulation (NCCS)  at  Goddard Space Flight Center. Charles J. R. Williams  acknowledges the financial support of the  UK Natural Environment Research Council (NERC) -funded  SWEET  project (research grant no.
04 Gt a -1 (see Supplemental Dataset 1), although future glacier changes (e. g. , due to surging  21  ; tidewater glacier cycle  22  ) will modify this list. To investigate which parts of the ocean are most affected by frontal ablation, and thus susceptible to related iceberg hazards and marine ecosystem impacts, we introduce an 'intensity index' that spatially aggregates frontal ablation from all glaciers occurring within 50 km of any point in the ocean on a 10 km grid (Fig. 3 ). The index depends on both the frontal ablation rates and the number of glaciers in the vicinity of each point, and thus varies in a pattern that is largely consistent with the spatial variation in frontal ablation of glaciers in each ocean region (Fig. 2 ). Several coastal regions show exceptionally high-intensity indices indicating hotspots of potential iceberg occurrence.
e. , no flow bands). For Svalbard and the Russian Arctic, where the optical velocity record from 2000 to 2010 was sparse in ITS_LIVE, we filled this gap using 1992-2012 observations from the JERS1, ERS1, ALOS PALSAR, and TerraSAR-X satellites. We used synthetic-aperture radar (SAR) offset-tracking to compute winter velocity from these scenes at 100 m resolution  43. We combined all available velocity observations over each decade to derive a mean decadal velocity along the flux gate for each glacier. We do not account for seasonal variability in ice motion due to a lack of data to quantify this variability. In the few instances for which no velocity data was available for the 2000-2010 period (3. 9% of flux gates), mean decadal values from 2010 to 2020 were used to fill the gaps.
Ice thickness Direct radar-derived measurements, primarily from the Glacier Thickness Database (GlaThiDa) 3. 0. 3  44 , were used to determine ice thickness along the flux gates and over deglaciated areas due to terminus retreat since 2000. For Arctic Canada North, Center for Remote Sensing of Ice Sheets 45 radar tomography observations acquired in ~3 km swaths from the 2014 NASA Operation IceBridge mission provided the main ice thickness source. Additional radio-echo sounding observations were incorporated for the Russian Arctic, where available  [46] [47] [48] [49] [50] [51] [52]. Whereas most of our observations were acquired during the study period, we used 1288 Svalbard thickness measurements in GlaThiDa from 1980 to 1995 and some observations in Russia that were collected in 1994  46  and 1997  47.
Glacier terminus front positions To determine 4S term (Eq. 3) for the periods 2000-2010 and 2010-2020, glacier terminus positions were digitized manually for every glacier in approximately 2000, 2010, and 2020, primarily from summer, cloud-free, true-color Landsat 5, 7, and 8 imagery (30 m spatial resolution). When Landsat images were not available, we first attempted to fill the gap with true-color ASTER data (15 m spatial resolution), particularly in northern Greenland and Arctic Canada North for 2000 and 2010. Where no cloud-free optical observations were available, we used Radarsat-1 fine beam mode (2000; 8 m pixel resolution; used in Northern Greenland) and ALOS PALSAR fine beam single polarization (2010; 10 m pixel resolution; used in Arctic Canada North) imagery downloaded from the Alaska Satellite Facility.
Climatic mass balance below the flux gate To account for mass changes unrelated to frontal ablation (e. g. , surface melting) between the flux gate and glacier front), we used the Python Glacier Evolution Model (PyGEM)  56  to model the climatic mass balance _ B clim for each glacier. _ B clim is defined as the mass changes due to snow accumulation, and melt at the surface and refreezing of melt and rain water below the surface  1. PyGEM assumes the basal mass balance and internal ablation is negligible. We calibrated PyGEM using geodetic mass balance estimates  9. We used monthly ERA5 air temperature and precipitation data  57  to compute glacier melt with a degree-day model, accumulation with a temperature threshold, and refreezing from annual air temperature. We extracted _ B clim for the lowest two 10-m elevation bins of each glacier.
For the remaining glaciers, we used ice discharge values neglecting _ B clim , but with 100% uncertainty. These 488 and 601 glaciers only made up 0. 70 ± 0. 44 Gt a -1 (2000-2010) and 0. 54 ± 0. 46 Gt a -1 (2010-2020) of ice discharge, accounting for 2% of total ice discharge during both time periods. See Dataset S1 for a list of glaciers that have this correction.
Similarly, we estimate the uncertainty in terminus mass change σ _ M term assuming each uncertainty component (Eq. 3) is independent, which yields: σ 2 _ M term = ðσ ΔS term HÞ 2 + ðσ H ΔS term Þ 2 + σ _ B clim ΔS term 2 2 + _ B clim σ ΔS term 2 2 ð6Þ where σ ΔS term is the uncertainty in the terminus area change and σ H is the uncertainty in the mean thickness over 4S term. Finally, the uncertainty in frontal ablation σ _ A is estimated by considering the ice discharge and terminus mass change uncertainties (Eq. 1) as independent: σ 2 _ A = σ 2 _ D + σ 2 _ M ð7Þ Uncertainty sources at flux gate point or glacier-scale We relied on the uncertainties reported by each dataset, whether at the scale of a flux gate point, flux gate, or glacier to propagate the uncertainties through Eqs. 5 and 6.
Acknowledgements The start of this project and early ideas that informed this work were a result of the "Workshop on the  Importance of Calving  for the  Mass Balance of Arctic Glaciers " held in Sopot, Poland in 2016, sponsored by the  Centre for Polar Studies, University of Silesia  and the  IASC Network on Arctic Glaciology. We thank  Ruitang Yang  for her help in improving the error analysis. W. K. acknowledges support from the  Vanier Graduate Scholarship. L. C. thanks the  Natural Sciences and Engineering Research Council of Canada, University of Ottawa and ArcticNet Network of Centres of Excellence Canada  for funding. T. S. acknowledges support from the  ESA Glaciers CCI  project. A. G. and I. L. thank  RG State  Contract  FMGE-2019-0004 , and H. J. the  Natural Sciences and Engineering Research Council of Canada.
Hydrogravimetric Estimate Although TWI provides clues on the location of potential preferential water storage areas, it does not allow to model WSC for comparison with observed gravity changes. For this purpose, we use NIHM which explicitly simulates surface and subsurface flow occurring in spatially heterogeneous subsurface compartments, providing a quantitative depiction of 3D water storage dynamics within the Strengbach catchment that we compare with observed gravity dynamics. We developed a Matlab R code to compute gravity temporal changes at gravity stations from the water content field θ NIHM x, y, z, t simulated by NIHM. The modeling code follows a 3-step workflow (Figure  9 ) inspired from  Pearson-Grant et al. (2018)  that consists in 3 functions: (i) cellcorners. m that discretizes the volume extending from the surface to the hydrological model substratum into a 3D gravity grid with prismatic cells based on a digital elevation model of the topography, (ii) density.
Gravity observations have been discussed by means of the topographic wetness index and by means of a physically based distributed hydrological model. This study demonstrates the ability of hybrid gravimetry to assess the water storage dynamics in the Strengbach mountain catchment, which contributes to a better understanding of the hydrological processes occurring in mountain hydro-systems in general. Last but not least, our gravity dataset provides observations not presumed by the applied physically based distributed hydrological model, indicating that this model is currently underconstrained, which likely result from the relatively large number of parameters required to capture the spatial heterogeneity of the catchment subsurface. Hybrid gravimetry therefore appears as a very promising measurement method to condition such distributed hydrological model.
ACKNOWLEDGMENTS We are very grateful for the considerable work done by the editor and three reviewers which greatly improved the quality of this paper. We thank OHGE (http://ohge. unistra. fr/)for providing all hydro-meteorological data used in this study. We thank  Julie Croisette  and  Célia Grün  for having performed the microgravity surveys during Summer 2019. We also thank  Sylvain Pasquet ,  Guillaume Modeste ,  Romain Pestourie ,  Matthias Oursin ,  Damian Kula , and  Sophie Chaffaut  for their generous support during the extensive fieldwork. Last but not least, we thank  Marianne Oudin  for her warm welcome at Les Grands Prés guest houses in Aubure.
FUNDING The installation of the iGrav30 superconducting gravimeter has been funded by EQUIPEX CRITEX (https://www. critex. fr). The doctoral school ED413 from the  Strasbourg university, as well as   ANR HYDROCRIZSTO ANR-15-CE01-0010-02  and the OZCAR network (http://ozcar-ri. prod. lamp. cnrs. fr/)provided the funding for this study.
g. ,  15, 28, 29  ). Independent shear-wavespeed models have been also retrieved from ambient noise data, and they are able to constrain V s with a uniform coverage, but techniques, dataset, and resolution are different (e. g. , at local/regional scale  [30] [31] [32] [33] [34] [35]  ). Full waveform adjoint tomography models are presently available for the European continent, offering large-scale hints on the Italian lithospheric structure inferred from wavespeed values in the shallow layers of their models  9, 11, 12, 36. In the last two decades, a large number of high-quality, threecomponent seismic instruments have been deployed in Italy  37  , and several moderate seismic sequences have been recorded. These data represent an invaluable source of information to construct, for the first time, a comprehensive 3D wavespeed model of the Italian lithosphere based on full waveform adjoint tomography ( 38 , and references therein).
Recent studies demonstrate the success of the adjoint-state method  39, 40  in imaging complex structures in many regions in the world (e. g. ,  11, 16, 36, [41] [42] [43] [44] [45] [46] [47]  ) and on a global scale  48, 49. The main advantages of the method are a) comparing entire seismic traces simulated for realistic 3D Earth models to corresponding recorded waveforms and b) performing inversions constrained by all the phases that appear in these seismograms. Here we present a new 3D tomographic model obtained by applying the method to the entire Italy. We start from a recent 3D traveltime tomographic model for the peninsula  29  and increase its resolution by taking advantage of a full waveform extended dataset (Fig. 1 ). We perform high-frequency wavefield simulations using a spectral-element method (SEM), implemented by the code SPECFEM3D_Cartesian  50  , and we exploit the SEM accuracy within the adjoint-state method to perform a 3D tomographic iterative inversion using full waveforms (Section "Forward and inverse modeling").
Some features of the seismograms are still not matched by our final model Im25, suggesting that future extensions of this work will need to add more data, perform some more iterations, and/or invert for additional model parameters, such as transverse anisotropy or attenuation. Nevertheless, the present model is able to reproduce many of the main features of the observed seismograms, which supports the robustness of the imaged structures. An independent test of the improved quality of fit for Im25 has been carried out by  51  , using the 2016, M w 6. 0 Amatrice earthquake, which was not included in our dataset. Tomographic model. In order to provide a general overview of model Im25, we analyze it in the same seven cross-sections chosen by  29  and one section through Mt. Etna volcano (Figs.
4, 5, Supplementary Figs. 6-10). This could also explain the uneven distribution of the Adriatic seismicity and the differences in the observed stress regimes. An evidence of the robustness and reliability of our Im25 model results from the agreement in the comparisons of V s profiles extracted at selected stations from our model and studies using different techniques, dataset, and regions (Supplementary Fig. 11 ). Horizontal cross-sections provide complementary views of our tomographic model (Supplementary Figs. 3g, h ). Variations in V p and V s are generally negative (i. e. , the new wavespeed is lower than the initial one) along the continent and positive along marine areas. In both sub-figures-but especially for V s -the footprint of the large sedimentary basin of the Po Valley is highlighted by a wide zone of lower wavespeeds in northern Italy.
Methods Dataset. We selected a total of 163 earthquakes that occurred in the chosen simulation volume between 2005 and 2014 and with a magnitude range 3. 0-5. 8 (red and blue dots in the inset of Fig. 1c ). The initial source mechanisms are obtained from Time Domain Moment Tensor (TDMT) solutions that are routinely calculated by INGV 72 by inverting full, three-component traces of regional, broadband stations using 1D wavespeed models (http://terremoti. ingv. it/). The selection of the events and collection of the corresponding source parameters were done using the e-Science seismological gateway, i. e. , the VERCE portal (European projects VERCE, www. verce. eu, and EPOS-IP, https://www. epos-ip. org). The computational cost of an adjoint tomographic inversion scales with the number of events, i.
For each event the corresponding recorded seismograms were downloaded from the INGV EIDA/Orfeus archive, and a preliminary data quality check was performed by discarding traces with gaps or low (<3) signal-to-noise ratio. We used only velocity traces recorded by permanent stations (see EIDA archive  73  https:// eida. ingv. it). Data recorded by temporary networks were neglected because these are short period (high frequency) instruments; these data could be useful for successive studies at higher resolution. Initial model and parameterization. The geographical volume considered for the tomographic inversion covers the entire Italian peninsula and neighboring areas in France to the west, in Switzerland and Austria to the north, and along the eastern Adriatic coast and Greece to the east (Fig.
1c ; see Section "Forward and inverse modeling"). We used the technique presented by  76  and implemented in their code cmt3d_FLEXWIN. The procedure requires an initial source solution for each event, which, in our case, is the TDMT solution based on a 1D wavespeed model and published in the INGV database (Section "Dataset"). SPECFEM3D_Cartesian (Section "Forward and inverse modeling") is used to simulate the synthetics for these source models and the starting 3D wavespeed structure, and to numerically calculate the Fréchet derivatives of a given waveform misfit function with respect to the source parameters. This misfit between data and synthetics is then minimized to invert for the six moment tensor components and the event location (a total of 9 source parameters). Supplementary Fig. 13a  shows the moment tensor solutions for each of the 123 events: the black beachball is the initial TDMT solution, the blue beachball is the solution of the source inversion performed using the starting 3D wavespeed structure (Supplementary Data 1).
This second source inversion allows us to refine the source descriptions based on a wavespeed model that was iteratively improved through the first 12 tomographic iterations (orange histograms in Supplementary Figs. 13b, c ). Forward and inverse modeling. Wavefield simulations for both forward and inverse modeling in our study were performed using the software package SPEC-FEM3D_Cartesian ( 50 , https://github. com/geodynamics/specfem3d), which is based on the spectral-element method (SEM;  77, 78  ). The SEM combines exponential spatial accuracy due to a high-degree polynomial representation of the wavefield, computational efficiency guaranteed by an explicit time marching scheme easily implemented on parallel computers (e. g. ,  50, 79  ), and geometrical flexibility. The solver SPECFEM3D_Cartesian is specifically designed for local and regional simulations and allows us to accommodate all complexities that affect seismic wave propagation, such as topography, lateral wavespeed variations, attenuation, anisotropy, absorbing conditions at model boundaries.
The wavefield is saved as synthetic seismograms at the chosen stations. For the first 12 iterations we used 123 events (red dots in the inset of Fig. 1c ) out of the 163 discussed in Section "Dataset". We added the remaining 40 events (blue dots in the inset of Fig. 1c ) to the dataset at iteration 13. Most of these events have lower magnitudes and can be useful to constrain higherfrequency features of the model. They also improve the spatial coverage for the adjoint tomography. For iterations 19 to 25 we removed 34 events located in the eastern part of the Adriatic Sea. This allows us to focus on resolving the volume under the Italian peninsula by reducing possible artifacts due to the inclusion of waves traveling through the Adriatic Sea and of earthquakes with large azimuthal gaps. Removing events also reduces the computational costs for the last iterations.
(ii) The synthetic seismograms and corresponding observed seismograms are processed in the same way, by removing the trend and mean, tapering the waveforms, deconvolving the instrument response from the data, and bandpass filtering. We initially considered the period range 10-100 s also based on the minimum resolved period of the mesh; then at iteration 16 we shifted to periods of 10-40 s to emphasize body-wave sensitivity over surface waves. (iii) Observed and simulated seismograms are compared using the code FLEXWIN  80  , which, based on specific tuning parameters, automatically selects time windows on waveforms where data and synthetics have quantifiable agreement and good signal-to-noise ratio. These windows are suitable for misfit measurements and are the only windows considered in the next steps of the tomographic inversion.
(iv) Misfit measurements between data and synthetics are calculated within each window using the code measure_adj 41 , which also calculates the corresponding adjoint sources based on the seismogram measurements. To quantify the discrepancy between data and synthetics, we performed multi-taper traveltime measurements: Δτ i ðω; mÞ ¼ τ obs i ðωÞ À τ i ðω; mÞ ð1Þ where m is the given wavespeed model and τ obs i and τ i are the frequencydependent traveltimes for observed and synthetic waveforms in the i th time window. (v) All the adjoint sources for each earthquake are considered as simultaneous sources located at the receivers. These adjoint sources generate the adjoint field of that event simulated using SPECFEM3D_Cartesian and the wavespeed model of iteration n. From the interaction between the forward and adjoint wavefields, the code constructs finite-frequency sensitivity kernels for each event and for the model parameters V p and V s , namely the isotropic event kernels K α and K β.
Acknowledgements This work was made possible thanks to the visionary computational geophysicist  Dr. Dimitri Komatitsch  (1970Komatitsch ( -2019)). Computational resources for the project IMAGINE_IT (PI  Dr. Dimitri Komatitsch ) have been provided by the  PRACE  consortium (40 million CPU-hours on the Tier-0 GENCI/TGCC CURIE supercomputer -2014 9 th PRACE call). Federica Magnoni  has been supported by the European  FP7-INFRASTRUCTURES project 'VERCE'  (n. 283543), the  Italian  project 'SHAKEnetworks' ( PRIN-Research Project of Relevant National Interest-MIUR ,  2012KMAEZF ), and the  H2020 projects 'EPOS-IP'  (n. 676564) and 'DARE' (n. 777413). We are grateful to  Giovanni Chiodini  and  Gianfilippo De Astis  for useful sharing of ideas, figures, and data. We thank  Ebru Bozdağ ,  Hejun Zhu , and  Daniel Peter  for fruitful discussions.
Data availability All the continuous seismic data are collected from the European Integrated Data Archives (EIDA; https://eida. ingv. it). The source parameters of the earthquakes considered in this study, for the initial TDMT solution and the solutions obtained with the starting 3D wavespeed model and model Im12, are available in the following public repository: https://doi. org/10. 5281/ zenodo. 5997587. A digital file of model Im25 in the study region of this manuscript is available in the following public repository: https://doi. org/10. 5281/zenodo. 6006660.