1 Introduction

1.1 Wind forecast and future scenarios

Prediction of wind speed is relevant in several crucial fields. It is vital for issuing of alerts, to prevent potential catastrophes represented by hurricanes (Pant and Cha 2019) and storm surge (Barbariol et al. 2022), and limit economic damages. Forecast of wind speed is useful in the air/sea transport sector, and one can e.g. reduce fuel consumption and CO2 emissions (e.g. Prats et al. 2022; Sun et al. 2022). In agriculture, wind forecast helps scheduling spraying, and dusting procedures (Zhang et al. 2017). In wind power production, a forecast with lead times of even few second may help to guarantee the stability of the electrical grid (Peng et al. 2016). From the times of Romans to WWII, forecasting wind patterns was even a key to success for military operations (e.g. Lowe 2000, Gross 2022). Thanks to the availability of large datasets, machine learning, or similarly conceived models could soon outdate physical modelling (Schultz et al. 2021). However, weather models are quite limited on the temporary scale, ca. 10 days (Zhang et al. 2019), making them unsuitable to predict long term wind variations, which is again a key variable.

In turn, feasibility studies of wind farms (Koletsis et al. 2016), insurances (Collier et al. 2021), and environmental policies (Zhao et al. 2021) require prediction of medium/long term wind fields, that can be affected by climate changes, cascading into changes in weather, and wind trajectories. For this reason, the IPCC panel uses Global Circulation Models (GCMs) that physically simulate multiple phenomena, including wind, and project climate scenarios based on various Social Economic Pathways (SSPs) until 2100.

1.2 GCMs downscaling

Because they operate at the global scale, GCMs have a large space resolution, ca. 200 km, which is mostly unsuitable for local studies. Indeed, climate variables, and particularly wind speed, can vary significantly, even in flat areas, within the span of few kilometres (Hubbard 1994). To tune the output of GCMs to local weather condition, one can implement downscaling, which provides weather outputs at a given/chosen resolution. Several techniques of downscaling were developed hitherto, and can be roughly folded into two main approaches, i.e. dynamic, and statistical (e.g. Keller et al. 2022).

Dynamical downscaling uses outputs from GCMs as input to regional/local climate models, in turn providing in output regional/local climate variables. Accordingly, dynamic downscaling considers the complexity, and topography of the study area, and may provide better results in terms of spatial and temporal variability (Vaittinada Ayar et al. 2016), at the cost of computationally extensive /expensive simulations.

Statistical downscaling on the other hand uses local climate series, to modify GCMs outputs, basing on statistical, data driven relationships. These statistical rules are then applied for future scenarios to produce local climate series from large scale projections.

While being faster, statistical downscaling requires long as possible historical data to provide robust parameterization. Both methodologies, and hybrid versions of the two, were used to assess future wind scenarios (e.g. Nolan et al. 2012; Gonzalez-Aparicio et al. 2017; Reyers et al. 2015). While dynamical downscaling models may provide consistent results to assess overall wind variations on a regional scale, at point scale, due to the extreme spatial variability of wind, a statistical model may provide better results.

Here we analysed daily wind speed, from 10 climate stations during 25 years in Lombardy region, Italy, and we projected future scenarios of wind using a statistical approach (namely, Stochastic Time Random Cascade STRC, properly introduced here, and tuned locally). Downscaled daily values were then aggregated in the three variables of interest: average wind speed, number of days per year of calm wind, and 95% quantile wind speed. To these variables we then applied test statistics (Linear Regression, Mann Kendall, Moving Window) to detect possible trends or non-stationary processes.

1.3 Study area

Lombardy is the most populated region of Italy, with almost 10 million people located in a 23,863 km2 area. It is also the largest industrial area of the country, and one of the largest in Europe. The region is almost equally divided between flatlands, and mountain regions, and the latter provides large water resources, used for industrial and agricultural purposes, in the Po Plain. Landlocked by the Alps in the North, and the Apennines in the South, wind recirculation in the Po Plain is mainly occurring from the Adriatic Sea on the East, leading to low wind velocity values. Wind power production in Lombardy is very modest, with only 12 wind farms, producing 0.06 MW (Terna 2023), but that is only a minor consequence of the scarce wind velocities in the region. The combination of the latter, with high urbanisation and industrialisation, makes Lombardy a most polluted area in Europe (Maranzano 2022), with several deaths yearly attributed to PM10 concentration (Baccini et al. 2011). Still, Po Plain is also the area of Italy most affected by tornadoes (Miglietta and Matsangouras 2018), as witnessed in July 2023, when a wind storm in Milan reached velocities over 100 km/h (Arpa 2023) causing severe damages.

