Parameterization-Driven Uncertainties in Single-Forcing, Single-Model Wave Climate Projections from a CMIP6-Derived Dynamic Ensemble

This study is focused on the impact of different parameterizations in the state-of-the art wave model WAVEWATCH3 (WW3) in describing the present climate and future wave climate projections. We have used a Coupled Model Intercomparison Project Phase 6 (CMIP6)-derived single-wind forcing (from EC-EARTH) to produce a dynamic wind-wave climate ensemble for its historic (1995–2014) and future (2081–2100) periods. We discuss the uncertainty due to the wave model (intra-model uncertainty) in simulating the present and future wave climate. The historical wave climate runs were compared against the ERA5 reanalysis and found to be in good agreement for the significant wave height. This gives a good degree of confidence to investigate the intra-model uncertainty in WW3 using the available physics packages such as ST2, ST3, ST4, and ST6. In general, for the historic period, ST3 and ST4 physics packages perform better in the tropics whereas ST6 performs better in the extratropics, based on M-Score performance assessment. The study also reveals that the extratropical South Indian Ocean and tropical eastern South Pacific areas exhibit a larger amount of uncertainty, mainly induced by the ST2 physics package. The results of this study shed new light on the impacts associated with the use of multiple physics parameterizations in wave climate ensembles, an issue that has not received the necessary attention in scientific literature.


