Our modelling framework consists of the Semi-implicit Cross-scale Hydroscience Integrated System Model (SCHISM, Zhang et al., 2016; Zhang et al., 2023) that is 2-way coupled with the 3^rd-generation spectral wave model Wind Wave Model (WWM-V, Roland et al., 2012; Sikirić et al., 2019).

SCHISM is an extension of the original SELFE code (Zhang and Baptista, 2008) that has undergone several enhancements and upgrades, including the inclusion of large-scale eddying and a seamless cross-scale capability from small bodies of water to the ocean. The model solves the Navier-Stokes equation, under hydrostatic and Boussinesq approximation, on unstructured grids, using a Finite-Element semi-implicit Galerkin scheme for horizontal motions, and a Finite-Volume scheme for vertical motions. The model accounts for the combined effects of wind and atmospheric pressure gradients. In this study, SCHISM was configured in its two-dimensional barotropic mode, which increases model efficiency and is a reasonable approximation for global applications with high computational and storage requirements (Fernández-Montblanc et al., 2019). In this work, SCHISM was used to simulate the generation and propagation of storm surges. The surface stress in SCHISM was calculated using a bulk formula with the drag coefficient computed estimated as a function of the wave motion. The bed shear stress was computed using the Manning approach, assuming a value of 0.02 for the Manning friction coefficient. A biharmonic viscosity filter with a diffusion number-like dimensionless constant of 0.025 is used to suppress grid-scale noises.

WWM-V is a 3^rd-generation state-of-the-art spectral wave model that solves the wave action equation on unstructured grids using the CRD-N-Scheme (e.g. Ricchiuto et al., 2005) for the advection part and discreitzing the spectral part using simple 1^st order implicit finite difference schemes, rending the whole approach monotone and convergent. For the integration, the source terms have been linearized (e.g. Patankar, 1980; Patankar, 1981) and if the linearization was not possible, we have underrelaxed the source terms as done for the triad-wave-wave interactions source term. For this study we employed the fully implicit RD-schemes with no splitting, using a Block-Gauss & Jacobi iterative solvers (e.g. Abdolali et al., 2020) on unstructured grids. The spectral domain used in this study consists of 36 frequencies and 24 directions. The frequencies are separated by a constant factor, from 0.04 Hz, corresponding to a period of 25 s, to 1 Hz, corresponding to a period of 1 second. The growth/dissipation of the waves was parameterized using (Ardhuin et al., 2010), setting the parameter BETAMAX to the value suggested for ECMWF data, while the non-linear interactions are reproduced by means of the Discrete Interaction Approximation (DIA) (Hasselmann and Hasselmann, 1985). With a resolution of 2 km along coasts, shallow-water interactions can play a significant role. In particular, the bottom friction was parameterized using the JONSWAP approximation (Hasselmann et al., 1973). Depth-limited breaking was represented using the approach by (Battjes and Janssen, 1978). In particular, the Unresolved Obstacles Source Term (UOST, Mentaschi et al., 2015b; Mentaschi et al., 2018; Mentaschi et al., 2019; Mentaschi et al., 2020) was introduced to parameterize the energy dampening of small unresolved islands. Furthermore, as ice plays an important role on wave dynamics at high latitudes, an approach based on UOST was implemented to dampen the energy as a function of the ice concentration. No treatment of ice was adopted in the circulation model. Finally, the methodology proposed by (Tracy et al., 2006) for spectral decomposition was implemented in the model. This approach allows for saving the mean parameter for the separate peaks of the wave spectrum, providing summarized information on the spectra at each node.

The models are 2-way coupled, in that the hydrodynamic component provides the wave model with current velocities and sea levels, while the sea wind stress is estimated by WWM-V within the growth-dissipation source term ST4 (Ardhuin et al., 2010) and passed back to SCHISM. In this parameterization the sea surface roughness is a function of the wave age according to (Janssen, 1989). The beneficial effect of such coupling on both waves and circulation was shown in several past studies, in particular with regards to sea levels and extreme conditions (Brown and Wolf, 2009; Bertin et al., 2012; Bertin et al., 2015; Staneva et al., 2016; Clementi et al., 2017).

The hindcast simulation was carried out on an unstructured mesh of 651,652 nodes covering the global ocean with a resolution ranging from about 50 km offshore to about 2-4 km nearshore globally. The mesh was generated using the OceanMesh2D tool (Roberts et al., 2019). The bathymetry was taken from the General Bathymetric Chart of the Oceans (GEBCO) dataset (Weatherall et al., 2015). The model was forced with hourly data of sea level pressure, wind speed and sea ice concentration from the fifth-generation reanalysis ERA5 (Hersbach et al., 2020). Apart from its resolution of 25 km, high for a global reanalysis, we selected ERA5 for its extended time span (from 1950 on), and its accuracy, compared to other products, in reproducing long-term trends (Erikson et al., 2022). The model time step was 400 s for both the circulation component and the wave component, with circulation-wave coupling occurring at each time step. The time span of the simulation is from 1950 to 2022. Model output was saved at each domain node with a time frequency of 3 hours for the variables listed in Table 1 .