2 Data and methods

2.1 Wind stations data

The wind speed (WS) observed data analysed in the study were measured by automatic weather stations operated by the Regional Environmental Protection Agency, Arpa Lombardia, which provides WS data every 10, or 60 min in more than 500 stations, since 1989. However, only 230 of such station operated for more than 10 years, and all of them have frequent data gaps. Here, we selected 10 stations, with the longest periods of data collection, and the fewer gaps. These stations display a wide range of altitude, from 96 to 2108 m a.s.l. (Fig. 1; Table 1). Each station was assigned to a specific elevation class based on criteria outlined by the Central Statistics Institute of Italy (ISTAT 2024), distinguishing between plain (below 300 m a.s.l.), hill (300 to 600 m a.s.l.), and mountain (above 600 m a.s.l.) regions. The elevation class is not strictly determined by the station’s elevation, but rather by the predominant elevations of the entire municipality. This decision was made to account for the topographic characteristics of the surrounding environment in our analysis, as these factors significantly influence the local wind patterns.

Fig. 1
figure 1

(a) Map of Lombardy region with considered wind station clustered for elevation class. (b) Data availability for each station

Full size image
Table 1 ID and features of the weather stations
Full size table

Starting from daily average data, we computed at monthly scale: the average WS, the number of days of calm wind (NCW), i.e. with WS < 1 m/s, and the 95th quantile (Fig. 2). We chose these additional indicators because Lombardy is both (i) the Italian region with lowest WSs, and thereby highest NCW, and (ii) one of the regions most affected by windstorms (Miglietta and Matsangouras 2018). The station with the highest WS values is CA, in the mountains, while lowest values are seen in BG and CP, on hilly, and flat turf, respectively. Seasonality of WS is evident only for DB station, while more relevant cyclical behaviour is detected for NCW, where 4 stations (AR, TR, VI, ED) similarly show lowest values in summer.

Fig. 2
figure 2

(a) Average, (b) 95th quantile and (c) number of days of calm wind computed at monthly scale for each station

Full size image

As a further preliminary analysis, we assessed spatial correlation of WS through Pearson coefficient. In Fig. 3 it is shown (a) the correlation of WS between stations, against horizontal distance, (b) the minimum, average and maximum correlation values for each station against others. As expected, correlation is strongly related to distance, and higher for stations within the Po Valley (BG, CP, AG, AR, TR). This seemingly indicates an overall homogeneity, and generally calm wind of the area. On the other hand, the stations placed in the mountains (CA, ED, DB) show lower correlation values, likely because wind fields depend upon local microclimate conditions. It seems clear here that WS is a complex variable, not only affected by topography, but likely by other important factors, playing a role in wind circulation.

Fig. 3
figure 3

(a) Pearson correlation values of WS for each station as a function of the distance. (b) Minimum, average and maximum correlation value for each station with all the others. Colour code represents station considered in the comparison

Full size image

2.2 Global circulation models

