Correction of ERA5 Wind for Regional Climate Projections of Sea Waves

Abstract

This paper proposes a method to infer the future change in the wind-wave climate using reanalysis wind corrected to statistically match data from a regional climate model (RCM). The method is applied to the sea surface wind speed of the reanalysis ERA5 from the European Centre for Medium-Range Weather Forecasts. The correction is determined from a quantile mapping between ERA5 and the RCM at any given point in the geographical space. The issues that need to be addressed to better understand and apply the method are discussed. Corrected ERA5 wind fields are eventually used to force a spectral wave numerical model to simulate the climate of significant wave height. The correction strategy is implemented over the Adriatic Sea (a semi-enclosed basin of the Mediterranean Sea) and includes the present-day period (1981–2010) and the near-future period (2021–2050) under the two IPCC RCP4.5 and RCP8.5 concentration scenarios. Evaluation against observations of wind and waves gives confidence in the reliability of the proposed approach. Results confirm the evolution toward an overall decrease in storm wave severity in the basin, especially under RCP8.5 and in its northern area. It is expected that the methodology may be applied to other reanalyses, RCMs (including multi-model ensembles), or seas with similar characteristics.

1. Introduction

Sea surface wind waves are an essential player in the coupled climate system [1]. The physical processes affected by waves include, but are not limited to, fluxes of energy and gas between sea and atmosphere (e.g., carbon dioxide), production of sea spray, Earth albedo, and its influence on the overall radiation budget. While wave climate analyses are relevant to scientific studies, they are fundamental to various applications including hazard assessments. For instance, wind waves have implications for coastal, near, and offshore structures’ design and operation [2] and contribute to the impacts of storm surge events [3]. For these reasons, ongoing and future changes in the wind-wave climate have received growing attention for global and regional impact studies [4]. The present climate of wind waves is commonly sourced from long-term qualified observations (in situ or satellite-borne) or model results, such as those provided by hindcast or reanalysis numerical simulations. These are complemented by climate models for projection studies extending from the historical to future periods. Despite the specific assessment requirements, reanalysis and climate model datasets can deliver valuable information on wind and wind waves.

Climate reanalyses combine past short-range weather forecasts with observations and provide a physically consistent source of parameters for the characterization of global climate and weather (e.g., [5]). Reanalyses include atmospheric, ocean, and land surface variables in a comprehensive coupled system and produce data at sub-daily interval spanning from the present to many decades back without any gap in time and space. Reanalysis data are being used in many applications, from monitoring ongoing climate change to research, policy-making, and business in renewable energy sectors. One of the most advanced reanalysis products is ERA5 [6], distributed by the Copernicus Climate Change Service at the European Centre for Medium-range Weather Forecasts ECMWF [7]. Many research centers are currently using ERA5 to release bulletins and studies on the climate’s current (present day) conditions at global and regional scales. (see, e.g., [8]). To maximize the benefits that ERA5 offers to the scientific community, ERA5 meteorological variables can be bias-adjusted using observations (e.g., [9]). For wind-wave studies, ERA5 products were used with different approaches, which include the analysis of the native wave products [10], the use of ERA5 wind to force a high-resolution wave model [11,12], or downscaling wind over a limited area [13].

Climate models provide an alternative numerical approach to modeling climate and weather systems at global and regional scales. The classical distinction is between fields produced by coarse resolution Global Climate Models (GCMs) and those by refined Regional Climate Models (RCMs), which were essentially developed to downscale over a limited area climate from GCMs [14]. GCM can describe the response of the climate to large-scale forcings, such as those due to greenhouse gases. At the same time, RCM can refine the large-scale information by accounting for the effects of sub-GCM grid-scale forcings and processes, such as those due to complex topography, land-sea boundaries, or dynamical processes occurring at the mesoscale (see the review by Giorgi [15]). GCMs and RCMs reproduce the climate over multi-decadal scales, from the past to the projected future, according to greenhouse gas emission scenarios. Often initialized at the onset of the 20th century, GCMs and RCMs models provide a possible realization of the climate system, consistent in terms of reproduced physical processes and setup of the numerical solutions. Since the GCMs providing boundary conditions to the RCMs do not employ assimilation of observed data, the simulated events are not in relation to the actual ones. Therefore, RCM evaluation is generally made through climatological statistics over mid to long temporal scales (e.g., months or seasons; [16]). The end-use compares statistics for future time slices with those for the reference period to assess greenhouse gas forced change signal at the local scale. For wind-wave projection studies, climate model wind is directly used to force spectral wave model simulations [17,18].

Intending to investigate the change in the wind-wave climate, we explore the benefit of correcting ERA5 wind to simulate the present day and the projected future climate of sea states. For this purpose, reanalysis winds are adjusted to match (in a statistical sense) those produced by a high-resolution RCM. For the wind correction, we rely on the quantile-quantile matching method [19,20,21]. This method has been commonly used in the literature to bias-correct the systematic errors of climate models using local observations or reanalysis data as a reference. In particular, this approach was adopted to gain knowledge in assessing the historical wind speed climate [22] and, under certain assumptions, the future wind speed climate [23]. Correction methods are also effective when applied to wave model data [24,25]. The correction strategy is herein used to produce a quantile-specific, model-to-model transfer function determined by the relationship between ERA5 and RCM cumulative distribution functions. The result is RCM-adjusted ERA5 wind fields for the present-day period (1981–2010) and the near-future period (2021–2050). The two IPCC RCP4.5 and RCP8.5 concentration scenarios are considered [26]. Since only model data are considered in the quantile-quantile matching, the advantages and disadvantages of this novel approach are discussed in the present study. Implicitly, this strategy also allows judging historical RCM fields using the reanalysis. Furthermore, corrected winds are used to predict the significant wave height fields and related climate by forcing a spectral wave model implementation based on WAVEWATCH III^®.

The climate projections are used to simulate wind and wind-wave future changes over the Adriatic Sea (Figure 1), as an implementation of the proposed methodology. This region is challenging for the modeling of weather and climate systems, given the nature of the atmospheric systems, combined with the complex orography [27]. The Adriatic Sea is a long and narrow semi-enclosed basin, extending in the North Mediterranean Sea for about 800 km along the major axis from the south-east to the north-west, with a width of about 200 km. Gale-force winds in the Adriatic Sea are primarily of two types [28,29]. On the one hand, southeasterly “Sirocco” winds blow over the Adriatic as a basin response to large-scale weather events in the Mediterranean Sea. On the other hand, strong northeasterly “Bora” flows are driven by the complex orography of the Dinaric Alps on the east of the Adriatic and create finely structured jets with strong sub-basin scale spatial gradients. Both types of winds generate energetic sea states (i.e., large wave heights). Still, Sirocco storms may force higher and longer waves, eventually impacting the northern littoral, close to Venice (Italy). Given this complexity, we adopt an implementation of the high-resolution (about 8 km) atmospheric model COSMO-CLM [30] as a reference RCM, which was proved to provide high accuracy in the climatological reproduction of all wind regimes and sea states in the Adriatic Sea [31,32]. Other regional climate models were used to infer the wind and wave future changes over the Adriatic Sea, such as the ones available at the highest resolution (about 11 km) from the European (EURO) Coordinated Regional Climate Downscaling Experiment (CORDEX) program [33,34] and the Adriatic Sea and Coast (AdriSC) modeling suite [35]. Projection studies based on different climate models, greenhouse gas emission scenarios, and future periods suggested that climate change will likely produce a general decrease of wind-generated stormy waves in the Adriatic Sea, from offshore [29,35,36,37] to the coast [32], also adopting an ensemble approach [38]. These studies report that the change has significance for Bora waves, while for Sirocco waves uncertainties remain.