Introduction
Surface gravity waves (hereafter called waves) play an important role in several aspects of ocean and coastal dynamics. Extreme waves due to intense storms are of great risk to coastal communities [1,2] particularly in the context of global sea level rise. Along with tides and storm surges, waves play a key role in coastal processes (e.g., [3,4]), being one of the main drivers of coastal erosion [5,6], flooding [7,8], sediment transport [9] and thereby changing the shape of the coast (IPCC report: [10][11][12][13]. Most importantly, waves influence the entire climate system due to their complex interactions with the atmosphere, sea ice, and the underlying ocean [14]. Understanding wave climate and its future evolution is a challenging task. Variations in surface wind fields, marine ice coverage, storm patterns (intensity and frequency) due to climate change will for sure have an impact on the future wave climate. In turn, relative sea levels, sediment dynamics, and shoreline processes will be affected. Based on the satellite altimeter measurements, Young et al. [15] observed an increase in the significant wave heights (H S ) throughout the global ocean, in the recent past, driven by the increasing wind speeds, and this trend is expected to continue under future greenhouse gas emission Forecast System (IFS) cycle CY41R2 release (ECMWF, 2016). The ERA5 was produced with an improved data assimilation technique (4D-Var scheme), carrying out hourly assimilation of altimeter-measured H S data from 1991 onwards. The horizontal resolutions of the atmospheric and wave components in ERA5 are about 31 km (0.25 • × 0.25 • ), and 40 km (0.36 • × 0.36 • ), respectively, both with an hourly time resolution. ERA5 wave data were generated using the WAM wave model, set up with a spectral resolution of 30 logarithmically spaced frequency bins (from 0.03453 Hz to 0.5478 Hz) and 24 directional bins, each with 15 • . The bottom topography used in ERA5 is based on the ETOPO2 (NGDC, 2006) dataset. Additional details about the ERA5 reanalysis can be found in Copernicus Climate Change Service (2017) and in Hersbach et al. [37].

The EC-Earth3 GCM
The EC-Earth is a widely used GCM in global and regional climate studies collaboratively developed by the European Consortium [38]. The current generation of the model, EC-Earth3, has been developed after CMIP5 and it is used in its version 3.3 for CMIP6 experiments [39]. EC-Earth3 comprises model components describing atmosphere, ocean, sea ice, land surface, dynamic vegetation, atmospheric composition, ocean biogeochemistry, and the Greenland ice sheet. The domains are covered by ECMWF's IFS cycle CY36R4 coupled with NEMO 3.6 (ocean), LIM3 (sea ice), and PISCES (Biogeochemical) models. Dynamical vegetation, land use, and terrestrial biogeochemistry are provided by the Lund-Potsdam-Jena General Ecosystem Simulator (LPJ-GUESS). For further details please refer to Döscher, et al. [38]. In this study, only one set of U 10 forcing data was considered, in order to minimize the uncertainty arising due to the use of different forcings, as different climate models (or even different runs from the same climate model) show different performances [40].
The EC-Earth was shown to be one of the best GCMs in simulating the historical U 10 among the CMIP5-based models. Particularly in the Southern Hemisphere, Casas-Prat et al. [40] showed that EC-Earth based wave climate simulations perform best among the remaining five GCMs used in their study. Such good performance was shown to be related to a good representation of the U 10 , with high correlation coefficients (above 0.90) and low biases. The performance of EC-Earth in simulating sea ice was also shown to be the best among other CMIP5 models [40,41] for both the Arctic and Antarctic.

The WW3 Wave Model
The WW3 is a third-generation spectral wave model that has been widely used for operational wave forecasting, research, and engineering applications. In this study, the latest public version of the WW3 (WAVEWATCH III 6.07; [36]) is used to simulate wave parameters across the global ocean. WW3 describes the evolution of the wave action density spectrum, N(k, θ, x, t), which is a function of wavenumber (k), direction (θ), time (t), and space (x). The action density (N) is conserved in the presence of ocean currents and the governing equation for the conservation of the wave action in its simplest form can be written as in Equation (1): where C g is the wave group velocity, and σ is the intrinsic frequency, obtained from the dispersion relation for ocean waves. The action density spectrum (N) is calculated by dividing the energy spectrum, F(k, θ, x, t), by the intrinsic frequency, σ. The S tot term on the right-hand side of Equation (1) represents all energy fluxes contributing to wind-wave evolution [42], and is comprised of parameterized source-sink terms (simply called source terms; STs) that represent several wave processes, mainly wave growth due to the action of wind (S in ), nonlinear wave-wave interactions (S nl ), scattering due to wave-bottom interactions (S sc ), triad wave-wave interactions (S tr ), and dissipation Climate 2022, 10, 51 4 of 23 due to wave breaking (S ds ), wave-bottom interactions (S bot ), wave-ice interactions (S ice ) etc. Thus, the general form of the total STs can be written as in Equation (2): S tot = S in + S nl + S ds + S sc + S tr + S ice + . . . + S user (2) where S user is the slot for additional user-defined STs.
Based on different combinations of S in and S ds Equation (2), various ST packages are available in WW3. The details of the different source term packages and methods used for this study are described below in Section 3.2.

Experimental Setup
The EC-Earth-WW3 wave climate simulations presented here were generated considering a geographical domain from 180 • W to 180 • E and from 80.64 • N to 80.36 • S, and a horizontal resolution of 0.703125 • × 0.7 • (longitude × latitude). The rectilinear grids were developed using ETOPO-1 bathymetry [43], together with v1.10 of the Global Self-Consistent Hierarchical High-Resolution Shoreline (GSHHS) Database. Three files were created: a bathymetry, a mask, and an obstruction grid which account for wave attenuation by unresolved islands using the gridgen software package [44,45]. The computational (global) time step in WW3 was set to three hours, using a spectral resolution of 29 frequencies ranging from 0.0350 Hz to 0.5047 Hz (logarithmic distribution with increment factor 1.1) and 24 directions (a constant increment of 15 • ).
In WW3, the physics and numerics are defined by switches during the compile stage of the wave model [46]. The common model settings used are described here with corresponding switches. The standard explicit fractional time-stepping scheme is used with a minimum source term time step of 30 s. We used the third-order Ultimate Quickest (UQ) propagation scheme along with the averaging technique (PR3) for garden sprinkler reduction [47]. Nonlinear wave-wave interactions are modeled using the discrete interaction approximation (DIA; [48]) by activating the switch NL1. A linear input source term (LN1), the parameterization of Cavaleri and Malanotte-Rizzoli [49] along with a filter for low-frequency energy from Tolman [50] are used for the consistent spin-up of the model from calm conditions, and to improve initial wave growth behavior. The bottom friction is modelled using JONSWAP (Joint North Sea Wave Project; [51]) parameterization (BT1). Depth-induced wave breaking is accounted for by using the Battjes and Janssen [52] formulation (DB1) with a Miche-style shallow water limiter (MLIM) for maximum wave height (H max ). Reflections by shorelines and icebergs are deactivated (REF0) and no bottom scattering (BS0) and sea ice scattering (IS0) are considered. We used a simple ice blocking scheme (IC0) for the wave damping due to sea ice with the default sea ice cut-off concentration of 50% (i.e., discontinuous treatment of ice, which means waves are generated/propagated when sea ice concentration (SIC) is below 50%). To generate each of the four ensemble members, four input-dissipation ST packages were used (by activating the relevant switches), namely ST2 [53], ST3 [54,55], ST4 [56], and ST6 [57][58][59]. A short description of each source term package and commonly used schemes are given below. Further details can be found in the WW3 user manual [36].

The ST2 Parameterization
The ST2 parameterization is based on Tolman and Chalikov [53], where the input source term is taken from Chalikov and Belevich [60] and Chalikov [61], the non-linear interaction term is based on Hasselmann et al. [48], and the dissipation term is set with two dissipation constituents, one for low frequency waves and the other for the high frequency tail of the spectrum. The strength of the non-linear term is reduced by tuning the numerical constant 'C' (from 2.8. to 1) in order to minimize the errors as it leads to an overestimation of the interaction at higher frequencies. Tolman [47] found, in a later testing of this ST at a global scale, that its swell dissipation is overestimated, and, as a result, this source term underestimates deep-ocean wave growth under stable atmospheric conditions. In order to Climate 2022, 10, 51 5 of 23 mitigate this problem, a stability correction can be activated with the STAB2 switch. In this study we are using the default parameter settings of the ST2 package implemented in the WW3 v6.07.

The ST3 Parameterization
This parameterization is based on WAM4 [62] as described in Janssen [55]. The windwave interaction source term described in this package was developed by Janssen [63], based on the theory of Miles [64], and the linear swell dissipation is based on Janssen [55]. The dissipation source term was developed by Bidlot et al. [54] and its later modification is known as "BJA". This ST is similar to the one in the EC-WAM model (the reader is referred to the 2019 WW3 manual [36]). Unfortunately, these parameterizations are sensitive to swell, and an increase in swell height typically reduces dissipation at the wind sea peak. Here, the default values of the namelist parameters corresponding to BJA [54] are used, which generally provide better results than the corresponding WAM4 parameterizations [36].

The ST4 Parameterization
Ardhuin et al. [56] proposed a package for wind input and dissipation, named ST4. The wind input is adapted from the Janssen [63] formulation, with an important reduction of input at high frequencies. A balance is required to achieve with the dissipation due to white-capping. The reduction of the input to the waves under high wind speeds is achieved by including a tunable sheltering factor that reduces the unrealistic large drag coefficients. The swell feedback on wind input is according to Ardhuin et al. [65]. We use the default model ST4 values with a β max value of 1.43 corresponding to the TEST471 value, which generally provides the best results at global scale when using ECMWF winds [36]. The term β max is a non-dimensional growth parameter that controls the wind to wave growth. Stopa [66] calibrated 12 different wind products by adjusting the β max parameter and found that modifying β max reduces errors in the wave field, thus the relevance of choosing a proper β max .

The ST6 Parameterization
It is the latest addition to the source term packages implemented in WW3 which include wind input and white-capping dissipation terms based on measurements taken at Lake George, Australia [67]. Additionally, a parameterization of swell decay based on field observations is included. The white-capping dissipation characteristics are similar to those of ST4 but are implemented differently [59,68]. A separate term is used to account for the losses resulting from the interaction of waves with opposing directions and is based on laboratory experiments of Donelan [69]. The ST6 package has been calibrated with flux parameterization and implemented as switch FLX4. Bulk adjustment to the wind field (for any uniform bias in the wind field) can be achieved by rescaling the drag parameterization FLX4 through a namelist parameter, CDFAC. This has a similar effect to the tuning of the β max parameter in ST4. One should also note that the computation of the friction velocity (u * ) is different in ST4 and ST6. In ST4, u * is calculated via a stress table, while in ST6, u * is computed using drag coefficients taken from Hwang et al. [70], also acting as a limit for the calculation of the wave-supported stress. The spectral wave growth γ parameter in ST6 (present in the S in term, Equation (2)) is a function of scaling wind speed following Rogers et al. [58] (U S = U 10 ≈ 28u * ) and thus the relevance of u * . Later, Rogers [71] found that using U S = 32u * could improve model skills and Liu et al. [72] implemented it in WW3 v6.07 along with a new set of parameters (Table 2.8; [36]) carrying out a thorough re-calibration of the ST6. The updated ST6 package performs well in simulating wave parameters (e.g., H S ) and showing an improved skill in the estimation of the high-frequency tail of the energy spectrum.
As described above, the four different input/dissipation source term packages (ST2/3/4/6) describe the wind generation and dissipation differently. The general form of wind input source term for ST3, ST4, and ST6 is similar, but ST4 and ST6 are more focused Climate 2022, 10, 51 6 of 23 on improving the dissipation source function in the model. A brief comparison of all the STs used for this study can also be found in Table S1.

Methodology
The ensemble is composed of four members, being the differentiating factor between each member the WW3 parameterization. Not only were the forcing data kept constant for all the ensemble members, but also their spatial and temporal resolutions, to minimize any uncertainties arising from interpolating the winds to a new model grid by the wave model.
The ensemble mean is computed considering a "democratic" approach: the unweighted mean of the four ensemble members. The ensemble's performance in representing the present wave climate is first evaluated by comparison with ERA5 during the 1995-2014 period (henceforth named PC20). The ensemble future wave climate projections and their uncertainty are then analyzed for the 2081-2100 period (henceforth named FC21). For convenience, the present climate historic ensemble is named PC20-E and future wave climate ensemble is named FC21-E, and the individual ensemble members are named PC20-ST2, PC20-ST3, PC20-ST4, PC20-ST6 and FC21-ST2, FC21-ST3, FC21-ST4, FC21-ST6 for present and future climate respectively. The 3-hourly U 10 and H S parameters have been processed for annual and seasonal means: DJF (December, January and February), MAM (March, April, May), JJA (June, July, August), and SON (September, October, November). The uncertainty analysis for the mean H S (H S ) and the 95% percentile H S (H 95 s ) was carried at global scale and for 13 different sub-areas. These regional sub-areas were chosen in accordance to Alves [73] and are detailed in Figure 1 and Table S2.
in simulating wave parameters (e.g., ) and showing an improved skill in the estimation of the high-frequency tail of the energy spectrum.
As described above, the four different input/dissipation source term packages (ST2/3/4/6) describe the wind generation and dissipation differently. The general form of wind input source term for ST3, ST4, and ST6 is similar, but ST4 and ST6 are more focused on improving the dissipation source function in the model. A brief comparison of all the STs used for this study can also be found in Table S1.

Methodology
The ensemble is composed of four members, being the differentiating factor between each member the WW3 parameterization. Not only were the forcing data kept constant for all the ensemble members, but also their spatial and temporal resolutions, to minimize any uncertainties arising from interpolating the winds to a new model grid by the wave model.
The ensemble mean is computed considering a "democratic" approach: the unweighted mean of the four ensemble members. The ensemble's performance in representing the present wave climate is first evaluated by comparison with ERA5 during the 1995-2014 period (henceforth named PC20). The ensemble future wave climate projections and their uncertainty are then analyzed for the 2081-2100 period (henceforth named FC21). For convenience, the present climate historic ensemble is named PC20-E and future wave climate ensemble is named FC21-E, and the individual ensemble members are named PC20-ST2, PC20-ST3, PC20-ST4, PC20-ST6 and FC21-ST2, FC21-ST3, FC21-ST4, FC21-ST6 for present and future climate respectively. The 3-hourly and parameters have been processed for annual and seasonal means: DJF (December, January and February), MAM (March, April, May), JJA (June, July, August), and SON (September, October, November). The uncertainty analysis for the mean ( ) and the 95% percentile ( ) was carried at global scale and for 13 different sub-areas. These regional sub-areas were chosen in accordance to Alves [73] and are detailed in Figure 1 and Table S2.  North Atlantic and Arabian Sea), showing a low bias (generally below 12%). The general trend is for EC-Earth to underestimate wind speeds in the equatorial ocean, and this agrees with previous study by Nolan and McKinstry [74]. Positive biases can however be found along the southern region of South China Sea (Maritime Continent). The EC-Earth winds show better performance along the mid-latitudes of the Southern Hemisphere areas (with biases below 4%). It is also observed that the performance of the EC-Earth is very good in the regions of active tropical storm areas, generally showing a low bias of less than 12%. The mean absolute difference between the EC-Earth winds and ERA5 winds also show similar pattern, exhibits a maximum negative bias (~2.2 m/s) along TSAO region ( Figure S1). Overall, the performance of EC-Earth in simulating the U 10 during the PC20 time-slice is good and well within the confidence level.

Surface Wind Speeds
The normalized differences (in %) between the EC-Earth and the ERA5 annual mean are shown in Figure 2 (normalized differences: EC-Earth minus ERA5 normalized by ERA5). In the Northern Hemisphere, the EC-Earth surface winds are comparable with the ERA5 in most of the extratropical areas and in some regions in the tropics (eastern tropical North Atlantic and Arabian Sea), showing a low bias (generally below 12%). The general trend is for EC-Earth to underestimate wind speeds in the equatorial ocean, and this agrees with previous study by Nolan and McKinstry [74]. Positive biases can however be found along the southern region of South China Sea (Maritime Continent). The EC-Earth winds show better performance along the mid-latitudes of the Southern Hemisphere areas (with biases below 4%). It is also observed that the performance of the EC-Earth is very good in the regions of active tropical storm areas, generally showing a low bias of less than 12%. The mean absolute difference between the EC-Earth winds and ERA5 winds also show similar pattern, exhibits a maximum negative bias (~2.2 m/s) along TSAO region ( Figure S1). Overall, the performance of EC-Earth in simulating the during the PC20 time-slice is good and well within the confidence level.

Wind-Waves
The present wave climate simulation ensemble (PC20-E) has been evaluated using ERA5 ( ). The normalized differences (in %) between the PC20-E and ERA5 (PC20-E minus ERA5 normalized by ERA5) for annual and annual are shown in Figure 3. In general, the PC20-E properly describes the wave climate at a global scale, being the performance of WW3 in simulating annual quite reasonable, with biases generally below 12% (Figure 3a). When compared to ERA5, over tropical areas, the PC20-E tends to underestimate the , with biases varying between −20% over the tropical Atlantic, −28% over eastern side of the Maritime Continent, and −4% to −12% over tropical western South Pacific (TWSP) which are consistent with similar negative biases in wind speeds ( Figure  2). Similarly, extratropical North Pacific (ETNP) region shows a positive bias (~5%) for ,

Wind-Waves
The present wave climate simulation ensemble (PC20-E) has been evaluated using ERA5 (H S ). The normalized differences (in %) between the PC20-E and ERA5 (PC20-E minus ERA5 normalized by ERA5) for annual H S and annual H 95 s are shown in Figure 3. In general, the PC20-E properly describes the wave climate at a global scale, being the performance of WW3 in simulating annual H S quite reasonable, with biases generally below 12% (Figure 3a). When compared to ERA5, over tropical areas, the PC20-E tends to underestimate the H S , with biases varying between −20% over the tropical Atlantic, −28% over eastern side of the Maritime Continent, and −4% to −12% over tropical western South Pacific (TWSP) which are consistent with similar negative biases in wind speeds ( Figure 2). Similarly, extratropical North Pacific (ETNP) region shows a positive bias (~5%) for H S , consistent with similar biases for the wind speeds. However, in some regions of the Indian Ocean, South China Sea, and Philippine Sea, a small negative bias (0.3 m, see Figure S2a) is seen for H S , together with positive biases in wind speeds~1 m/s ( Figure S1). These small differences could be partly due to the coarse resolution of the EC-Earth winds compared to ERA5 and partly due to the different source term packages used in this study.
Nevertheless, PC20-E shows good performance in simulating the present wave climate along the Southern Ocean. In general, the Northern Hemisphere has a relatively larger variability in wave climate conditions than the Southern Hemisphere, in agreement with the earlier wave climate studies (e.g., [75,76]). differences could be partly due to the coarse resolution of the EC-Earth winds com to ERA5 and partly due to the different source term packages used in this study. N theless, PC20-E shows good performance in simulating the present wave climate the Southern Ocean. In general, the Northern Hemisphere has a relatively larger var ity in wave climate conditions than the Southern Hemisphere, in agreement with th lier wave climate studies (e.g., [75,76]   Additionally, the treatment of sea ice in the ECWAM model (which is used for ERA5 wave simulations) is different from the WW3 setup used here; however, we intended to keep it as close as the ECWAM by using the simple ice blocking scheme (see Section 3.2). Enhanced positive biases are observed in the extratropical North Atlantic and North Pacific (ETNA and ETNP areas) compared to the H S (Figure 3a) illustrating a large difference between the PC20-E H 95 s and ERA5. This may be because ERA5 tends to underestimate the extreme wave events [77]. Overall, the PC20-E shows a good agreement in simulating the ERA5 present wave climate. To assess the global and regional (13 sub domain) performance of PC20-E and its ensemble members, a non-dimensional arcsin Mielke Measure or M-score [78,79] has been applied for H S and H 95 s . This metric has been used widely, for example, by Hemer and Trenham [80] and Semedo et al. [76] to measure the skill of wave climate simulations. For the present wave climate simulation (PC20) and reanalysis (ERA), the M-score is defined as where V is the spatial variance and G is the spatial mean of wave height and MSE is the mean-square error. The M-Score varies from zero (no skill) to a maximum hypothetical score of 1000 (i.e., MSE = 0). The skill assessment for H S is shown in Table 1. The highest M-Score value for PC20-E is found in the ETNP region (941), being the annual global average 930. The lowest annual M-Score (<700) are seen over the tropical South Atlantic (TSAO) and tropical eastern South Pacific (TESP). In general, PC20-E shows high M-Scores over the extratropical region of both the hemispheres. Similar performance is observed for all the ensemble members, except for PC20-ST2, with a low M-Score (695) in the ETSP region. Global performance of all the individual ensemble members shows that the PC20-ST4 (926) is the best among all. However, individual ensemble members show a better performance (best M-Score among all) in different sub regions. For example, PC20-ST6 performs well over most part of the extratropical areas (ETSA, ETSI and ETNP), while PC20-ST4 performs better in the tropics, mainly over the eastern Pacific and Indian ocean basins (TNEP, TESP, TNIO and TSIO), but also in the ETSP region. In the tropical Atlantic Ocean (TNAO and TSAO), ETNA and tropical western Pacific areas (TWSP and TWNP), PC20-ST3 performs better among the ensemble members. The skill assessment for the 95% percentile wave height also shows similar performance as that of the H S (Table S3).

Annual and Seasonal Mean H S
Annual ensemble H S and individual ensemble member anomalies for PC20 are shown in the left panels of Figure 4. The highest annual H S is found in the ETSI area (Figure 4a) and the seasonal annual maximum is observed during DJF and JJA along ETNA and ETSI area ( Figure S3b-d).  The PC20-ST2 shows a negative anomaly for the entire global ocean (Figure 4c) for annual as well as for seasonal ( Figure S3e-h) H S . The annual and seasonal anomalies are higher in the region of higher wave height (Figure 4a), which means PC20-ST2 largely underestimates the wave height during extreme wave height, an inherent drawback exist in ST2 (Table S1). PC20-ST2 has, in fact, the greatest global annual H S anomaly (−0.34 m) among all the ensemble members (Table 2). Seasonal global mean anomalies also show a similar trend (Figure S4a-d). PC20-ST2, in general, shows greater anomalies in the higher latitudes, down to −0.45 m in the ETSI area (Table 2). In the Northern Hemisphere (ETNA and ETNP), during DJF, PC20-ST2 peaks at −0.47 m ( Figure S4). Similar negative anomaly is observed in the Southern Hemisphere (ETSP and ETSI) during MAM and during SON at ETSI (Figure S4b,d). Here, PC20-ST2 exhibits its greatest seasonal anomaly of −0.51 m, in JJA ( Figure S4c). For PC20-ST3, the annual H S shows a positive anomaly for the entire global ocean. However, the core of the maximum positive anomaly (~0.55 m) resides over the tropical Eastern Pacific sector (TENP and TESP). The tropical eastern Indian Ocean, closer to Indonesia, also shows a positive annual anomaly of~0.45 m (Figure 4e). Seasonal mean anomaly shows that both DJF ( Figure S1i Figure S1k,l). The PC20-ST3 exhibits, however, the greatest regional annual mean anomaly in the tropical latitudes, along TESP, at 0.43 m ( Table 2). All seasons show similar behavior with the greatest seasonal mean anomaly visible during DJF, at 0.55 m (Figure S4a-d). In general, for the entire global ocean, PC20-ST4 shows lower anomalies: 0.15 m over the tropical oceans, 0.25 m over extratropical oceans, and the highest anomaly (0.35 m) occurring along the ETSI region (Figure 4g). Seasonal H S also shows a similar pattern to that of the annual H S (Figure S3m-p). Moreover, along the ETSI region, among all the seasons, JJA exhibits the highest H S anomaly (~0.75 m; Figure S3o). Similar to PC20-ST2, at ETSI, PC20-ST4 shows its maximum regional annual mean anomaly of 0.20 m (Table 2) and seasonal mean anomaly of 0.23 m during JJA ( Figure S4c). PC20-ST4 shows a greater positive anomaly in the extratropical Northern Hemisphere in DJF and in the Southern Hemisphere during MAM to SON. PC20-ST6, shows a negative anomaly in the tropics and a positive anomaly in the higher latitude areas (Figure 4i). During, DJF greater positive anomalies (~0.25 m) are seen over the extratropical latitudes in the Northern Hemisphere ( Figure S3q). On the other hand, during JJA, the greatest anomalies (~0.65 m) occur at ETSI. Similar to PC20-ST3, the maximum annual mean anomaly for PC20-ST6 is visible over the tropical eastern South Pacific Ocean (TESP) at −0.16 m (Table 2). TESP also shows the greatest anomaly for DJF and MAM (−0.22 m and −0.21 m, respectively; Figure S4a,b). During JJA, PC20-ST6 exhibits a positive anomaly of 0.13 m in ETSI. Conversely, during SON, a negative anomaly of 0.13 m is visible in TENP ( Figure S4d).
The future projected wave climate ensemble (FC21-E) shows similar pattern as that of PC20-E with a slight change over the higher latitudes for both hemispheres (Figure 4b and Table 3). An increase in the annual H S in the Southern Hemisphere, mostly along the ETSP whereas a slight reduction in the annual H S in the northern extratropical regions is observed (Figure 4b). Increased fetch due to the large ice-free region in the ETSP area and the projected increase in wind speed are the primary reason for this increase in the annual H S in this region. This slight reduction in the annual H S compared to PC20-E in the Northern hemisphere is mainly seen during DJF season ( Figure S5a  Annual H 95 s for PC20-E and the respective ensemble member anomalies are shown in the left panels of Figure 5. The annual H 95 s exhibits the highest values at ETSI and ETNA regions, up to 8 m (Figure 5a). Seasonal H 95 s shows that the Northern Hemisphere has its maximum H 95 s (about 10 m) during DJF over ETNA ( Figure S6a) and the Southern Hemisphere exhibits its maximum over ETSI during JJA ( Figure S6c). The PC20-ST2 shows an annual anomaly of about −1.1 m (Figure 5c) and a seasonal anomaly of −1.3 m in the higher latitudes of the Northern Hemisphere (ETSI and ETNA) during DJF ( Figure S6e). PC20-ST2 shows the highest global annual H 95 s anomaly (−0.54 m) among all the ensemble members (Table 4). In general, it exhibits higher anomaly in the higher latitudes with a maximum of −0.77 m over ETSI region. In the Northern Hemisphere, a negative anomaly of 0.87 m is noticed over ETNA during DJF and in the Southern Hemisphere, mainly over ETSI, a negative anomaly of more than 0.74 m is observed during MAM to SON, peaking at (−0.86 m) during JJA ( Figure S7a-d). PC20-ST3 exhibits a higher annual anomaly over the tropical areas of the eastern Pacific Ocean, but also in some parts of the Indian Ocean and ETNA (Figure 5e). The annual global mean anomaly for PC20-ST3 is 0.29 m with a regional maximum (0.48 m) at TESP (Table 4) In the Southern Hemisphere extratropical latitudes, PC20-ST6 exhibits a smaller anomaly relative to PC20-ST4, similar behavior is also noticed in all the seasons except during DJF ( Figure S6n-p). PC20-ST4 shows its highest regional annual and seasonal mean anomalies (0.30 m and 0.31 m, respectively) at the ETSI area (Table 4; Figure S7a-d). However, the highest seasonal mean anomaly (−0.25 m) for PC20-ST6 is observed over TESP during DJF, and slightly less (−0.24 m) during MAM.   The FC21-E H 95 s shows a similar pattern as that of PC20-E, despite a slight change over the higher latitudes for both hemispheres, where it exhibits the highest values (about 8 m) over ETSI (Figure 5b and Table 5). Compared to PC20-E, an increase in the annual H 95 s is visible over large portions of the ETSI area and slight reduction (about 1 m) is observed at ETNA. This slight reduction in the Northern Hemisphere is visible in all seasons except JJA and the increase in the Southern Hemisphere is visible in all seasons except DJF (Figure S8a-d). Similarly, the change in the annual and seasonal mean anomalies for FC21-E compared to PC20-E is only visible in FC21-ST2 and FC21-ST3 (Figures 5 and S8). FC21-ST2 in general shows greater annual mean anomalies in the higher latitudes reaching the maximum (−0.79 m) in the ETSI area (Figure 5d and Table 5). Seasonal mean anomalies also show the greatest values (−0.90 m) over the same region during JJA ( Figure S7g). For FC21-ST3, unlike FC21-ST2, the greatest positive seasonal mean anomalies (more than 0.40 m) are observed in the tropical latitudes (mainly in the Pacific Ocean) in all seasons, topping at 0.69 m along TESP during DJF ( Figure S7e-h). FC21-ST4 shows lower annual and seasonal mean anomalies compared to FC21-ST2 and FC21-ST3, but slightly higher seasonal anomalies are visible in the extratropical latitudes for all seasons, reaching a maximum of 0.34 m along ETSP during JJA ( Figure S7g). FC21-ST6 also exhibits the lowest annual and seasonal mean anomalies among all the ensemble members, similarly to PC20-ST6. Generally, negative anomalies are visible in the tropical areas, down to −0.27 m during DJF at TESP, whereas positive anomalies are visible in the higher latitudes, topping at 0.18 m during JJA at ETSP.

Uncertainty Due to Different STs
The global and regional inter-annual H S values for PC20-E, FC21-E, and each of the individual ensemble members along with the associated spread (the range between the minimum and maximum values) are shown in Figure 6. ST2 shows a larger negative deviation from the ensemble mean, globally and extra tropical Southern Ocean regions (ETSI, ETSP and ETSA). This larger negative deviation indicates a larger underestimation of wave height, which is partly associated with poor dissipation source term and partly due to the drag-dependent wind input source term present in the ST2 (Table S1). ST3 shows a positive deviation from the ensemble mean, globally and in all the sub regions, relatively less in magnitude compared to ST2. In general, ST3 shows larger deviation for the tropical ocean basins. As mentioned in Section 3.1.2, ST3 is sensitive to swell, an underestimation of swell dissipation will lead to an overestimation of wave heights (Table S1) and this could be the reason for a larger positive deviation in the tropics. ST3 and ST4 show almost similar anomalies in the extratropical ocean basins (ETNA, ETSA, ETNP, ETSP, and ETSI) and in some parts of the tropical ocean such as TWNP and TNAO. ST6 exhibits almost negligible anomaly in the Northern Hemisphere extratropical ocean (ETNA and ETNP); a slight positive anomaly in the Southern Hemisphere extratropical ocean (ETSP, ETSA and ETSI) and a slight negative anomaly in the tropics (TWNP and TNAO). TESP is the only region where ST4 has the least anomaly among all the members. The ensemble means and its spread of global and regional annual H 95 s for PC20 and FC21 also show similar behavior as that of H S (Figure 7). A reduction in the anomaly for ST4 particularly in the regions of North Atlantic (TNAO and ETNA) and most part of the Pacific Ocean (TENP, ETNP and TESP) is observed compared to H S whereas the anomaly is increased in some of the extratropical ocean basin (e.g., ETNA, ETSI, ETNP) for ST2 causing a larger uncertainty for H 95 s . In general, the uncertainty has increased for H 95 s ; in most of the regions ST3 shows similar behavior as that of H S ; small reduction in the anomaly for ST6 in the tropical ocean basins (TNAO, TSAO, TWNP, and TENP) and slight increase in the anomaly in the extratropical ocean basins such as ETNA, ETSA, ETNP, and ETSP are observed. Thus, the extratropical southern hemisphere ocean (ETSA, ETSP, and ETSI) showed larger intra-model uncertainty compared to northern hemisphere ocean (ETNA and ETNP) due to the larger negative anomalies contributed by ST2. Similarly, most part of the tropical ocean shows larger intra-model uncertainty, equally contributed by the positive anomalies due to ST3 and negative anomalies due to ST2. Among all the 13-sub regions, the tropical north Atlantic and Western North Pacific are the regions of low intra-model uncertainty.
Lastly, we assess the uncertainty in simulating the intra-annual variability of H S for both PC20 and FC21 due to different STs. Figure 8 displays the intra-annual variability of H S in the 13-sub regions for PC20-E, FC21-E, and each of the ensemble members. In general, the extratropical regions exhibit a larger intra-annual variability whereas the tropical ocean basins exhibit smaller intra-annual variability. In the tropical North Indian Ocean (TNIO) basin, the H S variation is very low throughout the year except for the monsoon season whereas TWSP shows almost no variability throughout the year. The intraannual variability for future climate show similar variability as that of the present climate, however a noticeable change from the present climate (PC20-E) is observed in the tropical North Indian Ocean (TNIO) basin during the southwest monsoon period ( Figure 8b). As mentioned before, the TNIO does not show any variability except in the monsoon season. It shows the similar increasing trend in June and reaches the maximum by the end of June but only continue until the last week of July, and started to decline nearly three weeks ago as compared to PC20. It is observed that PC20-E, FC21-E, and the ensemble members show similar behavior in all the 13-sub regions, however, the minimum and maximum values of H S differ, particularly for ST2. For example, ETNA and ETNP show a clear annual cycle of H S with minimum (~1.2 m) in summer (July) and gradually increases to its maximum (~3.8 m) in winter (Jan) for the ensemble mean ( Figure 8a) and PC20-ST3, PC20-ST4, and PC20-ST6 show similar variability in these regions with a minimum H S of 1.3 m to a maximum of 4 m, while PC20-ST2 ranges only from 1 m to 3.3 m. Thus, even the intra-annual variability is similar in all the ensemble members, there is a difference in the minimum and maximum values of H S in each member. In general, PC20-ST2 shows the least values for both minimum and maximum H S values among all, whereas the other three STs show almost similar ranges (min and max) of H S . Climate 2022, 10, x FOR PEER REVIEW 18 of 25 Figure 6. Time series of the global and regional (13 sub regions) annual mean Hs for the present, historical, and future ensemble mean, spread, and the deviation from the ensemble mean for individual ensemble members (ST2, ST3, ST4, and ST6).

Summary and Conclusions
An assessment of the intra-model uncertainty in the WW3 wave model, arising from the use of different source term packages to simulate the present and future projected wave climates has been carried out in this paper. The wind-wave climate ensemble data

Summary and Conclusions
An assessment of the intra-model uncertainty in the WW3 wave model, arising from the use of different source term packages to simulate the present and future projected wave climates has been carried out in this paper. The wind-wave climate ensemble data set was generated using a CMIP6-derived single-forcing (EC-EARTH) and single-future scenario (SSP5-8.5). We have analyzed the skill of the wind speeds (U 10 ) and the resulting wave heights (H S ) relative to the ERA5 reanalysis, for which a good agreement was found. However, an overestimation of extreme wave height is visible in the Southern Oceans, the smaller number of samples used for the H 95 s calculations could lead to these larger biases. The M-Score suggests that ST4 and ST3 are better than ST6 and ST2 in a global perspective. Regionally ST3, ST4, and ST6 have the maximum M-Score in different ocean basins. In general, ST3 and ST4 perform better in tropics whereas ST6 performs better in the extra tropics. The anomalies associated to each of the STs, comparing to the ensemble mean, show that ST2 tends to underestimate the H S and H 95 s in global and regional ocean basins, whereas ST3 and ST4 overestimate both. ST6 shows, in general, the lowest anomalies among all the STs, with positive values in the extratropical areas and negative ones in the tropical. The global and regional ensemble spread suggest that the uncertainty in the ensemble is higher in the extratropical regions for both historical and future climate, mainly contributed by ST2. As mentioned before, ST2 has a known weakness in estimating proper swell dissipation and resulting in an underestimation of wave height. But in the tropical ocean, the larger uncertainty is contributed by both ST2 and ST3, dissipation source term in both ST2 and ST3 are sensitive to swell. It is also observed that the uncertainties due to different STs are increased for the H 95 s . The uncertainty in simulating the intra-annual variability of H S for both PC20 and FC21 due to different STs shows that the intra-annual variability exhibits similar trend in all the STs for PC20 and FC21, however there is a difference in the minimum and maximum values of H S in each ensemble members showing the intra-model uncertainty present in the WW3 model.
Previous studies on wind-wave climate ensemble never addressed the intra-model uncertainty in the wave model (WW3). Here in this study, we focused on the impact of using multiple source term parameterizations for wave climate projections, even within the same wave model, to build associated ensembles. This is an often overlooked detail that different model physics present in the wave model can lead to larger uncertainty in the wave climate ensemble. Our study shows that uncertainty due to ST2 in the present climate (PC20) is more than the order of climate change signal and therefore does not make any significant contribution to the future climate projections. Therefore, proper attention needs to be taken while building wave climate ensembles considering different mode physics (STs) in the state-of-the art wave model, WW3.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cli10040051/s1. Author Contributions: R.K. wrote the manuscript with the support from G.L., A.S., R.K. and G.L. conceived the idea of this paper. A.S. supervised the whole project. G.L. and R.K. carried out the formal analysis. R.K., G.L., A.S. and F.A. contributed to writing-review and editing. All authors have read and agreed to the published version of the manuscript.