Six Global Circulation Models (GCMs), EC-EARTH3, CESM2, MIROC6, CMCC-CM2, MPI-ESM, and HadGEM3, were considered. These come from the Assessment Report 6 of the Intergovernmental Panel on Climate Change (IPCC). We assessed potential wind scenarios in Lombardy according to two Shared Socioeconomic Pathways: SSP1-2.6 (cut to net zero CO2 emissions around 2075 and estimated warming of + 1.4 °C by 2100) and SSP5-8.5 (CO2 emissions triple by 2075 and estimated warming of + 4.4 °C by 2100), giving somehow the expected extreme conditions in term of potential emission patterns, and climate evolution thereby. The GCMs considered have spatial resolution between 80 and 130 km. A first comparison between ground observations, and GCMs showed a low correspondence (Fig. 4). Specifically, this discrepancy is seen as:

  1. i)

    A discernible bias in the calculated averages (Fig. 4a), with most GCMs tending to overestimate the wind speed (MPIESM being the one with the highest bias). Conversely, HADGEM3 and MIROC6 display relatively smaller biases, but they exhibit a tendency towards underestimation, particularly pronounced when considering station CA, which ranks second in altitude and records the highest wind speed values (both in terms of average and the 95th quantile) and exhibits the lowest NCW values (Fig. 2).

  2. ii)

    Distinct differences of the coefficient of variation (CV). Most GCMs exhibit a negative bias in their CV values (Fig. 4b), indicating a smoother signal compared to the observed data. This disparity in wind variability is a somewhat common occurrence when comparing GCMs with local observational data, likely attributable to coarse spatial resolution of the former. However, it is noteworthy that two specific models, MIROC6 and MPI-ESM, exhibit a higher level of variability with respect to five, out of the ten monitoring stations. Station DB, in particular, stands out displaying lowest variability, despite being located in a mountain region, with high wind speed values, second only to the high-altitude station CA.

  3. iii)

    Substantial differences in the autocorrelation structure (Fig. 5). In the case of 8 out of the 10 monitoring stations, the GCMs, with the exception of MIROC6, display autocorrelation values with very rapid decline, when compared against observed data. GCMs thus fail to faithfully capture periodicity in wind speed variations.

To further examine the alignment between GCMs and observed data, we employed Pearson correlation coefficient as a quantitative measure. Correlation coefficients were calculated between observations, and GCMs data within the corresponding geographical tiles, at daily (Fig. 6a), and monthly scale (Fig. 6b). Correlation values are notably low, particularly at daily scale, failing to attain values above 0.1 for any combination of GCM and monitoring station. Considering monthly values, correlation is generally better, with a substantial improvement for HADGEM3, and yet not good enough to consider GCMs data as well representative of local (measured) values.

Collectively, the findings from our analysis thus far demonstrate that the GCMs under scrutiny cannot generate daily WS time series aligning closely (in a statistical fashion) to the observed data. Therefore, one has to adjust GCMs outputs, to better conform to local wind conditions in the study area, which is pivotal to understanding wind fields’ dynamics within the study area.

Fig. 4
figure 4

(a) Difference of average WS between GCMs estimates and observation. (b) Difference of the coefficient of variation (CV) between GCMs estimates and observation. Positive values point to overestimate of WS/variability of GCMs with respect to observed data

Full size image
Fig. 5
figure 5

Comparison between autocorrelation of observed WS (black line) and simulated WS by GCMs: uncorrected (red line) and downscaled (DS, green line), for the 10 monitoring stations

Full size image
Fig. 6
figure 6

Pearson linear correlation between (a) daily observed WS and daily simulated WS and (b) monthly observed WS and monthly simulated WS

Full size image

2.3 The Stochastic Time Random Cascade downscaling

To correct the discrepancies in the GCMs data, we employed a Stochastic Time Random Cascade (STRC) algorithm. This method has been widely used in the field of meteorology, primarily for modelling the spatial distribution of meteorological variables, with a notable focus on precipitation (Schertzer and Lovejoy 1987; Gupta and Waymire 1993; Over and Gupta 1994; Deidda 2000). Moreover, it has been specifically applied for spatial downscaling of GCMs precipitation outputs (Bocchiola 2007; Groppelli et al. 2011). Importantly, the STRC algorithm structure allows to extend its application to also model the temporal distribution of a variable, as we have done in this study with daily WS values.

The STRC algorithm is based on a branching tree structure (Fig. 7), where each layer of the tree represents the variable at a different scale (e.g. Bocchiola and Rosso 2006). The node positioned at the top of the tree, namely the root node, represents the coarsest resolution data (i.e. the input of the algorithm). Moving downwards, each subsequent layer is characterized by a progressively finer resolution, ultimately ending in the nodes at the bottom, which corresponds to the variable at the finest resolution (i.e. the output of the algorithm).

