dataesr/scientific-paragraphs-categorization huggingface dataset:
Such discrepancies have already been reported in several studies (e. g. Carmagnola et al. , 2013)  and attributed to uncertainties in the values of the ice complex refractive index from  Warren and Brandt (2008). The use of alternative databases, especially the one from Grundy and Schmitt (1998) noticeably improves the agreement in the range 1400-1500 nm as shown by Fig. 10. A more extensive comparison of the simulated and measured reflectance angular distribution at 30 • illumination angle shows that the BRF predicted by the model is in good agreement with the measured one. Spectral changes in the angular distribution are well reproduced. The impact of the CLD shape on the simulation is only detectable for absorptive wavelengths (1500 nm) and high viewing zenith angles. The anisotropy of the reflectance, quantified by the relative difference between the reflectances measured at a 70 • zenith angle in the forward and backward directions in the principal plane, is systematically overestimated in the simulations.
(  5 ). The last but most probable reason is the "skin effect". For small optical paths, the light penetration depth is really small and probably less than a millimetre (e. g. Mary et al. , 2013) , so the role of the skin layer, where the grains could be finer and have a preferable facet orientation, is crucial. This is also corroborated by the fact that S2 and S3 have only a small difference in SSA, measured in deeper layers, but a strongly different behaviour in reflectance. https://doi. org/10. 5194/tc-15-3921-2021 The Cryosphere, 15, 3921-3948, 2021   (1998)  and  Kou et al. (1993). For all the other simulations we used refractive index from  Warren and Brandt (2008). The effect of crystal shape and SSA on snow BRF is intricate (e. g. Dumont et al. , 2010; Stanton et al. , 2016). The dataset presented in this study is limited to only three snow samples, so no statistically robust conclusions on the effect of crystal shape on reflectance can be drawn.
The measurements, however, show that the faceted crystals exhibit more anisotropic reflectance than fragmented particles/rounded grains as observed in  Stanton et al. (2016)  and that the anisotropic behaviour is further reinforced if facets are orientated upward. More quantitative relationships between crystal shape and the BRF would require a much larger number of snow samples analysed. One way of addressing this problem is to take advantage of ray tracing models such as those developed by  Kaempfer et al. (2007 ), Picard et al. (2009 ), and Petrasch et al. (2007) , which can be run directly on 3D images of the microstructure. To conclude, this unique dataset combining X-ray tomography imaging of snow microstructure and high-accuracy measurements of snow BRF was used to demonstrate the following.
Code and data availability. The BRF dataset is available in the PANGAEA database along with the CLDs and the simulation results at https://doi. pangaea. de/10. 1594/PANGAEA. 935028 (Dumont et al. , 2021). All simulations are based on the equations from  Malinka (2014)  and Malinka et al. (2016). Other numerical data used in the study are available directly in the tables, in the main text or in the appendices.
Methods This manuscript is based on a full-plate tectonic reconstruction of the last billion years that has been developed on the open access software GPlates (www. gplates. org). All GPlates files needed to reconstruct the model are available here-https:// zenodo. org/record/3647901. The methodology for constructing full-plate tectonic reconstructions are detailed in Merdith et al. 10  and also in Domeier and Torsvik  56. Computation of latitudinal surface area was done using pyGplates (www. pygplates. org). To extract the latitudinal distribution of continent area, we first created a global equal-area mesh. We then used a grid-intersection between the nodes of the mesh and the polygons of the plate model that represent continental crust to estimate the area of crust at each latitude. These were summed for each timestep to create an array from 1000-520 Ma of the latitudinal distribution of crust on the earth.
Acknowledgements A. S. C. is supported by  Australian Research Council  grants  FT120100340  and with M. L. B. and J. D. F. are funded by  LP160101353 , with the support of the  Northern Territory Geological Survey, Origin Energy ,  Santos Ltd  and  Imperial Oil and Gas. A. S. M. is supported by the  Deep Energy Community  of the  Deep Carbon Observatory. We are grateful for the help of two reviewers who improved the manuscript.
We used altimetry measurements from the ENVISAT satellite, covering the 2002-2010 period. As each 0. 5 • TRIP cell is not fully covered by water, h(t) in Equation (  2 ) is an equiv-alent water height which cannot be directly compared to the radar altimetry measurements. We need to determine the variable surface water extent of each TRIP cell corresponding to the altimetric datasets. We adopt the methodology developed in  [54]  to derive lake water extent from 16-day MODIS/NDVI imagery  [53]. We define water when NDVI values are lower than 0. 1 and NIR (near infrared) reflectance is lower than 0. 2. The measured equivalent water height derived from radar altimetry and MODIS measurements over each 0. 5 • cell is defined as: h Measured (t) = 1 S tot h ALTI (t) h 0 S MODIS (h)dh (3) where h ALTI (t), h 0 , S MODIS and S tot are respectively the altimetric water height mea- surements, a reference altitude, the surface water extent over the TRIP cell and its total surface.
Validation of the River Models We validated our river models by comparing them to independent datasets: (1) mean annual river discharge gauges from GRDC, and (2) total mass variations over each individual basin with GRACE estimates. Figure  2  shows the comparison between the mean annual modeled discharge for a selection of six gauges in the Amazon, Orinoco and Parana River basins and gauge measurements from the GRDC. Although our river modeling is rather simple, there is a significant good agreement between our optimal river model and the discharge data, both in terms of amplitude and phase. A perfect match between the modeled and the measured discharge cannot be expected, as real discharge can vary significantly along the river whereas our model has a low resolution (0. 5 • ) and a constant flow velocity.
5 years, the gap ratio could not be more than one third of the time series. The observation period spans from January 2003 to February 2016 and 60% of GPS time series were longer than 6 years (148 sites). The formal error of the selected stations was lower than 1 mm and 3 mm for the horizontal and vertical components, respectively. We corrected the time series offset discontinuities related to equipment changes or earthquakes using the Hector software  [60]  with the potential steps epoch database provided by the NGL. Then, we adjusted and removed a linear trend of the GPS time series and outliers on all the three components at each station. We defined outliers as data points that exceed two times the sigma value. Before signal analysis, we filled the gap in GPS time series using an autoregressive (AR) model.
We can clearly notice the major improvement for the vertical component in the Amazon basin. We compute the number of sites in each case for different classes of correlation coefficient values: lower than 70%, between 70% and 80%, between 80% and 90%, and greater than 90%. Table  1  gives the distribution of the corresponding values for each component and each dataset. For the whole station set, the correlation coefficient values range from 0% to 100% for all components. The improvement between the initial (labelled nR in the table) and the new models (referred as wR in the table) can be highlighted for instance by the reduction of the percentage of sites with a correlation coefficient lower than 70%. Regarding the values in Table  1 , the improvement is not obvious considering this value over the entire network.
Data Availability Statement: GPS timeseries were obtained at NGL (http://geodesy. unr. edu/ accessed on 1 March 2021). Hector software is available at (http://segal. ubi. pt/hector/ accessed on 1 March 2021). Maps were obtained using GMT software (https://www. generic-mapping-tools. org/ accessed on 1 March 2021). We thank Hydroweb (data available at http://hydroweb. theia-land. fr/ accessed on 1 March 2021) research groups for providing water level data used in this study. Most of the loading estimates are available at the EOST loading service (http://loading. u-strasbg. fr accessed on 1 March 2021).
Two examples are presented: The first station pair, II. PFO-AZ. FRD represents an ideal scenario in which the train sources cover part of the stationary phase zone; for the second station pair, II. PFO-AZ. CRY this requirement is relaxed (see Fig. 6 ). For data selection and processing we follow a strategy similar to Pinzon-Rincon (2021). We use station CI. IDO located near the railway to time trains as they are passing by the study area  (Brenguier et al. 2019). A 12-min-long time window centred around each detection time is selected. We take the seismograms recorded at the target stations (vertical component only) and perform a cross-coherence measurement in the selected time windows. We perform the analysis over 10 yr of seismic recordings from 2010 to 2020 and average all the correlation functions obtained.
2018b). However, the interpretation of correlation observations in this case requires a different mindset and correlation modelling is a promising approach. Modelling body-wave interactions in correlation functions is a relatively expensive and tedious endeavour. If one were only interested in correlation waveforms, database approaches would be a good alternative  (Ermert et al. 2020 ) and the required modifications can be included easily. Our current framework provides insights into the reconstruction of body waves in correlation functions and is employed to study specific observations in the context of monitoring the San Jacinto Fault. In general, combining observed data and simulations has several potential benefits, for example facilitating/guiding interpretations, testing of different hypotheses that might explain the observations and initiating/accelerating algorithmic developments.
A C K N O W L E D G E M E N T S The authors would like to thank the Editor  Lapo Boschi ,  Jacopo Boaga  and one anonymous reviewer for their time and excellent suggestions that helped to improve the manuscript. In addition, we are grateful to  Yehuda Ben-Zion  and  Frank Vernon  for useful discussions. KS acknowledges support from the  Swiss National Science Foundation (SNSF)  under grant number  P2EZP2 184379.
DATA AVA I L A B I L I T Y Data used in this manuscript are publicly available at IRIS. The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms, related metadata, and/or derived products used in this study. IRIS Data Services are funded through the Facilities for the Advancement of Geoscience and EarthScope (SAGE) Proposal of the National Science Foundation under Cooperative Agreement EAR-1261681. The following networks were used in this study: (i) AZ: UC San Diego. (1982). ANZA Regional Network [Data set]. International Federation of Digital Seismograph Networks. https://doi. org/10. 7914/SN/AZ (ii) CI: California Institute of Technology and United States Geological Survey Pasadena. (1926). Southern California Seismic Network [Data set].
Space techniques based on GPS and SAR interferometry allow measuring millimetric ground deformations. Achieving such accuracy means removing atmospheric anomalies that frequently affect volcanic areas by modeling the tropospheric delays. Due to the prominent orography and the high spatial and temporal variability of weather conditions, the active volcano Mt. Etna (Italy) is particularly suitable to carry out research aimed at estimating and filtering atmospheric effects on GPS and DInSAR ground deformation measurements. The aim of this work is to improve the accuracy of the ground deformation measurements by modeling the tropospheric delays at Mt. Etna volcano. To this end, data from the monitoring network of 29 GPS permanent stations and MODIS multispectral satellite data series are used to reproduce the tropospheric delays affecting interferograms.
A tomography algorithm has been developed to reproduce the wet refractivity field over Mt. Etna in 3D, starting from the slant tropospheric delays calculated by GPS in all the stations of the network. The developed algorithm has been tested on a synthetic atmospheric anomaly. The test confirms the capability of the software to faithfully reconstruct the simulated anomaly. With the aim of applying this algorithm to real cases, we introduce the water vapor content measured by the MODIS instrument on board Terra and Aqua satellites. The use of such data, although limited by cloud cover, provides a two-fold benefit: it improves the tomographic resolution and adds feedback for the GPS wet delay measurements. A cross-comparison between GPS and MODIS water vapor measurements for the first time shows a fair agreement between those indirect measurements on an entire year of data (2015).
, 1998; Webley et al. , 2004) , 2) use of the relation between InSAR phase and elevation  (Beauducel et al. , 2000) , 3) use of the tropospheric delays derived from GPS analysis  (Wadge et al. , 2002; Aranzulla and Puglisi, 2015). In this work, for the first time, an entire year of GPS tropospheric delays has been combined with integrated water vapor contents measured with MODIS instruments on board Terra and Aqua satellites in order to calculate the tropospheric spatial and temporal variability delay. The use of MODIS data in the tomographic modeling of the troposphere has proved to be fruitful by  Benevides et al. , 2015a. However, in  this study, we expand this approach in very harsh experimental conditions dominated by the prominent topography of Mt. Etna. The Earth Observation (EO) data are the only ones that can be systematically used for water vapor content measurements over Mt.
The Global Positioning System Measurements Since 2000, a continuous GPS array has operated to monitor the ground deformations of Mt. Etna volcano  (Palano et al. , 2010). The actual array consists of 42 stations that provide a dense coverage of the volcano (Figure  2 ). Among the whole GPS array, we used a subset of 29 stations where a full data archive is available for the studied period  (2015). The stations are equipped with low-multipath choke-ring antennas, and we apply elevation-dependent corrections for the antenna phase centers  (Schmid et al. , 2005). We processed the GPS data using the GAMIT package developed by the Massachusetts Institute of Technology  (Herring et al. , 2010). This software uses double-differenced GPS phase observations to estimate for each observing span (in our case 24 h) a single set of station coordinates and orbital parameters together with piecewise linear models of zenith tropospheric delay (ZTD) and gradients at each station.
Standard models for precession, nutation, Earth rotation, and solid Earth and ocean tides are applied (IERS Conventions). The motions of the GPS satellites are taken from the Final Orbits of the International GNSS Service (https://www. igs. org/), which typically have an accuracy of ∼2 cm. We processed our data with the data from six surrounding IGS stations (Figure  2 ), used as reference stations to tie our network to the International Terrestrial Reference Frame (ITRF). Although GAMIT's output separately provides dry and wet tropospheric delay, caused by the high spatial and temporal variability, only the wet component is important for our purposes (Figure  1A ). Zenith wet delay (ZWD) is estimated during the GPS processing by assuming a 2-h interval and interpolated through time by a spline approach.
Moderate Resolution Imaging Spectroradiometer Data and Integration The MODIS (Moderate Resolution Imaging Spectroradiometer) instruments operate on board the Aqua and Terra satellites. The sensors collect data in 36 spectral bands within the range 0. 4-14. 4 µm and at different spatial resolutions (2 bands at 250 m, 5 bands at 500 m and 29 bands at 1 km). The satellites are operated by NASA (https://modis. gsfc. nasa. gov/) and cover the entire Earth every 1-2 days; in particular, the Mt. Etna area is covered two times per day, so these data are particularly suitable for the aim of this study. NASA web services disseminate MODIS data at different levels of processing as well as derived products for land, ocean, and atmospheric applications. The MODIS has the appropriate absorption bands to retrieve water vapor content in the atmospheric column  (Spinetti 2004 and Buongiorno, 2007).
Here, we use the IPWV (integrated precipitable water vapor) product of the multispectral EO data series that is stored and delivered at 5 km pixel resolution. Figure  3  shows one of MODIS IPWV related to a night passage of the Aqua satellite for the day August 22, 2015. In order to add these IPWV as observable in our tomographic model of Mt. Etna (Figure  1A ), we must first verify that they are consistent with the IPWV values estimated from the processing of our GPS data. The ZWD is directly proportional to Precipitable Water Vapor by a factor calculated during GPS processing  (Bevis et al. , 1994). The comparison has been performed for the entire year 2015 considering GPS and MODIS-retrieved water vapor data at the same time as satellites pass. In particular, 612 Terra and Aqua satellite acquisitions (both diurnal and nocturnal) were employed in the comparison, and 33 stations of the overall network were considered for GPS measurements.
The comparison was made by averaging the measurements within a considered common area (box in Figure  3 ). The averages result from all the pixels falling in the box for MODIS data and from all acquisitions by the GPS stations within the same box. Figure  4  shows the comparison results. MODIS data overestimate the GPS ones by ∼7%, and the resulting Pearson correlation is equal to 0. 91. Considering these differences, it can be stated that the IPWV measurements derived by MODIS data are usable as input for the tomographic model. Finally, it is worth to note that only the clear sky condition enables obtaining reliable IPWV measurements. We used the daily IPWV data from the MODIS as input for model of Eq. 6 together with SWD from GPS measurements gathered every epoch throughout the day. The SWD was retrieved every 2 h, and we considered the values interpolated by spline.
Sentinel-1A DInSAR Image We performed a DInSAR analysis of C-band Sentinel-1-A data referring to selected images acquired between 18 May and November 26, 2015. These data are acquired from the Sentinel-1A satellite operating under the Copernicus program of the European Space Agency (ESA) available from https://scihub. copernicus. eu/dhus/. The processed images were acquired in TopSAR (Terrain Observation with Progressive Scans SAR) Interferometric Wide (IW) mode (VV polarization; vertical transmit and receive polarization), along the descending orbit. The Sentinel-1A data were processed by GAMMA software  (Wegmüller et al. , 2015) , using a spectral diversity method, and a procedure able to co-register the image pairs with extremely high precision (<0. 01 pixel). The interferograms were produced by applying a two-pass DInSAR processing and multi-looking pixels (5 × 1 in range and in azimuth) in order to maintain the full ground resolution (11 × 13 m).
S. National Weather Service (NWS), we set the tomography volume and computed the A matrix of Eq. 6. Table  1  shows the selected test cases together with GPS at the time of Sentinel-1A passes, the IPWV MODIS data, and the available GPS receiver data. The tomographies were computed considering only the nocturnal acquisitions of the MODIS, as the closest in time to the Sentinel-1A acquisition time. The cloudy sky of 21 October and 26 November compromised the possibility to measure the water vapor content from MODIS data, excluding cloudy areas as observables. This is why the damping factor c of Eq. 7 is different for those two cases (see Table  1 ). As previously mentioned, the SF and c thresholds govern the selection of the representative voxels. The last column of Table  1  represents the percentage of the resolved voxels where the 21 October and 26 November cases show the effect of the lack of MODIS data clearly.
Besides the differences in the spatial resolution between the Sentinel-1A data and the tomography, due to the intrinsic differences in the pixel sizes in the two datasets, overall, the DInDelay shows the same fringes pattern of the SAR interferograms. It confirms that at the first order, the selected SAR interferograms contain significant atmospheric contributions as we supposed. The corrected interferograms (third row of Figure  10 ) show two notable characteristics. 1) The pattern of the fringes often mimics the boundary of the voxels ("boxiness"); this is particularly evident in the south-eastern (along the coast line), western and north-western flanks of the volcano. This effect, which prevents the analysis of subtle "small" features such as local faults, is somewhat expected due to the large dimension of the voxels, and might be overcome in the future by reducing the voxel dimensions.
Indeed, the same voxel-shaped anomaly is observed in the corrected interferograms of 18 May-17 July and 18 May-26 November. The latter shows a residual in the western flank, indicating an approach of the ground surface to the radar sensor which cannot be associated either to the inflation or to the spreading of the volcano, both phenomena often observed in this flank  (Bonforte et al. , 2008). Indeed, by independently processing the GPS time series for the same period, with GIPSY 6. 4 software (https://gipsy-oasis. jpl. nasa. gov/), we cannot observe any significant deformation consistent with these residuals (Figure  11 ). The maximum deformation in LOS among all the stations is in fact less than 1 cm and is recorded at the ECHR station. Thus, the residuals might be related to turbulent atmosphere that cannot be modeled by the tomography.
In the lower histogram, the relative difference between "clear sky" (GPS + MODIS) and "cloudy" (GPS) vs. "cloudy" (GPS) tomographies, for the different layers (1,000 m, 3,000 m, 5,000 m, 7,000 m, and 9,000 m) and the whole set of the tomographies. refractivity tomography results in about 37% of resolved voxels. Adding MODIS IPWV data increases on average the number of resolved voxels by about 8%. In the lower critical layer (1,000 m), the resolved voxels improve by about 80% in relation to the use of GPS data only. Introducing the GRIB file from NOAA global forecast analysis has allowed us to initialize the tomographic computation with reliable atmospheric conditions and assign the wet refractivity values in the remaining unresolved voxel. One of the main advantages of the availability of a tomographic model of the lower atmosphere is the possibility to estimate the LOS travel time of the SAR signals.
FUNDING This study has been carried out in the frame of the EC  FP7 Mediterranean Supersite Volcanoes project  ( ECGA 308665 ) funded by the  European Commission  and supported by  H2020   EVEREST  project funded by the  European Commission. Copernicus Sentinel-1 data (2015) are provided by  ESA  in the framework of  GEO-GSNL initiative  and available at the  Copernicus Open Access Hub  (https://scihub. copernicus. eu). The MODIS dataset were acquired by the  INGV-ONT  real-time KSG multi-mission satellite acquisition system funded by  Italian Presidenza del Consiglio dei Ministri -Dipartimento della Protezione Civile (DPC).