SSH and Hs exhibit different statistical properties that require different treatments of the signal in the validation. Specifically, H_s is always positive, allowing for an easy-to-read interpretation of indicators such as Bias and Root Mean Squared Error as fractions of the mean Hs, by means of a normalization. On the other hand, SSH can be negative and the mean value of SSH may be close to 0. Therefore, the normalization often used for H_s is not meaningful for SSH. Additionally, our model is unable to replicate long-term trends in SSH observed in the data, such as Sea Level Rise (SLR). As a result, it was necessary to remove these trends from both the observations and the model output prior to the validation, and an indicator such as the Bias, which for H_s is a proxy of systematic error, loses its meaning for SSH. In the reminder of this section, the different signal treatments and error indicators used for SSH and Hs are explained.

Modelled Sea Surface Height (SSH) was compared globally versus along-track altimetry data from the Copernicus Marine Service (CMS) (Pujol et al., 2022), which comprises several altimeter missions (Topex-Poseidon, Jason-1, OSTM/Jason-2, Sentinel-3A, ERS-1, ERS-2, Envisat, Geosat Follow On, Cryosat, SARAL/AltiKa, HY-2A) from 1993 ongoing. The dataset includes all the atmospheric forcing contributions to sea level variation, included as a linear component to the overall Sea Level Anomaly (SLA). The model output was space-time interpolated to the corresponding points along the tracks. The paired observation and modelled values were then binned on a validation grid with a spatial resolution of 0.5°x0.5° and the skill was evaluated for each cell.

In addition to the SSH validation at the open ocean, the performance of the nearshore SSH was assessed versus the tidal gauge data provided by the GESLA-2 database (Woodworth et al., 2016; Haigh et al., 2022). The astronomical tide was removed from the station data by means of tidal analysis on an annual basis, using the classical harmonic analysis tool developed by (Pawlowicz et al., 2002).

As our barotropic model is unable to reproduce low-frequency variability of sea levels and Sea Level Rise, both simulations and observations of SSH were detrended prior to comparison.

The statistical indicators used to quantify the skill both of the offshore and nearshore SSH are:

* Root Mean Squared Error (RMSE_SSH):

where N is the number of observation-simulation pairs, and

η i m and η i o are the modelled and observed SSH, respectively.

• Normalized Root Mean Squared Error, NRMSE_SSH, similar to RMSE but normalized by the difference between the 1^st and 99^th percentile of the observed water levels:

This normalization allows an easy-to-read expression of the error as a percentage of the observation range.

• Pearson correlation:

where ηm¯ and ηo¯ are the mean modelled and observed SSH, and σm and σo are the standard deviations of the 2 signals.

To assess the behavior of the model in normal and extreme water level conditions, the skill indicators were evaluated a) on the whole time series and b) on the data beyond the 95^th percentile of the model or observation of each cell/tidal gauge.

The modelled significant wave height H_s was validated with the along-track altimetry H_s data provided by CMS, in collaboration with the Climate Change Initiative (CCI) of the European Space Agency (ESA) (Dodet et al., 2020). Similar to SSH, the model data were space-time interpolated to the points along the satellite tracks, and binned on the same validation grid. H_s was also compared with buoy data from a dataset of wave buoys from different sources, including the National Buoy Data Center (NDBC) and the European Marine Observation and Data Network (EMODnet).

The statistical indicators used to quantify the skill of significant wave height are:

• Normalized Bias (NBI):

where

H s   i m and H s   i o are the modelled and observed significant wave height, respectively. This indicator is a proxy of the model’s systematic error expressed as a percentage of the mean significant wave height (e.g. Mentaschi et al., 2013).

* The Root Mean Squared Error (RMSE_Hs) with a definition similar to that of SSH:

• Normalized Root Mean Squared Error (NRMSE_Hs):

Similar to NRMSE_SSH, the normalization in eq. 6 allows an easy-to-read expression of the error as a percentage of the observational range.

Similar to SSH, the behavior of the model was assessed both in normal and extreme conditions of H_s.

Overall, our results show offshore skills in line with the ones of other contemporary large-scale models, for both storm surges (Muis et al., 2016; Vousdoukas et al., 2017; Fernández-Montblanc et al., 2020; Tadesse et al., 2020) and waves (e.g. Perez et al., 2017; Mentaschi et al., 2020; Smith et al., 2021; Morim et al., 2022).

