Guihuan Feng

Zhengxing Sun

Christian Viard-Gaudin
email: christian.viard-gaudin@univ-nantes.fr

Hand-drawn Electric Circuit Diagram Understanding Using 2D Dynamic Programming

Keywords: Sketch understanding, electric circuit diagram, two-dimensional dynamic programming

A difficult task of a sketch understanding system is that it should always try to balance between the drawing freedom and the complexity of recognition. Most online existing works are based on the assumption that people will not start to draw a new symbol before the last one has been finished. As obviously it is not always the case, we propose in this paper a method which relaxes this constraint. The proposed methodology relies on a twodimensional dynamic programming technique (2D-DP) allowing symbol hypothesis generation, which can correctly locate symbols even when they are drawn temporally overlapped with each other. A neural network classifier, embedded in the 2D-DP, is used to label each hypothesis. Besides, tolerant connectivity constraints are introduced to help increasing the performance as well as the efficiency. Experiments show that this novel approach achieves an accuracy of more than 90 percent.

Introduction

It is a quite well-known proverb "a picture is worth a thousand words", which means that sketches or graphics are more expressive compared with only texts or words. Besides, sketches are widely used in engineering and architecture fields. This is mainly due to the fact that a sketch is a convenient tool to catch rough idea, so that the designers can focus more on the critical issues rather than on the intricate details [START_REF] Calhoun | Recognizing multi-stroke symbols[END_REF].

The emergence of Tablet PCs, digital pens and papers, and electronic whiteboards allows to record such handdrawn sketches as online information. Again, with the development of human-computer interaction, artificial intelligence and some other technologies, there is an increased interest in sketch understanding. Lots of fields are concerned with sketches. In this paper, we focus especially on hand-drawn electric circuit diagrams.

A difficult task of a sketch understanding system is that it should always try to balance between the drawing freedom and the complexity of recognition. Generally, the more freely the system can endure, the more difficult sketch understanding will be. As for the drawing issues, they can be divided into three different levels, i.e. sketchlevel variance, symbol-level variance and stroke-level variance. With respect to sketch-level, we mean symbols can be drawn in different order. Consider circuit diagrams for example, some would like to draw in the order of "connector-component-connector", while others prefer to finish all the components at first, and then complement with the connectors. At the same time, symbols can even be drawn temporally overlapped with each other, which means that sometimes people start a new symbol before he or she finishes the last one. Furthermore, different user can draw the same symbol very differently; even for the same user, the drawing style will vary from time to time. At last, not only one stroke can contain several symbols, but also one symbol can be made up of many strokes, which we call the multi-symbol strokes and the multi-stroke symbols. All these uncertainties have made sketch understanding a great challenge.

In order to solve the problems mentioned above, a successful sketch understanding system should at least have the following parts. First, stroke segmentation is used to divide strokes into smaller but more meaningful units, namely primitive segments; second, symbol hypothesis generation is included to try to group the segments in some specific way, so that each group can match to a symbol in a given domain; at last, a robust symbol recognition is needed so as to adapt to the great variability of different user drawing styles.

The paper is organized as follows. We first discuss the related works in Section 2. In Section 3, we give a brief introduction of our approach, and some more detailed information is illustrated in Section 4. In Section 5, we discuss about the performance of our approach. Finally, the conclusion is proposed in Section 6.

Related works

At the early stage of sketch understanding, symbols are segmented manually by pressing a button or pausing for a long time before the start of a new symbol. More efforts are focused on the recognition of isolated symbols [START_REF] Rubine | Specifying gestures by example[END_REF] [START_REF] Fonseca | Cali-an online scribble recognizer for calligraphic interfaces[END_REF].

Right now, researchers are working on the automatic parsing of the continuous input, so that people can express their ideas more freely. Structural or syntactic methods are the most frequently used approaches in sketch understanding. Such methods represent objects in terms of their structure, such as constitutional components and their relationships. Hammond and Davis [START_REF] Hammond | LADDER, a sketching language for user interface developers[END_REF] propose a sketching language LADDER, which can automatically generate a sketch interface based on the user defined description of shapes in a given domain. The limitation lies on that their approach can only describe regular shapes without too much detail. Our system performs well even when symbols are drawn with over-traced strokes. Costagliola et al. [START_REF] Costagliola | A parsing technique for sketch recognition system[END_REF] propose a Left-Right based parsing approach. The rank value of each candidate is computed by combining the accuracy of both stroke types and relations. So it has the disadvantage of the dependence on preprocessing. Our approach is inspired from their work to introduce tolerance in relationship measurement. But our symbol recognizer is more robust to noises and errors produced by segmentation.

Sim-U-Sketch is a sketch-based interface for Simulink package [START_REF] Kara | Hierarchical parsing and recognition of hand-sketched diagrams[END_REF]. The system provides a hierarchical "markgroup-recognize" sketch understanding architecture, and a domain-independent, trainable symbol recognizer that can learn new definitions from a single prototype. Their approach is not suitable for electric circuit diagram in that it is guided from marker symbols.

Gennari et al. [START_REF] Gennari | Combining geometry and domain knowledge to interpret hand-drawn diagrams[END_REF] employ ink density to enumerate candidate symbols, and domain knowledge is introduced to prune away invalid symbols. But their work is based on the assumption that the user finishes drawing one symbol before drawing a wire or another component, while our approach can correctly locate symbols even when they are drawn temporally overlapped with each other.

Segzin and David [START_REF] Sezgin | Temporal sketch recognition in interspersed drawings[END_REF] [13] make use of the temporal information and present a parsing technique based on multi-scale models. Compared with other methods, its main advantage is the high efficiency. And the disadvantage is that it depends too much on the drawing style.

Saund et al. [START_REF] Saund | Perceptual organization as a foundation for intelligent sketch editing[END_REF] present a system that uses Gestalt laws to locate salient objects. Their approach concerns only the grouping of strokes and no symbol recognition is included. It seems to be more suitable for the recognition of texts instead of diagrams.

System overview

The Intuition

Sketch is generally made up of some basic elements, namely symbols in a given domain. The problem is that while scribing, all the strokes are input continuously. People will not pause between symbols, or when pausing, it does not always mean the start of a new symbol. Besides, there is no any prior knowledge about the constitution of the sketch, which is the main difference between the recognition of a sketch and that of the isolated symbols. So we can not simply match the input with some predefined templates.

In terms of strokes, we mean the sequence of points that is recorded between a pen-down and pen-up. Due to the noisy and inaccuracy inherence, before parsing we first divide strokes into perceptually salient segments. A symbol hypothesis is a set of segments grouped together under some hypothesis generation approach. For each of these hypotheses, symbol recognition is performed, and a cost as well as a label is assigned. The bigger the cost is, the less likely the hypothesis can be a real symbol. Therefore, if we consider Γ to be all the possible ways to do segment grouping, then sketch understanding aims at finding the one at the minimum expense. It can be formulized as Eq. 1, where L is one of the possible parsing solutions.

arg min( | , )

L cost Sketch L ∈Γ ∑ (1)

Flow Chart

There are two kinds of strategies to do sketch understanding. One is the immediate feedback, which means once a stroke is drawn, sketch understanding will start. The advantage of such strategy is that user can view the recognition results in real-time. But it will, to some extent distract the user during the design task [START_REF] Hong | Sketch recognizers from the end-user's, the designer's, and the programmer's perspective[END_REF]. Again, due to the lack of complete drawing context, such methods always need to add constraints to the drawing style of some specific symbol. In this paper, we adopt another strategy, i.e. lazy feedback, which means sketch understanding starts only after the whole sketch has been finished. We believe that the global context can help to increase the performance and decrease the ambiguity, while at the same time do not interfere with the users. Furthermore, this strategy is well suited for digital pen and paper solution, where no immediate feedback is available for the user.