The remainder of the paper will proceed as follows. Section 2 is dedicated to presenting numerical models (reanalysis, regional climate, and spectral wave) and measured datasets used to infer the wind and wind-wave climate. This section also provides the principles of the correction method adopted to adjust ERA5 reanalysis wind using present-day and future RCM data. Section 3 is dedicated to the model data assessment with a specific focus on the characteristic of the corrected ERA5 wind and related waves. Section 4 examines the use of the corrected reanalysis to assess the projected future changes in the Adriatic Sea wind and wave climate. An application on the storms offshore Venice is provided. A discussion and summary of the study’s main results are presented in Section 5.

2. Datasets and Methods

2.1. Wind and Wave Observations

To assess the present-day wave model’s wind forcing and output, time series of measured data are considered at three observatories in the Adriatic Sea (Figure 1). In the northern sub-basin, wind velocity and significant wave height (H_s) have been recorded on a routine basis from the Acqua Alta oceanographic research platform (from now on AA; see [39] for a detailed description of the station), which is located at 17 m depth, 15 km offshore the city of Venice (coordinates: 45°18′51.288″ N, 12°30′29.694″ E). Values of H_s are taken from the publicly available set after [40], which provides verified, three-hourly, recorded values spanning almost continuously the years 1979 through 2018. The horizontal components of wind speed are measured at 15-m height above mean sea level with a VT0705B SIAP anemometer and averaged every 5 min (i.e., twelve values every hour). Observed wind speed U was corrected to the standard 10-m height reference level (U₁₀) by assuming winds occurred in near-neutral conditions. The H_s climate has been inferred further south using wave buoy data recorded by the Italian national wave network at two different locations, namely Ortona (central Adriatic Sea: 42°24′54″ N, 14° 30′ 21″ E) and Monopoli (south Adriatic Sea: 40°58′30″ N, 17°22′36″ E). Buoys were operated from July 1989 to March 2008 and moored at a sea depth of approximately 100 m. For the sake of the present study, the original time series of observed or model data (see following sections for details) were re-sampled at a common temporal interval (either 1-h, 3-h, or daily) to allow intercomparison. The reference time for all sources of data is UTC. Skill metrics uses the following indicators for the comparison between time series of model data and observations: mean error (ME, model minus observation), mean absolute error (MAE), slope of the line of best fit in a least-square sense (Slope), and cross-correlation coefficient (R). Model quantiles and probability distributions are also assessed. We remark that unlike what is usually done in the bias-correction procedures (e.g., [20]), the set of observations is herein used with the only purpose of cross-validating model outputs. In other words, these will not be used for the statistical matching procedure proposed in this study.

2.2. ERA5 Reanalysis

For long-term climate studies, uniform model datasets are necessary. This requirement is fulfilled by the ECMWF ERA5 reanalysis. ERA5, the successor to ERA-Interim (the previous 0.75°-resolution ECMWF reanalysis), provides global, hourly estimates of atmospheric variables at a native lon/lat resolution of 0.28125° and 137 vertical levels from the surface to 0.01 hPa. ERA5 presently extends from the present day back to 1979. It includes the evolution of greenhouse gases, volcanic eruptions, sea surface temperature, wind waves, and sea ice (quality-assured monthly updates of ERA5 are published within three months of real-time). For storage issues, data have been gridded to a regular lat-lon grid of 0.25°. At the time of writing, ERA5 data are available from the Copernicus Climate Data Store [41]. For the present study, zonal and meridional components of the horizontal wind speed U₁₀ were downloaded for the 30 years spanning from 1981 to 2010 over a region that encloses the Adriatic Sea and a large portion of the Mediterranean Sea surrounding the Italian peninsula.

In the open ocean, ERA5 skill for near-surface wind speeds U₁₀ shows a 20% improvement relative to ERA-Interim and a performance similar to that of current operational ECMWF forecasts [42]. However, given the relatively coarse spatial resolution, ERA5 capacities are expected to degrade in semi-enclosed basins such as the Mediterranean Sea, especially for the most extreme conditions, for which the strong spatial gradients happen to be underestimated by modeled winds used to simulate wave fields [43]. At the 99th percentile, the U₁₀ model minus observation residual can be as low as −1.5 m/s [11]. This limitation might be further amplified in narrow basins such as the Adriatic Sea [44]. Indeed, ERA5 is inherently penalized compared to high-resolution models (10 to 5 km or less) that provide not only more highly detailed geographical patterns but significantly stronger and more accurate wind speeds (see, e.g., [13]).

2.3. COSMO-CLM Regional Climate Model

To understand how the global climate might respond to changing greenhouse gas concentrations in the atmosphere, the most valuable tools are the Global Circulation Models, with forcings prescribed according to, for instance, the Intergovernmental Panel on Climate Change (IPCC) scenarios [26]. However, typical horizontal resolutions of global GCMs are around 50–80 km to 500 km, making them unsuitable for studies on a regional scale where the local climate is affected by small-to-mid scale processes not fully resolved by global models [45]. This problem is especially true for enclosed and semi-enclosed basins with complex topography, such as the Adriatic Sea, where local meteorological events occur at scales of a few kilometers [27]. In this context, the dynamical downscaling technique allows a better description of the local climate, as it employs a Regional Climate Model nested in a GCM (used as initial and boundary conditions) to generate high spatial resolution information on a specific area of interest and time scales up to centuries [46].

For the present work, climatological wind forcings in the Adriatic Sea were produced using the RCM COSMO-CLM (CCLM, [47]), the climate version of the operational, non-hydrostatic mesoscale weather forecast model COSMO-LM [48], initially developed by the German Weather Service and currently developed by the Climate Limited-area Modeling Community in the framework of the Consortium for Small-scale Modeling (COSMO). The COSMO-LM model is a non-hydrostatic limited-area atmospheric prediction model. COSMO-LM is based on the primitive thermo-hydrodynamical equations describing compressible flow in a moist atmosphere. The model equations are formulated in rotated geographical coordinates and a generalized terrain-following height coordinate. A variety of physical processes are considered by employing parameterisation schemes. The advantage of the non-hydrostatic scheme is that it can be applied to any horizontal resolution, thus allowing a better description of convective phenomena [49]. The non-hydrostatic scheme enables a better representation of the local orography and the land-sea boundary, which control local meteorological features such as the development of fetch-limited winds [50,51]. For more details on the formulation of the model and the parameterisation settings, the reader is addressed to the COSMO model documentation page [49,52,53,54].