This method enables to increase the spatial or temporal resolution of a signal while concurrently preserving the spatial or temporal correlation structure observed in the measured data. In our specific application, it allowed us to downscale GCMs data from monthly scale, which we considered to be more consistent than daily data, to synthetic daily series, where it is preserved the time-wise autocorrelation structure identified in the observations.

2.4 STRC algorithm structure

Within the framework of the STRC algorithm, the first step requires to correct the GCMs data with a month specific correction factor \({\alpha }_{i}\) obtained by the following:

$${\alpha }_{i}=\frac{{\stackrel{-}{WS}}_{OBS,i}}{{\stackrel{-}{WS}}_{GCM,i}}$$
(1)

Where \({\stackrel{-}{WS}}_{\text{O}\text{B}\text{S},\text{i}}\) is the average value of the WS observations of the month \(i\), and \({\stackrel{-}{WS}}_{GCM,\text{i}}\) is the average value of the GCM WS data for the same month.

Then, the values in each layer are generated, by multiplying the corresponding value of the preceding layer by distinct multiplicative factors (MFs). The MFs are stochastically drawn from a probability distribution. Notably, the parameters of the distribution vary for each layer-to-layer downward step, and are estimated from the observations. In this study, we decided to employ the lognormal distribution to model the MF. This fits with established practice in the application of the STRC algorithm, basing on the following consideration:

  1. i)

    The product of two variables, each following a logarithmic distribution, is itself described by a logarithmic distribution. This property applies well in the field of meteorology, as numerous variables have been empirically observed to exhibit lognormal distribution characteristics (Baran and Lerch 2015; Bocchiola 2007; Tustison et al. 2002; Marsan et al. 1996; Tapley and Waylen 1990).

  2. ii)

    The two parameters of the lognormal distribution (i.e. m, s), are directly linked to the average value of the distribution. Consequently, if the average value is known, the two parameters are uniquely interrelated, reducing to one (i.e. variance) the number of calibration parameters.

The choice of the lognormal distribution was also tested against another typically adopted distribution, i.e. Weibull (Ozay and Celiktas 2016; Pishgar-Komleh et al. 2015), the former giving better performance (see Sect. 3.1).

The generic WS value i at the scale s can be computed starting from the value at the root node as follows:

$${WS}_{s}^{i}= {WS}_{0}\prod _{j=1}^{s}{MF}_{j}^{i}$$
(2)

Where WS0 is the input value (i.e. the data at the coarsest resolution) and \({MF}_{j}^{i}\) is the multiplicative factor for the value i at the layer j. The statistics of the lognormal distribution of the factor y are the following:

$$E\left[{MF}_{s}\right]=1$$
(3)
$$VAR\left[{MF}_{s}\right]= {\vartheta }_{{MF}_{s}}^{2}$$
(4)

Equation (2) yields the average conservation of the mass in the cascade and it is a direct consequence of the application of the correction factor in Eq. 1. Since the expected value of a lognormal distribution is given by:

$$E\left[{MF}_{s}\right]={e}^{{\mu }_{ln{MF}_{s}}+\frac{1}{2}{\vartheta }_{{lnMF}_{s}}^{2}}$$
(5)

From Eqs. (2) and (4) we can write the parameters of the lognormal distribution as follows:

$${MF}_{s}\sim LN(-\frac{1}{2}{\sigma }_{lnMF}^{2},{\sigma }_{lnMF})$$
(6)

As mentioned before, from Eq. (5) one only needs to estimate one parameter, \({{\upsigma }}_{\text{l}\text{n}\text{M}\text{F}}\). We employed an iterative estimation algorithm, identifying \({{\upsigma }}_{\text{l}\text{n}\text{M}\text{F}}\), by variance fitting, i.e. minimizing the difference of the standard deviation between the observation and the GCMs downscaled data, at each time scale (i.e. for every layer of the tree structure).