Figure 1 gives an illustration of the overall process of our sketch understanding technique. Raw strokes are firstly sent to the pre-processor to be divided into segments. Next, we use a two-dimensional dynamic programming (2D-DP) approach to manage segment grouping. But only those that can pass the validity test are valid hypotheses. The validity testifier is used to eliminate those groupings having other segments falling inside. This is based on the assumption that symbols will not be drawn spatially overlapped with each other. Later sketch understanding goes on two paths. We check first if the hypothesis can be recognized as a connector, and if the connectivity constraints of a connector are satisfied. For those that are not connectors, we will pass them to a Neural Network (NN) Classifier, and examine the connectivity according to the label given by the classifier. Besides recognition cost, there is also a connectivity cost to set to what extent the hypothesis fulfills the connectivity of a symbol. We use a tolerance to measure connectivity constraints, because we found the binary classification is prone to errors. Both the recognition cost and the connectivity cost are combined together to affect the decision of the 2D-DP.

Method in details 4.1. Stroke Segmentation

As the sequential coordinates are collected, they are segmented into lines and arcs using an HMM-based stroke fragmentation technique. We build a model that can represent an arbitrary stroke made up of lines and arcs, and segmentation is done through looking for the optimal path with the highest probability. As the model is domainindependent, it can be easily adapted to other areas. Besides, we take into account both the local and global information, so the approach can handle smooth curves properly. More details are available in [START_REF] Feng | HMM-based stroke fragmentation[END_REF].

Before generating the symbol hypotheses, we also introduced pre-grouping to help decreasing the searching space of 2D-DP. It is accomplished by merging small and overlapped segments.

Symbol Hypothesis Generation

We need to search on the segments to generate symbol hypotheses. If the sketch is of a large size, the searching space will be quite important, thus resulting in a great computational complexity. The characteristic of online sketch understanding is that besides geometrical information of the strokes, we also have the temporal information that can be utilized to simplify the problem. If symbols are not allowed to be drawn temporally overlapped with each other (as is stated in previous works [START_REF] Gennari | Combining geometry and domain knowledge to interpret hand-drawn diagrams[END_REF] [12] [START_REF] Sezgin | Sketch interpretation using multiscale models of temporal patterns[END_REF]), 1D-DP can solve the problem properly. But in practice, we found that it is not always the case. Therefore, we introduced a 2D-DP approach.

2D-DP is different from 1D-DP in that it can generate hypothesis with inconsecutive elements. With respect to the grouping of primitives, 1D-DP can only generate hypothesis like {S i , S i+1 , S i+2 ,…}, with i being a temporal index, while 2D-DP can also generate hypotheses like {S i , S i+k , S i+k+1 ,…} where k > 1. In this case, we say there is a time-jump. One main problem of using 2D-DP is the computational complexity. This is because the introduction of jumps will increase greatly the total number of hypotheses. Consider a sketch with N segments after preprocessing, the total number of hypotheses of 2D-DP will be 2 N -1. Hence, for each symbol hypothesis, we add the following restrictions:

(1) The maximum number of segments should be less than a fixed number. It has been set to 10 in the proposed experiments. This is based on the observation that most symbols in circuit diagrams are of simple structures. For example, many people draw a capacitor with only two lines. Resistor seems to be the one with the maximum number of segments, but we observed that 10 should be enough for most circumstances.

(2) Up to 2 jumps are allowed. A larger number of jumps would increase the searching space greatly. The example shown in Figure 2 is composed of a transistor symbol corresponding to the segment sequence (S 2 , S 3 , S 5 , S 7 , S 8 ), thus having two time jumps. With a pure 1D hypothesis generation scheme, it would not be possible to recover such a symbol. (3) The distance between segments should be less than a threshold. Since if two segments are too far away from one another, it is less likely that they belong to the same symbol. As for implementation, this threshold is dynamically computed according to the sketch. Consider a sketch with 50 segments, without any constraints, 2D-DP will generate O(10 15 ) hypotheses. Under constraint (1), the total number of hypothesis can be decreased to C(50,10)=50!/(10!*(50-10)!)=O(10 7 ); the same result can be achieved with constraint (2). In terms of constraint (3), it depends. If segments of a symbol are all close to each other in a given domain, this constraint will help a lot. Furthermore, we can also restrict the total number of hypotheses directly.

Symbol Recognition

Connectors and circuit components are handled differently. A group of segments is a connector if all the following constraints are fulfilled.

Number of segments should be less than 3. We found that people seldom draw a connector in many segments. Also, this requirement can help to avoid compressed resistors and inductors to be misrecognized as connectors . The minimum square error to fit a line should not exceed a threshold. There should not be a big gap inside the approximated result, to prevent that the two wires connected with both sides of a capacitor are merged together. The connectivity constraint of a connector should be satisfied.

Indications that a symbol hypothesis is a circuit component include the following: It does not meet the requirements of a connector. The NN classifier can match the group to a component with a relative high recognition score. In our system, only those with a normalized recognition score greater than 0.65 are taken as actual components, outputs of the NN being normalized with a softmax function between 0 and 1. The hypothesis has fulfilled the connectivity constraint of the component labeled by the NN classifier.

Our approach differs from [START_REF] Gennari | Combining geometry and domain knowledge to interpret hand-drawn diagrams[END_REF] in both the symbol hypothesis generation strategy and the introduction of tolerant connectivity constraint. Connectivity is always used accessorily as error correction measurements, where only some heuristic rules are presented. But here as most hypotheses are invalid groupings and NN classifier has poor performance for these outlier samples, so connectivity constraint is an important measurement to decrease errors, as well as increase the efficiency of the 2D-DP.

Connectivity Constraints

NN classifier can tell to how much extent the input can match with one of the predefined symbols, but has less capability to reject outliers. Also, the minimum square approximation can figure out whether the input can be a perfect line, but cannot distinguish a line, which would be a part of a symbol, from a real connector.

In this paper, we introduce the connectivity constraint measurement to help to solve this problem. In terms of connectors, we assume that there should be at least two symbols connected with both sides. These symbols can be either connectors or circuit components. While as for components, based on their types and orientations, we have defined 13 different connectivity constraints, shown in Table 1. Here, all the constraints are concerning connectors. If we consider for example constraint of type 3, we mean there should be at least two horizontal connectors connected at both left and right side of horizontal capacitors, where a horizontal capacitor is a symbol with 2 parallel vertical line segments. We separate type 11 from type 3, because as for capacitors, people generally draw two connectors connected at the middle of its bounding box. Hence during computation, we calculate the distance from the endpoint of the connector to the middle point of the bounding box of the symbol; while people usually draw resistors or inductors with connectors in a single stroke, therefore the connectivity constraint is not as strict as for those capacitors.

The computation of the connectivity cost is shown in Eq. 2, where n is the number of connectivity needed, d i is the i th minimum distance according to the rules, and we use a scalar s to normalize all these distances. This s value is computed during preprocessing of the global sketch from the histogram of the distances between consecutive extrema points in both x and y.

n i i d s Cost Connectivity n = = ⎛ ⎞ ∑ ⎜ ⎟ ⎝ ⎠ (2)
This connectivity cost will be combined with a recognition cost to define the global cost as introduced in Eq. 3. When dealing with connector hypothesis -left branch of flow chart presented in Figure 1, the recognition cost is based on the least square approximation error of the considered points with a straight line. Negative costs have been considered for perfect connectors, so that 2D-DP can get the solution with as many connectors as possible. Concerning component hypothesis -right branch of Figure 1, as the value of the winning output of the NN classifier is not so reliable, the corresponding recognition cost will be set to 0, meaning that only the connectivity constraints will guide the 2D-DP, however they are computed based on the label proposed by the symbol recognizer. The ultimate cost for each hypothesis is calculated based on Eq. 3. Here, w is a weighting factor, which has been set experimentally. In our experiment, w is selected as 1.9. _ _ Cost Cost Recognition w Cost Connectivity = + × (3)

System evaluation

NN Classifier Training

To train the classifier, several issues have been considered, so the recognition engine can take into account whether or not a shape can vary according to scale, flip and rotation. In the experiments reported here, we have defined 9 electrical circuit components shown in Table 2. Further, they are divided into 26 classes according to the orientations and different connectivity requirements. A training database of such isolated symbols has been collected. It involved 11 people, who have drawn each 4 to 5 symbols of each type, with more or less application. Finally, there are a total of 451 isolated symbols used to train the classifier.

Resistor

