Non-stationary GEV models for estimating design sea-states in a changing climate – applications to offshore wind farms along the French coasts

The rapid expansion of the French offshore wind sector requires a critical reassessment of structural durability in the face of evolving marine conditions driven by climate change. Traditional design methodologies assuming stationary conditions are inadequate for climate-resilient infrastructure, as they fail to account for evolving wave patterns, storm tracks, and extreme event frequencies (Amlashi, 2024; Barkanov et al., 2024).

This study introduces a non-stationary GEV framework to quantify future changes in significant wave height using monthly maxima, reducing uncertainty in design tools for offshore wind farms under climate change. While non-stationary GEV models are well established, their application to offshore wind design using climate projections remains limited. The main objective is to derive an equivalent stationary design level that accounts for time-varying extremes over the lifetime of structures. Based on CMIP6 climate models and reanalysis data, results reveal a projected trend towards a more pronounced seasonal contrast along the French Atlantic and English Channel coasts under future scenarios (SSP1-2.6 and SSP5-8.5), whereas the French Mediterranean Sea exhibits weaker increases in extremes and larger uncertainties (inter-model spread). Projections indicate more intense winters and calmer summers, along with a shift in the seasonal cycle. Overall, the multi-model ensemble suggests an increase in the design levels for extreme sea states.

The research concludes by defining a new methodology for calculating an equivalent design level over the structure's operational lifespan. This tool is deemed essential for ensuring the resilience and economic viability of future offshore wind farms in a changing climate.

1 Introduction

The French offshore wind sector is undergoing rapid expansion, with a projected cumulative capacity of 3.6 GW expected by the end of 2027, driven by the commissioning of seven major offshore projects. This growth occurs in a context that is both innovative and promising, yet increasingly competitive.

Long-term durability and reliability of offshore wind structures critically depend on their ability to withstand evolving marine conditions over operational lifespans that often exceed two decades. Climate change is modifying key metocean variables, including significant wave height, total sea level, and wind speed, as well as the frequency and intensity of extreme events (Amlashi, 2024; Barkanov et al., 2024). These evolving conditions introduce substantial uncertainty in design parameters, particularly in regions exposed to high wave energy and strong winds, which remain the most favorable for energy production (Susini et al., 2022). Neglecting these changes may result in under-designed structures, reduced operational lifespans, increased maintenance costs, or even catastrophic failure (Chellaa et al., 2012; Schloer, 2013).

Metocean studies are essential for understanding and forecasting the environmental conditions that offshore structures will face throughout their service life. In particular, extremes in significant wave height, total sea level, and wind speed are critical design parameters that directly influence structural performance, safety, and operational efficiency (Slater et al., 2021). Significant wave height, representing the average height of the highest third of waves, is a key factor in determining wave loading on structures and plays a central role in fatigue analysis, foundation stability, and the overall design of wind turbine components (Zhang et al., 2019).

Traditionally, offshore wind design assumes stationarity, meaning that the statistical properties of environmental conditions (mean, variance, extreme values…) are constant over time. However, the non-stationary nature of climate change challenges this approach. Sea level rise, shifting storm tracks, and intensifying wave conditions can significantly affect wave loading, fatigue life, and structural integrity (Chellaa et al., 2012). Capturing non-stationary behavior is essential for robust design and accurate estimation of future extreme conditions (Chavez-Demoulin and Davison, 2004). Seasonality is also crucial, as more intense winter extremes and calmer summer conditions may affect additional fatigue (Susini et al., 2022). Correctly characterizing these contrasts is critical for fatigue analysis, foundation stability, and wave loading calculations (Slater et al., 2021; Zhang et al., 2019). Although non-stationary GEV models are widely studied in the statistical literature (Coles, 2001; Chavez-Demoulin and Davison, 2004), their application to CMIP6 climate projections for French offshore sites is limited. And projecting return levels up to 2100 raises additional challenges of data extrapolation and model uncertainty, as estimating long period return values in a changing climate is known to be associated with large uncertainties related to model variability and methodological choices (Ewans and Jonathan, 2023). So, existing methods relying solely on historical stationary data may underestimate future extremes, leading to inadequate design loads.

Neglecting these evolving trends may lead to under-designed structures, reduced operational lifespans, and increased maintenance costs, ultimately compromising the economic viability and safety of offshore wind projects. It is therefore imperative to assess the impact of climate change on sea state conditions to inform the design criteria of offshore wind farms and other coastal infrastructure. This research aims to address this challenge by suggesting a new methodology to better investigate projected changes in significant wave height. The final objective is to provide engineers, designers, and stakeholders with updated tools that account for non-stationary environmental conditions, ensuring that offshore wind farms remain resilient and cost-effective in a changing climate.