Specifically, in the present work, we adopted a specific configuration of CCLM developed for the Italian and surrounding areas by the Euro-Mediterranean Centre on Climate Change (CMCC, Centro Euro-Mediterraneo sui Cambiamenti Climatici); more information on this CCLM configuration (here and after CCLM-CMCC) are available in [30,55]. This regional climate simulation is forced by the global climate model CMCC-CM [56], a coupled atmosphere-ocean GCM, whose ocean component is the global ocean model OPA 8.2 [57] with a horizontal resolution of 2° × 2° on 31 vertical levels. The atmospheric component is the global atmosphere model ECHAM5 [58,59] with a T159 horizontal resolution, corresponding to a Gaussian lon/lat 0.75° grid on 31 hybrid sigma-pressure levels. The coupling between the atmospheric and the ocean models is performed via the Ocean-Atmosphere-Sea Ice-Soil version 3 (OASIS3) coupler [60].

The CCLM-CMCC wind field (zonal and meridional components of U₁₀, available with a daily temporal resolution) spans the Italian peninsula and a large portion of the surrounding Mediterranean Sea (longitude range: [2° E–22.3° E]; latitude range: [35.7° N–47.6° N]; see Figure 2 for an example of two simulated sea storms), with a high spatial resolution of 0.0715° (approximately 8 km). The Adriatic Sea is wholly included, and the computational domain extends further south of the Otranto Strait to cover the North Ionian Sea. CCLM-CMCC was run with three forcing setups: (i) control run (CTR), from 1981 to 2010 to represent the present-day climate condition; (ii) representative concentration pathway 4.5 run (RCP4.5), which is described by IPCC [26] as an intermediate scenario for future climate; (iii) representative concentration pathway 8.5 run (RCP8.5), which prescribes emissions continuing to rise throughout the 21st century and is taken as the basis for worst-case climate change scenarios. In our analysis, the near-future period covering the years 2021–2050 is considered for projected changes compared to the present climate.

2.4. Spectral Wave Modeling

Wind-wave simulations were carried out with the latest version (version 6.07.1) of the state-of-the-art spectral model WAVEWATCH III^® (from now on WW3), released in April 2019. WW3 [61,62] is a third-generation wave model developed and maintained at the National Center for Environmental Prediction (NOAA/NCEP) in the spirit of the WAM model [63,64]. Since model version 3.14, WW3 evolved from a wave model into a wave modeling framework (open-source repository available [65]), which allows for easy development of additional physical and numerical approaches to wave modeling. WW3 solves the random phase spectral action density balance equation for wavenumber-direction spectra. The implicit assumption of this equation is that properties of medium (water depth and current) and the wave field itself vary on time and space scales that are much larger than the variation scales of a single wave. The key steps of the implementation phase are described as follows.

Spectral wave models require the definition of both a physical (x, y) and a spectral grid (k, θ). For the physical grid, we have used a high-resolution grid of the Adriatic Sea and the Ionian Sea (north of 39° N in Figure 1). The grid is structured (regular in lon-lat) with 361 00647 280 nodes and square cells of 0.025° side length (corresponding to a linear distance of 2 km in longitude and 2.8 km in latitude in the northern Adriatic). The model bathymetry was obtained by interpolating a high-resolution bathymetric dataset of 0.0625° (80 m in the north of Adriatic [66]) on the grid. The bathymetric dataset itself was obtained by harmonizing a very high-resolution survey (personal communication, CNR-ISMAR) and the GEBCO relief (General Bathymetric Chart of the Oceans, 330 m resolution). For the future climate runs, the water level was uniformly raised by 0.12 m, by extrapolating, without distinction between RCPs, the present observed rate of Mediterranean sea level rise (2.8 mm/year; [67]) to the 2021–2050 period. This choice does not introduce further uncertainties and it is supported by the proximity of the future time period and by the consideration that the rise is a small fraction of the water depth (always larger than 10 m, except in the more coastal regions) and has a negligible impact on wave dynamics. Furthermore, no waves were prescribed to enter the computational domain from the Mediterranean Sea. This assumption has no impact on the cross-basin Bora waves. In contrast, it necessarily leaves to the growing waves a buffer of about 150–200 km further south when the southeasterly Sirocco blows [37]. This choice guarantees that all the Adriatic Sea states north of the Otranto Strait are correctly reproduced.

The spectral grid in the wavenumber and direction domains comprises 32 wavenumbers and 36 directions. The wavenumbers are related to intrinsic frequencies via linear dispersion relation, which are geometrically distributed with an increment of 1.1 between one frequency and the next and cover a range between 0.07 and 1.34 Hz. Such a frequency range is convenient for the Adriatic Sea environment since it prioritizes a good representation of the short-wave components. The discrete directions are evenly distributed between 0° N and 360° N, with an interval of 10° and the first direction at 5° N. For the stability of the numerical method adopted, the time discretization (four timesteps in WW3) is determined automatically from the physical grid characteristics. In particular, the following time steps Dt were used: Dt(global) = 600 s; Dt(x, y) = 150 s; Dt(k, q) = 300 s, Dt_min = 10 s. The wave model has been implemented by setting, with default coefficients, the following principal parametrization [68]: wind-wave interaction and dissipation due to deep-water breaking as in WAM cycle four and updates to make the WW3 implementation as close as possible to the one in the ECMWF wave model [69]; discrete nonlinear interaction approximation; JONSWAP bottom friction formulation; higher-order propagation schemes with Tolman’s averaging technique [70]. Although possible, we decided not to calibrate the model parameters since the interest mainly focuses on comparing projected future and present-day wind sea states. However, wave model data were assessed using in situ observations for the historical period. We analyzed the hourly values of the significant wave height H_s from the 2D spectral representation in the lon-lat geographical space and temporal domain. The value of H_s has been computed diagnostically as four times the square root of the variance of the wave field and represents an integral measure of the energy of waves. The peak direction of wave propagation was also considered to identify selected storms.

2.5. Climatological Constraint Condition for ERA5 Sea Wind

