The sensitivity of present-time electricity demand on past climate change: a case study for Italy

<jats:title>Abstract</jats:title>           <jats:p>A methodology for estimating secular daily minimum, mean and maximum (T<jats:sub>n</jats:sub>, T<jats:sub>m</jats:sub> and T<jats:sub>x</jats:sub>) temperature records for any urbanised point of a 30-arc-second-resolution grid covering Italy is presented. It is based on the superimposition of 1961–1990 climatologies and departures from them (anomalies). The anomalies are obtained by applying inverse distance weighting to 143 Italian high-quality records, whereas the climatologies are based on a larger dataset and on the application of local weighted linear regression of temperature versus elevation.</jats:p>           <jats:p>The grid-point T<jats:sub>n</jats:sub>, T<jats:sub>m</jats:sub> and T<jats:sub>x</jats:sub> records are then used to set up secular records (period 1801–2013) of temperature-derived variables that influence Italy present-time national electricity demand. They are national averages over Italian urbanised areas of cooling degree-days (CDD), heating degree-days (HDD) and solar radiation deficit with respect to a defined threshold (S), with solar radiation estimated using daily temperature range as a proxy.</jats:p>           <jats:p>The monthly and yearly sums of the daily CDD, HDD and S records are then used, alongside with a model allowing to link these variables to present-time Italy electricity demand, in order to understand the impact of climate variability and change on present-time Italian electricity demand. We find that temperature changes as the ones observed in the last two centuries are capable of altering significantly the present-time monthly profile of the electricity demand, raising (lowering) summer (winter) months contributions. The impact is higher in summer months where it exceeds 5 % of present-time Italy average monthly electricity demand, whereas the decrease of the winter demand is rather low because of a very limited use of electricity for heating. The summer and winter opposite-sign changes result globally in an increase of the yearly demand of about 5 TWh, corresponding to about 1.5-2.0 % of present-time Italy yearly electricity demand.</jats:p>

section_content: Spatial climate datasets in digital form are currently in great demand and gridded estimates of 30 years climatological normals are requested by a variety of models and decision support tools, such as those used in agriculture, engineering, hydrology, ecology, energy management and natural resource conservation (Daly et al. 2002; Daly 2006 ). These datasets have to be set up providing a realistic representation of the major forcing factors that affect spatial climate patterns, in order to give reasonable estimates also for areas with poor station coverage, such as high elevation sites in mountain areas (Daly et al. 2008). Beside the spatial distribution of the climatological normals, it is also important to describe the spatio-temporal behaviour of climate variability and change. It is therefore necessary to set up methodologies that allow for estimating both the climatological normals of any point of the territory and the corresponding long-term records. As the desired resolution turns often out to be rather high (e. g., 1 km 2 ), such spatialisation methodologies present a number of non-trivial issues. 

In this context, this paper discusses the estimation of daily secular records of minimum, mean and maximum temperatures (T n, T m, and T x ) for any cell of a 30-arcsecond resolution grid covering Italy classified, according to the EC JRC Global Land Cover 2000 (European Commission 2003), as "artificial surfaces and associated areas" (these grid-cells are hereinafter called "urbanized"). These records are then used to investigate past variability and change of temperature-derived variables that influence present-time national electricity demand and to estimate the impact of past variability and change on present-time Italy electricity demand. Such an estimation may be useful within the estimation of the costs of climate change. 

The basic assumption of the methodology we use to set up the grid-cell daily temperature records is that the spatio-temporal structure of the signal of a meteorological variable over a given area can be described by the superimposition of two fields: the climatological normals over a given reference period (i. e. the climatologies) and the departures from them (i. e. the anomalies) (New et al. 2000; Mitchell & Jones 2005; Brunetti et al. 2009; Brunetti et al. 2012). The climatologies, which we produce with the procedure described in Brunetti et al. (Brunetti et al. 2014), are basically linked to the geographical features of the territory, therefore they exhibit remarkable spatial gradients. On the contrary, the anomalies, which we produce by inverse distance weighting, are linked to climate variability and change and are characterized by higher spatial coherence. We reconstruct the two fields in a completely independent way from each other, using different datasets: for the anomalies the priorities are data quality and the availability of long records, whereas for the climatologies the most relevant aspect is the availability of a large number of stations. In this case, in fact, we consider the issue of record length as less relevant, as a 30 year period is enough to estimate climate normals. 