While RMSE_SSH appears generally low both at high and low latitudes ( Figure 1A ), the offshore skill of SSH expressed in terms of normalized indicators such as NRMSE_SSH and correlation shows a clear dependence on latitude ( Figures 2 , 3 ), with values of NRMSE_SSH and correlation significantly better at mid-high latitudes with respect to tropical and equatorial areas. This is a common feature in storm surge models and is due to a range of factors. The variability of sea levels at mid-high latitudes is higher and dominated by the eddy activity and baroclinic instability in the atmosphere, which result in short-term inverse barometric effect and wind setup that our model can capture. The lower atmospheric eddy activity in tropical/equatorial regions is a reason for the poorer skills for SSH in these areas (Hoskins and Valdes, 1990). At low latitudes the short-term variability of SSH is low with respect to middle-high latitudes, and storm surges are typically caused by tropical cyclones, hurricanes, and typhoons, which are relatively rare events, with a limited contribution to the mean signal. Though ERA5 can reproduce such events to some extent, limitations still exist, partially related to the still relatively low spatial resolution (Hersbach et al., 2020), with results usually characterized by negative biases of wind. Another reason is the complex ocean currents and coastal geography in the equatorial region, which is characterized by weak Coriolis acceleration. The equatorial region is home to the Equatorial Undercurrent (EUC) and the Equatorial Counter Current (ECC), which are difficult to model due to their strong variability and complexity (e.g. Stellema et al., 2022). These currents can also cause upwelling and vertical circulation, which can have an impact on SSH, but the ability of barotropic models to reproduce it is limited. Related to this is the low-frequency variability that characterizes the sea levels in tropical and equatorial areas, due to the Tropical Instability Waves and the seasonal cycle, and the associated thermostatic expansion. Another source of uncertainty is the oscillating intensities of the monsoons and of the Walker cells, and the associated patterns of variability such as the El Nino Southern Oscillation (Nidheesh et al., 2013). Other areas characterized by low-frequency variability are the western boundary currents, where the water levels depend on the intensity of the currents. The better skill in relatively sheltered tropical/equatorial basins such as the Red Sea, the Persian Gulf, the East China Sea, part of the South China Sea, the Gulf of Thailand and the Arafura Sea, is likely due to the relatively smaller low-frequency variability. In this respect, a 3D circulation model could be able to better capture the low-frequency component of SLA, resulting in better overall skill at low latitudes.

That low-frequency variability is a major contribution to the error is confirmed by a comparison of our model’s SSH against the data of Dynamic Atmospheric Correction (DAC) from Pujol et al. (2022), which better identifies the short-term variability than SLA ( Figure 10 ). The comparison, done for the mean signal, shows a significantly improved model skill, with a mean correlation of 0.75 versus 0.55 against SLA ( Figure 1B ), and better performances everywhere, with the exception of the Red Sea and the coasts of Vietnam.

Wave results show overall good skill, though with a tendentially negative bias, especially low in areas characterized by seasonal presence of ice, e.g. in the southern portion of the Arctic Ocean, in the Labrador Sea, and in the Southern hemisphere in the Weddell and Ross seas, and in semi-enclosed basins, with consequent higher values of NRMSE_Hs ( Figure 5 ). The error in H_s in these areas is related to the high uncertainty in the sea ice concentration datum provided by ERA5, especially in the Marginal Ice Zones (Renfrew et al., 2021), together with the fact that the presence of ice invalidates the altimetric observations of H_s (Dodet et al., 2020). On the other hand, introducing a parameterization of ice-induced wave dampening comes with benefits, as simulations neglecting it result in significant positive biases of H_s in the Southern Ocean and in the Northern Pacific Ocean and Atlantic Ocean that are now absent. The negative bias and NRMSE_Hs increase beyond the 95^th percentile ( Figure 6 ), reflecting a general tendency of our wave models in underestimating the peaks, partly due to an underestimation of extreme winds in ERA5 (e.g. Campos et al., 2022) as well as possible limitations of 3^rd generation models in extreme conditions (Cavaleri, 2009). At the same time, in enclosed and semi enclosed basins (the Mediterranean Sea, the Red Sea, the Persian Gulf, the inner seas of Indonesia and others) the negative bias of H_s beyond the 95^th percentile improves with respect to the whole signal. This suggests that the higher resolution of ERA5 compared to previous reanalyses and that of our model allow for a better reproduction of mesoscale events relevant for the statistics of the extremes in short fetches (e.g. Mentaschi et al., 2015a).

The advantage of the high resolution used in this study becomes clear when examining the model skill vs tidal gauges and buoys located close to the coasts. Here our model generally outperforms lower resolution ones, in particular the ones from previous studies by the authors (Mentaschi et al., 2017; Vousdoukas et al., 2017; Vousdoukas et al., 2018). At global scale, the performance of SSH versus tidal gauges is remarkably in line with the one versus altimeters. The global mean NRMSE_SSH/correlation is 17.09%/0.58 versus altimeters, and 14.49%/0.54 versus tidal gauges ( Figures 1 , 3 ). This picture is consistent for SSH beyond the 95^th percentile, with a global mean NRMSE_SSH/correlation of 27.41%/0.74 versus altimeters, and 20.07%/0.73 versus tidal gauges ( Figures 2 , 4 ).