This study develops a statistical framework using non-stationary GEV models applied to monthly maxima to better quantify projected changes in significant wave height. Its main contribution lies in applying these models to CMIP6 projections along the French coasts to assess future extreme sea states, while also deriving an equivalent design level that accounts for non-stationarity over the operational lifespan of offshore structures, following the approach recently formalized by Barbaux et al. (2025). This framework provides engineers, designers, and stakeholders with updated tools to ensure resilience, safety, and economic viability of offshore wind farms under evolving climate conditions. In this paper, the statistical methodology is thoroughly presented in Sect. 3, focusing on the application of non-stationary GEV models to monthly maxima. Subsequently, Sect. 4 highlights the results, emphasizing projected changes in return levels and extreme sea states extending up to 2100. Finally, Sect. 5 discusses the implications of these findings for offshore wind design and adaptation strategies, while critically examining the uncertainties associated with long-term extrapolation and model assumptions.

2 Sea-state data

2.1 Sea-state data

2.1.1 CMIP6 Climate models

Climate models are an essential tool for assessing future environmental responses, with the Coupled Model Intercomparison Project Phase 6 (CMIP6) representing the latest generation available from the IPCC (Intergovernmental Panel on Climate Change – Sixth Assessment Report, AR6, IPCC, 2022). The wave data used in this study are specifically derived from a numerical wave model (WaveWatch III) that was forced by eight General Circulation Models (GCMs) from the CMIP6 ensemble (listed in Table 1). More details about this data can be found in the article of Meucci et al. (2024). The dataset covers a historical period from 1985 to 2014 and a future projection period from 2071 to 2100, with 30-year duration for each period. This choice ensures a consistent basis for comparison between past and future conditions, avoiding bias that could arise from periods of unequal length. The selected periods also align with AR6 recommendations and standard climatological practice for assessing baseline and future climate scenarios, specifically SSP1-2.6 and SSP5-8.5, which correspond to low and high greenhouse gas emission pathways.

Bi et al. (2020)Semmler et al. (2020)Lovato and Peano (2020)EC-Earth Consortium (2019)Boucher et al. (2018)Pak et al. (2021)Wieners et al. (2019)Yukimoto et al. (2019)

Table 1 List of the wind/waves climate models used in this study.

Data prior to 1985 were not included due to limited availability and reduced reliability of reanalysis-forced wave data. Similarly, the historical and future periods are “discontinuous” because the projections are intended to represent long-term climate scenarios rather than a continuous time series. The wave data have a 3 h temporal resolution and a 0.5° spatial resolution, providing sufficient detail for extreme wave analysis.

2.1.2 Reanalyses