Inductor Capacitor

AC voltage DC voltage Current

Earth ground Chassis ground Transistor

Diagram collection

To test the performance of the global system, we asked 10 subjects to draw a series of electric circuits. Note that only one of them is also involved in sample collection of the NN classifier training. Samples are collected using the Anoto digital pen and paper, so that they can draw as freely as usual. Each of them is asked to copy 10 diagrams and the information is stored as a series of 2D coordinates. All the schemes are selected from two electrical textbooks, and the subjects are not told about how the system works. As the samples are drawn on the paper, modification is not allowed. However, if they are aware of any big mistakes, they will have to redraw a new one.

It is very interesting to notice that people draw quite differently. Some would like to start from left to right, and others prefer to draw from top to bottom. Although most participants draw symbol by symbol, they do jump from time to time. Among the 100 diagrams we have collected, 31 sketches have time-jumps inside, where 15 are caused by transistors. Sometimes, people even jump more then 2 times when scribing a symbol.

During sample collection, we found that due to the problem of the sampling devices partial information has been lost. Also, there are some noisy strokes while people are pausing. All these make it a real big challenge to traditional approaches. In addition to the data collection process, we have manually segmented and labeled the sketches to define the ground truth.

Measurements

The accuracy of a specific class is defined as the average of the precisions of all instances, so as to make it to be size independent. Similarly, the accuracy of each user can be defined in the same way. The precision of one symbol instance is computed at point-level, namely the number of points that are correctly labeled divided by the number of points labeled in the ground truth, shown in Eq. 4. Figure 3 provides an illustration. Here, some part of a resistor is mis-classified as vertical connector (strokes in yellow). Therefore, the precision of the resistor is computed as the number of points located in red divided by the number of points located both in yellow and in red. After getting precisions of all the instances, precision of the class is achieved as Eq.

Performance Evaluation

The recognition results are presented in Figure 4. We have achieved an average symbol-level accuracy of more than 90% (See Figure 4a). Only for capacitors, transistors and ac-voltages, the precisions are less than 90%. The performance varies between users (See Figure 4b), and we found this is also due to the mis-recognition of capacitors and ac-voltages. In terms of ac-voltage, errors are mostly included by the NN classifier, as it is likely to misrecognize ac-voltages to currents. As for capacitors, if the user tends to draw the two poles too close, our connectivity constraints will fail to distinguish it from connectors. Out of the 100 samples, 17 are left unrecognized, among which 2 are due to stroke segmentation errors. Since if we mis-segment a connector, there will be no required connectivity constraints, and the 2D-DP will fail. Most failures are owing to the scribing errors, such as forgetting to draw a stroke of a symbol and etc. In such circumstance, there will be an error result provided, but they are not included in the performance evaluation. Also, there are 2 schemes where the users do not respect the standard connectivity constraint of transistors, which is another cause of the recognition failure. We argue that this is based on the consideration of the efficiency instead of the limitation of our approach.

Conclusion

In this paper, we present a sketch understanding technique that can automatically extract symbols from the continuous stream of strokes. Unlike the previous works, we do not insist that symbols should be drawn one by one. The 2D-DP based symbol hypothesis generation approach can correctly locate symbols even when they are drawn temporally overlapped with each other. The introduction of tolerant connectivity measurement helps to increase the performance, and avoid errors that might be included due to the lack of reject capability of NN classifier to the outliers. The preliminary experiment has achieved an accuracy of more than 90 percent.

In the future, we plan to solve problems when there is no result provide by the 2D-DP and conduct a more convictive evaluation of the proposed approach.

Figure 1 .

 1 Figure 1. Flow chart of the understanding system.

Figure 2 .

 2 Figure 2. An example of how people draw a transistor with jumps.

Figure 3 .

 3 Figure 3. An illustration of the computation of symbollevel precision.

  capa tran acvo dcvo curr eart chas (a) Symbol-level precision. User-level precision.

Figure 4 .

 4 Figure 4. Recognition results.

Table 1 .

 1 Type of component connectivity constraints.

	No. Component type	Constraints
	1	groundings to the left	left only
	2	groundings to the right	right only
	3	h_capacitors	left & right (s)
	4	groundings to the top	top only
	5	groundings to the bottom	bottom only
	6	v_capacitors	top & bottom (s)
	7	transistors to the top	top, left & right
	8	transistors to the bottom	bottom, left & right
	9	transistors to the left	left, top & bottom
	10	transistors to the right	right, top & bottom
	11	h_resistors or inductors	left & right (l)
	12	v_resistors or inductors	top & bottom (l)
	13	alternating current	(left & right) or
		voltages	(top & bottom)

Table 2 .

 2 Electric circuit components used in our experiment.

Acknowledgement

We thank the anonymous reviewers for their valuable feedback.

This research is jointly supported by French ANR grant CIEL 06-TLOG-009, and the Program for New Century Excellent Talents in University of China (Project No.NCET-04-0460)


Nazzareno Diodato

Dr Iñigo Gómara
email: i.gomara@ucm.es

Gianni Bellocchi

Modelling snowfall in southern Italy: A historical perspective in the Benevento Valley (1645-2018)

Keywords: Climate Variability and Change, Historical reconstruction, Modelling, Snowfall, Southern Italy