The technique we consider to capture the dependence of present-time Italy electricity demand on meteorological variables is that presented by Scapin et al. (Scapin et al. 2015). They model the demand of ordinary days by means of a linear regression model made up of i) a timedependent term capturing long-term trends, ii) a term accounting for the differences among the different days of the week and iii) a set of terms depending on temperature-derived variables that can be obtained from high-resolution T n, T m, and T x fields. The latter are national averages over Italian urbanised areas of cooling degree-days (CDD), heating degree-days (HDD) and solar radiation deficit with respect to a defined threshold (S), with solar radiation estimated using daily temperature range as a proxy. Applying this technique, Scapin et al. (Scapin et al. 2014) quantified the present-time dependence of the Italian electricity demand of an ordinary day on CDD, HDD and S in, respectively, 24. 6 GWh degree-day -1, 6. 6 GWh degree-day -1 and 2. 9 GWh day -1 MJ -1 m 2. These values give evidence of a relevant contribution of meteorological conditions on the electricity demand, especially in the summer period: in this season, the fraction of the electricity demand driven by meteorological conditions can exceed 20 %. Temperature is indeed the most important meteorological variable influencing the electricity demand: it regulates the request for conditioning in summer and for heating in winter (Pardo et al. 2002; Moral-Carcedo & Vicéns-Otero 2005; Hor et al. 2005; Bessec & Fouquau 2008; Hekkenberg et al. 2009; Apadula et al. 2012). The latter request is however rather low in Italy as the use of electricity for heating is limited. 

Variability and change of CDD, HDD and S from the beginning of the 19 th century to present can be used, alongside with the model proposed by Scapin et al. (Scapin et al. 2015), to get a rough estimation of how the present-time electricity demand would be if the climate were still in the situation prior to global warming. To this purpose, we apply this model keeping, on the one hand, the time-dependent term and the term accounting for the differences among the different days of the week as they are at present-time and considering, on the other hand, the secular CDD, HDD and S records we get from the corresponding high-resolution T n, T m, and T x fields over the 1801-2013 period. In other words, we investigate what would happen in present-time conditions if the meteorological variables changed as they did in the past. The shortcoming of this approach is that presenttime meteorology-electricity demand relationships are not completely independent from past climate variability and change as global warming may, e. g., have caused a wider diffusion of air conditioning devices. We assume, however, that this dependence is weak and that meteorology-demand relationships are mainly influenced by socio-economic factors such as the time evolution of the Gross Domestic Product (GDP). The CDD, HDD and S secular records and the virtual electricity demand record we get using them in the model proposed by Scapin et al. (Scapin et al. 2015) allow therefore giving a rough assessment of the impact of past climate variability and change on Italy present-time electricity demand. 

As far as we know, this assessment was still missing for Italy. In fact, although some papers have discussed the relationships between meteorological variables and the national electricity demand (Bessec & Fouquau 2008; Apadula et al. 2012; Manera & Marzullo 2005; Lee & Chiu 2011), none of them studied how much of present-time demand can be ascribed to past temperature changes. Therefore this paper represents a first step toward the wider objective of assessing the impact of past temperature changes on present-time Italy energy demand, which has to be estimated considering, beside the electricity demand, also the use of natural gas, fuel oil and other fuels. They are particularly relevant for heating energy use. 

The paper is organised in six parts. After the introduction, we focus on the climatologies; then we present the anomaly records and the procedure we adopt in order to superimpose the climatologies and the anomalies. After these parts on the meteorological data, we discuss the method we use to set up national average records allowing to capture the dependence of present-time Italy electricity demand on meteorological variables; then we present past variability and change of these records and discuss the relevance of their changes on the light of the sensibility of Italy present-time electricity demand on CDD, HDD and S. In the end, we summarize the results of the paper and discus some open issues. (Hiebl et al. 2009) for the Greater Alpine Region, regression kriging (RK), widely used in the literature (see e. g. (Hengl 2009 )), and local weighted linear regression (LWLR) of temperature versus elevation, based on the PRISM (parameter elevation regression on independent slopes model) conceptual framework (Daly et al. 2002; Daly 2006; Daly et al. 2008; Daly et al. 1994). The performances of these methods were evaluated by estimating, with a leave-one-out approach, the climatologies at the station locations and then comparing the predicted values with the corresponding observed station normals. All three approaches led to quite reasonable estimates of the station normals, with the lowest errors in spring and autumn and the highest errors in winter. However the LWLR approach showed slightly lower errors than the other two approaches (root mean square errors (RMSEs) range from 0. 74 °C (April and May) to 1. 03 °C (December)). The better performance of LWLR was even more evident when selected station clusters were considered, giving evidence of a greater reliability of a local approach in modelling the behaviour of the temperature-elevation relationship in Italy's complex territory.
section_name: Introduction

section_content: For the reasons discussed above, we use the same LWLR approach for T n and T x as well. For this purpose, we consider a dataset of 1109 stations distributed over the entire Italian territory (Fig. 1a ) with monthly 1961-1990 T n and T x normals that have been subjected to the same quality-control procedure carried out for T m in Brunetti et al. (Brunetti et al. 2014 ) and we apply LWLR as described in that paper. 

Specifically, we use a weighted linear regression of the data from nearby stations to predict the monthly temperature normal at any GTOPO30 grid-cell as a function of its elevation. The weights of the stations involved in the regression depend on their geographic similarity with the grid-cell. They are obtained by the product of five Gaussian weights (w) of the form:

cvar where (λ g, ϕ g ) is the position of the grid-cell, var is the specific geographical variable which is considered (position, elevation, distance from the sea, slope steepness and slope orientation), Δ i var is the difference between the value of this variable at the grid-cell and that at the i-th station location and c var is a coefficient which regulates the decrease of the weight for increasing values of Δ i var. The stations selected for the regression are the 35 with the highest weights among those within 200 km from the given grid-cell. Instead of progressively extending the search range until a sufficient number of stations is available, we consider a very large search range (200 km) and choose only the stations with the highest weights. This is because the radial distance from the grid-cell is not always the leading discriminant. 

As the position weight is the most important one for the performance of the procedure, we optimized its c var for each month by minimizing the average RMSEs over all stations, with errors estimated with a leave-one-out approach. This procedure was performed independently for T n and T x. For the other geographical variables, we simply used the c var values proposed by Brunetti et al. (Brunetti et al. 2014). 

After such optimization of the parameters, we estimated the errors of the climatologies at the station locations. They are reported in Table 1 in terms of mean error (BIAS), mean absolute error (MAE) and RMSE. Table 1 gives evidence of higher errors for T n ; moreover, we observe that summer and winter months have in general higher errors than spring and autumn months. 

Once the monthly T n, T m, and T x 1961-1990 climate normals are available, the next step of our procedure consists in obtaining corresponding daily values by fitting them by means of the first two harmonics of a Fourier series. We assume then that the error due to this step is negligible and that Table 1 represents the errors of daily climatologies too.
section_name: High-resolution climatologies

section_content: Beside the climatologies, our methodology requires the estimation of secular daily temperature anomaly fields that have to be superimposed to the grid-cell temperature normals.
section_name: Grid-cell anomaly and temperature records

section_content: A key issue for the application of our methodology is the availability of a high-quality database of long-term secular records. Here we use updated and improved versions of the datasets presented by Brunetti et al. (Brunetti et al. 2006a) and Simolo et al. (Simolo et al. 2010). Updating mainly concerns the last 10 years, whereas the activities for improving the records mainly concern data homogenisation. 

The Brunetti et al. (Brunetti et al. 2006a ) dataset was used to assess Italian temperature trends up to 2003. It is mainly based on the data collected at observatories established in the 19 th century: most of them ended their observations in about the last 30 years. The Simolo et al. (Simolo et al. 2010 ) dataset concerns the stations of the network of the Italian Air Force (Aeronautica Militare, hereinafter AM): a significant fraction of them are currently managed by the Italian Agency for Civil Aviation (ENAV). We update most of the stations of this network in near real time by means of the Global Surface Summary of the Day, (GSOD) managed by US National Climatic Data Center (NCDC). 

The dataset we use in this paper includes a significant fraction of composite records. In some cases it was just necessary to merge data from different sources for the same station (e. g. AM/ENAV and GSOD), in other cases the merging concerned also data from different stations, such as some of the oldest observations from the Brunetti et al. (Brunetti et al. 2006a ) dataset merged in the last decades with the AM/ENAV records from the Simolo et al. (Simolo et al. 2010) dataset and updated for the last years with the GSOD records. A detailed homogenisation was therefore mandatory. We performed it subjecting all monthly temperature records to a relative homogeneity test based on the procedure described by Brunetti et al. (Brunetti et al. 2006a). In this procedure, each series is tested against 10 other series by means of a multiple application of the Craddock test (Craddock 1979). When a break is identified in the test series, some reference series are chosen among those series that prove to be homogeneous in a sufficiently large period centred on the break and that correlate well with the test series. Several series are used in order to better identify the break and to get more reliable adjustments. When a break is homogenised, the preceding portion of the series is corrected, leaving the most recent portion of the series unchanged. This allows for updating the records without considering any adjustment. Homogenisation is performed at daily resolution, with the adjustments obtained fitting the first two harmonics of a Fourier series to the adjustments we get for the monthly records. 

The dataset we use in this paper to get the grid-cell anomalies consists of 143 daily T n, T m and T x records. Stations sites are shown in Fig. 1b, whereas Fig. 2 shows data availability versus time. 

Table 1 Errors of T n, T m and T x climatologies at station locations. T m errors are retrieved from Brunetti et al. (Brunetti et al. 2014 ). All errors have been estimated with the leave-one-out approach Tn Tm Tx
section_name: Station anomaly records