This section presents the statistical method used to correct the ERA5 wind speed U₁₀ using CCLM-CMCC data as a reference. It is worth reminding the reader that by combining model data with observations, ERA5 has the advantage of reproducing at a high-temporal resolution (1 h) the realistic sequence of historical meteorological events. Therefore, the correction strategy aims to inherit the reanalysis sequence and temporal resolution of sea storms as the basis for the historical and projected future wind and wave climate characterization. The 1-h resolution allows for accurate reproduction of the wave storm peak and shape, which is not guaranteed when employing a wind forcing with a resolution of 6 h or more [37,43].

The method we used to climate-correct ERA5 wind speeds stems from the Quantile-Quantile Matching approach (from now on QQM; [19,20,21,23]), which is a statistical downscaling technique initially developed to assess the information about the regional climate. QQM is a prior probabilistic downscaling strategy with bias correction to deal with and provide local-scale cumulative distribution function (CDF) of random climate variables through a transformation applied to large-scale CDFs. QQM was used in climate change studies as a probabilistic method to downscale GCM values using local observations or model reanalysis of the present climate (for a review, see [71]). Unlike this traditional approach, the QQM method was here used with a novel strategy, aiming to correct the reanalysis to ensure that its CDF matches (in a stochastic sense) that of the climatological model datum for both the historical and the projected future periods. In the following paragraphs, the datasets resulting from the climate-corrected ERA5 are labeled as ERA5c.

The advantage of QQM is that it potentially corrects all biases in the model CDF since the transfer function is quantile-specific. This choice allows us to correct the percentiles of average and extreme wind speeds separately [22,72]. Climate wave models forced by bias-corrected winds significantly improved the simulated climate at different scales [18]. A disadvantage of posterior bias correction using the QQM method may be the physical inconsistency of the corrected variables (e.g., wind speed and air temperature), as well as the loss of spatio-temporal correlation of a given variable since each grid point is corrected independently (i.e., the integrity of the forcing is not guaranteed). While the former does not apply in our analysis (only one variable is constrained), the latter will be evaluated by assessing the historical corrected winds and related waves against measured data.

In the historical climate experiment, time series of reanalysis wind speed U₁₀ are corrected using the RCM data (constraint condition), and relying on the QQM method, for each grid point of the reanalysis, in order to account for the local wind weather and the distance from the land-sea boundary. Assuming that both the RCM and the reanalysis data are temporally stationary, this method constructs the respective CDFs and scales the reanalysis wind in order to ensure that its CDF matches that of the RCM.

Let F_E(U_10(i)), defined over the interval [0, 1], stand for the cumulative histogram value from the reanalysis order statistics U_10(i), and F_R(U_10(i)) the corresponding one for the RCM dataset bi-linearly interpolated on the reanalysis points during the same period. F_E(z) and F_R(z) are nonlinear and give the probability that the random variable Z is below or equal to a given speed value z, i.e., Pr(Z ≤ z). The basic idea to adjust ERA5 speed $U_{10}^{\mathrm{E}}$ is to find a corrected value $U_{10}^{\mathrm{EC}}$ based on the assumption that

$${\mathrm{F}}_{\mathrm{R}}\left(U_{10}^{\mathrm{EC}}\right)={\mathrm{F}}_{\mathrm{E}}\left(U_{10}^{\mathrm{E}}\right).$$

(1)

The constraint to the RCM in the present climate simulation is therefore given by

$$U_{10}^{\mathrm{EC}}={{\mathrm{F}}_{\mathrm{R}}}^{-1}\left[{\mathrm{F}}_{\mathrm{E}}\left(U_{10}^{\mathrm{E}}\right)\right],$$

(2)

where F_R⁻¹ is the inverse function of the CDF (see Figure 1 of [71] for a graphical interpretation of the procedure).

For the future period, ERA5 winds are corrected assuming we know the constraint with respect to the future, which is given by the climate model signal, as in the CDF-transform method (CDF-t) proposed by Michelangeli et al. [23]. It is assumed that there exists a transformation T allowing correcting the CDF of RCM winds into the CDF representing the reanalysis variable. Let F_Eh stand for the CDF of reanalysis for the historical period and F_Rh for the CDF of RCM outputs for the same period. F_Ef and F_Rf are the CDFs equivalent to F_Eh and F_Rh but for a future period. The transformation T: [0, 1] → [0, 1] takes the following form

T[F_Rh(z)] = F_Eh(z),

(3)

and it can be then modeled as

T(u) = F_Eh[F_Rh⁻¹(u)],

(4)

where u belongs to [0, 1]. Then, assuming that relationship in Equation (3) will remain valid in the future, i.e., that T[F_Rf(z)] = F_Ef(z), and that we know F_Rf, the desired correction of the CDF of ERA5 is written as follows:

F_Ef(z) = F_Eh{F_Rh⁻¹ [F_Rf(z)]}.

(5)

One can interpret Equation (5) as the ranked category of ERA5 speeds, which accounts for the wind speed changes in the RCM between the historical period and the future period.

From the algorithmic point of view, the QQM method has been used in the following way. At first, we assume that the empirical estimates of CDFs from ERA5 and RCM data are good estimators of the actual distributions. Moreover, since the transformation that matches hourly ERA5 with hourly COSMO-CLM wind speeds is unknown, and we instead relied on daily average quantities, we assume that the functions that constrain daily average and hourly ERA5 winds are identical. This allows retaining the correction for daily averages to scale ERA5 hourly data. Since QQM is applied to continuous model outputs, no uncertainties will arise from sparse data.

The CDF is estimated at each ERA5 grid node over the COSMO-CLM computational domain, based on the whole time series of daily average wind speeds $U_{10}^{\mathrm{E}}$. In practice, CCLM-CMCC daily average wind speed meridional and zonal components (30 years from 1981 to 2010 and 30 years from 2021 to 2050, under RCP4.5 and RCP8.5 scenarios) have been bi-linearly interpolated on the ERA5 0.25° grid. To reduce gaps in the matching, the histogram bins were set at 0.1% step, and linear interpolation is used to approximate values between percentiles. The nonlinear correction factor CF takes the functional form given by

$$CF=CF(U_{10}^{\mathrm{E}}, \mathrm{lon}, \mathrm{lat}).$$

(6)

Hence, we obtain that

$$U_{10}^{\mathrm{EC}}(\mathrm{lon}, \mathrm{lat}, \mathrm{time})=\mathrm{CF}·U_{10}^{\mathrm{E}}\left(\mathrm{lon}, \mathrm{lat}, \mathrm{time}\right).$$

(7)

Given that for the largest wind speed values the variability of percentile estimates is relatively high [73], the CF has been estimated in the range between the 1st and 99th percentile values, beyond which the CF was extrapolated (nearest neighbor interpolation). Although the choice of keeping a constant CF for very high wind speeds can constrain ERA5 values below an upper limit, it should not provide highly biased values for the rarest and more extreme events [21]. Given this uncertainty, we did not consider fitting the highest quantiles with a parametric distribution. However, the error committed is small since the CF varies little for higher percentiles. Moreover, no smoothing was applied to the empirical curves since the small-scale variability was generally within 1%.