Reanalyses are necessary to evaluate the quality of climate models over the past period, called “historical period”, and to make bias-adjustment for future projections when necessary. In the English Channel and the Atlantic Ocean, the HYWAT (https://doi.org/10.5150/jngcgc.2024.094, Michaud et al., 2024) reanalysis was used to make this bias correction. The HYWAT reanalysis extends from latitudes between 43 and 52° N and longitudes between 7° W and 5° E in a non-structured grid. Data is hourly with a spatial resolution from approximately 500 m on the coast to a few kilometers offshore. The reanalysis is based on the configurations of HYCOM and Wavewatch III with a parametrization corresponding to TEST 471 (Michaud et al., 2024). HYWAT reanalysis is forced by the ERA5 atmospheric reanalysis for the wind and atmospheric pressure fields (Hersbach et al., 2020). HYCOM currents, water levels and surges (Jourdan et al., 2020) are also provided to the Wavewatch III model every 10 min. In the Mediterranean Sea, the MED-WAV reanalysis (https://doi.org/10.48670/MDS-00376, European Union-Copernicus Marine Service, 2025) was used to make the bias correction. The MED-WAV reanalysis has two grids. The coarse grid covers the North Atlantic Ocean from 75° W to 10° E ($\mathrm{1}/\mathrm{6}$° resolution) and from 70° N to 10° S while the fine grid covers the Mediterranean Sea from 18° W to 36° E and from 30 to 46° N ($\mathrm{1}/\mathrm{24}$° resolution). Data is hourly with $\mathrm{1}/\mathrm{24}$° resolution. The reanalysis is based on the WAM 4.6.2 model. MED-WAV reanalysis is forced with daily averaged currents from Med-Physics and with 1 h, 0.25° horizontal-resolution ERA5 reanalysis 10 m-above-sea-surface winds from ECMWF.

2.2 Bias correction

We applied the Cumulative Distribution Function transform CDF-t method (Michelangeli et al., 2009) to both correct the bias and to downscale large-scale GCMs data to the specific locations, using local-scale reanalysis as the reference. This correction was applied at six offshore locations in this study (situation in Fig. 1 and coordinates in Table 2), selected as representative sites, representative of different climates, for future offshore wind farms across the French seafronts. They are selected to be sufficiently spaced (see Fig. 1) to fall in different GCM grid cells, and in grids over the open waters, as explained in Kervella et al. (2025).

Figure 1 Locations of the Selected Representative Point (SRP) along the French coasts.

Table 2 Seafront representative points coordinates.

The CDF-t method (Michelangeli et al., 2009) assumes that a mapping exists between the Cumulative Distribution Functions at the global scale (corresponding to the GCM data) and the local one (corresponding to the reanalysis data) for the current period, which is defined by:

$$\begin{array}{}\text{(1)}& T\left(F_{m,h}\right(x\left)\right)=F_{o,h}\left(x\right)\iff T\left(u\right)=F_{o,h}\left(F_{m,h}^{-\mathrm{1}}\right(u\left)\right),\end{array}$$

with F_m,h the CDF of a GCM and F_o,h the CDF of the reanalysis over the historical period, h. Moreover the authors assume that the bias correction between the past and the future remains the same, meaning that the difference between the GCM and the reanalysis is supposed to be stationary. Thus, it is possible to estimate the CDF of the future local data through the transformation of the GCM projections. The final equation, to correct the cumulative density of the projection period of the GCMs is given by:

$$\begin{array}{}\text{(2)}& F_{o,f}\left(x\right)=F_{o,h}\left(F_{m,h}^{-\mathrm{1}}\right(F_{m,f}\left(x\right)\left)\right),\end{array}$$

with F_m,f the CDF of GCM on the future period, f, and F_o,f(x) its corrected cumulative density function. The CDF-t has been applied independently for each GCM, using the reanalysis data (HYWAT or MED-WAV depending on the sea-front), and for each month. For data that are beyond the rage of the observations, a constant correction is applied as in Déqué (2007), Michelangeli et al. (2009), and Kervella et al. (2025).

Figure 2 illustrates the impact of this correction through the annual average, the cumulative distribution function (CDF), and the monthly average. The correction significantly reduces the dispersion of the models' annual averages relative to the reanalyses, and aligns GCM distributions with reanalysis data. Although the CDF-t preserves the trend across all quantiles, adjusting only the rank and distribution properties (not the mean), the dispersion in the seasonal cycle is also notably reduced. This latter improvement is achieved because the bias correction was applied separately to each month. More details about the results of the CDF-t method for this application can be found in Kervella et al. (2025).

Figure 2 Validation of the CDF-t bias correction method at the Eastern English Channel representative point. (Upper row) Comparison of annual average significant wave height (H_s) before (raw GCM outputs, left) and after (right) correction, with the HYWAT reanalysis (black solid line) as reference; x-axis: time (years); y-axis: H_s (m). (Middle row) Cumulative Distribution Functions (CDFs) of H_s for raw GCMs (left), corrected GCMs (right), and HYWAT reanalysis (black); (Bottom row) Monthly average H_s before (left) and after (right) CDF-t, with HYWAT reanalysis (black solid line); x-axis: months (January–December); y-axis: H_s (m). Colors represent individual CMIP6 models (see Table 1).

3 Modeling of extreme values

3.1 Models for annual maxima

Following Coles (2001), we model the yearly maxima using a GEV distribution:

$$\begin{array}{}\text{(3)}& M_n\sim \text{GEV}(\mathit{\mu },\mathit{\sigma },\mathit{\xi }),\end{array}$$

where μ is the location parameter, σ a scale parameter and ξ a shape parameter.

From a practical point of view, the data is split by year, model and scenario, and the GEV distribution is fitted on the yearly maxima of each group, using a maximum-likelihood approach.

Once the model is fitted, we are able to go beyond the range of observed values, and estimate any quantile of the distribution of the yearly maxima, what is referred to as return level in the literature. In Fig. 8, we show the estimations of such quantities from the fitted model, for each location and under each scenario. We have also checked individually that the fitted model fit well the data by looking at Q–Q-plots (not shown here).

3.2 Models for monthly maxima

3.2.1 Methodology

As described before, using monthly maxima allows to discard less data that using the yearly maxima, at the added complexity to account for seasonality. Generalized Additive Models are very convenient to account for such effect, as already introduced in Chavez-Demoulin and Davison (2004), and further applied to sea-state data in Reinert et al. (2021), Roustan et al. (2022), and Cheynel et al. (2025). This models are briefly described hereafter. When taking the maxima on yearly blocs, a lot of data is not used, which results in large uncertainties in the fitted parameters due to the paucity of data. However, when using smaller blocs size, e.g. monthly blocs, we have to take into account the seasonality of the data, because the monthly maxima cannot be assumed to be identically distributed (see e.g. Trasch et al., 2023). This study adopts as its starting point the following equation, which assumes that the distribution of the monthly maxima adheres to a GEV distribution:

$$\begin{array}{}\text{(4)}& M_{m,n}\sim \text{GEV}(\mathit{\mu }=f_{\mathrm{1}}(m),\mathit{\sigma }=f_{\mathrm{2}}(m),\mathit{\xi }),\end{array}$$

where m is the month of the year. Because estimating ξ is notably difficult (see Coles, 2001), it is assumed to be constant for all the months.

3.2.2 Computation of return levels

The N-year return levels are defined (see Coles, 2001), as the value which is expected to be exceeded on average once every N years. From Eq. (4), it can bee seen that this quantity is not accessible directly because we are dealing with monthly maxima. However, the distribution of yearly maxima can be computed back, by taking the product of the monthly CDF, using:

$$\begin{array}{}\text{(5)}& \mathbb{P}(M_n\le x)=\prod _{m=\mathrm{1}}^{\mathrm{12}}\mathbb{P}(M_{n,m}\le x),\end{array}$$

which is justified if the monthly maxima are independent. We tested this hypothesis by analyzing the auto-correlation and partial auto-correlation functions. Our analysis revealed no remaining time dependence in the residuals of the models, which strongly supports the initial hypothesis. The corresponding figures for the representative point in the Eastern English Channel are shown in Fig. 3. This finding is consistent with results observed in other studies (e.g. Cheynel et al., 2025; Youngman, 2022).

Figure 3 Comparison of auto-correlation function (ACF), left columns of each figure, and partial auto-correlation functions (PACF), right columns of each figure; The Location used here is the East English Channel, before modeling the non-stationarity of monthly maxima (left plot) and after (right plot). It can be observed from this plot that the seasonality is removed from the Monthly Maxima, and that, conditionally on the month, they can be considered independent.

3.3 Long-term non-stationary models

3.3.1 Methodology

The 100-year return value, which is the value having an annual exceedance rate of $\frac{\mathrm{1}}{\mathrm{100}}$, is in widespread use in design studies (Leach et al., 2025). Such a low probability of occurrence requires to base the study on sufficient amount of data, in order to reduce the statistical uncertainty associated with the extrapolation, at least 20 to 30 years of data. However, as the climate is evolving at a quick pace (Calvin et al., 2023), the long-term stationary assumption may no longer be valid, and extensions of such models are needed. In this study, we incorporated a non-linear, non-parametric trend in addition to the non-parametric cycles. More specifically, starting from Eq. (4), we propose a non-stationary GEV model where seasonal cycles evolve over time via spline functions on the tensor space of month and year (Wood, 2006):

$$\begin{array}{}\text{(6)}& M_{m,y}\sim \text{GEV}(\mathit{\mu }=f_{\mathrm{1}}(m,y),\mathit{\sigma }=f_{\mathrm{2}}(m,y),\mathit{\xi }),\end{array}$$

where m is the month and y is the year. We will assume here that for $i=\mathrm{1},\mathrm{2},f_i$ is a spline function on the tensor space of month and year (Wood, 2006), allowing to take into account both the seasonality and the change thereof as the time varies, due to the influence of a non-stationary climate. More precisely, the seasonal effect is assumed to be a cyclic P-spline function, with cycle equals to one year. For the time effect, we used thin plate regression splines. Each spline has its own maximal degree of freedom, which is estimated using the method described in Wood et al. (2016), and allows for the spline to shrink to a null value, leading to a constant term. This model has been fitted thanks to the {evgam} package (Youngman, 2022) of the R programming language (R Core Team, 2024). This model has been fitted independently on each GCM output and for each scenario described in Sect. 2 after the downscaling and bias correction step. We also compared models with only the location parameter as stationary, others with non-stationary shape and scale parameters, and several combinations of tensor products of month and year. We ultimately selected the model with the lowest AIC, as in e.g. Rohmer et al. (2020).

3.3.2 Computation of design sea-state in a non-stationary climate

We cannot assume, using model described in Eq. (6), that the annual maxima are identically distributed once the month is accounted for, contrary to the derivations made in Sect. 3.2.2. Indeed, the distribution of the annual maxima depends now on the year considered. The common practice in the design of structure, is to assume that the failure probability is directly linked to a quantile of the maximal value observed over the whole lifetime of the structure, e.g. 50 years. This value can be estimated from the annual maxima, assuming that these maximal are independent. To define an equivalent return level in a changing climate, we decided to take into account the whole lifetime of the structure and to compute a quantile of distribution of the maximal value observed during the whole installation time of the structure at its location. In the following paragraph, we define which quantile of this distribution we should consider, and how we can compute this distribution.

Say, without loss of generality, that the offshore wind farm is intended to stay for D_L=30 years, and is design to have an annual failure rate of $P_{\text{annual}}=\frac{\mathrm{1}}{\mathrm{100}}$. Then, the overall failure probability over the duration D_L is, assuming that the failure rate is independent of years:

$$P_{\text{survive over}{D}_L\text{years}}=(\mathrm{1}-P_{\text{annual}}{)}^{D_L}\approx \mathrm{0.74}.$$

In a non-stationary climate, the annual probability of default can no longer be considered constant, as it is influenced either by the internal variability of the Earth system or by the long-term climate response to anthropogenic forcing. In that case, we can write the distribution of the maxima over the whole lifetime of the structure, D_L years:

$$\begin{array}{}\text{(7)}& P(M_{D_L}\le x)=\prod _{y=\mathrm{1}}^{D_L}\prod _{m=\mathrm{1}}^{\mathrm{12}}\mathbb{P}(M_{m,y}\le x),\end{array}$$

assuming only that the maximal values across years are independent. Building on Eq. (7), we propose to define an equivalent design condition, obtained by taking the appropriate quantile of the distribution of the maximal value over the lifetime of the structure to obtain the same default probability over the lifetime, in a stationary climate, e.g. ≈ 0.74 in the example above.

3.4 Uncertainty quantification

The methods described in Sect. 3.2.2 and 3.3.2 allow for obtaining point estimates of the return levels derived from the different models presented in this paper. However, from a statistical perspective, it is also critical to quantify the uncertainties associated with these estimates and to assess whether the proposed methodology can outperform a stationary approach, in which climate change is not accounted for.

To quantify the uncertainties, we adopted a Monte Carlo approach, in which the parameters of the models are simulated from their posterior distributions B times using the fitted model. For each simulated set of parameters $(\mathit{\mu },\mathit{\sigma },\mathit{\xi })$ (which may be functional, as in Sect. 3.3.2 – though the functional notation is omitted here for simplicity), we compute the associated return level or the equivalent design condition, as applicable. This yields to B independent estimates of the quantity of interest, from which we derive uncertainty intervals by computing the empirical quantiles of the sample. In the Results section, we set B=1000.

4 Results

We now present the results obtained from fitting the non-stationary models to the wave data. For clarity, we report detailed findings only for a representative location in the Eastern English Channel; results for all other sites are provided in the Supplement.

Figure 4 displays the 99th quantile, calculated using Eq. (6), and estimated for each GCM at the East English Channel location. As stated in Sect. 3.3, this quantity is dependent on the month and on the year, as a function of both the position and scale parameters. The years are represented on the x-axis, while the y-axis contains the months, from January (bottom of each subplot), to December. It is important to notice in this figure that the time period is discontinuous between the historical and future period.

Figure 4 Comparison of monthly 99th quantile, at the East English Channel Representative point, obtained from the completely non-stationary model of Eq. (6). The different models are represented in each column of this figure, while the three row corresponds to the three scenarios available. For each sub-figure, the x-axis corresponds to the year, while the y-axis is mapped to the month of the year. This figure shows that the model is able to capture very different time-evolution patterns, and can also reduce to a stationary model (e.g. for EC-EARTH model under SSP1-2.6). See the text for a detailed interpretation.

The figure clearly illustrates the pronounced seasonality of significant wave height (H_s), with consistently higher values during winter (yellow) compared to summer (blue) across all models and time periods, suggesting a shift in seasonal storm patterns. A seasonal contrast in H_s emerges between the historical period and scenario SSP1-2.6, and becomes even more pronounced under SSP5-8.5. Historical H_s values generally range from 3 to 7 m, with gradual transitions represented by thick bands in Fig. 4. Under SSP1-2.6, the range remains similar (3 to 7 m), but seasonal variability increases, as indicated by narrower bands. In contrast, SSP5-8.5 exhibits a wider range (2 to 8 m) and sharper seasonal fluctuations, with H_s shifting from about 2.5 m in August to 5 m in September, for example. These results suggest a trend towards more intense seasonal contrasts, characterized by higher winter values (up to 8 m) and lower summer values (around 2 m), particularly evident when comparing the historical period to SSP5-8.5.

Furthermore, the figure reveals a time-evolution of the seasonality: for most models, the summer period tends to become longer and to peak later in the year, though exceptions exist (AWI, KIOST, and MRI models). Finally, the models successfully reduce to the stationary case when the data shows no evidence of shifting seasonality (e.g., KIOST model under SSP1-2.6).

An additional remark can be made regarding the results of MRI-ESM2-0 under scenario SSP1-2.6, for which the seasonality is very pronounced in the projection period, with a very large shift in the seasonal cycle (see Fig. 5). This behaviour is also observed at the other locations (see the Supplement). As this behaviour is unexpected and cannot be explained by any physical mechanism, this model should perhaps be disregarded when drawing a precise assessment of the evolution of extreme sea states under climate change. However, since the focus of this paper is on the methodology, we decided to retain this model in the remainder of the analysis.

Figure 5 Comparison of monthly 99th quantiles, at the East English Channel Representative point, obtained from the completely non-stationary model of Eq. (6) at different years (2010 for Historical run; 2100 for Scenarios SP1-2.6 and SSP5-8.5). Each figure is a transect of each sub-figure of Fig. 4 at its rightmost end. This figure shows, for each model, the seasonality of monthly maxima at the end of the century (red and orange curves) compared to 2010 (grey curve). The rise in the seasonal amplitude and the shift of the calmer period toward the end of the year can be clearly seen for each model. It is also noticeable that the effect is stronger under the SSP5-8.5 scenario.

Figure 5 illustrates the shift in seasonality discussed earlier by showing, for each model, the monthly extreme quantiles estimated with the non-stationary approach for two reference years: 2010, corresponding to the end of the historical period, and 2100, representing the end of the future period. The most notable feature highlighted by this figure is the seasonal shift observed in most models when comparing the grey curve (2010) with the colored curves (2100 under two scenarios). The latter exhibit lower values during summer and higher values in winter. Another striking aspect is the asymmetry between the beginning and end of the year in future scenarios, with January values exceeding December by approximately 0.2 to 0.5 m, depending on the model. Furthermore, the minimum H_s values tend to occur later in the year, and the adjacent slopes become steeper, indicating a more pronounced and extended summer season compared to the historical period. Overall, all models suggest a shift of the seasonal cycle toward the end of the year, implying a potential redefinition of conventional seasons. For instance, the lowest H_s values now occur in July rather than June. Finally, the comparison between SSP1-2.6 and SSP5-8.5 reveals that the latter amplifies these effects. While all models generally agree on this trend, the KIOST model shows the smallest changes, whereas the MPI model exhibits the largest.

Because of the shifts identified both in the summer and winter values, it may be difficult to assess if the annual return levels also evolve with time, because the two effects can compensate each other. This can be investigated in Fig. 6. To compute the annual return level from the fitted model, we first compute the distribution of the annual maxima, using Eq. (5) of Sect. 3.2.2, by taking the product of the CDF of monthly maxima. This leads to the CDF of the annual maxima, that is then inverted numerically to obtain the return levels, which still depends on the year. Uncertainties are estimated using a Monte-Carlo approach, by simulating from the fitted model.

Figure 6 Comparison of yearly return levels, defined as a quantile of the yearly maxima distribution from Eq. (5), at the East English Channel Representative point. The layout of this plot is the same as for Fig. 4, the y-axis of each sub-figure being now the return level of significant wave height. Each colors corresponds to one projection scenario and the shaded bands are the 95 % confidence intervals computed by parametric bootstrap details in Sect. 3.4.

We see that the shift of the summer and intensification of the winter may lead to an increase of the annual return level for AWI, MRI and also MPI but for SSP5-8.5 only. While some models as ACCESS-CM2, EC-EARTH, IPSL and KIOST can also indicate that there are no long-term evolution of the return levels. And finally, the CMCC model shows a contrasted result with an increase of the return level under SSP1-2.6 and a decrease under SPP5-8.5.

We can also notice a more complex pattern in the MPI model during the historical period, which would not have been captured by a classical linear approach. Looking back at Fig. 4, this behavior results from a shift of the summer season towards the beginning of the year, followed by an intensified winter that fails to compensate, producing an “inverted bell curve”.

In order to compare the estimated values to the ones obtained using a classical stationary model, we estimated the equivalent design condition using Eq. (7). In a stationary climate, this definition would coincide with the 100-year return level. These quantities are compared on Fig. 7, for the different models and scenarios. We included the 100-year return level estimated from annual maxima (points), monthly maxima (triangles) and the proposed non-stationary approach (squares).

Figure 7 Comparison of yearly return levels, defined as a quantile of the yearly maxima at the East English Channel Representative point. Each panel corresponds to a scenario, the x-axis of each sub-figure corresponds to one GCM. The dots corresponds to the return levels computed from a stationary model as in Eq. (3), the Triangles corresponds to the time-stationary model from Eq. (4), and the Squares to the equivalent life-time design point from the full non-stationary models using Eq. (6) and definition in Eq. (7). The errors bars corresponds to the 95 % confidence interval, computed using a parametric bootstrap approach described in Sect. 3.4. The model mean and confidence intervals are computed from each method, by using the values, their spreading and associated uncertainties, assuming they are i.i.d and normally distributed. Some values are missing (e.g. for CMCC-CM2-SR5 model due to a numerical failure of the GEV fit).

The main insight from this figure is that using monthly maxima instead of annual maxima significantly reduces the uncertainty in estimating 100-year return levels. This improvement is largely due to the model being fitted to twelve times more data, despite its increased complexity. The figure also shows strong consistency across methods, with estimates generally falling within each other's confidence intervals. In contrast, the annual maxima approach may yield spurious results because of the limited data available for fitting the GEV model in this case.

For stationary cases as identified before (e.g. KIOST model under scenario SSP1-2.6), all three approaches produce similar estimates, with the fully non-stationary model providing the narrowest confidence intervals. When averaging return levels across all models, we obtain approximately 6.15 m for the historical period, 6.75 m for SSP1-2.6, and 7 m for SSP5-8.5, indicating that future scenarios predict higher return levels.

Lastly, Fig. 8 contains the estimation of the equivalent design condition for each seafront representative point detailed previously (Table 2), and using the methods used throughout the paper. All trends shown are progressive, first between the historical period and SSP1-2.6, then with SSP5-8.5, even though the return periods are very similar for the future scenarios. Only the East Channel location shows a slight decrease, the other representative points tending toward an increase. West Channel location has a slight increase, North Atlantic location has a clear increase of approximately 2.5 m, South Atlantic location has an increase of approximately 1 m, West Mediterranean location has a significant increase of 5 m for the historical scenario and 10 m for future scenarios, and finally East Mediterranean location has an increase of 1 m for the historical scenario and SSP5-8.5 and 2 m for SSP1-2.6. Globally, we can derive the same comments as the previous figure for points located in the English Chanel and in the Atlantic. The non-stationary monthly method reduces uncertainties and has values close to those of the monthly maxima method. However, for representative points in the Mediterranean Sea, the annual maxima method gives the lowest margin of uncertainty. It should be noted that West Mediterranean location has very high H_s values compared to all other representative points, reaching up to 20 m and uncertainties of up to 30 m. For the locations in the Mediterranean sea, the results are more mitigated compared to the English Channel and the Atlantic, which can be explained by a less important seasonality (see Figs. S9 and S10 in the Supplement) and to a less-convincing representation of the sea-states from the GCMs in this area. In fact, the coastal processes are very important and not well represented in the very large-scale climate models used to force the numerical wave model. This has a particular impact on the Mediterranean Sea, which is surrounded by coastlines, compared to the Atlantic Ocean or the English Channel.

Figure 8 Comparison of average 100-year return level over the different GCMs, at the different representative points and the different scenarios, with the different methods proposed. Each line corresponds to one representative point, each column a scenario. The errors bars corresponds to the 95 % confidence interval, computed using a parametric bootstrap approach described in Sect. 3.4.

To conclude, the lifetime design level represents the wave height to not exceed within the structure lifetime of the wind farms. In Fig. 9 the lifetime design level is represented for the English Channel. For all models the lifetime design level is smaller for the historical scenario from 1985–2014 (7.4 m) than for future scenarios from 2081 to 2100, the SSP1-2.6 (8 m) and SSP5-8.5 (8.2 m) are close to each other. The largest uncertainties are for SSP5-8.5. The models with the highest uncertainties are CMCC and MPI for the historical period and SSP5-8.5 and it is KIOST for SSP1-2.6. The investigation of the lifetime design level is interesting through all maritime seafronts (see Supplement – Figs. S13 to S16). The Atlantic coast presents similarities with the English Channel even if values are higher with 11 to 13 m for the historical period and 13 to 15 m for both SSPs. The South part of the Atlantic has the lowest values compared to the North part. And for the Mediterranean Sea SSP1-2.6 seems higher than SSP5-8.5, showing really high uncertainties even if compared to all other seafronts. For the eastern part of the Mediterranean Sea, values are low around 10 m for the historical period and 12 m for the two SSPs, whereas for the West part values are high, around 20 m for all scenarios. The Mediterranean Sea highlights great uncertainties from 4 to 14 m so around 10 m difference between all models for SSP scenarios while uncertainties are lower in the Atlantic and the English Channel with a maximum of 3 m difference between all models.

Figure 9 Comparison of lifetime design level at West English Channel location,, estimated from . This figure shows a Comparison of lifetime design level over the GCMs and the different scenarios, as Fig. 7, but for West English Channel location. Only the the fully non-stationary model of Eq. (6) is represented here for a better comparison of the results.

5 Discussion

The three French maritime seafronts – Atlantic, English Channel, and Mediterranean – display fundamentally different projection behaviors and uncertainty envelopes in climate-driven wave simulations, with significant implications for offshore wind development. Ensemble analyses and regional studies show that the Atlantic and Channel coasts generally produce more coherent signals across models, notably clearer changes in storm-track activity and a pronounced intensification of winter extremes, whereas the Mediterranean sea emerges as both a climate hotspot and as a basin with much larger inter-model spread and several physically implausible outliers in some model-wave ensembles (e.g., extremely large simulated H_s in a few runs) which are likely linked to coarse GCM resolution, poor representation of local cyclogenesis, and unresolved coastal/fetch effects (Hemer et al., 2013; Morim et al., 2019; Cos et al., 2022). Several authors therefore caution that Mediterranean projections must be interpreted with care: ensemble medians often differ from individual model extremes, and values that appear unrealistically high typically arise from biased wind forcing or from insufficiently resolved dynamical drivers rather than from robust physical signals (Giorgi and Lionello, 2008; Davison et al., 2024). By contrast, the Atlantic basin frequently shows a distinct “deepening” of seasonality in extremes, a stronger winter peak and starker summer contrast in H_s, a signal that is reproduced more consistently across multiple modelling studies and wave projections (Hemer et al., 2013; Meucci et al., 2024). Finally, CMIP6 (and earlier CMIP5) model families exhibit substantial inter-model differences in the representation of the large-scale drivers of wave climate (storm tracks, surface wind fields and seasonal variability), with specific models (e.g., KIOST-ESM, some CMCC configurations, MPI-ESM variants) producing divergent wind-forcing fingerprints and therefore divergent wave projections; this heterogeneity has been documented in model evaluation and ranking studies and argues for multi-model, bias-corrected and high-resolution downscaling strategies when deriving local extreme H_s estimates (Meucci et al., 2024; Lorenzo et al., 2025). In summary, while Atlantic and Channel projections provide relatively coherent seasonal signals – particularly a clearer intensification of winter extremes – Mediterranean projections remain more uncertain and prone to isolated, physically questionable extremes. Risk assessments in the Mediterranean sea therefore demand rigorous model selection, downscaling, and plausibility checks rather than blind reliance on raw GCM-driven wave outputs (Hemer et al., 2013; Cos et al., 2022; Davison et al., 2024).

The seasonal variability of extreme significant wave heights shows marked contrasts between France's three maritime seafronts, reflecting both atmospheric circulation patterns and basin morphology. Along the Atlantic coast and in the English Channel, studies by Bricheno and Wolf (2018) and Bulteau et al. (2013) highlight a clear winter dominance of extreme events, associated with intensified westerly storms and North Atlantic depressions, whereas summer remains generally calmer, confirming a strongly seasonal regime. These works suggest that winter maxima are projected to intensify consistent with our results suggesting the increase of severe winters and really low values of H_s in summer. In fact, Bricheno and Wolf (2018) anticipate an increase of annual maxima H_s in the three French seafronts for mid-century and end of century with RCP4.5 and RCP8.5. Overall, all studies agree that the cold season remains dominant for extreme H_s, even though future projections point to a more complex seasonal cycle, with longer and later summers. Chirosca and Rusu (2022), using ERA5 from 2001 to 2020, found that H_s 95 percentile values are increasing in the North Atlantic coast and in the South of the English channel. As Vanem (2015) discerned a shift towards higher extremes in a future wave climate similar to our results in the North Atlantic with RCP4.5 and RCP8.5 comparing 1970–1999 to 2071–2100. An intensification of the 100 return level is found in the Atlantic and the English Channel by Bulteau et al. (2013) using the POT and GPD using numerical wave hindcast BOBWA database from 1958 to 2002 in accordance with our results. Within this context, our results showing a tendency towards more intense winters and calmer yet delayed and prolonged summers, together with an increase in the return period of 100-year H_s are consistent with these previous findings, suggesting a future regime characterized by stronger seasonal contrasts and a temporal redistribution of calm conditions towards early autumn.

The adoption of non-stationary extreme-value models has emerged as a critical advance in assessing changes in wave extremes, offering superior performance compared to stationary assumptions under evolving climate conditions. Traditional stationary models presume that the statistical characteristics of wave climate remain constant over time, an assumption increasingly invalid in the presence of long-term trends, changing storm regimes or evolving wind and fetch patterns as there are important in French seafronts. Vanem (2015) showed that including time-varying location and scale parameters in extreme-value modelling substantially improves the estimation of return levels (e.g., centennial H_s) under non-stationary forcing and reduces bias in projections. For RCP4.5, Vanem (2015) reported an uncertainty of 13.15 m for the 100-year return level, which is substantially larger than the uncertainty obtained with our method (maximum of 10 m). This reduction in uncertainty is likely attributable to the larger amount of data used in our analysis, as the increased sample size across all models leads to more robust extreme value estimates. Turki et al. (2020) highlighted that non-stationary approaches can account for local, seasonal and directional variability in wave extremes across basins (e.g., North Atlantic, Mediterranean) where stationary assumptions fail to reproduce observed shifts in intensity, timing and frequency of events. More recently, Ewans and Jonathan (2023) emphasized that non-stationary extreme-value analyses are indispensable for climate adaptation and coastal engineering, since they quantify how the probability of rare, extreme wave events evolves under future scenarios rather than assuming historical stability. Complementing these findings, earlier studies (Mendez et al., 2008; Menendez et al., 2009) illustrated how modelling non-stationarity in seasonal parameters yields more realistic return value estimates for marine extremes by capturing seasonal shifts in storminess and fetch, thereby supporting the use of non-stationary parameterization for extreme H_s. Altogether, these works converge on the conclusion that non-stationary models not only offer a more realistic statistical representation of changing wave climates but are also essential for deriving credible return periods and extreme-value projections in a non-equilibrium climate system. Our study on wave extremes therefore provides a new, effective method for monthly maxima with reduced bias and taking climate change into account. Using monthly maxima instead of annual maxima provides a finer temporal resolution that better captures intra-annual variability in wave extremes and their seasonal drivers. This approach allows the detection of shifts in the timing and intensity of extreme events, an information that is lost when only annual maxima are considered. Moreover, monthly maxima improve the robustness of non-stationary analyses, as they provide a larger sample size and make it possible to explicitly model seasonal covariates (e.g., storm seasonality, climate indices), leading to more accurate and physically meaningful estimates of extreme return levels.

6 Conclusion and perspectives

This study introduces a non-parametric approach based on monthly maxima to estimate return levels for monthly and annual periods from 8 GCMs, explicitly accounting for climate change effects. A key contribution lies in proposing a new method for calculating design levels – equivalent to return periods – under non-stationary conditions driven by a changing climate. The consistency of our results with existing literature reinforces the credibility and the relevance of the approach.

By leveraging monthly maxima, we substantially reduce uncertainty in extreme-value estimation. The parsimonious, data-driven model captures seasonality and its potential shifts, as well as changes in extremes themselves, making it particularly suited for assessing rare events critical to offshore wind farm design.

Our findings reveal coherent signals of stronger seasonal contrasts under future climates, with intensified winters and calmer yet prolonged summers, especially in the Atlantic and English Channel. Despite significant inter-model uncertainty – most pronounced in the Mediterranean sea – the ensemble generally indicates higher extreme sea states under both SSP1-2.6 and SSP5-8.5 scenarios.

This work opens several avenues for future research: addressing biases in extreme-value estimation to improve historical consistency; refining wave models along the French coastline for better representation of local dynamics; developing parametric approaches to seasonality and its evolution; and exploring joint distributions of significant wave height (H_s) and peak period (T_p) under future climate scenarios to enhance the robustness of offshore wind design.