section_content: 1 -0,04 1,06 1,35 -0,04 0,77 1,01 -0,01 0,81 1,10 2 -0,04 1,02 1,30 -0,04 0,69 0,90 -0,01 0,76 1,02 3 -0,04 0,94 1,19 -0,03 0,60 0,78 -0,01 0,75 0,98 4 -0,04 0,90 1,14 -0,03 0,58 0,74 -0,01 0,76 0,99 5 -0,03 0,93 1,17 -0,02 0,58 0,74 -0,01 0,80 1,02 6 -0,02 1,00 1,26 -0,01 0,62 0,79 -0,01 0,87 1,11 7 -0,02 1,09 1,37 -0,02 0,66 0,85 -0,01 0,94 1,19 8 -0,02 1,07 1,35 -0,02 0,65 0,84 -0,01 0,91 1,15 9 -0,03 1,01 1,27 -0,02 0,62 0,79 -0,01 0,86 1,08 10 -0,03 0,96 1,20 -0,02 0,63 0,81 -0,01 0,78 1,03 11 -0,03 0,95 1,19 -0,03 0,67 0,86 -0,01 0,74 1,00 12 -0,04 1,02 1,29 -0,04 0,78 1,03 -0,01 0,83 1,13

Errors of the climatologies

After the homogenisation of the station records, the next step of our methodology consists in filling the gaps in the monthly-homogenised records over the period by means of the procedure described by Brunetti et al. (Brunetti et al. 2006a ). The completed records are then used to calculate monthly 1961-1990 station normals from which corresponding daily normals are estimated by fitting them with the first two harmonics of a Fourier series. We prefer this approach with respect to directly calculating daily station normals because i) it is the same approach we use to get the daily grid-cell normals on the basis of monthly high-resolution climatologies, ii) the monthly records are easier to complete than the daily records and iii) direct calculation of the daily normals from the daily data would produce rather noisy values. The homogenised daily temperature records are finally transformed into anomaly records with respect to the 1961-1990 normals, by simply subtracting from each data the corresponding daily normal.
section_name: Month BIAS MAE RMSE BIAS MAE RMSE BIAS MAE RMSE

section_content: Once homogenised station anomaly records are available, the next step of our methodology consists in using them to evaluate, by means of a leave-one-out approach, the performance of the interpolation technique. Specifically, we predict the anomaly record A(t) at each station site (λ i, ϕ i ) by means of a weighted average of the anomaly records of the other stations and we compare it with the corresponding observed record. The basic assumption of this procedure is that the comparison of the predicted and the observed records at the station-sites, gives also a reasonable picture of the accuracy of the predicted records at the grid-cells. 

The predicted record at the i-th station is constructed by first building weighting terms accounting for distance and elevation difference from each of the other stations:

if j ≠ i and if station j data are available at time t; otherwise w

where d j,i and Δz j,i are, respectively, the distance and the elevation difference between the station under analysis and the j-th station and τ d and τ Δz regulate the extent to which a station is weighted with respect to distance and elevation difference: exponentials decrease by half when distance reaches τ d and, in analogy, when elevation difference equals τ Δz. For these parameters we use here the values proposed by Scapin et al. (Scapin et al. 2015). They are 80 km for τ d and 400 m for τ Δz. We use then these weighting terms for the projection of the observed anomalies on the site of the i-th station:

With the data availability of our dataset -before about 1870 the distribution of the stations is so inhomogeneous that they can not be considered as representative of the grid-cells -this simple procedure can however not be applied to assess the accuracy of the grid-cell records. Moreover, as the urbanised grid-cells we consider in this paper are mainly located at low elevation (99. 5 % of the urban grid-cells are below 750 m) we prefer restricting the comparison to the station-sites below 750 m. Therefore, we first divided the 1801-2010 interval in 42 consecutive 5 year periods. Then we looked for the stations with available data in each of these periods and we used only the records of these stations to estimate, by means of (1) and ( 2 ), the station-site anomalies in the 1974-2003 period (i. e. the period with the best data availability). In case a record used for the estimation was not available in the 1974-2003 period, it was simply substituted with a neighbouring one, in order to prevent missing data to bias the estimation process. Finally, we compared the estimated and the observed station anomalies in the 1974-2003 period. The comparison was therefore performed for a subset of stations (the 122 stations with available data in the 1974-2013 period) covering the entire Italian territory, whereas the estimation of the anomalies at the station sites was performed with a subset of stations defined on the basis of the data availability in each of the 5 year periods. 

The results of this comparison are shown in Fig. 3a by means of box-plots, which give evidence of the distributions of the RMSEs of the estimated station anomalies different station subsets we use for the estimation of the 1974-2003 anomalies give equivalent results. In this period the median of the RMSEs of the station T m anomalies ranges between 0. 97 °C (stations with available data in the 1956-1960 period) and 1. 09 °C (stations with available data in the 1996-2000 period), whereas for T n and T x the values range, respectively in the intervals 1. 33-1. 40 °C and 1. 35-1. 51 °C. On the contrary, for the station subsets corresponding to data availability before 1950 the errors are larger, showing that the smaller is the subset of the stations used for the estimation of the anomalies, the larger is the error of the reconstructed anomalies. With the station subsets corresponding to data availability before 1880, and even more with those corresponding to data availability before 1870, there is an additional difficulty in the reconstruction of the station anomalies. In fact, in this period, almost all the stations are in northern Italy. Due to this problem, for a significant fraction of the stations the errors become so high that the RMSE of the estimated anomalies becomes higher than the standard deviations of the observed anomalies. 