In this study, we apply the algorithm for each station, by starting from the mean monthly value of the corresponding GCM cell, which is then downscaled at daily scale. The tree structure consists of four layers (branches), then split into two, three and five new branches, respectively (see Fig. 7). The upper layer represents the monthly average data from the GCM, while the layers below represent 15 days, 5 days and daily averages. Every value at each layer is obtained multiplying one value of the layer above by a MF. Due to the inherent structure of the algorithm, it follows that values derived from the same “parent” branch exhibit a higher degree of correlation with each other, compared to values obtained from a different “parent” branch. For the same reason, values originating from neighbouring “parents” display stronger correlations than values originating from “parents” farther apart in the tree structure. This structure aids production of downscaled data characterized by autocorrelation patterns that diminish as the time lag increases, with a rate that closely approximates the observed autocorrelation behaviour. Note that the calibration of \({{\upsigma }}_{\text{l}\text{n}\text{M}\text{F}}\) is done for every station individually, since ground based WS data showed a marked difference in the autocorrelation structure.

Fig. 7
figure 7

Tree structure of the STRC algorithm employed to downscale monthly average WS data from the GCMs to a daily time scale

Full size image

2.5 Trend analysis tests

To pursue a comprehensive trend analysis on the corrected GCMs data, we employed three distinct trend assessment methods:

  1. i)

    The Linear Regression F-test was utilized to investigate the presence of a linear trend within the time series data.

  2. ii)

    The Mann-Kendall test (Kendall 1975), a non-parametric test, was applied to examine whether there exists a monotonic trend in the data, without making any assumptions about the rate of the trend.

  3. iii)

    The Moving Window Average test (De Michele et al. 1998) was employed to assess whether the data series exhibit stationarity or exhibit significant deviations from the mean. This was accomplished by comparing the moving average over 30-year periods (WMO, 2017).

These trend tests were pursued on the same variables previously analysed (i.e. yearly averages, 95th quantiles and NCW). The data series subject to testing pertained to the future projections of WS under the SSP1-2.6 and SSP5-8.5 scenarios, during 2015–2100. These projections were downscaled for each monitoring station. All tests were executed with a 95% confidence level.

3 Results and discussion

3.1 Distribution of observed multiplicative factors

The initial assumption requiring evaluation pertains to the distribution of the MF, derived computing the ratio between daily WS observations and corrected GCMs data (i.e. multiplied by the factor \({\alpha }_{i}\)). Quantile-Quantile (QQ) plots in Fig. 8 (Marden 2004), compare quantiles of the natural logarithm of observed MF values (Y axis), against those of a normal distribution (X axis). For most stations, MF values are well approximated by a lognormal distribution. However, there are exceptions, particularly evident in the higher MF values of specific stations (e.g. CP, VA, and ED). Notably, the QQ plot shapes suggest that these stations exhibit MF value distributions with a fatter right tail, as compared to a lognormal distribution. Consequently, for such stations, we expect a decreased reliability of the STRC algorithm for high WS values. The insights from the QQ plots find confirmation in the Goodness of Fit tests conducted on the MF (Fig. 9). A substantial number of stations pass the Chi-Squared test with the majority of the GCMs. However, fewer GCMs pass the Kolmogorov-Smirnov and Anderson-Darling tests, primarily because these tests place greater emphasis on the tails of the distribution, which corresponds to a vulnerability in the MF of specific stations. In conclusion, we believe these findings provide enough justification for employing a lognormal distribution to characterize the bulk of the MF distribution, while showing potential limitations in representing high-extreme WS events.

Fig. 8
figure 8

Quantile-Quantile plots of the MF computed between WS observations and GCMs data, for each station

Full size image
Fig. 9
figure 9

Number of GCMs that passes the goodness of fitness tests, divided by test (Chi squared, Kolmogorov-Smirnov and Anderson-Darling) and by station

Full size image

3.2 STRC downscale calibration

The STRC algorithm is calibrated for each weather station and every GCM. \({{\upsigma }}_{\text{l}\text{n}\text{M}\text{F}}\) is computed for every layer of the tree structures with an iterative algorithm, which minimizes the difference on the standard deviation between the observations and the downscaled GCMs data, at every time scale (i.e. every layer of the tree). In this study \({{\upsigma }}_{\text{l}\text{n}\text{M}\text{F}}\) was computed specifically for individual months. This decision was made based on the observation of pronounced seasonality in WS data across multiple weather stations, that may imply different distributions parameters of the MF. As expected, the estimated \({{\upsigma }}_{\text{l}\text{n}\text{M}\text{F}}\) increases as time scale decreases (Fig. 10). Notably, a clear time dependence is visible for 5-days average and, more evidently, for daily averages, showing lower variability of the MF during summer months, and higher variability during winter months, which suggests a corresponding behaviour of the variability of the WS (confirmed by observations, Fig. 2).