An illustration of how QQM was applied helps elucidate the correction process. Figure 3A shows the U₁₀ quantiles at AA from CCLM-CMCC (blue markers) and ERA5 (red markers) model data against observations (OBS, years 2004–2010). One can note that CCLM-CMCC quantiles for the historical period 1981–2010 over the whole range of speeds (daily average from 0 m/s to about 19 m/s) are in good agreement with observations. This skill indicates an overall good performance of the model for all regimes. On the other hand, a large ERA5 underestimation of measured data are present at all quantiles. For the correction, once the quantile of the reanalysis is obtained, the value of the RCM data corresponding to the quantile is specified. Figure 3B, the shape of the empirical CF for the historical period (CTR) as derived from Equation (2) is shown. After matching with CCLM-CMCC, adjusted ERA5 speeds turn out to be enhanced (CF > 1). Except for the lowest percentiles, the correction is the highest for the strongest winds (where it exceeds +25%), compensating for the systematic decrease of the ERA5 performance for the more intense storms. We further note that the minimum of CF is for speeds close to the 50th percentile (about 2.5 m/s), for which we have no explanation. In the same panel, the shapes of CF at AA for the two RCP4.5 and RCP8.5 future projections are also shown. These curves and their use will be described later in this paper, but they anticipate that RCP8.5 is corrected less than CTR and RCP4.5.

3. Model Assessment

3.1. ERA5 Wind

An evaluation of the ERA5 wind speed U₁₀ against observations in the North Adriatic Sea (AA platform) is shown in Figure 4A. ERA5 performance (i.e., underestimation) heavily depends on the adopted threshold. The model minus observation error is smaller (i.e., larger in absolute value) than −5 m/s above 20 m/s. Overall, the mean error ME is −0.59 m/s, the mean absolute error MAE is 1.49 m/s, and the slope of the line of best fit is 0.79. At the highest percentiles, the data scatter presents a unimodal shape compared to observations, which one can interpret as evidence that Bora and Sirocco storms are underestimated by a similar amount. Although they are representative of a single location, these poor scores suggest that without a correction, ERA5 data cannot be used for a reliable characterization of the most severe wind and wave regime over the Adriatic Sea. However, the high level of dependence between ERA5 and observed speeds (the correlation coefficient R is 0.81) and the good representation of the wind directions (Figure 4B) suggest that conveniently scaled, ERA5 is appropriate to provide realistic forcings to a spectral-wave model. Regarding the wind direction, there seems to be a clockwise rotation of about 10° of ERA5 winds compared to measurements, for which a misalignment with respect to the true north might be present. Both model and measured data highlight that orography has a substantial impact on local winds, being clear the dominance of northeasterly cross-basin Bora and southeasterly along-basin Sirocco regimes.

3.2. CCLM-CMCC Wind

The simulated climatic scenario for the historical period 1981–2010 has been evaluated by comparing the wind statistics (climatology and probabilities) with observations. Regarding the climatology of the model, Figure 5 shows the annual cycle (months from 1 to 12) of the monthly average and monthly 99th percentile of U₁₀ computed by CCLM-CMCC (blue line) and derived from observations (black line). In agreement with previous studies based on RCMs in the Adriatic Sea [31,36,37], the annual cycle shows great year-to-year variability, and the model well reproduces it. Overall, no systematic model bias to be considered appears evident in ordinary and stormy conditions. Given the relative proximity of the AA platform to the coast, this result gives confidence in the capability of CCLM-CMCC in reproducing the winds just offshore the coastal regions. Performance of the analogous CCLM-CMCC implementation over the Adriatic Sea for the present-day period was also assessed in [31], which verified significant achievements in terms of wind directionality, with unprecedented performance of this regional climate model in reproducing the bimodal dominance of Bora and Sirocco weathers.

The frequency distribution of wind speed (Figure 6A) features very good agreement between CCLM-CMCC and observations over the full range of speeds. On the other hand, as shown above, ERA5 substantially underestimates the occurrence probabilities above 6–7 m/s and it achieves lower speeds than CCLM-CMCC. From the histogram of wind direction in Figure 6B, one may notice that prevailing CCLM-CMCC winds are triggered by northeasterly Bora storms. The peak’s observed and model flow direction (around 230° N) are consistent, accounting for the mentioned angular difference. We note the good match between historical CCLM-CMCC and ERA5 directions, and between CCLM-CMCC historical and near-future (RCP4.5 and RCP8.5) directions. The directional agreement can also be found at the basin scale (Figure 7), where there is an evident control of northeasterly winds in the north and southeasterly winds in the central Adriatic, marking the regions where Bora and Sirocco dominate, respectively. The effect of the native resolution of ERA5 is visible (Figure 7A), whose 0.25° fields are less able to reproduce the small-scale features of surface winds at the land-sea boundary, near the coasts, though consistent with those of CCLM-CMCC over all periods (Figure 7B–D). It is worth noting that for future projections to the end of the 21st century, the steadiness of directions for the extreme Bora sea states does not seem to be guaranteed [35].

3.3. ERA5c Wind

Before going into details of the corrected ERA5 wind fields, it is convenient to sum up the attributes of the correction method in order to help interpret the model performance. First, the ERA5c dataset may answer the question as to how past storms would behave under climate change. The present-day climate fluctuations also represent the future year-to-year variability of the system, and only the external climate forcing effects are extracted from the model estimates of the scenario. In other words, the assessment of the regional climate responses is not affected by the changes in internal variability. This condition applies not only to long-term climate simulations but also to short-term event simulations. Moreover, by scaling wind speeds, information on future climate is available only for the intensity of storms, while their temporal structure and directional properties are not handled. This assumption is reasonable in a semi-enclosed basin, such as the Adriatic Sea, where the wind field structure is mainly controlled by the geographical land-sea distribution [27,74]. The steering is reproduced by numerical models with fine grid spacing, which allow an accurate description of the orography and a better representation of small-scale physics. This condition is fulfilled by the regional climate model CCLM-CMCC and the reanalysis ERA5, as shown in the previous sections. However, the invariance of the wind patterns may be no longer valid for projections to the distant future (toward the end of the 21st century, for instance), when, for a specific emission scenario, a marked projected change in the atmospheric properties at the synoptic scale can be found [33]. This issue is common to bias correction techniques applied to future climates very different from the present [15].

