dataesr/scientific-paragraphs-categorization huggingface dataset:
e. warm sea surface temperatures enhancing convective snow precipitation) play a role in increasing cold-spell hazards? In this paper, we focus specifically on Italy: recent cold and snowy spells in this country have caused casualties in the population, strongly affected ground and air transport, and caused disruptions in services. Our strategy to tackle these questions is to analyse simulations produced in a global circulation model (GCM) under different emission scenarios. We first validate the cold and snowy spells produced by a simplified GCM of intermediate complexity with historical forcing, i. e. the Planet Simulator (PlaSim)  (Fraedrich et al. , 2005a, b) , against those detected in a reanalysis dataset. Then we analyse dynamic analogues of cold and snowy spell events under different climate change scenarios.
Sources and dataset Our study is based on the detection of synoptic meteorological configurations leading to cold spells over Italy in PlaSim, considering a control run based on the recent historical climate and a set of three increased emission scenarios at steady state. In order to do so, we will proceed with the following steps: 1. identify large-scale, high-impact winter cold spells over Italy; 2. describe the dynamic and thermodynamic conditions associated with such cold spells; 3. detect cold-spell analogues in a historical climate dataset; 4. detect cold-spell analogues in PlaSim runs and evaluate whether climate change can significantly modify their frequency and in which direction; 5. characterize the PlaSim cold-spell analogues by analogy with point 2 to assess the potential of the considered dynamic configurations in producing relevant winter phenomena in a sensibly warmer climate.
In order to identify relevant cold spells over Italy, we consider documented events that have produced at least a record low temperature and/or a record snowfall amount (or snow at locations where snowfall has never been previously reported) at one or more locations in Italy. We combine official sources and both professional and avocational websites dedicated to weather and climate, where collections of weather event reports are available, and we countercheck their validity with station data and trusted documentary sources  (Bailey, 1994; Payne and Payne, 2004). Our documentary sources include local networks, newspapers and periodicals (see Appendix A); news and commercial meteorological websites  (ansa. it, 2020; 3bmeteo. com, 2022a; meteociel. fr, 2020; meteogiornale. it, 2022a) ; and temperature and hydrological records (evalmet.
Observed cold-spell dynamics Besides the qualitative analysis involving the cities of Bologna and Campobasso briefly presented in Sect. 2. 1, we now aim to characterize the dynamic and thermodynamic features of the considered cold spells at the synoptic scale. To this purpose, we rely on the National Centers for Envihttps://doi. org/10. 5194/esd-13-961-2022 Earth Syst. Dynam. , 13, 961-992, 2022 ronmental Prediction (NCEPv2) Reanalysis dataset. In particular, we consider geopotential height at 500 hPa (Z 500 ) and sea-level pressure (SLP) as dynamical fingerprints and to compute the analogues  (Jézéquel et al. , 2018) ; temperature at 850 hPa (T 850 ) to track cold air advection without surface disturbances  (Grazzini, 2013) ; 2 m temperature (T 2 m ) to characterize near-surface conditions; and daily precipitation rate (PRP).
In order to choose the optimal number of clusters, we first performed a scree plot (not shown), obtained by plotting the within-groups sum of squared differences from the cluster centroids. This analysis did not give clear indications about the ideal number of clusters. Therefore, we compared clustering results at different values of k, finding that for k = 3 two of the three clusters displayed very similar spatial features, and with larger k the resulting clusters can always be reduced to two main patterns: we then choose k = 2. We remark that the k-means algorithm and other clustering techniques are based on assumptions such as equal size and sphericity of the clusters, which can be met only in coarse approximation in real-world high-dimensional datasets. In particular, the poor indications from the scree plot may be due to the different number of events assigned to each cluster (respectively, 22 and 10 for k = 2).
The main differences concern the UK, Iceland and Scandinavia, due to the different direction of the high-pressure axis. In cluster 1, the tilted high-pressure drives warmer air towards the British Isles and Scandinavia, while the cold anomaly is confined to more southern latitudes. In cluster 2, the meridionally oriented axis brings warmer air towards Greenland and Iceland, while the core of the cold air is located between Scandinavia and central Europe. The T 2 m fields (Fig. 3c  and d ) generally show above-zero temperatures over central southern Italy during these events: this is partially due to the coarse resolution of the NCEPv2 dataset that regards southern Italy and Sicily as sea grid points. Italy is also a country with a complex geography and mountain ranges (the Alps and the Apennines) extending from the north to the south.
Bias correction Climate models, even those with higher complexity than PlaSim, are characterized by a finite resolution, thus leaving smaller scales unresolved, and contain several physical and mathematical simplifications that make climate simulations computationally feasible, while also introducing a certain level of approximation. This results in statistical biases https://doi. org/10. 5194/esd-13-961-2022 Earth Syst. Dynam. , 13, 961-992, 2022 that can be easily observed when comparing control runs to observations or reanalysis datasets. In order to mitigate the effects of these biases, a bias correction step can be performed. Bias correction usually consists of adjusting specific statistical properties of the simulated climate variables to a validated reference dataset in the historical period.
A total of 8 storm parameters were selected from a literature review to be used as proxies. Two different proxy types emerged from this literature review: microphysical and dynamical proxies. The proxies are calculated from 1h AROME forecasts from the assimilation cycle. Machine learning regression models are used to relate observed FEA and the simulated proxies. The training of the observation operator is performed on a dataset of 44 days and 3 additional days are used for the validation. The data are processed as a climatology over the whole domain (i. e. France) and time period. The performances of each proxy are evaluated by computing Fraction Skill Scores (FSS) between observed FEA and proxy-based FEA. The present study suggests that microphysical proxies seem to be more suited than the dynamical ones to model satellite lightning observations with the AROME NWP system.
Satellite lightning data present numerous advantages including their large spatial coverage, providing observations where ground-based radar data are scarce or nonexistent and making them well-suited for convective-scale data assimilation in limited-area NWP systems. By 2023, lightning observations from the Lightning Imager (LI) onboard the Meteosat Third Generation satellite (MTG) will be available over Europe, Africa, the Atlantic Ocean and a part of South America. The LI will be able to detect both IC and CG lightning, providing total lightning information, day and night. One of the products that will be provided by the LI is a flash count per pixel accumulated over time and hereafter referred to as flash extent accumulation (FEA)  1. In modern data assimilation systems like variational or ensemble based systems, an observation operator is required to establish a link between observations and the NWP system background.
For lightning observations, developing an observation operator is not trivial since most operational NWP systems do not include an electrification and lightning scheme. Indeed, the equations to relate the NWP system prognostic variables and the microphysics to produce an electric field are complex, nonlinear and computationally expensive (e. g. Barthe and Pinty, 2007; Mansell et al. , 2002). In consequence, the observation operator developed in this study is based on empirical relationships between the lightning observations and a set of proxies derived from the NWP system variables, in a similar approach as used in the assimilation papers mentioned above. The objective of the present study is to prepare the assimilation of lightning satellite data in storm-scale NWP systems by designing a suitable observation operator in the scope of development of a new EnVar assimilation system for the French regional NWP system AROME-France.
The lightning data employed here are generated from ground-based lightning data from the Météorage network  (Erdmann et al. , 2022)  to mimic the future MTG-LI data. This paper is organized as follows. In sections 2. 1 and 2. 2 the lightning generator to mimic the MTG-LI data as well as the NWP system configuration are described. The list of the tested proxies is detailed in section 2. 3, along with the way they are retrieved from the prognostic variables of AROME-France. Section 2. 4 depicts the method to link the lightning pseudo-observations to the selected proxies. The results for 10 min FEA are presented in section 3. In section 4, both lightning threats introduced by  McCaul et al. (2009)  are also investigated and the McCaul lightning calibration is reproduced to be compared with the results of the method described here.
Data and methodology In the context of the present study, 1 h forecasts from the AROME-France assimilation cycle for 47 stormy days in 2018, described in Table  1 , were used to compute a set of 8 storm variables to be tested as proxies for MTG-LI observations. Those storm variables are hereafter directly referred to as "the proxies". Synthetic MTG-LI observations were generated for the same 47 days as targets to train observation operators first from individual proxies and then from combinations of proxies based on their individual performances and their inter-correlations built on machine learning (ML) regression algorithms. In the following, a brief description of the lightning generator and the NWP system configuration is provided, as well as a detailed description of each proxy and how they are retrieved from the AROME-France forecasts.
MTG-LI synthetic observations As mentioned in the introduction, the geostationary LI will provide observations by 2023. The geostationary LI is an optical sensor that detects the cloud top illumination due to lightning with a resolution of 7 km at European latitudes and 4. 5 km at its nadir. Total lightning activity will be detected both day and night but LI detection efficiency will most likely vary with the time of the day, and possibly with the geographical position of the lightning activity within LI field-of-view. This assumption is based on related studies on the LIS instrument onboard the International Space Station (ISS-LIS) and the GLM  (Bateman et al. , 2021; Peterson et al. , 2017). Erdmann et al. (2022)  developed a method to generate GLM pseudo-observations using lightning data from the National Lightning Detection Network (NLDN), consisting of ground sensors in the contiguous US.
As the MTG-LI and GLM instruments are expected to provide similar observation and data structure, and an intercomparison study between the NLDN and The simulation of synthetic observations is performed in two steps. First, GLM (or MTG-LI in our case) flash characteristics, such as flash duration or flash extent, are simulated using ML supervised models  (Erdmann et al. , 2022). The ML models have been trained with flash characteristics of coincident NLDN and GLM flashes from a database of 10 lightning active days in the South-East of the USA, spread over a 6 months period. The second step consists in generating the pseudo FEA gridded observations from the simulated flash characteristics. The method was evaluated by comparing the GLM synthetic observations with operational GLM data. Erdmann et al.
NWP system configuration The AROME-France NWP system  (Seity et al. , 2011; Brousseau et al. , 2016)  resolves deep convection with a horizontal resolution of 1. 3 km. It is operational since the end of 2008 with a major update in 2015 concerning its resolution and a reduction of the period of the data assimilation cycle from 3 to 1 h. It is a limited-area NWP system with a geographical domain shown in Figure  2. AROME-France is computed on 90 vertical levels with a maximum altitude of 10 hPa but the levels are mostly concentrated in the troposphere. The model simulates twelve prognostic variables including two components of the horizontal wind (U and V ), the temperature T , specific content of water vapour q v and the surface pressure p. The ICE3 microphysical scheme used in AROME-France predicts 5 out of the twelve prognostic variables: specific contents of rain q r , snow q s , graupel q g (including different types of large rimed crystals like graupel, frozen drops and hail), cloud droplets q c and ice crystals q i.
Rimed Particle Column The rimed particle column (hereafter called RPC, in m) is described by Figueras i  Ventura et al. (2019)  as the difference between the upper limit of the highest level where rimed particles are predominant and the lower limit of the lowest level where those species are predominant. The data they used are from a lightning measurement campaign that took place in the Alps in 2017 where a LMA was deployed for the occasion. The RPC was retrieved out of a hydrometeor classification from radar data. They noticed an increase in lightning activity from a RPC thickness of 3 km onward and a high CG lightning activity when the RPC was larger than 8 km. For the purpose of our study, the levels where rimed particles are predominant are the levels where the specific content of graupel q g is higher than each specific content of the other hydrometeor variables.
One day (26 May) is used for a sensitivity study to the accumulation period (section 5). In the context of this sensitivity study, FEA with accumulation periods of 5, 10, 15, 20, 30 and 60 minutes were investigated. However, it is the 10 min FEA that was used to compare the proxies' performances (following section), pixel-to-pixel comparison is not performed and the data are stacked and rank-ordered over the whole time period and domain (i. e. dashed rectangle in Figure  2 ) so that they are treated as a global sorted distribution without time and space dependence. centered Then, all data points with proxy values equal to zero are disregarded. This classification allows to identify a proxy threshold X th corresponding to the first non-zero FEA value, that can be interpreted as the minimum required amount of the proxy quantity to observe lightning.
Finally, the data are normalized to the [0,1] range to perform the regressions. All the results presented in this paper are re-scaled to physical units. To relate each proxy individually to the FEA, five regression models were tested: an ordinary least squares linear regression (LinReg), a cubic polynomial regression (PolyReg3), a linear support vector regressor (LinSVR), a multi-layer perceptron (MLP) with 20 hidden layers and a random forest regressor with 20 decision trees (RF20). The depth of the trees was unbounded. All the regression models were implemented in Python using the Scikit-Learn package. For the linear regressions (LinReg and LinSVR), the established function f (proxy) is piecewise-defined so that the regression is performed only with the non-zero points where a and b are the coefficients of the linear regression curve.
This score is also computed on the validation set in order to measure the model's predicting ability on an independent dataset. The Fraction Skill Score (FSS;  Roberts and Lean, 2008 ) is also computed using the fast calculation method introduced by  Faggian et al. (2015)  to evaluate the displacement error of the modelled FEA values compared to the observed ones for the validation set. The FSS is a neighbourhood verification score where the spatial distribution of events is treated probabilistically which is particularly valuable for data with a small horizontal grid spacing. For fractions p o and p f of observed and forecast events, defined as the ratio of the number of pixels with values higher than a fixed threshold and the total number of pixels in a defined neighbourhood size (behaving as a sliding window), the FSS can be written as FSS = 1 - 1 N N i=1 (p f -p o ) 2 1 N N i=1 p 2 f + 1 N N i=1 p 2 o ( 12 ) where N is the total number of windows in the domain, depending on the size of this window.
Univariate models The results of the fitting on the training set for all the regression models are presented in Figure  6  and Figure  7  for the graupel mass and wmax, respectively (other proxies not shown here). Because data from the all the training days are concatenated and 290 sorted, scatter plots of FEA versus each proxy are monotonic increasing curves. The lower FEA values, below 100 flashes, are well fitted with all the regression models for the graupel mass and wmax, but their predominance in the dataset (see colour shades in Figures  6  and 7 ) weighs the regression curves of the LinReg and LinSVR models too heavily. In consequence, the higher FEA values are not well fitted by the LinReg and LinSVR models, resulting in lower R 2 scores. All the scores for  F1 2. 92 × 10 -2 8. 17 × 10 -4 1.
989. A good fit on the validation set is not enough to conclude on the ability of a proxy to predict lightning at the right position and time. In consequence, the FEA is computed for each hour for a given proxy of the validation dataset and compared to the observed FEA using FSS. As stated in section 2. 2, the proxies are calculated using 1-hour AROME-France forecasts from the assimilation cycle. Figure  9  shows a typical example of the observed and computed 10 minutes FEA fields using the RF20 regression model with the graupel mass, a microphysical proxy, for 9 August 2018 at 13 UTC. The minimum observed FEA value is 1 flash and the colormap is in log scale. The lightning activity observed in the North-East of Spain is well reproduced by the regression model, as well as the storms above South-East of France, although slightly shifted northwards because of the displacement error of the graupel mass forecast from the NWP system.
Multivariate models The objective here is to take advantage of the information contained in several proxies by combining them into a multivariate regression model. A problem arose here because of the high correlation between some of the proxies, called features in ML, as derived from the same fields from the AROME forecast. The multicollinearity among the features is demonstrated when trying to assess the importance of each feature using the Permutation Feature Importance method  (Breiman, 2001). The goal of this method is to determine how much the performance of a model relies on a single feature by calculating the difference between the R 2 score from the original dataset for a given regression model and the score when one feature of the dataset is randomly shuffled. If the shuffled feature is important, the score difference should be large.
However, when the Permutation Feature Importance is calculated on our training dataset, the score differences for the random forest regression model are all smaller than 0. 1 (not shown) which suggests that none of them are important. This is in contradiction with the results presented in section 3. 1 with high skill for most proxies, meaning that when a feature is corrupted, very similar information can be found in another one. Multicollinearity is a problem because it decreases model accuracy and robustness and increases overfitting  (Cohen et al. , 1983; Neter et al. , 1996; Chatterjee and Hadi, 2006). To avoid these issues, we went through a process of features selection using two different methods: i) hierarchical clustering on the Spearman correlation-rank ordered features to study their correlation and drop redundant variables, and ii) a principal component analysis (PCA).
e. on 9 August 2018 at 13 UTC. The general areal coverage is very similar to what was predicted with the graupel mass alone, especially in South of France and North-East of Spain, with some additional isolated patches that can be attributed to the contribution of wmax. However, they do not improve significantly the overall prediction compared to the graupel mass alone and this is also demonstrated with the FSS for the whole validation period shown in red in Figure  14. Other combinations were tested, for instance using the updraft volume and the IWP, for which the FSS is plotted as well in Figure  14 , but the performances were all quite similar. The PCA is a dimension reduction technique that creates new uncorrelated variables with linear combinations of each feature from the original dataset  (Jolliffe and Cadima, 2016).
The linear coefficients are called the principal component (PC) loadings and are calculated such as the PCs have successively a maximum variance. The proportion of total variance explained by each PC is often expressed as a percentage and only the first PCs that cumulatively explain 70 % of the total variance are kept (common cut-off point according to  Jolliffe and Cadima, 2016). In our case, the PCA is applied to the non-sorted standardized (mean equals to 0 and standard deviation to 1) array of stacked column-wise features. Analysis of the cumulative explained variance ratios (not shown) reveals that only the first two PCs are required to describe at least 70 % of the variance. The random forest regression algorithm is then fitted with the selected PCs, that were sorted beforehand. The same transformation to a dimensionally reduced dataset is applied to the validation dataset.
The FEA is then calculated for each hour of the validation. Figure  13b  depicts the FEA modelled using the first two PCs for 9 August 2018 13 UTC. It presents little to no difference with the FEA modelled using simultaneously graupel mass and wmax, highlighting the contribution from both microphysical and dynamical proxies. The FSS for multivariate models using either the first two PCs or the first five PCs are plotted in Figure  14  as well. Overall, they do not present a clear improvement compared to microphysical proxies alone. Ultimately, the F3 function presents a simple empirical way to forecast lightning density using both microphysical and dynamical properties of the cloud and was able to reproduce the areal coverage quite well overall for the validation period. However, the predicted F3-FEA values have too low amplitude whereas our method, where all data from the training dataset are used to fit the regression models, produces FEA with values closer to the observed ones.
The sensitivity of the relationship between a proxy and observed FEA to the accumulation period of the FEA has been evaluated. In the following, the graupel mass is the proxy selected to perform the study because it is one of the proxies presenting the strongest correlation with the FEA (R 2 > 0. 9 for all regression models on the validation set). Only the cubic polynomial regression model is considered for this sensitivity study. The method to calibrate the regression functions is quite similar to the one described in section 2. 4: graupel mass is extracted from the 1 h AROME forecasts from the assimilation cycle and FEA accumulated successively on 5, 10, 15, 20, 30 and 60 minutes are generated for every day of the training dataset (44 days). The graupel mass is not accumulated, it is an instantaneous data field, valid at time t.