In order to better highlight the influence of stations in central and southern Italy on the errors before 1880, we applied the previous procedure also evaluating for the station subsets concerning the years before 1880 only the errors of northern Italy stations (i. e. the stations with latitude above 44°). The results are shown in Fig. 3b. They give evidence of much lower errors than in Fig. 3a, with almost all T m 1974-2003 RMSEs within 2 °C.
section_name: Assessment of the anomaly record interpolation technique

section_content: Once the reliability of the technique for the projection of the anomalies has been verified, it has been applied to any urbanized grid-cell of the DEM. The construction of the grid-cell records was naturally limited to the areas and the periods for which the analyses performed for the station-sites highlighted reasonable errors. In particular, we did not reconstruct central and southern Italy anomalies before 1880. 

Once daily T n, T m, T x anomaly records were estimated for all the relevant grid-cells, we obtained temperature records by simply superimposing the anomaly records to the corresponding grid-cell climate normals. The construction of such daily grid-cell records had not yet been performed for Italy, whereas corresponding monthly records had already been presented by Brunetti et al. (Brunetti et al. 2009; Brunetti et al. 2012).
section_name: Grid-cell temperature records

section_content: A number of recent papers have shown that temperature has a strong impact on electricity demand in many areas of the world (see e. g. (Bessec & Fouquau 2008; Apadula et al. 2012; Feinberg & Genethliou 2005) ). 

In this context, Scapin et al. (Scapin et al. 2015) have developed a linear model linking the Italian daily aggregate electricity demand in the 1990-2013 period to temperature-derived variables. The model is based on the superimposition of deterministic components related to the weekly cyclical demand pattern and to long-term demand changes and on weather sensitive components. It assumes that the Italian daily electricity demand (D(t)), can be described by means of the following relation:

where:

The first term consists of a third order polynomial which aims to describe the temporal evolution of electricity demand caused by the economic conjuncture and by long-term changes in consumption habits. 

The second term consists of dummy variables I j (t), introduced in order to account for the strong weekly pattern of electricity demand. In particular I 1 (t) is set to 1 on Monday and 0 otherwise; in analogy, I 2 (t), I 3 (t) and I 4 (t) are used to model the behaviour of Friday, Saturday and Sunday. No dummy variables are used for central weekdays (Tuesday, Wednesday, Thursday) which can be grouped together since they exhibit similar behaviour. In this way, the dummy variable term accounts for the differences between central weekdays and the rest of the week. The third term consists of a summation over three exogenous variables (CDD, HDD and S), describing the influence of weather factors on the Italian electricity demand. The last term represents the error of the model (i. e. the difference between the actual and the estimated demand). 

The coefficients α i, β j and γ k of relation (3) are obtained by means of least squares regression. 

The model was applied to the Italian aggregated National electricity demand of ordinary days (holidays and special events are excluded) considering twelve two-year periods, starting from 1990-1991 and ending in 2012-2013 (Scapin et al. 2015). 

The CDD, HDD and S records considered in relation (3) were obtained with a bottom-up approach: first, local CDD, HDD and S records were estimated for all Italian grid-cells of the GTOPO30 DEM that are classified as urbanised according to GLC2000 land cover (European Commission 2003), then an average was computed over them. 

The first step consisted in estimating T n, T m and T x records for each Italian urbanised grid-cell. Then effective temperature (T*), a delayed signal obtained through exponential smoothing of T m series (it allows taking into account temporal inertia of buildings), was considered and cooling and heating degree-days were defined according to:

where T* is the grid-cell daily mean effective temperature, T S 1 is 20 °C and T S 2 is 15 °C (Scapin et al. 2015). 

Finally, the grid-cell daily T n and T x records were used to estimate grid-cell global solar radiation records (H) by means of the following formula (Hunt et al. 1998 ):

where a 0 and a 1 are site-dependent empirical coefficients, ΔT is the daily temperature range (T x -T n ) and H 0 is the exo-atmospheric radiation on the horizontal plane, i. e. the daily integral of solar irradiance that would be observed on a horizontally oriented surface placed at the top of the atmosphere. H 0 can be easily determined by standard computation (see e. g. Iqbal (Iqbal 1983 )), whereas a 0 and a 1 can be recovered from previous studies. Scapin et al. (Scapin et al. 2015) used for the entire Italian territory the values (a 0 = 0. 190 K -0. 5, a 1 = -2. 041 MJ m -2 ) that Abraha and Savage (Abraha & Savage 2008) proposed for Padua (northern Italy). Term S, which was introduced to take into account the effect of lighting on electricity demand, was then defined from solar radiation as (Scapin et al. 2015) :