Furthermore, over the Adriatic Sea, dominant northeasterly Bora and southeasterly Sirocco winds were quantile-corrected with the same set of coefficients; that is, CF was not fashioned depending on the wind direction. This shape of CF might produce a misestimation in regions where the wind climate is not unimodal in direction. In the North Adriatic (AA platform), where the wind direction distribution is bimodal, an analysis we made has shown that directional differences are small. With respect to the corrected winds computed from all directions, we have found that northeasterly would experience a −0.2 m/s correction at 50th and 99th percentiles, while for southeasterly the difference would be −0.2 m/s and +1.5 m/s at 50th and 99th percentiles, respectively. The overall good performance of ERA5c confirms the validity of such an approach for all regimes. However, we expect an even better matching by further partitioning to account for the wind direction and allowing a smooth transition among regimes. Separation is also possible for the seasons using separated CFs, even though in the Adriatic Sea, the most energetic Bora and Sirocco storms typically occur in winter. All this said, we must keep in mind that by increasing the dimensionality of the statistical correction, the physical consistency of the resulting fields is degraded. As the goal is to use the wind fields to force a physically based spectral wave model, the inconsistencies might be alleviated. For the future period, the independence of CF from the wind direction is allowed if the different regimes undergo an even change. We verified this condition in the northern Adriatic, where previous studies showed, for the distant future (end of the 21st century), a possible shift in the interplay of Bora and Sirocco waves [32]. At AA, the CF values from the overall and northeasterly winds are consistent under RCP4.5 and RCP8.5 compared to CTR. This agreement is also found between overall and southeasterly winds under RCP4.5. On the other hand, under RCP8.5, the southerly storm conditions (speed greater than about 10 m/s) might be underestimated by the adopted correction strategy.

The effect of the statistical correction of the ERA5 wind speed fields over the Adriatic Sea for the present-day period is shown in Figure 8. The geographical patterns of ERA5 (Figure 8A) and climate-corrected ERA5 (ERA5c, Figure 8B) average winds are very similar, but the latter provides higher speeds over every sub-basin; indeed, as mentioned before, CF is always greater than unity and shows a substantial regional dependence (Figure 8C). The increase is greatest close to the land-sea boundary (up to +50%) near the Italian coasts (West Adriatic), where the high-resolution CCLM-CMCC data allows for improved reproduction of the offshore and onshore wind conditions. High values of CF (from 1.2 to 1.3) also characterize the regions with strong winds, for instance, along the northeasterly Bora wind jets in the north and central basins. Areas with dominant Sirocco (see Figure 2) and between jets are corrected less (CF~1.1 to 1.2). This correction is similar to the one (CF = 1.11, obtained with a least-square fitting) applied to ECMWF operational forecast winds (9-km resolution) for the simulation of the October 29, 2018, exceptional Sirocco event in the Adriatic Sea [75]. We finally note that despite the fact that the QQM statistical method is technically blind to the source of the difference between the datasets that are quantile mapped, we recognize that this is the spatial resolution that penalizes ERA5 in the semi-enclosed Adriatic Sea.

Using the QQM method may raise the question as to whether the climate signal from CCLM-CMCC is preserved in the adjustment of ERA5. Figure 9 shows the comparison between the annual mean 50th and 99th percentiles of U₁₀ from CCLM-CMCC (Figure 9A,D), daily average ERA5c (Figure 9B,E), and hourly ERA5c (Figure 9C,F). For CCLM-CMCC, the annual mean 50th percentile shows a clear northwest-to-southeast negative gradient, with Bora-like conditions in the northern and central Adriatic never exceeding 4 m/s. In contrast, maxima are dominated by long-fetched, persistent wind sea conditions in the south, with values above 5 m/s. On the other hand, the annual mean 99th percentile is dominated in the north by short fetched, orography-dependent winds funneled through gaps in the Dinaric Alps of the eastern coast, where the strongest Bora jets reach values around 17 m/s in the northern Adriatic and 16 m/s in the central Adriatic, while Sirocco wind speeds in the south never exceed speeds of 16 m/s. There is a substantial similarity between CCLM-CMCC and ERA5c. The latter preserves the climate signal of U₁₀, even though there are local ERA5c minus CCLM-CMCC differences, which are about +0.5 m/s at most for the 50th percentile. The absolute difference at the 99th percentile is smaller and below 0.3 m/s. The likely reason for this mismatch at both percentiles is that the correction defined on the averages and then applied to the hourly data does not preserve the mean values since the CF is not constant over speeds. Therefore, the daily average of corrected winds is slightly different from the correction of daily average winds. This difference may be used to quantify the effect of the assumption we made on the similarity between daily average and hourly matching functions. Finally, Figure 9C,F show the same statistics from ERA5c hourly data; they indeed give higher values than their daily counterpart and provide peak values (99th percentile) of 18 m/s and 17 m/s in the northern and central Adriatic and around 17 m/s in the southern Adriatic.

To further evaluate the effect of QQM, we used the set of observed data to assess the winds after correction and the corresponding wind-generated waves. Figure 10B shows that in the northern Adriatic Sea (AA platform), the time series comparison between observed, ERA5, and ERA5c sequence of storms reveals an improved model performance at all speeds and wind directions (both Bora and Sirocco weathers) for the climate-corrected ERA5c dataset, which largely compensates for the systematic underestimation of ERA5. Figure 10A depicts the scatter diagram and the quantile plot between observations and ERA5c, which mark the very good performance of the present-day ERA5c at all speeds (ME = 0.15 m/s and Slope = 0.96). The scatter spread does not seem to be reduced by the correction procedure, probably because it is associated with the timing of events and with the small-scale atmospheric patterns not well reproduced by the native ERA5.

3.4. Waves

Wave model results display similar skills to those seen for ERA5c wind. Figure 11 shows the comparison (scatter and quantile plots) of measured and modeled values H_s at AA and the two wave buoys of Ortona and Monopoli. In all three locations, the wave model provides a slight overestimation of the average wave energy, being the ME for H_s < +0.04 m, and the MAE about 0.2 m on average; the Slope is at least 0.98 at AA and 1.08 at most in the southern part of the Adriatic Sea (Monopoli buoy). The performance of this wave model is similar to that obtained in the Adriatic Sea by a dedicated high-resolution hindcast purposely developed to reproduce the sea states over the entire Mediterranean Sea [76].

4. Application to Wave Climate Scenarios over the Adriatic Sea

4.1. Basin Scale Analysis

Projected future evolutions of the wind and wave climate over the Adriatic Sea by the middle of the twenty-first century (2021–2050) are assessed from a single GCM-RCM and for two representative concentration pathways, namely a medium (RCP4.5) and a high-emission scenario (RCP8.5). As before, we focus on the 50th (median value) and 99th percentiles and consider them to be representative of the typical and the extreme climates, respectively. For the selected quantities, statistics was computed for each year and averaged over the entire period to provide an annual average value (30 samples for each simulation). The impact of the projected climate is expressed by considering the relative change for the period 2021–2050 compared to the period 1981–2010. We provide an assessment of the U₁₀ change first; then, we focus on the H_s fields obtained by forcing the high-resolution WW3 model with the climate-corrected ERA5c winds. The stochastic significance of the projected change is assessed using the Wilcoxon–Mann–Whitney U-test with sample size N = 30 and significance level α = 0.1 [77].

