dataesr/scientific-paragraphs-categorization huggingface dataset:
Computational efficiency Besides the accuracy and precision of the different packages that compute carbonate system variables, some users with large data sets may be concerned with computational efficiency. To assess differences in computation time between packages, we chose to use the global gridded data set described in Sect. 3. 1. With nearly 1 million ocean grid points, the computational time needed to compute all carbonate system variables varies by more than a factor of 1800 between packages (Table  2 ). The slowest package (swco2-Excel) required more than half a day, while the fastest (mocsy) needed 30 s. Packages based on spreadsheets are generally slower than those run by directly calling routines with programming languages (Fortran 95, MATLAB, Visual Basic). The latter are usually coded so that the equilibrium calculations are made one time for each set of input data, whereas spreadsheets often repeat the same set of calculations for each computed variable (each cell).
5. (6)  Millero (2010)  provides this set of the six α i coefficients for each of K 1 and K 2 on the original seawater pH scale but also on the free and total scales, by making analogous fits after converting the measured constants to each of those other scales  (Millero, 2010, Eqs. 9 to 12). With the  Millero (2010)  formulation, the four packages with internally consistent results are CO2SYS-MATLAB, CO2SYS-Excel-Pierrot, CO2SYS-Excel-Pelletier, and mocsy (Fig. 14 ). All use the coefficients for his formulation on the seawater scale. In seacarb, however, it is more complicated. At the surface, seacarb uses Millero's set of coefficients on the total scale to compute K 1 and K 2 ; conversely, below the surface seacarb uses Millero's coefficients for the seawater scale and then converts the resulting constants to the total scale after making the pressure correction (as is appropriate).
Although CO2calc's main code was taken from CO2SYS, the CO2calc developers included the Millero (2010) K 1 and K 2 formulations themselves, before they were available in CO2SYS. Lacking the CO2calc source code, we studied discrepancy patterns and made sensitivity tests. These differences appear to come from a different number of significant figures in one of the coefficients. More precisely, Fig. 15  compares the constants computed from the published set of coefficients  (Millero, 2010, Tables   2  and 3 ) to those computed with the unpublished yet more precise set of coefficients used by Millero, i. e. , his spreadsheet for the same publication (F. J. Millero, personal communication 2013). Only the α 5 coeffi-   Millero (2010)  over the same ratio for the reference. We refer to that ratio of ratios as the K 1 /K 2 quotient.
The reference is chosen arbitrarily as Millero's set of coefficients on SWS scale from his spreadsheet. Curve colors and line patterns are as in Fig. 15. cient differs, being given with one additional decimal place in the spreadsheet (Table  13 ). The difference between constants computed with the published coefficients and those computed with unpublished coefficients, i. e. , the spreadsheet coefficients on the seawater scale, match the pattern and magnitude of the differences between CO2calc and CO2SYS-MATLAB (compare Figs. 15 and Fig. 14 ). Hence CO2calc developers appear to have used the set of seawater-scale coefficients from Millero's spreadsheet. To further confirm that the different sets of coefficients yield fundamentally different results, we made two additional comparisons. First, we compared their K 1 /K 2 ratios (Fig.
16 ). By definition, that ratio should be independent of the pH scale, yet with different sets of coefficients it differs by up to 1. 5% over the observed range of T and S. Secondly, we made a sensitivity test with CO2SYS-Excel-Pelletier's two options for the Millero (2010) formulation (two sets of coefficients), both of which became available after our discussion paper was published. With that package, users choose to use either the published seawater-scale coefficients (as in the CO2SYS-MATLAB) or the unpublished seawater-scale coefficients from the spreadsheet (as used by CO2calc). With the former, CO2SYS-Excel-Pelletier agrees with the reference; with the latter, results are like those for CO2calc. Therefore, all three simple analyses indicate that it is the different sets of coefficients that cause substantial differences in computed variables.
K B and K W : principal non-carbonate alkalinity constants Previous comparison revealed significant discrepancies in subsurface variables computed from the swco2 package (Fig. 5 6 7 8 9 ), which are not due to K 1 and K 2 (Fig. 13 ). These discrepancies occur only when A T is a member of the input pair or when A T is computed, suggesting that they stem from the need to correct from total to carbonate alkalinity. The largest factor to correct for is borate alkalinity; hence it is also the most likely cause. Indeed, comparison of K B computed by the different packages does reveal discrepancies for swco2 (Fig. 17 ). Furthermore, swco2's divergence from CO2SYS-MATLAB increases linearly with depth, consistent with discrepancies in computed variables. Although we do not have access to the code for swco2, this discrepancy is consistent with a sign error in its a 2 pressure-correction coefficient for K B , as verified by sensitivity tests in seacarb.
K F and K S : constants to change pH scales For K F , all packages agree with the check value to its third and final digit after the decimal. The K F from the swco2 package does diverge from the others as salinity increases, but it still agrees with CO2SYS-MATLAB within pK F = 0. 00006, and its discrepancies do not change with temperature or pressure (Fig. 18 ). For K S , all packages agree with the check value, again to its third and final digit after the decimal. Beyond that, packages agree even more tightly, with the largest divergence reaching pK S = 0. 0002 for swco2. There are visible differences for K S computed from the different packages that merit further investigation (e. g. , positive excursions under surface conditions for swco2), but they remain quite small. Given the negligible consequences, we leave their resolution to future work.
Conclusions To assess the consistency of carbonate chemistry software packages, we have compared computed variables from 10 publicly available distributions, identifying discrepancies and causes. This comparison has led to improved agreement. Since our discussion paper was published  (Orr et al. , 2014) , there has been a fivefold reduction in discrepancies in pCO 2 and CO 2- 3 when all packages use the set of constants recommended for best practices  (Dickson et al. , 2007). The small discrepancies that currently exist remain insignificant even after pressure adjustments are made for the cold waters that pervade the deep ocean. Only in warm deep waters, such as found in the Mediterranean Sea, are there significant discrepancies and only for one package (swco2), e. g. , when pCO 2 and CO 2- 3 are computed from the A T -C T input pair.
, 1973, Eq. 8). For reproducibility, users should cite not only the package name but also its version number and all chosen equilibrium constants as well as any other selected package options. Developers would facilitate future comparison and reproducibility by providing computed constants as output, releasing their source code, and offering access to older versions, ideally through a public revision control system. Developers are also encouraged to consider providing results for potential and in situ pCO 2 , which differ greatly in deep waters. Although our focus has been on public packages, it is just as necessary to validate privately developed code. To facilitate such validation, we provide a supplement containing a subset of the data produced for this comparison. The Supplement related to this article is available online at doi:10.
Acknowledgements. This comparison was conceived by the  SOLAS-IMBER Ocean Acidification Working Group (SIOA). Design of the comparison, data production, and preliminary analysis was originally supported by the  International Ocean Carbon Coordination Project (IOCCP)  with funds from the  U. S. National Science Foundation  (grant  OCE-1243377 ) to the Scientific Committee on Oceanic Research (SCOR). Follow-up support came from the  Ocean Acidification International Coordination Centre (OA-ICC ), a project of the  Peaceful Uses Initiative of the International Atomic Energy Agency (IAEA). Final analysis and the publication of this article were supported by the  French ANR  project  MACROES (MACRoscope for Oceanic Earth System   ANR-09-CEP-003 ) and two  EU FP7 projects, CARBOCHANGE  (grant  264879 ) and  MedSeA  (grant  265103 ).
Geographic Area The area of study is a band crossing the Canadian polar shelf at about 76 N from Byam Martin Channel in the west to Lady Ann Strait in the east (Figure  1 ). Since ice within this band is typically immobile (fast) between December and July, its wintertime development can be conveniently monitored. Byam Martin Channel is dominated by cold water flowing south from the Arctic Ocean. It has a thick stable ice cover in winter and thereby is a suitable reference site against which to compare ice and ocean characteristics at locations further east-tidal polynyas in Penny and Cardigan Straits, fast ice in Jones Sound and fjords of southern Ellesmere Island, and a wind-forced polynya in Lady Ann Strait (see Figure  1  for locations). The data central to this study were collected during the winters of 2011-2012 and 2012-2013 by one author (EB) in northern Jones Sound.
During the first winter, EB was based on his yacht Vagabond in South Cape Fjord; the focus at this time was acquiring time series of air temperature, ice and snow thickness, and ocean temperature and salinity at a few locations in the fjord. In addition, several ice-thickness surveys were completed within and beyond the fjord. During the second winter, Vagabond was moved to the village of Grise Fjord; the emphasis for this winter was on wide-ranging surveys via snow-mobile of ice and seawater within Jones Sound and the fjords on its northern side; a time series was maintained near Grise Fjord but less intensively than in South Cape Fjord during the preceding winter. Additional data have also been utilized including: regional CTD and ice thickness surveys in late winter in  1983, 1984, 1999, and 2001 : T, S, and ice thickness measurements from instruments moored in Byam Martin Channel  (2011) (2012) (2013)  and Penny  Strait (2009 Strait ( -2010)).
Ice thickness and ocean properties on the Canadian polar shelf have been surveyed during winter since the 1970s  [Birch et al. , 1983]. Data from surveys in  April 1999 April , 2001 April , and 2011   where ice was thinnest. The association of thin winter ice with elevated under-ice temperature (Figure  3 ) suggests strongly that oceanic heat flux is retarding ice growth in these areas. However, since snow depth was not measured, its possible variation across the area could plausibly be implicated. Although not shown here, conditions were quite different during a 5 h interval before mooring retrieval in Byam Martin Channel in April 2012: here, the temperature just below the ice was steady at 1 m C below freezing (61 m C) while that at 31 m was 320 6 5 m C above freezing. It seems that near-surface water here was ready to freeze, a condition consistent with the 1.
If the inflow at the seabed occupied the bottom 5 m, its average speed would be 2 cm/s. These are weak currents. Because energy is not available from the wind or an estuarine circulation in winter, tidal energy is the only plausible source to drive upwelling. The steadiness of upwelling in South Cape Fjord is consistent with tidal forcing, as also are the differences among the fjords of southern Ellesmere Island in the impact of upwelling on ice growth. Differences are to be expected because each fjord is topographically unique and interacts differently with the tide in dissipating energy [e. g. ,  Farmer and Smith, 1980]. Although the lack of good bathymetric data precludes informed discussion of tidal influence, a continuous series of temperature-salinity data recorded in South Cape Fjord from 29 April to 10 July 2012 provides some insight.
CTD data were analyzed to estimate the phase speed of the first-mode internal wave; it decreased from 29 cm/s in November 2011 to 15 cm/s in March 2012. Froude number was calculated following  Farmer and Smith [1980] , using the speed of tidal flow at the sill (not that on approach because the height-above-surrounds of fjord sills typical almost blocks inflow below sill depth). The mode-1 Froude number at the sill in South Cape Fjord ranged between 0. 7 (subcritical) in November and 1. 3 (supercritical) in March. These estimates are with tidal generation and release of lee waves at the South Cape sill, as documented elsewhere by  Farmer and Smith [1980]. The observation of apparent tidally generated wave packets in South Cape Fjord raises the possibility that a range of hydraulic phenomena occurs here as tidal flow and stratification vary at the sill.
The resolving kernel approach allows one to quantitatively describe the improvement of the resolution of the density models achieved by merging gravity data and muon radiographies. The method developed in this paper may be used to optimally design the geometry of the field measurements to be performed in order to obtain a given spatial resolution pattern of the density model to be constructed. The resolving kernels derived in the joined muon-gravimetry case indicate that gravity data are almost useless for constraining the density structure in regions sampled by more than two muon tomography acquisitions. Interestingly, the resolution in deeper regions not sampled by muon tomography is significantly improved by joining the two techniques. The method is illustrated with examples for the La Soufrière volcano of Guadeloupe.
, 2010 Lesparre et al. , , 2012c;; Shinohara and Tanaka, 2012; Eppelbaum and Khesin, 2012; Portal et al. , 2013; Carbone et al. , 2014). Muon radiography is a straight-ray transmission method involving a Radon transform that markedly differs from the 3-D integrative gravity method. As will be discussed in the next section, muon radiography uses the flux of muons "coming from above", and it is limited to the imaging of shallow structures located above the particle detector. Also, muon radiography only provides information on parts of the density structure that are crossed by the rays, contrarily to the gravity method that brings information on the whole density distribution. These differences between both methods motivate the joining of both types of data to elaborate density models of complex geological structures.
Studies combining muon data and gravity measurements remain scarce, and we emphasize the early study by  Caffau et al. (1997) , who compared muon tomography with gravity measurements. More recently,  Davis and Oldenburg (2012)  and  Nishiyama et al. (2014)  presented joined inversions of gravity data and muon tomography using a straightforward linear regularized inversion based on block models. In the present study, we develop a quantitative methodology to examine how information brought by gravimetry data and by muon radiography may be joined to improve the resolution of density models of highly heterogeneous structures like altered active volcanoes. Fusion of information is studied through an approach based on the resolving kernels which provide a way to quantitatively evaluate the resolution of the resulting density model independently of any particular parametrization (e.
g. block discretization). Resolving kernels only depend on the geometrical properties of the data acquisition (i. e. locations of measurement points and telescope acceptance functions), and they allow one to perform prior analysis to evaluate the model improvements that may be expected by joining additional gravimetric data or muon radiographies. In this way, field measurements may be optimized with respect to the characteristics aimed at for the resulting density model. We begin by establishing the relationships between the density structure and both muon tomography and gravity data. Next, we derive the resolving kernels respectively corresponding to individual gravity and muon inversions and to joint gravity-muon inversion. The resolving kernels translate the information contained in the data into information concerning the density structure.
Basic principles of cosmic muon radiography The small cross section of muons in ordinary matter  (Barrett et al. , 1952)  allows the hard component of the muon spectrum  (Dorman, 2004; Tang et al. , 2006; Gaisser and Stanev, 2008)  to cross hectometres, and even kilometres, of rock. Most muons crossing the rock volume have a negligible scattering relative to the instrument angular resolution and travel along straight trajectories, ranking muon tomography among the class of straight-ray scanning imaging methods (2) muon tomography data coverage. The lines represent the observation axes of the muon telescope when located at Rocher Fendu (red), Ravine Sud (yellow) and Savane à Mulets (green). (Malmqvist et al. , 1979; Marteau et al. , 2011). In practice, muon tomography is performed by using a series of pixelated particle detectors that allow one to determine the trajectories of the muons passing through the rock body.
Muon tomography The primary information used in muon tomography consists in cosmic muon flux attenuation measurements resulting from the screening produced by the geological volume to be scanned. Attenuation is measured by counting the number of muons emerging from the volume for each observation axis, s m = (r m , P m (ϕ, θ )), of the telescope (Fig. 1. 2). r m represents the position of the telescope, and P m the observation axis acceptance pattern which depends on (ϕ, θ ), the azimuth and zenith angles referenced at r m (see Fig. 2 ). Note that r m is the same for all the observation axes on a given site. P m depends on the telescope geometry and angular orientation on the site. Our standard field telescopes count 31 × 31 observation axes and, in a field experiment where the telescope successively occupies several places around the target, the number, M, of data may easily reach several hundredths.
, 1985). This implicitly assumes that the ζ k are linearly independent with respect to the retained inner product; i. e. no acquisition kernel can be written as a linear combination of the other kernels. The noticeable instances where this important assumption is not satisfied correspond to situations where several data have been acquired identically, i. e. either at the same location for gravity measurements or with the same position and orientation of the telescope for muon tomography. In such cases the dimension of X K is reduced since the redundant data may be merged (i. e. averaged) into a single one. The best density distribution that can be recovered through the inversion process (it is the best because it takes all the information contained in the data and makes the fewest hypotheses about the X K complementary subspace) is a linear combination of the generating functions, ρ(r) -ρ 0 (r) = K k=1 a k × ζ k (r).
This is achieved through the projection γ of (r, r ) onto a δ(rr ), γ (r) = (r, r ), δ(r -r ) X = V w(r, r ) × (r, r )dr. (26) Here the X scalar product is defined by Eq. (  15 ). Here, the Hamming function is used as the weight function w(r, r ) = H L (||r -r ||) K w × 1 + cos 2π × ||r -r || L , (27) Geosci. Instrum. Method. Data Syst. , 4, 177-188, 2015 www. geosci-instrum-method-data-syst. net/4/177/2015/ where K w is a normalizing constant, H L is a rectangular pulse that restricts w to |r -r | ∈ [-L/2; L/2], and L = 25 m. If r or r is located outside the volcano, we take w(r, r ) = 0. This choice is explained in Sect. 5. 4. We now display resolving kernels corresponding to the data acquisition shown in Fig. 1. 1 for muon tomography and in Fig. 1. 2 for gravimetry. Muon radiographies are taken from three sites equidistantly located along the southern edge of the volcano.
Gravity measurements are assumed to be made on a regular grid over the entire lava dome. Accounting for the fact that the acquisition kernels ζ k are either for gravity or for muon tomography (Eq. 16), we successively consider the case of resolving kernels obtained for muon tomography alone, for gravity data alone, and for a combination of muon tomography and gravity data. We compute (r, r ) for two positions r = {r 1 ; r 2 } located along a vertical line that goes through the centre of the dome (Fig. 5 ). Points r 1 and r 2 are respectively inside and below the volume of the lava dome spanned by the lines of sight of the telescopes (Fig. 1. 1). The parameter γ is computed and plotted on the four characteristic slices represented in Fig. 5.
The important heterogeneities inside the muon tomography acquisition kernels forbid us from using the Radon transform mathematical corpus. For an equivalent resolution and scanning, the kernels can be regularized by taking the telescope away from the volcano and reducing the angular aperture. We then get into the typical experimental conditions of a medical X-ray tomography. But, the consequences are a weaker particle flux (a longer acquisition time) and a greater sensitivity to potential noises. So, a compromise has to be found, but (2) Ravine Sud tomography data (1) Gravimetry data (4) Gravimetry and tomography data (3) All sites tomography data 1 0 -1 -2 -3 -4 Resolving kernel value normalized at point 1 in logarithmic scale Figure  6. Resolving kernel at point r 1 (Fig. 5):  (1) the gravity data alone; (2) Ravine Sud muon radiography alone; (3) combination of three muon radiographies; and (4) joined muon and gravity data sets.
We draw the reader's attention to the fact that, despite their compact support, the acquisition kernels overlap each other for neighbouring main lines of sight (see Fig. 3 ). As will be seen below, this characteristic is fundamental to understanding the shape of the resolving kernels. A muon tomography resolving kernel is a linear combination of muon tomography acquisition kernels M (Eq. 23), and the inversion process optimizes the b i coefficients to obtain the best density model (Eq. 24). Figure  6. 2 and 6. 3 show the resolving kernel for point r 1 for different combinations of tomography data sets. When using data acquired from a single place located at the southernmost edge of the volcano (Ravine Sud; see Fig. 1. 2), the resolving kernel (Fig. 6. 2) encompasses lines of sight spanning a limited range of azimuths.
Joined muon tomography and gravimetry kernels We now consider resolving kernels computed by using both muon M and gravity G acquisition kernels. Gravimetry does not improve significantly the inversion process at r 1 , and the resolving kernel (r 1 , r ) (Fig. 6. 4) looks very similar to the one obtained for the muon radiographies alone (Fig. 6. 3). The information provided by muon tomography is dominant relative to gravimetry, except at the immediate vicinity of the gravity measurement points. The situation is very different for point r 2 , where the resolving kernels obtained by joining muon radiographies and gravity data (Fig. 7. 2) appear very different from the gravity kernel (Fig. 7. 1). The most conspicuous effect is that muon data compensates for the great sensitivity of gravimetry at near-surface locations by shifting the centre of the mass of the resolving kernel downward.