where H S was set to 17 MJ m -2. Relation (3) explains from 97. 7 % (2008-2009) to 99. 4 % (1996-1997) of the variance of the Italian daily demand record (ordinary days only), with a mean absolute percentage error (MAPE) of about 1 % (Scapin et al. 2015). 

The CDD term shows a strong positive trend in the 1990-2007 period, followed by a tendency toward stationarity in the following 6 year period. The data of the last 3 two year periods considered by Scapin et al. (Scapin et al. 2015) allow therefore for quantifying the present-time dependence of the Italian electricity demand of an ordinary day on CDD in 24. 6 GWh degree-day -1 (Scapin et al. 2014). For HDD and S we prefer quantifying this dependence considering the average values over the latest 6 twoyear-period coefficients. This because these terms show a tendency toward stationarity in the last 12 years. The values are, respectively, 6. 6 GWh degree-day -1 and 2. 9 GWh day -1 MJ -1 m 2 (Scapin et al. 2014 ). This paper focuses on the impact of climate change on CDD, HDD and S and on the corresponding impact on present-time electricity demand. We are therefore much more interested on yearly, seasonal or, at least, monthly values, rather than on daily values. We sum then equation (3) over all days of each year, season or month (it is actually a virtual year, season or month as it is composed only of ordinary days), getting cumulated contributions obtained multiplying the γ k terms in equation (3) (i. e. 24. 6 GWh degree-day -1, 6. 6 GWh degree-day -1 and 2. 9 GWh day -1 MJ -1 m 2 ) by the yearly, seasonal or monthly sums of the daily CDD, HDD and S values. We consider therefore hereinafter only these cumulated CDD, HDD and S data. 

We know we are able to construct Italy CDD, HDD and S records covering the 1880-2013 period, whereas we are not able to construct them before 1880, as the estimation of grid-cell temperatures can be performed only for the northern part of Italy. We observe however that, even though northern and southern Italy temperature anomalies may have strong differences at daily resolution, the agreement increases when monthly, seasonal or yearly periods are considered (see e. g. (Brunetti et al. 2006a) ). The same behaviour concerns naturally the temperature-derived variables we use in this paper. We used therefore the period 1880-2013 to check whether yearly sums of CDD, HDD and S obtained only from northern Italy urbanised grid-cells can be used to estimate the corresponding national CDD, HDD and S values. The results of this analysis are shown in Fig. 4. They give evidence that northern Italy CDD, HDD and S yearly records capture a large fraction of Italy corresponding records. The latter records can therefore be estimated from the former. 

Based on these results, we extended the estimation of the yearly Italian CDD, HDD and S records also to the years before 1880. 

In order to further check the errors of the yearly CDD, HDD and S data in the first years of the 19 th century, when only Milan, Padua and Turin were available, we first estimated northern Italy urbanised grid-cell temperatures of the 1974-2003 period using only these 3 stations. Then we used these temperatures to get northern Italy yearly CDD, HDD and S records from which we estimated corresponding national records by means of the linear regressions shown in Fig. 4. The comparison of these records with those estimated by means of the full dataset gives evidence of a bias of -1. 7 % for CDD, of 0. 4 % for HDD and of 0. 5 % for S. We conclude therefore that the yearly CDD, HDD and S records we get by means of only Milan, Padua and Turin data have rather low bias. It is also interesting to check the RMSEs of the CDD, HDD and S yearly records we get by means of these 3 stations. They turn out to be 7. 1 % for CDD, 2. 5 % for HDD and 1. 5 % for S, where the percentages are expressed in relation to the average values in the 1974-2003 period, obtained from the full dataset.
section_name: Local secular records and electricity demand

section_content: The yearly Italian CDD, HDD and S records are shown in Fig. 5, together with a 10 year standard deviation Gaussian low-pass filter. 

These curves give evidence of a strong rise of CDD across the last 2 centuries, showing a more pronounced slope in the last 30 years. Consequently, the low-pass filter curve increases from about 200 degree-days at the beginning of the 19 th century to over 400 degree-days in the last years. Year 2003 has by far the highest CDD value of the whole series (699 degree-days): summer 2003 was indeed characterised by exceptionally hot weather, which was a matter of concern for many reasons (Grazzini et al. 2003; Schär et al. 2004 ). The minimum CDD value is found in year 1816 (88 degree-days), which is widely renowned as the year without a summer for its abnormally low summer temperatures (Stommel & Stommel 1979 ). On the contrary, heating degree-days have been decreasing in the examined period. The decrease, which starts at about 1865, is more regular than the increase of CDD. As for CDD, the trend has strengthened in the last 30 years. Solar radiation instead does not show a significant trend in the considered period. 