Figure 12 shows the results for the sea surface wind speed U₁₀ from CCLM-CMCC. The present-day statistics are shown in Figure 12A (annual average of the 50th percentile) and Figure 12D (annual average of the 99th percentile). Relative changes under RCP4.5 (Figure 12B,E) and RCP8.5 (Figure 12C,F) exhibit characteristic spatial patterns, but the latter shows generally larger maximum changes. A weak signal of projected changes (~1%) of the annual average 50th percentile of U₁₀ is depicted for both scenarios, resulting in decreased speeds under RCP4.5 and increased speeds under RCP8.5. Except for a few sporadic places, these changes are statistically not significant. A more robust and spread decrease at the basin scale is seen for the stormy events (99th percentile). The change is biggest in the north (up to −8% and −5% for RCP8.5 and RCP4.5, respectively), while it is smaller (~−3%) in the mid-Adriatic region, where along-basin winds eventually dominate the weather.

These wind variations are similar to those found by Bonaldo et al. [31] under RCP8.5 for the end of the 21st century (2071–2100). The authors found a −10% decrease across the Adriatic Sea, although mainly with a stronger reduction along the Bora jets patterns. For the same scenario (RCP8.5) but a different future period (2041–2070), Belušić Vozila et al. [33] used an ensemble with nineteen members to find that the number of Bora events in the northern Adriatic may be reduced, though intensities may be higher. However, a spread between the members was observed, suggesting the need for multi-model simulations to assess the significance of the projected wind/wave simulations, particularly for the distant future [17]. For completeness, we report that the multi-decadal historical observations from AA (1979–2015; [78]) showed, on the one hand, a tendency toward a reduction of the wave activity in the northern Adriatic Sea and, on the other that the decrease of intensity is mainly due to Bora regimes (particularly in the upper percentiles). Later in this section, we will focus on separating the north-east (Bora) and south-east (Sirocco) storm events by characterizing the long-term wave statistics in front of the Venice littoral.

Since the local winds determine surface waves in the Adriatic Sea, a similar geographical pattern for the variations has been obtained from the wave model simulations (Figure 13B–F). This similarity is not necessarily valid for the annual mean values of 50th and 99th percentiles for the present-day period (CTR run, Figure 13A,D), which show a statistical dominance of along-basin sea states (SE-NW directions) compared to northeasterly Bora ones due to the effect of the longer fetch available for waves to grow [79]. On an annual average base, the roughest stormy conditions occur in the southern Adriatic (H_s up to about 4 m for the 99th percentile), while in the central and northern parts of the basin, conditions are milder (H_s < 3.5 m and H_s < 3 m, respectively). Similar patterns, but with less intense waves, from ERA5-based hindcast simulations (1981–2019) were found in the Adriatic Sea by Barbariol et al. [12] who extended the analysis to the whole Mediterranean Sea. The geographical distribution of stormy H_s does not necessarily match the one where the maximum H_s may occur, e.g., offshore the northeastern Italian coast during Bora storms or offshore Croatia during Sirocco storms. For this specific topic in a projected climate, the reader is directed to Denamiel et al. [35].

Projected to the near future (2021–2050), typical significant wave heights (50th percentile) show minimal changes (<3%) for both scenarios (Figure 13B,C), although RCP4.5 indicates slightly negative variations in the central/north part of the basin. A decrease of extreme significant wave height (99th percentile; Figure 13E,F) up to −10% is visible in the northern part of the Adriatic. Weaker changes (around −5%) are reproduced for the central and southern parts of the basin. RCP8.5 produces generally greater changes than RCP4.5, albeit locally, the opposite also holds (e.g., along the Croatian coast). This widespread decrease in wave energy confirms previous scenario analyses conducted in the same area with different models and projection strategies [32,35,37,38,80].

These results can be placed in a broader perspective. The connection between the Adriatic Sea wind and wave regimes and the Mediterranean area behavior is straightforward [81], but less simple is the relationship on the larger scale, for instance, with the North Atlantic Ocean [78], and this is beyond the scope of the present study. However, we point out that projected changes in Sirocco weather may indicate a variation of the cyclogenesis in the central-western Mediterranean Sea and the movements of cyclones southeastward. These low-pressure systems also generate Bora winds, driven by high pressure over central Europe [82]. The frame is further complicated by pressure lows and resulting local winds passing over the Adriatic Sea [44,83]. This complexity and its relative changes in the mid-to-long-term future have practical applications such as for instance, the impact on the coasts [84]. Indeed, wind-wave climate changes, together with sea level rise, are the main hazard to coastal stability, especially in the case of sandy coasts [85]. This is the case for the central and northern Italian coast facing the Adriatic Sea. Particular relevance can be given to the sea area surrounding the Venice lagoon (Figure 1), where the Mo.S.E. multiple barriers system at the lagoon’s three inlets ([86]) has been operational since 2020 to protect the historical city of Venice and its lagoon from the extreme surges associated with Sirocco events [87].

4.2. Extreme Value Analysis Offshore Venice (Italy)

We may want to use the projected wave storms to assess what might be the sea state conditions just offshore the Mo.S.E. barriers during their first 30 years of operation. Propagated shoreward and eventually inside the inlets, the wave heights would help define the loads on the structures. For this analysis, we have selected a location offshore the middle inlet (called “Malamocco” by the locals), where the depth is about 14 m (Figure 14). For the present-day simulation covering 1981–2010, the 30-year long time series of hourly values of H_s has been analyzed to infer the local wave climate. The statistics of extreme values of the significant wave height has been estimated following the peak-over-threshold approach (POT; [88]), where the maximum value of H_s in each of a large number of individual storms is considered. A storm is defined as an uninterrupted sequence of H_s values exceeding a specific threshold. Following [89], the threshold for storm identification was set equal to 1.5 times the 30-year average H_s. Storms with successive peaks that occur within less than 10 h were aggregated together to ensure stochastic independence. The peak of H_s for each storm was defined as an extreme value for the characterization of the wave climate. Finally, to isolate the more energetic storms, only H_s values above 2 m were considered for the determination of the empirical probability curves.