The lack of long-term, homogeneous snowfall records is a limitation in environmental studies. Statistical modelling holds potential to extend back in time snowfall records with a limited set of predictors: snow severity and winter-spring temperatures (with their variability) to reflect elevation influences. The annual number of snow days (SDY) in the Benevento Valley (southern Italy) was detailed for the period 1870-2018. Calibrated in the period 1870-1968 (R 2 =0.85) and validated in the period 1969-2018 (R 2 =0.67), the model developed here enabled the reconstruction of a time-series of SDY between 1645 and 2018. This unique series (the longest in southern Italy) shows that SDY peaked during the Little Ice Age (until ~1850), dominated by cold air masses or characterized by winter seasons extending until May (1655May ( , 1684May ( , 1763May ( and 1830) ) or June (1620). After the change-point detected in 1866, the modelled SDY time-series declined rapidly (Modern Warming Period, 1867-2018). The atmospheric conditions that favoured snowfall in the Benevento Valley throughout the study period were generally associated with an anomalous highpressure system located over northern-north-western Europe and a low in the eastern Mediterranean. This configuration allowed the incursion of cold continental air from the east-northeast into southern Italy. Our results are consistent with similar studies of snowfall in other European and mid-latitude regions of the northern hemisphere.

INTRODUCTION

In hydrological studies, spatial and temporal heterogeneity is a main challenge for observing, understanding and modelling snowfall, snow accumulation and snowmelt, which are among the processes that have the greatest impact on the global water cycle [START_REF] Webster | Snow in the changing seaice systems[END_REF]. Snowfall occurrence in the northern Hemisphere indicates that winter periods in western Europe (including Italy) are characterized by only a few days of snow per year [START_REF] Jennings | Spatial variation of the rainsnow temperature threshold across the Northern Hemisphere[END_REF]. Different aspects of extreme climate (e.g., extreme cold and snow) are of interest in understanding regional climate characteristics because their changes have impacts for water resources, ecosystems and society, with far-reaching environmental and socio-economic implications [START_REF] Braconnot | Past and future contexts for climate and water-cycle variability, and consequences for the biosphere[END_REF].

To support global water-cycle modelling, the number of snow days per year (days with snow depth equal or greater than 1 cm) is an important climatological indicator (WMO 2009) of winter-snow storage, especially in areas where snow falls infrequently. Snowfall delays the transfer of precipitation to surface runoff, infiltration and streamflow generation.

In this way, it affects the timing and magnitude of peak flows [START_REF] Wang | Estimating Snow Mass and Peak River Flows for the Mackenzie River Basin Using GRACE Satellite Observations[END_REF], the recession of hydrographs [START_REF] Yarnell | Ecology and Management of the Spring Snowmelt Recession[END_REF], the onset of snowmelt-driven streamflow [START_REF] Grundstein | Trends in Average Snow Depth Across the Western United States[END_REF], and the magnitude and duration of summer baseflow [START_REF] Godsey | Effects of changes in winter snowpacks on summer low flows: case studies in the Sierra Nevada, California, USA[END_REF]. However, defining these impacts is difficult and they have only been explored on a large scale. For instance, the variability of autumn snowpack in Eurasia can drive the atmospheric circulation in the northern Hemisphere [START_REF] Henderson | Snow-atmosphere coupling in the Northern Hemisphere[END_REF].

Ongoing snow measuring networks are not yet well established at regional and local scales, especially in remote regions [START_REF] Dewalle | Principles of Snow Hydrology[END_REF]. In the absence of continuous and reliable data, consultation of historical press reports can help to extend and improve the time-series of snowfall events, but characterisation of the temporal and spatial distribution of the snowfall and the associated synoptic patterns remains limited to a few decades [START_REF] Martínez-Ibarra | Reconstruction and analysis of 1900-2017 snowfall events on the southeast coast of Spain[END_REF]. Regional investigations using satellite data focus mainly on the extent of snow cover on a global and continental scale, and only refer to the last few decades, during which several satellite missions have been launched or planned [START_REF] Capozzi | Retrieval of snow precipitation rate from polarimetric X-band radar measurements in Southern Italy Apennine mountains[END_REF]. Thus, the understanding of snowfall patterns in several world regions remains limited due to the brevity of the record and poor knowledge of the hydroclimatic mechanisms that control snowfall dynamics.

In Europe, the occurrence of snowfall and the number of snowy days are mainly dependent on latitude and the large-scale atmospheric circulation [START_REF] Croce | The snow load in Europe and the climate change[END_REF], even if snowfall frequency and intensity do not show clear spatial patterns [START_REF] Navarro-Serrano | Spatio-temporal analysis of snowfall events in the spanish Pyrenees and their relationship to atmospheric circulation[END_REF]. In Italy, snow is common in all mountains, occasionally falling at low elevations and not unfrequently in the northern plains, along the Adriatic coast (east side) and in hinterland areas of the central and southern regions. However, during the centuries of the Little Ice Age (LIA, ~1300-1850 CE; [START_REF] Miller | Abrupt onset of the Little Ice Age triggered by volcanism and sustained by sea-ice/ocean feedbacks[END_REF], snowfall was more abundant throughout Italy. During the LIA, the climate was highly variable in Europe, and winters were characterized by very cold and snowy weather. This period coincides with an increase in atmospheric radiocarbon [START_REF] Stuiver | Sun, ocean, climate and atmospheric 14CO2 : an evaluation of causal and spectral relationships[END_REF], several intense volcanic eruptions [START_REF] Briffa | Influence of volcanic eruptions on Northern Hemisphere summer temperature over the past 600 years[END_REF]) and a decrease in solar activity with a low sunspot number [START_REF] Lean | Reconstruction of solar irradiance since 1610: Implications for climate change[END_REF]. The snow-capped volcano Vesuvius (1281 m a.s.l.), depicted in the 19 th century anonymous painting in Fig. 1, shows that in the past there were cold snaps and severe snowfall also in places, such as southern Italy, where the climate is generally mild and temperate today.

Studies aimed at collecting, reconstructing and analysing snowfall episodes began in

Italy in 1681 with the publication of La Figura della neve ("The shape of snow") by the Italian natural philosopher and Catholic priest Donato Rossetti (1633-1686), who described the configuration of snow crystals seen through a microscope. Italy is home of the longest series of daily snow records in the world, starting in Rome in 1741 [START_REF] Mangianti | La neve a Roma dal 1741 al 1990[END_REF] and in Turin in 1788 [START_REF] Mercalli | Il clima di Torino. Società Meteorologica Subalpina[END_REF]. [START_REF] Enzi | Snowfall time-series reconstruction in Italy over the last 300 years[END_REF] provided a historical perspective of snowfall occurrence in Italy back to the LIA. More recent studies showed that distinct climatic patterns played an important role in driving snowfall fluctuations over centennial time-scales [START_REF] Diodato | Mediterranean winter snowfall variability over the past millennium[END_REF], and that temperature and precipitation changes are dominant controlling factors of hydrological changes (Diodato & Bellocchi 2020). However, no specific assessment of long-term snowfall dynamics in southern Italy (below the 42 nd parallel) was reported in previous studies.

Here, we present an annually resolved reconstruction of snow days per year (SDY) in the period 1645-2018 for the Benevento Valley (BNV). Extending southwest of Benevento (41° 08' N, 14° 47' E), the BNV hosts the longest time-series of observed snowfall in southern Italy. Benevento city observatory (BNOBS) and the Met European Research Observatory (MetEROBS), 5-km apart within the BNV, offer a unique heritage of snowfall data to understand and interpret the snowfall regime of the prevalently Mediterranean climate that characterizes this inland valley. This was done with a physically consistent statistical model, hereafter referred to as NLMRM (Non-Linear Multivariate Regression Model). First, we developed a parsimonious model for the reconstruction of the SDY data from climate anomalies related to snowfall severity, mean winter and spring air temperatures, and their standard deviations (Section 2). Then, we used the model to capture climate variability at multiple timescales and to identify patterns of change in snowfall in the BNV over the period 1645 to 2018 and sub-periods (Section 3). We discussed the analysis of historical snowfall in the BNV (Section 4) and concluded the article with a summary of results and future research directions (Section 5).

MATERIALS AND METHODS

Environmental setting

The southern peninsular part of Italy stretches from the area south of Rome (below 42° N) to the most southerly provinces of Apulia (around the 40 th parallel) and Calabria (around the 38 th parallel). The Benevento Valley is located in the Campania Region, within the southern Apennine range (Figs. 2A,B). Centred around 41° 11' N and 14° 27' E, and with elevations ranging from over about 100 to about 300 m a.s.l., the BNV has an area of about 100 km 2 (Fig. 2C) at the transition between the central and southern Apennines. The Benevento Valley has a typical hot summer Mediterranean climate (Köppen Csa) with more continental features than the maritime (sub-tropical) climate of the Tyrrhenian coast (e.g., around the volcano Vesuvius). According to De [START_REF] Renzi | Osservazioni sulla topografia-medica del Regno di Napoli, P. II. Tipografia Fratelli Criscuolo[END_REF], The western and southern winds have little influence over this province, less so than the mountains that surround it... The winds of the north-east, and of the north-west, are those which blow there with the greatest impetus. ... But there are places where the west and south winds blow furiously, after crossing the gorges of the mountains, where they have gained strength.

Synoptic large-scale configurations favourable for snowfall and snow persistence in the BNV are generally associated with the negative phases of the North Atlantic Oscillation (NAO), a see-saw pattern of mean sea-level pressure between the sub-polar and sub-tropical latitudes of the North Atlantic [START_REF] Pinto | Past and recent changes in the North Atlantic oscillation[END_REF]. Under a negative NAO phase (NAO-), a meandering and relatively weak extratropical jet stream favours the incursion of cold Siberian or Arctic air over central and southern Europe. For the former to occur, a low-pressure system is typically located over southern Italy, associated with an anticyclone over Scandinavia [START_REF] Fazzini | Snowfall analysis over peninsular Italy in relationship to the different types of synoptic circulation: First results[END_REF]. For the incursion of Arctic air, the anticyclone must be located further upstream in the North Atlantic, favouring the advection of air from the north instead of the north-east (D' Errico et al. 2020). Another recurrent configuration is the placement of a low-pressure system south of Italy (over the Ionian Sea) and an anticyclone over northwest Europe [START_REF] Dafis | The atmospheric circulation characteristics favouring snowfall in an area with complex relief in Northwestern Greece[END_REF], bringing over the BNV cold air from the Balkan Peninsula that moistens after crossing the Adriatic Sea. Local orographic features also favour occasional snowfall and snow cover in the BNV, as cold air can be trapped in the area due to vertical thermal inversion. Seasonal and monthly patterns in the BNV (Fig. 3) indicate that snow falls mainly from December to March (rarely in November and April).

While mean snow depths in the winter months are rather contained to a few millimetres, the monthly 90 th percentiles are expected to be over 40 cm.

Observatories and data resources

Weather data in the BNV were mostly derived from the observatories of Benevento (BNOBS) and Monte Pino (MetEROBS), at 180 and 170 m a.s.l., respectively (Fig. 2C).

We refer here to the BNV homogenized time-series of observed SDY, covering the period from the winter 1869-1870 to 2017-2018 (Diodato & Bellocchi 2020). This time-series was created from the BNOBS and MetEROBS data, and additional data from the former Servizio Idrografico e Mareografico Nazionale (Italian Hydrographic and Tidal Service;SIMN, 1991SIMN, -1997) ) observatory, which was about 400-m away from the BNOBS. The latter started its records in 1870, under the guidance of Nicola Orazio Albino, and stopped them in 1968. The SDY series at MetEROBS is continuous from its start in 1985 until today [START_REF] Diodato | Paesaggi d'Inverno: Aspetti naturalistici e climatologici delle nevicate sulla Campania interna[END_REF]. Before 1985, the BNV time-series of observed SDY data is divided into five main periods: (i) 1870-1886, including the observations made by Ambrogio Di Renzo;

(ii) 1887-1910, with the observations of Venanzio Vari;(iii) 1911-1949, reconstructed by Diodato and Bellocchi (2020);(iv) 1950-1968, with the measurements of Nazario [START_REF] Doretti | L'Osservatorio Meteorologico di Benevento nel suo ottantennio[END_REF]and(v) 1969-1984, which embraces the observations of the former SIMN. SIMN snow observations before 1969 are available with discontinuities and were not used. The SIMN time-series (1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984) is the one that presents a relative lower reliability than BNOBS and MetEROBS because the station was not always manned, and snow collection did not systematically take place. BNOBS and MetEROBS were monitored 24 hours a day and the snowfall data were validated at the observation time. The BNOBS series would have been the longest if some data had not been lost during the World War II, when regular measurements stopped until 1949 [START_REF] Doretti | L'Osservatorio Meteorologico di Benevento nel suo ottantennio[END_REF]. For the period of missing data , we relied on the data modelled by Diodato and Bellocchi (2020) based on SDY data from neighbouring sites. The SDY observations of the most recent period have been taken directly by the authors at MetEROBS, which is the only operational weather observatory delivering continuous snow measurements since 1985 in the BNV. In fact, the operational monitoring networks of the Centro Funzionale Multirischi Protezione Civile della Regione Campania (Multirisk Functional Centre of Civil Protection-Campania Region; http://centrofunzionale.regione.campania.it/#/pages/dashboard), which has replaced the SIMN, were not set up to perform snow observations. The Mann-Whitney-Pettitt test [START_REF] Pettitt | A Non-Parametric Approach to the Change-Point Problem[END_REF], widely used to assess the homogeneity of the distribution of hydro-meteorological variables [START_REF] Xie | Exploring the ability of the Pettitt method for detecting change point by Monte Carlo simulation[END_REF], did not detect any abrupt changes (p > 0.05) in the distribution of the homogenized SDY data. This supports our assumption that limited changes in geographic location and elevation (and observer changes) do not critically bias the unique realization of the joint time-series (with no overlapping periods for observation site data).

Gridded climatic datasets

The following gridded climatic datasets were used in this study:

• The European seasonal reconstruction of air temperatures: Temperature data for the period 1645-2002, with a horizontal resolution of 0.5° [START_REF] Luterbacher | European Seasonal and Annual Temperature Variability, Trends, and Extremes Since 1500[END_REF], were considered to reconstruct SDY back to 1645 (interpolated to the nearest BNV grid-point). We used seasonal winter (December to February) and spring (March to May) temperatures. For the exceptionally cold conditions of 1684 (with an expected return period of ~1000 years), we used the seasonally adjusted mean temperatures derived by [START_REF] Diodato | Discovering the anomalously cold Mediterranean winters during the Maunder minimum[END_REF] for southern Italy, i.e., 0.1 and 7 °C for winter and spring respectively (instead of 4.1 and 11 °C from the regional reconstruction).

Once the SDY time-series was reconstructed, its trends and variability where also analysed based on this gridded dataset (Section 3.3).

• The European monthly reconstruction of mean sea-level pressure (mslp):

The available period 1659-1999 of surface pressure reduced to mean sea-level data (hPa) at 5º horizontal resolution [START_REF] Luterbacher | Reconstruction of sea level pressure fields over the Eastern North Atlantic and Europe back to 1500[END_REF] was considered to analyse trends and variability of SDY.

• The Climatic Research Unit (CRU) time-series of variations in climate with

variations in other phenomena v3: Both seasonal temperature (°C) and mslp (hPa) abovementioned reconstruction datasets [START_REF] Jones | Climatic Research Unit (CRU): Time-series (TS) datasets of variations in climate with variations in other phenomena v3[END_REF] were updated from the Climatic Research Unit (CRU) global climate data until 2018.

• The Last Millennium Reanalysis (LMR) Project, Global Climate Reconstructions Version 2: Sea-surface temperatures (°C) at 2º horizontal resolution were retrieved globally as annual ensemble means to evaluate SDY trends and variability. The period considered was 1645-2000 [START_REF] Tardif | Last Millennium Reanalysis with an expanded proxy database and seasonal proxy modeling[END_REF]).

• The sea-surface temperature (SST) and sea ice data from the Met Office Hadley Centre (HadISST): The global monthly 1° sea-surface temperature (°C) dataset from 1870 onwards [START_REF] Rayner | Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century[END_REF] was also analysed for the period 1870-2018 to evaluate SDY trends and variability.

Climate indices

The following climate indices were considered in this study:

• Snow-Severity Index: As a covariate for SDY modelling, we referred to the snowseverity index (SSI) developed by [START_REF] Diodato | Mediterranean winter snowfall variability over the past millennium[END_REF]. It establishes five classes of snow severity, discernible in narrative texts drawn from historical documentation, from SSI = 0, which identifies winters without comments on their severity or their impacts on the society and economy to SSI = 4, which records extraordinarily snowy winter seasons with a low recurrence rate over many years.

• The Pacific Decadal Oscillation (PDO) index: The PDO is a recurring sea surface temperature (SST) anomaly pattern in the North Pacific with ⁓40-yr periodicity (warm/cold phases) and reported global impacts in climate and marine ecosystems [START_REF] Mantua | The Pacific Decadal Oscillation[END_REF]. The PDO index available from the LMR reanalysis (1645-2000;[START_REF] Tardif | Last Millennium Reanalysis with an expanded proxy database and seasonal proxy modeling[END_REF]) was utilized to assess SDY variability.

• The Atlantic Multidecadal Oscillation (AMO) index: The AMO refers to a warming/cooling pattern of multidecadal North Atlantic SST anomalies with ⁓60-yr periodicity [START_REF] Knight | Climate impacts of the Atlantic Multidecadal Oscillation[END_REF]. The AMO index also available from LMR was considered for SDY variability analysis.

• The North Atlantic Oscillation (NAO) index: As described in the introduction, the NAO is a see-saw pattern of mean sea-level pressure that characterises the strength and location of the westerly flow in the North Atlantic. Additional to the NAO index available from LMR , the winter NAO index reconstruction from [START_REF] Trouet | Persistent Positive North Atlantic Oscillation Mode Dominated the Medieval Climate Anomaly[END_REF] for the period 1645-1995 and the NAO Gibraltar-Stykkisholmur index from [START_REF] Jones | Extension to the North Atlantic oscillation using early instrumental pressure observations from Gibraltar and south-west Iceland[END_REF] were also considered for comparison on SDY variability assessment.

Development and assessment of the statistical model

In order to minimize uncertainties arising from SDY data estimation, we developed a non-linear multivariate regression model with three predictors and five parameters.

Following the principle of parsimony [START_REF] Mulligan | Modelling and Model Building[END_REF], for which parsimonious modelling has greater explanatory prediction power and lower complexity, we constructed a parsimonious model that could be easily parameterized and validated.

The model was sufficiently reliable to be used to reconstruct long-term historical SDY. To increase skill prediction and reduce uncertainties in the modelling estimation of annually resolved SDY data, we performed a multivariate nonlinear regression with four input variables and five parameters considering the original observational period 1870-2018. This original SDY data (149 years) was segregated into two groups: about two-thirds (99 years) were used for calibration ) and the rest (50 years) for validation . In this way, the calibration work was brought into line with the reality of past snowfall to allow a reasonably accurate reconstruction of historical snow days. As a first step, to minimize the number of inputs, we evaluated the effect of individual variables, or groups of variables, on SDY during climatically inhomogeneous but effective periods to take into account different situations. An iterative process (trial-and-error to compose the relevant drivers) subsequently allowed us to explain the dynamics of SDY in relatively simple way. A multi-stage selection logic was used to alternate between term addition and deletion and derive the following non-linear multivariate regression model to estimate SDY (d yr -1 ) at BNV:

SDY(BNV) = A • [(φ + SSI) η + (β -Tw -Ts) σ • VC(Tw → s) 𝑡=-1 𝑡=0 ] (1)
where: A (d yr -1 ) is a scaling parameter relative to the number of SDY when the term in brackets is equal to one, φ and β (°C) are position parameters; η and σ are shape parameters.

The first term of Eq. ( 1) includes the snow-severity index (SSI), which contains a qualitative information relating snowfall anomalies at local scale. In addition to this qualitative term, winter (Tw) and spring (Ts) temperatures indicate that changes in SDY are temperature-driven, mostly associated with long-distance transport of warm or cold air masses over the European continent. In fact, a rise in temperature will lead to more precipitation in the form of rain and more snowmelt [START_REF] Croce | The snow load in Europe and the climate change[END_REF]). The last term in Eq. ( 1) represents the temperature variability, expressed as the coefficient of variation (VC) calculated as the ratio between the standard deviation and the mean of the temperature values from December to May (including winter and spring temperatures, Tw and Ts, of the previous and current year, t=-1 and t=0). This term is a key factor in inferring seasonal temperature fluctuations, as it is higher in the presence of snow (Diodato et al. 2020a).