Fig. 10
figure 10

\({\sigma }_{lnMF}\) computed for each month, for the different layers of the STRC algorithm. The value of the parameter and its monthly changes increase with the temporal scale

Full size image

3.3 Downscaling output

The adjusted WS time series were compared with the observed data, revealing a significantly improved alignment in contrast to the unadjusted ones. Specifically, three key observations can be made:

  1. i)

    Discrepancies in the total average values become negligible. The use of a correction factor based upon average values results in deviations from observed averages not exceeding 10%.

  2. ii)

    Enhanced representativity is evident when considering the coefficient of variation (CV) of WS (Fig. 11). Deviations from the observed CV values do not exceed 20% for 8 out of 10 stations. Conversely, station VA exhibits a pronounced deviation from the observed CV, meaning likely that the downscaled GCMs data display greater variability, compared against the observed data. This outcome could be caused by little adherence of MF in this station to a lognormal distribution (especially on the right tail, see Figs. 8 and 9).

  3. iii)

    Enhanced fit is observed in the autocorrelation patterns between the observed and the downscaled GCMs data (Fig. 5). This finding is particularly relevant within the scope of WS analysis, as no downscaling methods that we know of deals with for this facet. We suggest that this circumstance is noteworthy when carrying out WS analysis (e.g. proper representation of autocorrelation in time may be key for short term forecast when using the STRC, e.g. Bocchiola and Rosso 2006).

Fig. 11
figure 11

Error on the coefficient of variation (CV) computed subtracting the CV of the observed WS to the CV of the corrected GCMs data. Compared to the uncorrected GCMs data the error on the CV is lower for almost all stations

Full size image

3.4 Trend tests results

A comprehensive assessment of the time series downscaled data was conducted for two distinct scenarios, SSP1-2.6 and SSP5-8.5. The results of the trend analyses show as follows:

  1. i)

    Absence of statistically significant trends in the annual average wind speeds (Table 2). However, it is noteworthy that a considerable portion of the models exhibited non-stationary behaviour in certain weather stations. This phenomenon was more pronounced in the SSP5-8.5 scenario, maybe indicating an amplification of non-stationary behaviour of wind patterns in a higher emissions scenario.

  2. ii)

    Substantial positive trends were identified in the 95th quantile of wind speed across all models (Table 3). Specifically, these trends indicated an average increase along the century, ranging from + 0.24 to + 0.39 m/s (from + 10% to + 17%), for the SSP1-2.6 scenario. Remarkably, one model displays a negative trend under SSP1-2.6. For the SSP5-8.5 scenario, the trends in the 95th quantile varied between + 0.15 and + 0.46 m/s per 100 years (from + 6% to + 20%), with all models displaying positive trends.

  3. iii)

    Consistent non-stationary behaviour of NCW values across all models and nearly all weather-stations (Table 4). Furthermore, the application of the Mann-Kendall test and the linear regression test identified significant trends in NCW for specific stations. However, the direction of these trends remains ambiguous: it is not clear whether the NCW is projected to either increase (as suggested by 5 models for the scenario SSP1-2.6 and 2 models for the SSP5-8.5), or decrease (as suggested by one model for the SSP1-2.6 and 2 models for the SSP5-8.5).

Table 2 Results of the trend tests (Linear Regression, LR, Mann-Kendall, MK and Moving-Window, MW) on the annual average of WS for the two scenarios (SSP1-2.6 and SSP5-8.5). For every model, the number of stations passing the test is indicated. The coefficient of the linear regression is indicated when the corresponding test is passed
Full size table
Table 3 Results of the trend tests (Linear Regression, LR, Mann-Kendall, MK and Moving-Window, MW) on the 95th quantile of WS for the two scenarios (SSP1-2.6 and SSP5-8.5). For every model, the number of stations passing the test is indicated. The coefficient of the linear regression is indicated when the corresponding test is passed
Full size table
Table 4 Results of the trend tests (Linear Regression, LR, Mann-Kendall, MK and Moving-Window, MW) on the annual number of calm wind for the two scenarios (SSP1-2.6 and SSP5-8.5). For every model, the number of stations passing the test is indicated. The coefficient of the linear regression is indicated when the corresponding test is passed
Full size table