The extreme value analysis was made of both the entire, omnidirectional time series and the storms coming from the two angular sectors (0° N, 110° N) and (111° N, 180° N) for the peak direction of wave propagation (Figure 14A). These are preferentially driven by Bora and Sirocco winds, respectively. We note that the incoming direction of northeasterly waves is rotated clockwise relative to the forcing wind due to the wave refraction in the relatively shallow environment close to the inlets. Sirocco wave directions, on the contrary, are less affected by the interaction with the bottom since their incoming direction is almost orthogonal to the coastline. The POT analysis for the period 1981–2010 led to the identification of 246 storms (i.e., 8.2 per year on average) in the entire series, and 154 (5.1 per year) and 92 (3.1 per year) storms for Bora and Sirocco, respectively (red markers in Figure 14B–D). Maximum values of H_s were found at 3.47 m for Bora seas and 4.30 m for Sirocco seas. The return period T_R of the largest value is about 30 years, according to Weibull [90]. Furthermore, empirical POT data were fitted to a generalized Pareto distribution (GPD; see, e.g., [91]), whose parameters were determined with a maximum likelihood estimator. Results from the GPD fitting are shown in Figure 14B–D for the omnidirectional, Bora, and Sirocco storms, respectively. From the GPD fitting, the H_s values at T_R = 30 years are estimated at 3.30 m and 4.14 m for Bora and Sirocco, respectively. As expected, the longest fetch available for Sirocco waves to develop impacts the more energetic storms that may reach the area. The good agreement we found between model and measured H_s at AA (a few kilometers offshore of the Malamocco inlet) gives confidence in the reliability of these estimates.

To assess what the wave conditions just offshore the Mo.S.E. barriers might be for the projected future, the same extreme value analysis has been applied to the H_s results of the RCP4.5 and RCP8.5 simulations. Figure 15 shows the comparison of the peak H_s at different T_R (2, 5, 10, and 30 years) between the present-day period (CTR, 1981–2010) and the two scenarios RCP4.5 and RCP8.5 for the years 2021–2050. It appears that under both RCP4.5 and RCP8.5, sea states at all return periods will undergo a decrease of the wave energy, which will be larger under RCP8.5 (about −9%) than under RCP4.5 (about −4%), and with minor differences (within 1.5%, on average) between Bora and Sirocco. Changes seem to be stochastically significant, with stronger evidence under RCP8.5. The similarity between the two regimes is the integral result of the two factors that most influence the wave growth: the wind speed (Bora speeds are generally stronger) and the fetch (for Sirocco, almost the whole Adriatic Sea is available). Especially for the southeasterly Sirocco sea states, the cumulative effect along the Adriatic Sea SE-NW axis is at a maximum near the north coast. On a slightly different subject, the same trend can be found in the analysis of Sirocco-driven storm surge projections in the same area [29,35]. However, since additional factors such as subsidence and sea level rise will impact the effectiveness of the Mo.S.E. barriers, a future reduction of the wind-supported sea energy does not necessarily imply increased safety for the coastal structures and the lagoon of Venice.

5. Concluding Remarks

The present study investigated the use of a statistical method to climate-correct the sea surface wind speed from the ECMWF reanalysis ERA5. We relied on a quantile-quantile matching method to correct the reanalysis winds in order to match the present-day (1981–2010) and near-future (2021–2050) states of the climate system. The procedure is applied over the Adriatic Sea and makes use of the fields provided by the high-resolution (0.0715° in lon-lat) regional climate model COSMO-CLM with the configuration over Italy developed by CMCC under two representative concentration pathways: a medium (RCP4.5) and a high-emission scenario (RCP8.5). Therefore, this study used matching to achieve a statistical correspondence between the reanalysis and the climate model data. To assess the climatological fields of the significant wave height (a proxy of the total wave energy), the climate-corrected ERA5c winds were used to drive spectral wave model simulations at 0.025° of resolution over the present and future periods. No observation of wind and waves was incorporated in the correction procedure, but the present-day simulations of wind and waves were evaluated against measured data at three independent locations along the basin to assess uncertainties. The main results of the study can be summarized as follows.

The model-to-model quantile correction method proved to effectively adjust the ERA5 bias over the Adriatic Sea, still maintaining a high degree of physical consistency. The correction was considered wind speed-dependent and geographically variable. The reanalysis ERA5 has the advantage of reproducing at high-temporal resolution (1 h) the realistic sequence of the historical events. Both features are required to characterize the wind-wave storm and climate properly. However, over the Adriatic Sea, the effect of the relatively coarse, native spatial resolution of ERA5 fields was to smooth the sharp gradients of the evolving wind field. The resulting bias could be then transmitted to the wave field. On the other hand, the high-resolution CCLM-CMCC was particularly skilled in the reproduction of wind magnitudes and directions of the dominant regimes. After matching, corrected ERA5 data have proven helpful for the wind and wave models evaluation since, for the historical period, a storm-based analysis against observations can be performed. Moreover, constraining the sequence of storms in the present-day and future climate experiments has advantages because the assessment of the regional climate responses is not affected by the changes in internal variability. It is expected that the technique may be applied to other reanalyses and basins with similar characteristics. The analysis may be further extended by adopting, for the correction, multi-model ensembles of GCM-RCM to account for the climate model uncertainties.

It is worth noting that the statistical correction of the reanalysis wind speed might introduce undesired secondary effects on the resulting fields [92]. Known deficiencies of bias correction methods, such as the inability to correct the storm track and timing, are largely alleviated using reanalysis variables as input. For the historical period, additional bias in the wind pattern can occur in those regions where the climate is at least bimodal in direction, and the correction of ERA5 differs among regimes. The frame is further complicated for regimes expected to undergo a different change in intensity and trend in the future. For all those situations where they do not strictly apply (e.g., for simulations covering the distant future or running over a large sea basin with weak or null orographic constraint), an improved strategy should be adopted, for instance, using a direction-dependent correction factor for ERA5. A possible method to correct the wind vector bias (i.e., both directional components) can use a bivariate quantile adjustment as in [72]. We note that this additional complexity might reduce the consistency of the corrected variables.

The proposed approach was implemented to simulate the wind-wave climate of the Adriatic Sea. Although in the absence of an assessment of the uncertainty associated with the choice of an individual climate model, our results show that for the near future (2021–2050), there is evidence of a statistically robust and consistent reduction of storm intensity and significant wave height (99th percentile) in the northern part of the Adriatic Sea, where Bora dominates the climate, with a larger change (up to −10% of H_s) under the high-emission scenario RCP8.5 than under RCP4.5. This result is consistent with previous climate studies and is compatible with an expected future shift of the winter cyclones over the Mediterranean Sea, consolidating the established relationship between large-scale changes induced by carbon dioxide increase and local effects. For the central and southern regions of the Adriatic Sea, a pattern of uniform change is hardly visible (not significant), albeit a spread decrease in the sea state severity has been found. For the milder wind regimes, the change we have found is small and not statistically significant. When focusing on the storms just in front of the Venice lagoon (where mobile barriers are protecting the city from floods), we demonstrated that H_s values at 2-, 5-, 10-, and 30-year return periods will undergo a decrease (−4% to −9%) in the future (2021–2050) under both RCP4.5 and RCP8.5 scenarios. This result does not necessarily guarantee increased safety for the coastal structure, as wave impacts must be combined with the local subsidence, sea level rise, and storm surge. All these issues deserve further in-depth study.