For significant wave height the skill versus buoys is significantly worse than that versus altimeters, with a more pronounced negative bias, especially beyond the 95^th percentile. When considering the whole signal, the global mean NRMSE_Hs/NBI is 16.58%/-3.73% versus altimeters, and 25.79%/-11.8% versus buoys ( Figures 5 , 7 ). Beyond the 95^th percentile the figures are 15.94%/-7.85% versus altimeters, and 26.67%/-16.26% versus buoys ( Figures 6 , 8 ). The NRMSE_Hs versus buoys can be decomposed into its random and systematic components according to (Mentaschi et al., 2013), revealing that beyond the 95^th percentile the bias account for roughly 20% of the mean squared deviation. The larger bias of H_s nearshore is, beside the underresolved coastal zone, likely due to multiple factors, such as limitations in reproducing growth-dissipation in short fetches, and uncertainties in the parameterizations of shallow-water dynamics.

The significant improvement over lower resolution, non-coupled models for both SSH and H_s relates to the ability to better reproduce the dynamic effect of small coastal features, and to considering the interactions between waves and currents, and the effect of waves on the sea surface roughness. A comparison of the new hindcast versus previous ones (Mentaschi et al., 2017; Vousdoukas et al., 2018, Figure 11 ) shows how the present model better captures local short-term variability. For the 20% of tidal gauges with the best RMSE_SSH in the present study, the old hindcast showed nowhere a better skill ( Figure 12A ). Furthermore, RMSE_SSH/NRMSE_SSH was larger than the 80^th percentile of RMSE_SSH/NRMSE_SSH of the present study in >50% of the tidal gauges. If SSH beyond the 95^th percentile is considered, this percentage exceeds 60% ( Figures 12A, B ). The increase in skill for SSH >95th percentile can be observed especially at middle-high latitudes, while at low latitudes the improvement is less pronounced ( Figure 12E ). For H_s, in normal conditions the improvement of skill is in line with that of SSH, while if we consider H_s beyond the 95^th percentile the improvement of skill with respect to the old hindcast is still clear, but less marked, due to the significant negative bias of H_s ( Figures 12C, D ). For Hs >95th percentile the improvement in skill can be observed especially for buoys in close proximity to the shore, while it is less strong for offshore buoys ( Figure 12F ).

It is worth mentioning that the setup of the physics in WW3 and D-FLOW in (Mentaschi et al., 2017; Vousdoukas et al., 2018) is similar to that of the present study, and the offshore skill comparable. Therefore, the nearshore improvement of skill is due to the increase in resolution and to the coupling.

Our model represents a clear improvement on previous lower resolution hindcasts. However, it comes with a number of limitations that will be listed in this section. Tidal forcing was not included in the setup, as we aimed to isolate the contribution of the storm surges. While tide elevation can be added a-posteriori, from tidal databases such as FES2014 (Lyard et al., 2021) or TPXO (Dushaw et al., 1997), recent studies pointed to the importance on non-linear interactions between tides and storm surges nearshore (Krien et al., 2017; Fernández-Montblanc et al., 2019; Arns et al., 2020; Tausía et al., 2023). Such non-linear interactions are not explicitly accounted for in this study.

As a resolution of 2-4 km is not enough to resolve exhaustively everywhere coastal processes such as nearshore wave dissipation, we opted not to introduce the nearshore wave-circulation coupling described by (Bennis et al., 2011) and implemented in SCHISM by (Martins et al., 2022). Therefore, the water levels at the coast do not include the wave setup. This is a common approach for large-scale models, but it should be mentioned because it is important for certain applications.

The ability of our model to capture tropical cyclones (TCs) is still limited by the offshore grid resolution (~ 50 km) and by the ability of the forcing fields from ERA5 to reproduce TCs. This is another known limitation of large-scale models (e.g. Bloemendaal et al., 2019) and requires separate special modelling to be addressed.

As mentioned in the methods section, ERA5 forcing was selected for its advantages over other products (high spatial resolution for a global reanalysis, long time extent, ability to detect long-term trends). However, ERA5 comes with significant underestimation of extreme winds (e.g. Campos et al., 2022), which reflects in negative bias especially of extreme H_s.

Last but not least, this study was conducted using a 2D barotropic model for storm surges, which comes with a limited ability to reproduce 3D ocean circulation and currents. Though an accurate representation of 3D circulation is not an objective of this study, this could have an impact on storm surges, e.g., in relation to phenomena such as the Coriolis setup. Furthermore, this can have consequences on waves, related to the current-induced Doppler effect.