dataesr/scientific-paragraphs-categorization huggingface dataset:
First, we integrate in time the 6-hourly OPA-nex output to obtain an output with the same time frequency as the observations (i. e. daily for sea ice concentration, 2 d displacement for the drift, and weekly for the thickness/volume in CS2SMOS). We then interpolate OPA-nex outputs onto the same grid as the observations. There is an additional step for the ice thickness and volume: as PIOMAS is only available monthly, we do a monthly average of OPA-nex outputs and CS2SMOS data and compare the three datasets with PIOMAS data interpolated on the CS2SMOS 25 km grid.
Sea ice volume and area export through Fram Strait We also evaluate the ice volume export through Fram Strait, as it is an important contributor to the Arctic sea ice mass balance  (Spreen et al. , 2009; Ricker et al. , 2017; Spreen et al. , 2020). We use the time series produced by  Spreen et al. (2020) , which covers the period 1992-2014, overlapping with most of our simulation. This dataset is based on sea ice thickness measurements derived from one to four upward-looking sonar measurements installed on moorings across the strait  (Vinje et al. , 1998)  (from which the section of ice thickness is extrapolated), and sea ice drift is retrieved from the JPL sea ice motion dataset  (Kwok et al. , 1998). We also compare OPA-nex results with the time series of  Spreen et al. (2009 Spreen et al.
( ), spanning 2003 Spreen et al. ( -2009 , that is obtained by combining sea ice thickness from ICESat altimeter observations and sea ice area and drift retrieved from AMSR-E 89 GHz. To compare our export to these datasets in a consistent way, we first estimate the sea ice transport (i. e. motion vectors) at the output frequency of OPA-nex (6 h), and then we interpolate these transports onto the same grid as  Spreen et al. (2020). The gate used for the computation is located at 79 • N (see Fig. 2 ). We also compare the simulated sea ice area flux through Fram Strait to the time series estimated by  Smedsrud et al. (2017). They combined a blended historical and modern record of sea ice concentration from the National Snow and Ice Data Center (NSIDC,  Walsh et al. , 2017)  with sea level pressure observations across Fram Strait to retrieve the sea ice area flux over the period 1935-2015.
Observed contributions to sea ice mass balance In order to evaluate the ice mass balance and its spatiotemporal variations, we make use of the ice volume change dataset from  Ricker et al. (2021). The originality of this dataset is that it separates the dynamic contribution (i. e. the import and export of ice in each region) from the thermodynamic contribution (i. e. the net sea ice growth in winter) to ice volume change in the freezing season (November to March) for six Arctic seas (Barents, Kara, Laptev, East Siberian, Chukchi, and Beaufort; see Fig. 2 ). To estimate these contributions,  Ricker et al. (2021)  combine sea ice volume  (Hendricks et al. , 2018)  and motion information  (Girard-Ardhuin and Ezraty, 2012)  to retrieve the dynamic volume change and then subtract it from the total volume change to estimate the net sea ice growth.
Therefore, OPA-nex successfully captures the seasonal cycle of the ice extent, and its interannual variability, particularly in the summer, with the exception of 2016 and 2017 where the ice loss is overestimated. We then examine sea ice volume (Fig. 4a ). Agreement with PIOMAS is generally good, although sea ice volume in OPA-nex is generally lower than in PIOMAS in the early 2000s. After 2008, the agreement becomes very good in both winter and summer, and the two models show a similar interannual variability. A lot of factors could explain the discrepancies between OPA-nex and PIOMAS (differences in atmospheric forcings, differences in the dynamics and thermodynamics of the models, and the use of data assimilation in PIOMAS), and it is difficult to attribute these differences to one of these factors or the other.
OPA-nex agrees well with CS2SMOS, with an average RMSE of 0. 34 m and an average bias of 0. 03 m for the whole period when observations are available (from October to April each year from 2011). This is also true for the ice thickness distribution during the ice growth season (Fig. 4b, c, d ). Biases in the distribution in OPA-nex compared to CS2SMOS are quite similar to those in PIOMAS and also found in most sea ice models: thick ice is not as thick as the observations in the central Arctic, while thin ice is often too thick, particularly on the western side of the Arctic Basin  (Stroeve et al. , 2014; Watts et al. , 2021). The simulated drift generally shows a good agreement with the OSI-SAF data (Fig. 5 ), with a low negative bias (-0. 35 km d -1 on average from 2010 to 2018) and a low RMSE (3.
Regional winter ice mass balance We now investigate the dynamic component of the ice mass balance (the net balance between import and export of sea ice) within the Arctic Basin, in a similar way to Fig. 3a  of https://doi. org/10. 5194/tc-17-617-2023 The Cryosphere, 17, 617-638, 2023 Ricker  et al. (2021). The ice transport contribution to the mass balance for the regions where data are available is very well estimated in OPA-nex (Fig. 7 ). The variability is well captured, with determination coefficients generally exceeding 0. 50 between the detrended OPA-nex results and estimates from observations. Similarly to  Ricker et al. (2021) , we do not find any significant trend over the period 2002-2018 (note that they also include 2019 in their study) for any of the regions analysed here. This is also true for the central Arctic regions (west and east) that are not included in  Ricker et al.
(2021). Figure  8  shows the same analysis but for thermodynamic processes, comparable to Fig. 3b  of  Ricker et al. (2021). The thermodynamic processes included in OPA-nex are the production and growth of young ice, the basal growth of (thicker) ice, and ice formed due to the flooding of snow, as well as basal and surface melt. The magnitude of the net winter growth is estimated well in general. The main discrepancies between OPA-nex and  Ricker et al. (2021)  are found in the Kara and East Siberian seas, where OPA-nex overesti-mates the amount of ice formed every winter. This overestimation of ice growth in these seas is also visible in the data from PIOMAS, shown as a reference in  Ricker et al. (2021). As in  Ricker et al. (2021) , we find a small but statistically significant (i. e. with a p-test result lower than 0.
org/10. 5194/tc-17-617-2023 The Cryosphere, 17, 617-638, 2023 case of divergence in Fig. A1 ). We found that P = 3 kPa and c lab = 2 MPa was a good compromise between thickness distribution (Fig. 4 ), sea ice drift (Fig. 5 ), and deformation patterns (Fig. A1 ). We also modified some values associated with the thermodynamics. The maximum thickness value of the young ice category has a strong effect on the slope of the sea ice volume evolution during the autumn growth. We chose a value that gave a slope similar to ice volume growth as estimated in CS2MOS. The values for the snow and sea ice albedos were chosen, in their physical range, to give a reasonable match between modelled and observed sea ice extent and thickness. Financial support. This research has been supported by the Norges Forskningsråd  (grant nos.
Acknowledgements. This research has been funded by the  Norwegian Research Council  ( Nansen Legacy  grant no. 276730 ,  FRASIL  grant no. 263044 , and  ARIA  grant no. 302934 ),  JPI Climate  and  JPI Oceans  ( MEDLEY  project, under agreement with the  Norwegian Research Council , grant no  316730 ), the  SASIP  project, and by  Copernicus Marine Environment Monitoring Service (CMEMS)   WIzARd  project. SASIP is supported by  Schmidt Futures , a philanthropic initiative that seeks to improve  societal outcomes through the development of emerging science and technologies. CMEMS is implemented by  Mercator Ocean  in the framework of a delegation agreement with the  European Union Copernicus Marine Environment Monitoring Services  (contract no. 69) and the  European Space Agency  through the  Cryosphere Virtual Laboratory (CVL , grant no.
Data availability. The OSI-450 sea ice concentration product is available at ftp://OSISAF. met. no/reprocessed/ice/conc/ v2p0 (last access: September 2021, Lavergne et al. , 2019). The OSI-430-b sea ice concentration product is available at ftp://OSI-SAF. met. no/reprocessed/ice/conc-cont-reproc/v2p0 (last access: September 2021, Lavergne et al. , 2019). PIOMAS outputs are available at http://psc. apl. uw. edu/research/projects/ arctic-sea-ice-volume-anomaly/data/model_grid (last access: August 2021,  Zhang and Rothrock, 2003). The CS2SMOS sea ice thickness product is available at ftp://ftp. awi. de/sea_ice/product/ (last access: September 2021, Ricker et al. , 2017). The lowresolution daily sea ice drift product from OSI-SAF can found at ftp://osisaf. met. no/archive/ice/drift_lr/merged (last access: September 2021,  Lavergne et al.
, 2010). Monthly outputs of all quantities discussed in the manuscript are available on Zenodo as NetCDF files. We also share the data used for each figure, also as NetCDF files (https://doi. org/10. 5281/Zenodo. 7277523,  Boutin et al. , 2022). The neXtSIM code is still in development and will be made open source in the coming months (in a dedicated publication).
A simulation study shows that the CP generates patterns in comparative data that are distinct from any Gaussian process, and that Restricted Maximum Likelihood (REML) parameter estimates and root trait reconstruction are unbiased and accurate for trees with 200 tips or less. The CP has only two parameters but is rich enough to capture complex pulsed evolution. It can reconstruct posterior ancestral trait distributions that are multimodal, reflecting the uncertainty associated with the inference of the evolutionary history of a trait from extant taxa only. Applied on empirical datasets taken from the Evolutionary Ecology and Virology literature, the CP suggests nuanced scenarios for the body size evolution of Greater Antilles Lizards and for the geographical spread of the West Nile Virus epidemics in North America, both consistent with previous studies using more complex models.
In order to model pulsed evolution, we consider here the Cauchy process (CP) on trees, that is a simple and (now) tractable pure-jump process in which the trait increment along each branch follows an independent centered Cauchy distribution of dispersion proportional to its length. The Cauchy process is comparable to the Brownian motion in many aspects. Both have only two parameters, namely the root trait value and the process variance (for the BM) or dispersion (for the CP). In addition, both are limiting cases of a Lévy-stable process: where the BM can be seen as a purely gradual model, the CP is a pure-jump process, with no diffusion component, and can be seen as a purely pulsed model. Both models are hence useful to have in our analysis toolbox, and the relative fit of one model against the other can quickly give us an idea of the general pattern of trait evolution on a specific dataset, either gradual or pulsed.
We provide here some guidelines and diagnostic tools to asses whether computations succeeded on a particular tree. An important feature of the CP is that, contrary to Gaussian models that always give rise to Gaussian posterior ancestral state distributions  (Royer-Carenzi and Didier, 2016) , it can reconstruct multimodal trait distributions at ancestral nodes (such a behavior was observed for other Lévy processes in  Elliot and Mooers 2014). The different reconstructed modes are well suited to represent the various evolutionary scenarios that are compatible with the observed dataset on extent species, and are an elegant way to handle the identifiability issues that arise in shifted Gaussian models  (Bastide et al. , 2017). The article is organized as follows. We first outline the mathematical and modeling properties of the CP, and present the algorithms for likelihood computation and ancestral state reconstruction.
To illustrate the behavior of the CP in such a situation, we simulated a dataset using one of the shifted BM models presented on Figure  2 , and then performed ancestral and increment reconstructions using our CP framework. We observe that these reconstructions reflect all three scenarios (a), (b) and (c). In particular, the posterior density of the "dot-node" state is multimodal, with three modes around to 0, -1 and 2 (Fig. 2-e ), which are indeed the ancestral expectations of this node in scenarios (a), (b) and (c), respectively. Similarly, the modes of the posterior increment density on the two daughter branches are consistent with the possible shift values on these branches (Fig. 2-f,g ). As we will see in our empirical datasets analyses below, the CP reconstruction hence provides us with at least two important pieces of information.
Model Selection Setting. To test our ability to distinguish gradual from pulsed evolution using the CP, we simulated datasets using various models of evolution, and tried to recover the true underlying model. In all simulations, we used the empirical tree from the Lizard dataset  (Mahler et al. , 2013 , see below), normalized to unit height. We simulated a first dataset according to a simple BM, with root fixed to 0, and variance σ 2 BM = 1. As the CP does not have a variance, we used the population median absolute deviation (MAD) to produce datasets with similar characteristics. It is well known that the MAD of a standard Gaussian is equal to its third quartile Φ -1 ( 3 4 ) = 0. 6745, and that the MAD of a centered Cauchy variable is equal to its dispersion parameter (see e. g. Rousseeuw and Croux, 1993).
We hence simulated a second dataset according to a CP with dispersion parameter σ = 0. 6745, so that it had the same MAD as the standard Brownian dataset. To see if stabilizing selection or explosive change could lead to the spurious selection of a CP, we simulated two datasets using an Ornstein-Uhlenbeck (OU;  Hansen, 1997)  or and Early Burst (EB;  Harmon et al. , 2010)  process. For the OU, we took the phylogenetic half-life of the process to be 20% of the tree-height, keeping a stationary variance of γ 2 OU = 1, so that the selection strength was equal to α OU = 3. 47 and the process variance to σ 2 OU = 6. 93. Exploiting the similarities between the OU and the EB  (Uyeda et al. , 2015; Aristide et al. , 2018) , we took and Early Burst with the same variance σ 2 EB = σ 2 OU , and opposite rate r EB = -α OU.
This α N IG parameter can hence be used to produce datasets that have patterns of variation in between the BM and CP models. We simulated datasets according to the NIG process, all with the same MAD value equal to Φ -1 ( 3 4 ) as before, and with an excess kurtosis varying from 1 to 10 (larger kurtosis being associated with lower α N IG values). All the parameters used in this simulation study are summarized in Table  1. Table  1 : Parameters used for the simulation of the processes on the Lizard tree. σ is the standard error parameter for Gaussian processes, and the dispersion parameter for the CP. µ is the starting value. The population MAD is the same for all processes, except for the EB which parameters were taken from the OU. The excess kurtosis for Gaussian processes is 0, varies from 1 to 5 for the NIG processes, and is not defined for the CP.
Each dataset was simulated 100 times. On each dataset, we used function phylolm  (Ho and Ané, 2014a)  to fit the the Gaussian processes (BM, OU with fixed root, EB), both with or without error, and function fitCauchy from our cauphy package to fit the CP. We used the maximum likelihood framework (with fixed root) for each fit, and selected the best model using the simple AIC score  (Akaike, 1974). Gaussian Processes are not Mistaken for Cauchy Processes. When the dataset was simulated under a Gaussian process (BM, OU or EB), the AIC always selected a Gaussian process, with the right model being chosen in the large majority of cases, confirming that there was a good recovery power also among Gaussian processes (see Fig  3 ). The other way around, when the dataset was simulated according to a CP, the AIC always selected a CP.
These results show that the signature of the CP is very different from the signature of a Gaussian process, and that there is a very high segregation power between these two modes of evolution. The Cauchy Process Captures Pulsed Evolution. When the data was simulated according to a NIG process, the true simulation model was not included in the fitting models. From Figure  3 , we can see that the Cauchy process is largely preferred over Gaussian processes for the NIG datasets, as soon as its excess kurtosis is large enough, that is, when the NIG goes away from the BM. As the NIG is a pure jump process that can also model non-gradual evolution  (Landis and Schraiber, 2017) , this result shows that the CP is able to pick up a different kind of pulsed evolution.
Parameter Inference Setting. We assessed the quality of the ML and REML estimators of the dispersion and root position parameters of a CP on a tree using simulations. We simulated pure birth trees conditioned on having from 10 to 250 tips using package ape  (Stadler, 2011; Paradis, 2015) , that we re-normalized to unit height. We then simulated trait data on each tree using a Cauchy process with dispersion 1 and root value 0. Each configuration was replicated 500 times. We then fitted a CP on each dataset, using either ML (with fixed root) or REML, and compared the estimated values of the parameters with the true simulation values. For the REML, we used ancestral trait reconstruction to get the posterior trait density at the root on a grid of values between -10 and 10 with step 0. 01, and constructed a 95% highest posterior density interval (HDPI) using R package HDInterval  (Meredith and Kruschke, 2022).
A small fraction of 5. 5% of the datasets over all scenarios produced point estimates that had an absolute value larger than 2, far away from the true value of 0. These could be cases where, because of the jumping nature of the Cauchy process, the root value gets "forgotten" along the tree. However, Figure  6  shows that, in the REML case, the empirical coverage of the HDPI (fraction of simulation where the true value lies in the estimated interval) is close to its norminal rate of 95%, for all sizes of tree. 4 Empirical Examples.
Greater Antilles Lizards Data. The adaptive radiation of anoles in the Greater Antilles is well documented  (Losos, 2009) , and has been hypothesized to follow a Simpsonian "evolution by jumps"  (Mahler et al. , 2013). Previous analyses of this dataset includes models based on Brownian Motion or Ornstein-Uhlenbeck with shifts  (Mahler et al. , 2013; Eastman et al. , 2013; Bastide et al. , 2018) , or Lévy processes  (Duchen et al. , 2017). We considered here the dataset from  Mahler et al. (2013) , that consists of 100 lizard species on a timecalibrated maximum clade credibility tree, with each species associated with an ecomorph. As in  Duchen et al. (2017) , we focused our analysis on the log-transformed snout-to-vent length (SVL), which is highly correlated with all the other morphological measurements  (Bastide et al.
, 2018). Methods. We fitted a Cauchy process on this dataset, using the ML or REML method. Using the AIC criterion  (Akaike, 1974) , we compared these results with fits obtained using traditional Gaussian processes, such as the Brownian Motion  (Felsenstein, 1973) , the Ornstein-Uhlenbeck  (Hansen and Martins, 1996)  and the Early Burst  (Blomberg et al. , 2003; Harmon et al. , 2010) , all fitted with phylolm (Ho and Ané, 2014a) (ML) or Rphylopars  (Goolsby et al. , 2017 ) (REML). We then carried an ancestral state reconstruction using the Cauchy REML fit, on a grid of 100 values between 3 and 5. 5 for all trait values at internal nodes, and between -1. 5 and 1. 5 for all increments on branches. Model fit. Compared to Gaussian processes, the Cauchy process increased the log-likelihood by about 8 units, and had a much smaller AIC (Table  2 ).
Consistent with a pulsed evolution hypothesis, the CP was hence clearly favored over any Gaussian process. Table  2 : Fits of Gaussian (BM, OU and EB) and CP models on the lizard dataset from  Mahler et al. (2013) , with 100 species, and using the snout-to-vent length (SVL) trait measurement. The first two lines are the (restricted) log likelihood and AIC scores for the ML (left) and REML (right) fits. The parameters are µ for the ancestral root value (ML), σ for the Gaussian standard error or Cauchy dispersion, and θ represents the selection strength α for the OU, and the rate r for the EB. The computation time for each fit on a standard laptop is given in seconds. Ancestral Reconstruction. As in  Duchen et al. (2017) , we found evidence for large positive increments on branches leading to clades with species from the "crown giant" ecomorph, that have particularly large body sizes (Fig.
7 ). The clade highlighted with bold edges on Figure  7  is particularly interesting, as it exhibits bi-modal posterior trait and increment distributions, that can support two possible hypotheses about body size evolution of these species. One possible scenario is consistent with the one exhibited in  Duchen et al. (2017) , with medium sized ancestors ("orange" mode on ancestral nodes), and shifts towards large values on branches leading to species from the "crown giant" ecomorph, and to the clade of "false chameleons", namely A. guamuhaya, A. porcus, A. chamaeleonides and A. barbatus, from the former genus Chameleolis, that are known to have unique characteristics. However, an other scenario that is supported by this dataset would see a shift toward large body sizes at the crown of the clade leading to large ancestors ("purple" mode on ancestral nodes), and shifts towards smaller values on branches leading to the two smaller species of A.
West Nile Virus Data. The epidemic of WNV in America, from its first detection in New York city in 1999 to its spread in California and Mexico in the early 2000s, caused more than 1,200 deaths in the United States and was closely monitored by the health authorities. Pybus et al. (2012)  studied the geographical spread of the virus in an integrated framework using BEAST  (Suchard et al. , 2018) , with a dataset of 104 whole genomes collected between 1999 and 2007 along with the sampling coordinates (latitude and longitude). We reproduced the analysis using the best fitting model they selected, which consisted for the nucleotides of a general time-reversible (GTR) substitution model with a gamma distribution of among-site rate heterogeneity, an uncorrelated lognormal relaxed molecular clock and a Bayesian skyline coalescent model prior, and for the continuous locations of a relaxed random walk (RRW) with a gamma prior for the mixing rates, and an added uniform jitter with window size 0.
5 applied to duplicated locations. We followed the instruction from the associated BEAST workshop tutorial (available at https://beast. community/workshop_continuous_diffusion_wnv) to set up the analysis, and ran the MCMC chain for 250 million states, sampled every 50,000 states, using BEAST v1. 10. 5 (prerelease #dc4f12c). The convergence of the chain was assessed with Tracer  (Rambaut et al. , 2018) , and we obtained a maximum clade credibility (MCC) tree using TreeAnnotator, with times summarized with the common ancestor heights option to avoid negative branch lengths. In all subsequent analyses, we used the MCC tree and the jittered continuous locations as the primary dataset. Methods. We fitted a CP using our framework on both the latitude and longitude independently, and reconstructed ancestral nodes an increments on all branches using a grid of 1000 values.