The is mainly controlled by winter precipitation, temperature control is dominant in transient snow regions [START_REF] Clark | Atmospheric controls on Eurasian snow extent[END_REF].

The parameters from Eq. ( 1) were calibrated with the SDY observed at BNV according to the criteria below:

{ R 2 = max MAE = min |𝑏 -1| = min (2)
The first condition maximizes the determination coefficient (0 ≤ R 2 ≤1, optimum), i.e., the variance explained by the model. The second condition minimizes the distance between modelled versus observed SDY data, reducing the mean absolute error (optimum, 0 ≤ MAE < ∞), in d yr -1 . The last condition approximates the unit slope of the straight line of the linear regression of the actual versus modelled data (b = 1, optimum). In addition, the Kling-Gupta index (-∞ < KGE ≤ 1) was considered as a measure of efficiency [START_REF] Kling | Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios[END_REF], where KGE > -0.41 indicates that a model performs better than the mean of the observations as a benchmark predictor. To select the set of relevant covariates for the parsimonious modelling of SDY, we iteratively added in predictors (one at a time) until small MAE and large R 2 values were obtained. The third criterion then determined the final selection. Convergence was obtained by repositioning each predictor in > 50 iterations.

Analysis of variance (ANOVA) was then performed to find out whether predictors were all necessary (and not redundant) for the modelling approach. The Durbin-Watson (DW) test was used to assess autocorrelation in model residuals.

RESULTS

Model parameterization and evaluation

With the calibrated parameters of Eq. ( 1 In the validation phase, R 2 (with a = -0.50 ± 0.49 and b = 1.12 ± 0.16) indicates that the model explains 67% of the total variance of the observations. MAE equal to 0.92 d yr - 1 is lower than the standard error of the estimates (1.16 d yr -1 ), while the Kling-Gupta index is -0.18. Again, there is no significant serial autocorrelation in the residuals (DW = 2.27, p > 0.05). The model's predictive capability may have decreased with the generally more sporadic snowfall events occurred in recent decades, which appears to be a southern Europe-wide process [START_REF] Diodato | New insights into the world's longest series of monthly snowfall (Parma, Northern Italy, 1777-2018)[END_REF]). In fact, there is a mismatch from 1978 onwards in the validation period timelines showing the co-evolution of SDY between observations and estimates (Fig. 6A). However, the associated distribution of the model residuals (Fig. 6B) does not significantly differ (p > 0.05) from what would be expected for a normally distributed sample.

In determining whether the SDY (BNV) model can be simplified, we have fitted a multiple linear regression model between SDY and three independent inputs: (0.25 + SSI) 1.4 , (28.5 -Tw -Ts) 2 and VC(T(w → s) t=-1 t=0 ) of Eq. ( 1). All of the variance components in the model are significant (the highest p-value being 0.0061 for the term VC) and the model cannot be simplified further.

Reconstruction of snow days per year in the Benevento Valley

Eq. ( 1) was used to reconstruct the evolution of SDY over the period 1645-2018 (Fig. 7A). Each year includes the winter season started on December of the previous year.

The reconstructed data were then analysed to find out possible trends and/or climatic oscillations explaining most of the observed variance in the long-term time-series of SDY, and to compare contemporary snowfall with historical snowfall. Overall, the estimated SDY values present a significant decreasing trend (p < 0.05) for the total period (Fig. 7A).

However, this decrease was not gradual and dominated during the 20 th and early 21 st centuries. The wavelet power spectrum of the time-series also shows temporal inhomogeneity in the SDY (Fig. 7B), with significant periodicities detected mainly at the beginning of the study period (90% confidence interval).

Consequently, we distinctly analysed SDY during three climatic sub-periods: (i) the Maunder Minimum (MM, 1645-1715), as it marks the end of a period of rapid cooling in 1600s; (ii) the Late Little Ice Age (LLIA, 1716(LLIA, -1866)); and (iii) the Modern Warming Period (MWP, 1867(MWP, -2018)). The application of test statistics suggests the existence of a change-point in 1866 [START_REF] Buishand | Some methods for testing the homogeneity of rainfall records[END_REF] or 1910 [START_REF] Pettitt | A Non-Parametric Approach to the Change-Point Problem[END_REF], which can be seen as a transition phase indicator from the LLIA to the MWP. Change detection suggests however that 1866 is a main breakpoint as it detects a structural change owing to the exit from the LLIA (~1850).

Table 1 shows that the mean and median values of SDY (modelled series) undergo a jump between MM and LLIA (6.0 and 4.7 versus 4.5 and 3.6 d yr -1 ) and between LLIA and MWP (4.5 and 3.6 versus 3.3 and 2.8 d yr -1 ). For the three time-slices, a higher mean than median indicates that in each sub-period few years have much higher SDY values than most of the rest. The percentage of years with snow days above the 95 th and 98 th percentiles (9.8 and 11.3 d yr -1 , respectively) is also markedly different in the three sub-periods considered, with a greater difference between the first and the following two sub-periods.

In particular, the higher percentage of years exceeding the 95 th (14%) and the 98 th (7%) percentiles confirms that the Maunder Minimum corresponded to an extreme climatic phase in southern Italy, implying the occurrence of recurring highly impacting snowy events. We know from [START_REF] Moio | Diario di quanto successe in Catanzaro dal 1710 al 1769[END_REF] manuscript (p. 244) that winter seasons during the Maunder Minimum were often hard, and winter conditions extended over spring months: In the year 1658, on April 18, the day of Holy Thursday, a lot of snow fell which raised a palm followed by an ice; all the vines were dried which were tender with great shortage of wine which for two years was sold at a very high price. Similar conditions also occurred in the autumn 1680, springs 1684 and 1694, and in May 1705 [START_REF] Diodato | Mediterranean winter snowfall variability over the past millennium[END_REF].

The recent trend towards warmer conditions indicates that rising temperatures are the main factor that triggered the decline of both SDY and cold spells during the 20 th and early 21 st century. The trend test returns a statistically significant decrease in SDY during the MWP (p < 0.05), whereas no trend (p > 0.05) is detected in the MM and the LLIA (Table 1). Observed records indicate (Fig. 7A) that the 95 th percentile was exceeded four times after the change-point year: in 1935, 1940, 1956 and 1963 with 10, 13, 10 and 10 d yr -1 , respectively. In 1940, with 13 snow days, the 98 th percentile was also exceeded. It is important to consider that whether the 95 th percentile was caught in 1940, the modelled peak values are all below that value after the change-point. According to the overall tendency of the observed records, which is representative of snowfall changes in the BNV, the model estimates do not always fit the snowfall peaks. Thus, it can also be assumed that the high SDY values occurred in the past could not be fully captured in the modelled timeseries.

Influence of atmospheric-oceanic circulation patterns

We analysed the links between large-scale atmospheric-oceanic circulation patterns and snowfall in different climatic sub-periods (MM, LLIA and MWP) through linear correlations between low-pass filtered standardized SDY anomalies (>11 years cut-off period) and gridded reconstructions of: (1) air surface temperature (T; Luterbacher et al.

2004);

(2) mean sea-level pressure (mslp; [START_REF] Luterbacher | Reconstruction of sea level pressure fields over the Eastern North Atlantic and Europe back to 1500[END_REF]; and (3) sea-surface temperature (SST; LMR and HadISST).

The increase of SDY at the BNV during the three sub-climatic periods is generally associated with an anomalous high-pressure system located over northern (Figs. 8A,B) or north-western (Fig. 8C) Europe and an anomalous low over the eastern Mediterranean.

These large-scale configurations (consistent with those described in Section 2.1) favour the advection of cold continental air from the east-northeast towards Southern Italy [START_REF] Fazzini | Snowfall analysis over peninsular Italy in relationship to the different types of synoptic circulation: First results[END_REF][START_REF] Dafis | The atmospheric circulation characteristics favouring snowfall in an area with complex relief in Northwestern Greece[END_REF], which can moisten as it crosses the Adriatic Sea in a ways analogous to lake-effect snow in the Great Lakes region [START_REF] Shi | Impact of Lake Surface Temperature Variations on Lake Effect Snow Over the Great Lakes Region[END_REF]. As expected, colder-than-usual surface air temperatures are observed over large parts of Europe (including the BNV) for all three sub-climatic periods, associated with the SDY increase (Figs. 8A-C). Differences in the position of the anomalous high over northern European latitudes could be attributed to the modulating role of global average temperatures (minimum during the MM) for the advection of cold and moist air towards the BNV, but also to the uncertainties associated with our reconstructed time-series (larger during the MM; [START_REF] Luterbacher | European Seasonal and Annual Temperature Variability, Trends, and Extremes Since 1500[END_REF]). Despite correlation between SDY and low-frequency (> 11yr) NAO reconstructed time-series [START_REF] Jones | Extension to the North Atlantic oscillation using early instrumental pressure observations from Gibraltar and south-west Iceland[END_REF][START_REF] Trouet | Persistent Positive North Atlantic Oscillation Mode Dominated the Medieval Climate Anomaly[END_REF][START_REF] Tardif | Last Millennium Reanalysis with an expanded proxy database and seasonal proxy modeling[END_REF] returns negative values for their temporal overlapping periods, these are not statistically significant (p > 0.05). The fact that mslp anomalies associated with SDY appear displaced downstream in the North Atlantic (Figs. (PDO). Linear correlation between SDY and the PDO index (from LMR) during the MM returns a statistically significant value (0.46, p < 0.05) and is consistent with the significant periodicities detected in SDY in Fig. 7B (around 20 to 60 years). However, this is not evident during the LLIA and MWP periods (not shown). We refrain at this point from proposing any causal relationship between PDO and SDY during the MM, which would require a much more detailed analysis supported by additional climate reconstruction databases (beyond the scope of this paper). For the AMO, no robust correlation was found with SDY in any of the sub-periods studied.

Finally, the SST correlation results in Fig. 8F for the MWP are also worth mentioning.

As this period is clearly dominated by the impact of global warming, the figure shows a robust global anti-correlation signal with SDY. When the linear trend was removed from the time-series, no correlation pattern with SST anomalies (from both LMR and HadISST datasets) was observed (not shown).

DISCUSSION

According to [START_REF] Camuffo | 500-year temperature reconstruction in the Mediterranean Basin by means of documentary data and instrumental observations[END_REF], for the Mediterranean area, a major effort to transform early observations into homogenized series through rigorous quality controls, validation and correction became possible only after the 17 th century. We thus used 1645 data onwards because weather data at earlier times were reconstructed from noninstrumental proxies only and, as such, were affected by larger uncertainties (Diodato et al. 2014). Thanks to our reconstruction, we can now comment and discuss on single, exceptional winter seasons. Observing Fig. 7A, it stands out the surprising number of snowy days (32) occurred in 1683-1684 during the MM. This is the estimated absolute maximum value of the modelled time-series, which exceeds by over 10 times the standard deviation of the entire series (2.96 d yr -1 ). The extraordinary cold winter 1683-1684 was particularly long and harsh, with unusual low mean temperatures, passed into the historical chronicles because of its bitter coldness and snowstorms [START_REF] Diodato | Discovering the anomalously cold Mediterranean winters during the Maunder minimum[END_REF]. We refer to the Moio and Sussana (1977) manuscript (p. 247), reporting a long period of winter weather from January to April 1684: with heavy snowfalls and continuing frost in southern Italy, with winter that was severe almost always covered with snow and ice and so continued and large that froze rivers and perished all the citrus trees. In the Corradi (1865-1890) annals (vol. II, p. 257), 1684 is also identified as a remarkable year, during which snow covered the ground until after Easter (occurred on April 2 nd ). In Fig. 9A, the January to March 1684 anomalies of reconstructed surface air temperature [START_REF] Luterbacher | European Seasonal and Annual Temperature Variability, Trends, and Extremes Since 1500[END_REF] and mslp [START_REF] Luterbacher | Reconstruction of sea level pressure fields over the Eastern North Atlantic and Europe back to 1500[END_REF] compared to the long-term seasonal mean are shown.

Consistent with the historical chronicles, colder than usual conditions were present over the BNV and the rest of Europe, with regional departures of up to -6 ºC from the long-term seasonal mean (north-eastern Europe). The anomalous large-scale mslp conditions of that year are also consistent with the findings in Fig. 8, with positive anomalies over Scandinavia and negative over the Mediterranean, promoting the advection of cold continental air from the east into the BNV.

To find a similar atmospheric configuration throughout the winter season, one has to go to the winters of 1829-1830 and 1939-1940. In the winter of 1830 (Fig. 9B), atmospheric conditions were similar to those of 1684. The most notable difference is the location of the coldest temperature anomalies over central and eastern Europe. That winter was one of the earliest and longest of those examined in the BNV, with abundant snow from mid-November to the end of March, which was lethal for agriculture in many parts of Europe (Corradi 1865(Corradi -1890)). Copious snow also fell in Rome and, in January 1830, the press was suspended due to heavy snowfall, which prevented the mails from arriving [START_REF] Serbati | Studi sul pensiero filosofico e religioso dei secoli XIX e XX[END_REF]).

In the winter of 1939-1940 (Fig. 9C), the cold anomalies appear less pronounced over southern Europe and the anomalous mslp conditions differ slightly, with the anomalous high centred near Iceland (as in Figs. 8C-D) and the low west of Iberia. This configuration also favoured the advection of cold continental air from the northeast into Europe. For that winter, almost 10 days of snow were estimated in the BNV. In Naples, eight days of snow were recorded, with 35 cm depth and trains at a standstill [START_REF] Majo | Le nevicate a Napoli (1866-1957)[END_REF]. Copious snowfall was also recorded in Rome, with 15-30 cm of snow depth [START_REF] Mangianti | La neve a Roma dal 1741 al 1990[END_REF], while 11 days of snow were registered in Florence [START_REF] Borchi | La neve a Firenze[END_REF].

After the winter 1955-1956, which exceeded the 95 th percentile with its 12 days of snowfall in February, winters never again reached that extreme value. These results are consistent with data indicating that temperature-induced precipitation has shifted from snow to rain during the MWP (Fig. 7A), with a general trend of decreasing snowfall in the central Apennines (Perugia station, [START_REF] Pandolfi | Due secoli di nevicate a Perugia: Analisi preliminare[END_REF] and the Swiss alpine region (Diodato et al. 2020b). This is not surprising, as the reconstruction approach used here relies mainly on parameters related to air temperature. Although this variable is crucial for snowfall, additional factors may play a determinant role in the context of global warming (D' Errico et al. 2020), which makes future projections of snowfall still a difficult task.

CONCLUSIONS

This study provides insights into the long-term fluctuations of snow days per year

  BNV, situated at a mean latitudinal position of Italy and surrounded by mountains, experiences a wide range of climatic characteristics common to those occurring in most of peninsular Italy (Diodato & Bellocchi 2020). The concept of the model (Fig. 4) summarizes the mechanisms mainly driving the common patterns of change of snow-cover extent and duration, and of snow days. It shows that in Italy, ~10-15 % of the territory (central tendency) is covered by snow over winter. In this season, percentiles above the 50 th percentile of the fraction of snow-cover extent deviate more from the central values than in other seasons. This indicates that SSI > 0, which captures extreme snowfall years (roughly above the median), is an important driver of snowfall, as it reflects the winter snowfall regime. Snowfall then responds to local temperature forcing during the winter and spring months. In contrast to the Italian mountain chains, where the increase in snow depth

  ) -A = 20, φ = 0.25, η = 1.4, β = 28.5 °C and σ = 2the model approximated the actual data. Fig.5Ashows a satisfactory agreement between observations and estimates, as the points tend to line up around the 1:1 line, with only one outlier (year 1906, black square in the scatterplot, with 5.5 modelled snowy days against one observed) in the 99-year dataset (~1%), which was not considered for model calibration. The R 2 -value indicates that the NLMRM SDY(BNV) explains 85% of the variance in actual data. The regression line has a y-intercept a = -0.36 ± 0.34 and slope b = 1.08 ± 0.07 near the optimum values (a = 0 and b = 1). Although a few zero-values on the x-axis indicate that the predictive ability of the model is lower for SDY values close to zero, the intercept is not statistically different from zero (Student-t p > 0.05). The MAE, equal to 1.02 d yr -1 , depicts a lower value than the standard error of the residuals (1.29 d yr -1 ), while the Kling-Gupta index is 0.54. There is no indication of serial autocorrelation in the residuals (DW = 1.78, p > 0.05). The model residuals return a normal distribution shape (Fig.5B; Jarque-Bera test p > 0.05) and the quantile-quantile (Q-Q) plot (Fig.5C) shows a distribution of sample-quantiles around the theoretical line, indicating only a few skewed high SDY values.

  Figs. 8E, F (only significant patterns for climatic sub-periods shown). During the MM, the

Fig. 1 :

 1 Fig. 1: The slopes of Mount Vesuvius under the snow, seen from Torre del Greco (40° 47′ N, 14° 22′ E) in an image of an anonymous of the 19 th century

Fig. 2 :

 2 Fig. 2: (A) Geographical setting in Europe, with central-southern Italy. (B) location of the

Fig. 3 :

 3 Fig. 3: Monthly mean distribution of snow depth (thick black line) and seasonal

Fig. 4 :

 4 Fig. 4: Monthly snow-cover extent fraction in Italy (50 th percentiles in red, 80 th percentiles in green, 98 th percentiles in blue), with the associated driving factors. Percentiles (prcs)

Fig. 5 :

 5 Fig. 5: (A) Scatterplot of observed and modelled SDY for the calibration sub-set (1870-

Fig. 6 :Fig. 7 :

 67 Fig. 6: (A) Timeline co-evolution of observed and modelled SDY (see legend) for the

  Modelling construction and evaluation was carried out in this work with the statistical and graphic software STATGRAPHICS (https://www.statgraphics.com) and WESSA (http://www.wessa.net). Exploratory and time-series analyses (breakpoint detection tests,

	wavelet	power	spectrum)	were	performed	using	AnClim
	(http://www.climahom.eu/software-solution/anclim), PAST routines (https://palaeo-
	electronica.org/2001_1/past/issue1_01.htm)	and	MATLAB	functions
	(https://www.mathworks.com/products/matlab.html).

Table 1 :

 1 (SDY) over a 374-year period (from 1645 to 2018) in the Benevento Valley in southern Italy. A detailed dataset was compiled in the region based on records available from two main stations (Benevento city observatory and the Met European Research Observatory) since 1870. Then, they were reconstructed back to 1645 using a parsimonious model sensitive to winter and spring temperatures and a storm-severity index based on historical documentary sources. It is important to highlight the relevance of this unique long-term dataset for studying changes in snowfall in an area where the literature on this subject is scarce. Indeed, this variable is crucial to better understand climate-change impacts on the energy balance and the water cycle, with important implications for sectors such as agriculture and tourism.The reconstructed SDY time-series shows a change-point in 1866, after which a significant decrease is observed. Overall, the trend of decreasing snow days since the outbreak of the LIA is in line with previous literature, not only in relatively close areas of Italy and Europe, but also elsewhere. Before the change-point year, no trends were detected during the Maunder Minimum (1645-1715) and the LLIA (1716-1866). The statistically significant link observed between SDY and the Pacific Decadal Oscillation during the Maunder Minimum, which is not sufficient to uncover the underlying physical mechanisms of this relationship given the uncertainties present in past data, is a clue for further analyses.What instead emerges in all sub-periods is that the large-scale circulation conditions fostering snowfall in the BNV were generally characterized by an anomalous high-pressure system located over northern-north-western Europe and an anomalous low over the eastern Mediterranean. These conditions favour the advection of cold continental air from the Balkan Peninsula and Siberia, which moistens as it crosses the Adriatic Sea towards the BNV.Although the results presented in this study are specific for the BNV, our statistical modeling approach could potentially be applied to adjacent Mediterranean or additional regions where reliable long-term snowfall records were available and surface temperature advection from poleward regions (determined by atmospheric circulation patterns; Croce et al. 2018) was the dominant driver for snowfall occurrence. Nevertheless, further studies on extended spatial scales are necessary to disentangle the mechanisms behind climatic drivers and snow cover responses, especially in the context of future climate projections (D'Errico et al. 2020). Descriptive statistics of the modelled SDY time-series for three climatic subperiods (MM: Maunder Minimum; LLIA: Late Little Ice Age; MWP: Modern Warming Period). Trends detected through Mann-Kendall test (p < 0.05).

		Central tendencies	Percentiles
	Climatic sub-	Mean				Trend
	period	(standard	Median	>95 th	>98 th
	(CE years)	deviation)
		d yr -1		% of years	d 10 -2 yr -1
	MM (1645-1715)	6.0 (4.9)	4.7	14.0	7.0	No
	LLIA (1716-1866)	4.5 (2.6)	3.6	5.9	2.0	No
	MWP (1867-2018)	3.3 (1.9)	2.8	<1.0	0.0	-1.1

Acknowledgements

This was an investigators-driven study without financial support.

Data availability

Long-term seasonal temperature and monthly sea level pressure reconstruction datasets over Europe from [START_REF] Luterbacher | Reconstruction of sea level pressure fields over the Eastern North Atlantic and Europe back to 1500[END_REF][START_REF] Luterbacher | European Seasonal and Annual Temperature Variability, Trends, and Extremes Since 1500[END_REF] and reconstructed North Atlantic Oscillation indices from [START_REF] Trouet | Persistent Positive North Atlantic Oscillation Mode Dominated the Medieval Climate Anomaly[END_REF] and [START_REF] Jones | Extension to the North Atlantic oscillation using early instrumental pressure observations from Gibraltar and south-west Iceland[END_REF]: Climate Explorer (monthly and seasonal historical reconstructions; monthly and annual climate indices): https://climexp.knmi.nl; CRU Global Climate Dataset: Time-series of variations in climate with variations in other phenomena v3: https://catalogue.ceda.ac.uk/uuid/3f8944800cc48e1cbc29a5ee12d8542d; Last Millennium Reanalysis (LMR) Project Global Climate Reconstructions Version 2: https://www.ncdc.noaa.gov/paleo-search/study/27850; Met Office Hadley centre for the

Competing interests

The authors declare no competing interests.