In order to investigate better the increasing and decreasing tendencies highlighted by Fig. 5, we subjected the records to running trend analysis (Brunetti et al. 2006b ). Specifically, we estimated the slopes of the time series over all time windows with minimum length of 30 years. The slopes are computed using the Theil-Sen method (Sen 1968; Theil 1983), which estimates the slope as the median of the slopes of lines crossing all possible pairs of points: such methodology is particularly suitable in presence of outliers as it significantly reduces their influence on the results. The results are shown in Fig. 6, where window widths and the starting years of the windows that the trends refer to are represented on y and x axes, respectively. Slopes are represented by means of a color scale. Significances are evaluated by the Mann-Kendall non-parametric test (Sneyers 1990). This type of running trend analysis, first introduced by Brunetti et al. (Brunetti et al. 2006b), is an instrument to investigate trends in depth and to produce plots that allow for visualizing trends on a wide range of timescales. Therefore, Fig. 6 captures the whole spectrum of significant trends in the CDD, HDD and S records and provides quantitative description of the peculiarities observed in Fig. 5. 

The results in Fig. 6 give evidence of the statistical significance of the positive (negative) trends of CDD (HDD): the trend is significant for all windows with length of at least 100 years. For shorter periods the picture is more complex. Specifically, for CDD the positive trend pixels reflect the increase in Fig. 5 between 1820 and 1870, 1920 and 1950 and from the beginning of the 1980s to present-time, whereas the negative trend pixels reflect the decrease in the last decades of the 19 th century and in the 1960s. The strongest trends are those of the last decades, with the highest value corresponding to 7. 0 degree-days/year in the 1976-2006 period. For HDD the variability is lower and only very few pixels correspond to positive trends. As well as for CDD the strongest trends are those of the last decades, with the strongest trend corresponding to -9. 4 degree-days/year (1978-2008 period). Finally, for S the trends of the longest windows are very low. For shorter time windows the strongest trends are those starting in the first half of the 20 th century: they are associated to the clear maximum in Fig. 5 between the beginning of the 1950s and the end of the 1970s. This maximum in solar radiation deficit seems to be associated to solar radiation dimming, probably caused by atmospheric aerosol linked to air pollutant emissions, even though in Italy the reduction of solar radiation in the global dimming period (Wild 2009) has been partially masked by a reduction of total cloud cover (Maugeri et al. 2001; Manara et al. 2015). 

Finally we investigated what would happen to Italy present-time electricity demand if CDD, HDD and S changed as they did in the past. Specifically, we first computed the virtual yearly series of the weatherrelated terms in relation (3), using for γ k values which are representative of present-time conditions. Then we constructed a summer contribution (driven by CDD) and a winter one (given by the sum of the HDD and S contributions), considering all days as ordinary. These contributions are shown in Fig. 7. This figure gives evidence that the winter contribution is much greater than the summer one with the temperatures we had at the beginning of the 19 th century, whereas with present-time temperatures the winter and the summer contributions are rather similar. The net effect of these changes (total curve in Fig. 7 ) is a slight rise of the weather-related contribution to electricity demand, driven by the recent increase of CDD. In fact, with the temperatures corresponding to the beginning of the 19 th century, the weather related contribution to present-time electricity demand would be nearly 20 TWh, whereas with nowadays temperatures it is around 25 TWh. We estimate therefore the contribution of past climate change on Italy presenttime electricity demand in roughly 5 TWh, which corresponds to about 1. 5 %-2. 0 % of present-time Italy electricity demand. 

As we did for the CDD, HDD and S yearly data, we subjected also the total curve we show in Fig. 7 to running trend analysis (Brunetti et al. 2006b ). The results are shown in Fig. 8. They give evidence that the trends of this curve are much more influenced by CDD trends than by HDD trends. The net effect of the weather related terms in relation (3) reflects therefore very closely the evolution of the CDD record, with slightly lower positive trends due to the fact that the HDD decrease causes a lower electricity demand. 

Beside the analysis performed at yearly resolution, it is also interesting to study how climate change affects the present-time monthly profile of the electricity demand. For this reason, we plot the monthly average values of the CDD, HDD and S terms of relation (3) for 3 different 40 year periods and we add to the plot also the corresponding values for the 2001-2013 period (Fig. 9 ). As well as for the yearly analysis, all values were calculated with γ k values that are representative of present-time conditions. It is evident that, as expected, changes in air temperature, as those that occurred in the past, substantially modify the yearly profile of the demand, lowering winter month contributions and raising summer month contributions. The results also indicate that summer months are those with the highest impact of climate change on present-time energy demand, whereas in May and October present-time electricity demand seems to respond in a very limited way to temperature changes as those that occurred in the past.
section_name: Long-term evolution of CDD, HDD and S and corresponding impact on Italy present-time electricity demand