3.5 Wind scenarios in Mediterranean area

Overall, our trend analysis results highlight the considerable uncertainty, and variability of projected wind speed scenarios in the Mediterranean area, where Lombardy region is laid. Trend’s magnitude, direction and location are dependent upon the area, and influenced by the chosen GCM and downscale method, often with opposite results. Notably, the main source of uncertainty typically arises from variance between GCMs (Koletsis et al. 2016).

Future scenarios in the Mediterranean region show a predominantly negative trend in WS. Donat et al. (2011), using Climate Models prior to 2007, found loss of wind storms (98th quantile) for the Mediterranean region. Similarly, Tobin et al. (2015), and Hueging et al. (2013) found a negative trend for wind energy production in the Mediterranean area. In another recent study (Bonanno et al., 2023), specific for Italy, a negative trend of wind power production was found, under the sole RCP 8.5 scenario. More specifically, the trend was detected in coastal regions, where wind power plants are located, while in inner areas, and particularly in the North, the trend is negligible. Some other scholars (Bloom et al. 2008; Kjellström et al. 2011) reported future negative trends, limited to the winter season, whereas in summer positive trends are found. This may be explained by alteration in large scale circulation patterns, e.g. resulting from a rise in anticyclonic impact, and a fall in cyclonic impact, while summer positives trend can be attributed to growth of air temperature.

However, numerous studies have reported highly uncertain results, also in regions close to case study area. Medugorac et al. (2021), when studying with RCMs storm surges from the Adriatic Sea, the main source of air recirculation in Po Valley, found small probability of change in wind fields. Often times, significant differences are found between GCMs, making difficult to detect a consistent trend (e.g. Rockel and Woth 2007; Outten and Esau 2013). Pryor et al. (2006) studied the 90th percentile trends in future scenarios, finding fluctuations of about +/- 15% at the end of 2100. However, the magnitude of these trends was no bigger than the scatter between GCMs. The divergence between GCMs in the same study resulted in an inability to predict future storminess.

4 Conclusions and outlook

In this study we applied a statistical method to correct the bias of WS data from GCMs. We kept the monthly signal of the GCMs, and through a downscaling technique we obtained daily WS data, resembling statistical properties of the observations (i.e. average, variability and autocorrelation). The downscaled timeseries demonstrated highly accurate replication of the observations, with no significant variations observed among the GCMS, but with some variations noted among the weather stations.

The downscaling algorithm was used to project future WS scenarios (namely SSP1-2.6 and SSP5-8.5) during 2015–2100. To detect potential implications on wind circulation resulting from climate change, trend tests were pursued using the projected daily WS data. Even though average annual wind speed exhibits non-stationary dynamics, particularly in the context of the highest emissions scenario (SSP5-8.5), no significant trends were evidenced. The 95th percentile of wind speed exhibited consistent, and significantly positive trends in most models for both scenarios, maybe calling for a future rise in extreme wind speeds. NCW values exhibited non-stationary behaviour and significant trends in select stations. Nevertheless, future NCW indicator remains uncertain.

The findings of this study suggest that climate change is expected to affect the circulation of wind in Lombardia region. However, it remains challenging to predict the evolution of WS quantitatively, by the end of the 21st century. Due to the contrasting behaviour of the variables, it appears that the probability distribution of WS is likely to be altered in the future. A more comprehensive investigation, incorporating a greater number of GCMs and weather stations, may result in less uncertain findings, and may be pursued henceforth. Furthermore, estimation of the STRC algorithm parameters may be tested for robustness/dependence against other variables affecting wind circulation (e.g. pressure gradient, vorticity). Finally, the examination of seasonal or monthly data separately should be undertaken, because WS clearly displays strong seasonal dependence.