section_content: A methodology was developed for estimating secular daily temperature records (T n, T m and T x ) for any urbanised grid-point of a high-resolution DEM covering Italy. It is based on the superimposition of two fields: The root mean square error of the climatologies turns out to be between 0. 7 °C and 1. 4 °C, with higher errors for T n and lower errors for T m. The error of the anomalies is generally higher, with values strongly depending on data availability of the considered period. In order to keep the root mean square error of the T m anomalies within 2 °C, we decided to avoid reconstructing the anomalies for central and southern Italy before 1880, whereas for northern Italy we reconstructed them from 1801. It is also worth noticing that the anomalies of high-elevation points have rather high errors. This problem is not relevant for the application presented in this paper as a very high fraction of Italian urbanised areas is at low elevation. It has however to be addressed in case of study areas which include a significant fraction of high-elevation grid-points. A better representation at those grid-points is likely to be obtained by using both a higher density network of anomaly records and a 1920 1921-1960 1961-2000 2001-2013 Fig. 9 Monthly averages of weather-driven contributions to presenttime Italy electricity demand evaluated using CDD, HDD and S of 3 consecutive 40 year (1881-1921, 1921-1960, 1961-2000). In addition, we also plot CDD, HDD and S contributions relative to the years 2001-2013 in order to highlight the impact on present-time demand of the recent temperature rise As in Fig. 5, we show also a low-pass Gaussian filter methodology that allows to take into account the elevation-dependence of the temperature anomalies in a more sophisticated way than by means of the elevation weight used in equation (1) (see e. g. (Frei 2014) ). 

The estimated daily temperature records allowed investigating past variability and change of temperature-derived variables that influence Italy present-time national electricity demand. They are national averages over Italian urbanised grid-points of cooling degree-days (CDD), heating degree-days (HDD) and solar radiation deficit with respect to a defined threshold (S). CDD and HDD yearly totals give evidence of strong and highly significant trends. Specifically, CDD increase across the last 2 centuries, with 30 year normal values increasing from about 200 degree-days at the beginning of the 19 th century to over 400 degree-days in the last years. On the contrary, HDD 30 year normals decrease in the examined period from about 1500 degree-days to about 1100 degree-days. 

The strong trends of these variables makes it relevant to compare the present-time Italy electricity demand with the virtual demand Italy would have at present if the climate were still in the situation prior to global warming. Assuming that the present-time dependence of the electricity demand from meteorological variables is independent from climate change, this comparison can be used to assess the impact of climate variability and change on present-time Italian electricity demand. The results, which were obtained by means of an electricity demand model discussed in Scapin et al. (Scapin et al. 2015), give evidence, at a yearly scale, of an impact of global warming of about 5 TWh, which correspond to 1. 5-2. 0 % of Italy electricity demand. It results from the increase of CDD which is partially counterbalanced by the decrease of HDD. The impact is highest in summer months where it exceeds 5 % of Italy average monthly electricity demand, whereas the decrease of the winter demand is rather low because of a very limited use of electricity for heating. 

Even though the impact of climate change we found is rather small with respect to long-term electricity demand changes driven by socio-economic factors, assessing it may give an interesting contribution in the context of the evaluation of the impact of climate change on the energy sector. This assessment, together with a corresponding assessment based on future climate scenarios (see Scapin et al. (Scapin et al. 2014), may be considered within climate change costs estimation. 

In the future, it may be interesting to extend the methodology presented in this paper to investigate the impact of HDD variability and change on total energy requested for heating by considering, in addition to electricity, natural gas, fuel oil and other fuels. Moreover, it may be interesting considering the different climatic areas of Italy, using specific relations between the meteorological variables and the electricity demand for each area. Such an approach is however not yet possible for us, as the demand data at our disposal are aggregate over the entire Italian territory.
section_name: Conclusions

section_name: Acknowledgments

section_content: We acknowledge all the institutions that provided the meteorological data used in this work. Specifically, we refer to those listed in Brunetti et al. (Brunetti et al. 2006a) and in Simolo et al. (Simolo et al. 2010), and we mention, for recent updates, US National Oceanic and Atmospheric Administration -National Climatic Data Center, NESDIS -Global Surface Summary of the Day, Asheville, NC (Available at https://data. noaa. gov/dataset/global-surface-summary-of-the-day-gsod). The present work was partly funded by the Research Fund for the Italian Electrical System under the Contract Agreement between RSE S. p. A. and the Ministry of Economic Development -General Directorate for Nuclear Energy, Renewable Energy and Energy Efficiency, stipulated on July 29, 2009 in compliance with the Decree of March 19, 2009 and partly by the EU FP7 project ECLISE ( 265240 ).
section_name: Acknowledgments

section_name: S

section_content: Fig. 5 Italian CDD, HDD and S yearly series (thin lines), plotted together with 10-y standard deviation Gaussian low-pass filters (thick lines) y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 y = 62. 64 + 1. 09 ⋅ x, r 2 = 0. 96 HDD (North) HDD (Italy) y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90 y = 280. 12 + 0. 78 ⋅ x, r 2 = 0. 90
section_name: S

section_content: The authors declare that they have no competing interests.
section_name: Competing interests

section_content: All authors participated in the conceptualization of the work, contributed in drafting parts of the manuscript and approved its final version.
section_name: Authors' contributions