Oke, P. R., Brassington, G. B., Griffin, D. A., and Schiller, A. (2008). The Bluelink ocean data assimilation system (BODAS). Ocean Modelling, 21(1-2), 46-70, https://doi.org/10.1016/j.ocemod.2007.11.002

Ollinaho, P., Lock, S., Leutbecher, M., Bechtold, P., Beljaars, A., Bozzo, A., Forbes, R.M., Haiden, T., Hogan, R.J., Sandu, I. (2017). Towards process-level representation of model uncertainties: stochastically perturbed parametrizations in the ECMWF ensemble. Quarterly Journal of the Royal Meteorological Society, 143, 408-422, http://dx.doi.org/10.1002/qj.2931

Palmer, T. (2018). The ECMWF ensemble prediction system: Looking back (more than) 25 years and projecting forward 25 years. Quarterly Journal of the Royal Meteorological Society, 145, 12-24, https://doi.org/10.1002/qj.3383

Penduff, T., Barnier, B., Terray, L., Sérazin, G., Gregorio, S., Brankart, J.-M., Moine, M.-P., Molines, J.-M., Brasseur, P. (2014). Ensembles of eddying ocean simulations for climate. In: CLIVAR Exchanges, 65(19), 26-29. Available at: https://www.clivar.org/sites/default/files/documents/exchanges65_0.pdf

Pham, D. T., Verron, J., and Roubaud, M.C. (1998). A singular evolutive Kalman filters for data assimilation in oceanography. Journal of Marine Systems, 16(3-4), 323-340, https://doi.org/10.1016/S0924-7963(97)00109-7

Pinardi, N., Lermusiaux, P.F.J., Brink, K.H., Preller, R. H. (2017). The Sea: The science of ocean predictions. Journal of Marine Research, 75(3), 101-102, https://doi.org/10.1357/002224017821836833

Quattrocchi, G., De Mey, P., Ayoub, N., Vervatis, V., Testut, C.-E., Reffray, G., Chanut, J., Drillet, Y., (2014). Characterisation of errors of a regional model of the bay of biscay in response to wind uncertainties: a first step toward a data assimilation system suitable for coastal sea domains. Journal of Operational Oceanography, 7(2), 25-34, https://doi.org/10.1080/1755876X.2014.11020156

Ren, S., Zhu, X., Drevillon, M., Wang, H., Zhang, Y., Zu, Z., Li, A. (2021). Detection of SST Fronts from a High-Resolution Model and Its Preliminary Results in the South China Sea. Journal of Atmospheric and Oceanic Technology, 38(2), 387-403, https://doi.org/10.1175/JTECH-D-20-0118.1

Ryan, A. G., Regnier, C., Divakaran, P., Spindler, T., Mehra, A., Smith, G. C., et al. (2015). GODAE OceanView Class 4 forecast verification framework: global ocean inter-comparison. Journal of Operational Oceanography, 8(sup1), S112-S126, https://doi.org/10.1080/1755876X.2015.1022330

Sakov, P., Counillon, F., Bertino, L., Lisæter, K.A., Oke, P.R., Korablev, A., (2012). TOPAZ4: an ocean-sea ice data assimilation system for the north atlantic and arctic. Ocean Science, 8 (4), 633-656, http://dx.doi.org/10.5194/os-8-633-2012

Sandery, P.A., and Sakov, P. (2017). Ocean forecasting of mesoscale features can deteriorate by increasing model resolution towards the submesoscale. Nature Communications, 8, 1566, http://dx.doi.org/10.1038/s41467-017-01595-0

Santana-Falcón, Y., Brasseur, P., Brankart, J.M., and Garnier, F. (2020). Assimilation of chlorophyll data into a stochastic ensemble simulation for the North Atlantic Ocean. Ocean Science, 16, 1297-1315, https://doi.org/10.5194/os-16-1297-2020

Sasaki, Y. (1970). Some basic formalisms in numerical variational analysis. Monthly Weather Review, 98(12), 875-883.

Sein, D.V., Koldunov, N.K., Danilov, S., Wang, Q., Sidorenko, D., Fast, I., Rackow, T., Cabos, W., Jung, T. (2017). Ocean Modelling on a Mesh With Resolution Following the Local Rossby Radius. Journal of Advances in Modelling Earth Systems, 9:7, 2601-2614, https://doi.org/10.1002/2017MS001099

<figure><img src="image_233.png" /></figure>
Simon, E., and Bertino, L. (2009). Application of the Gaussian anamorphosis to assimilation in a 3-D coupled physical-ecosystem model of the North Atlantic with the EnKF: a twin experiment. *Ocean Science*, 5, 495-510, 2009, https://doi.org/10.5194/os-5-495-2009

Smith, G.M., Roy, F., Reszka, M., Surcel Colan, D., He, Z., Deacu, D., Belanger, J.-M., Skachko, S., Liu, Y., Dupont, F., Lemieux, J.-F., Beaudoin, C., Tranchant, B., Drévillon, M., Garric, G., Testut, G.-E., Lellouche, J.-M., Pellerin, P., Ritchie, H., Lu, Y., Davidson, F., Buehner, M., Caya, A., Lajoie, M. (2016). Sea ice forecast verification in the Canadian Global Ice Ocean Prediction System. Quarterly Journal of the Royal Meteorological Society, 142(695), 659-671, https://doi.org/10.1002/qj.2555

Storto, A., and Andriopoulos, P., (2021). A new stochastic ocean physics package and its application to hybrid-covariance data assimilation. Quarterly Journal of the Royal Meteorological Society, 147, 1691-1725, https://doi.org/10.1002/qj.3990

Tchonang, B.C., Benkiran, M., Le Traon P.-Y., Van Gennip, S.J., Lellocuhe, J.M., Ruggiero, G. (2021). Assessing the impact of the assimilation of SWOT observations in a global high-resolution analysis and forecasting system. Part 2: Results. Frontiers in Marine Science, 8:687414, https://doi.org/10.3389/fmars.2021.687414

Thacker, W.C., Srinivasan, A., Iskandarani, M., Knio, O.M., Le Hénaff, M., (2012). Propagating boundary uncertainties using polynomial expansions. Ocean Modelling, 43-44, 52-63, http://dx.doi.org/10.1016/j.ocemod.2011.11.011

Thoppil, P.G., Frolov, S., Rowley, C.D. et al. (2021). Ensemble forecasting greatly expands the prediction horizon for ocean mesoscale variability. Communications Earth & Environment, 2, 89, https://doi.org/10.1038/s43247-021-00151-5

Toublanc, F., Ayoub, N.K., Lyard, F., Marsaleix, P., Allain, D.J., 2018. Tidal downscaling from the open ocean to the coast: a new approach applied to the Bay of Biscay. Ocean Modelling, 214, 16-32. http://dx.doi.org/10.1016/j.ocemod.2018.02.001

Usui, N., Ishizaki, S., Fujii, Y., Tsujino, H., Yasuda, T., Kamachi, M. (2006). Meteorological Research Institute multivariate ocean variational estimation (MOVE) system: Some early results. Advances in Space Research, 37(4), 806-822, https://doi.org/10.1016/j.asr.2005.09.022

Vervatis, V. D., Testut, C.E., De Mey, P., Ayoub, N., Chanut, J., Quattrocchi, G. (2016). Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing. Ocean Modelling, 100, 1-19, https://doi.org/10.1016/j.ocemod.2016.01.003

Vervatis, V.D., De Mey-Frémaux, P., Ayoub, N., Karagiorgos, J., Ciavatta, S., Brewin, R., Sofianos, S., (2021a). Assessment of a regional physical-biogeochemical stochastic ocean model. Part 2: empirical consistency. Ocean Modelling, 160, 101770, http://dx.doi.org/10.1016/j.ocemod.2021.101770

Vervatis, V. D., De Mey-Frémaux, P., Ayoub, N., Karagiorgos, J., Ghantous, M., Kailas, M., Testut, C.-E., and Sofianos, S., (2021b). Assessment of a regional physical-biogeochemical stochastic ocean model. Part 1: ensemble generation. Ocean Modelling, 160, 101781, https://doi.org/10.1016/j.ocemod.2021.101781

Waters, J., Lea, D.L., Martin, M.J., Mirouze, I., Weaver, A., While, J. (2014). Implementing a variational data assimilation system in an operational 1/4 degree global ocean model. Quarterly Journal of the Royal Meteorological Society, 141(687), 333-349, https://doi.org/10.1002/qj.2388

Zaron, E.D. (2011). Introduction to Ocean Data Assimilation. In: Schiller, A., Brassington, G. (eds) "Operational Oceanography in the 21st Century". Springer, Dordrecht. https://doi.org/10.1007/978-94-007-0332-2_13

<figure><img src="image_234.png" /></figure>
# 6. Sea Ice modelling

<figure><img src="image_235.png" /></figure>

CHAPTER COORDINATOR

Laurent Bertino

CHAPTER AUTHORS (in alphabetical order)

Ed Blockley, Johnny A. Johannessen, and Einar Órn Ólason

<figure><img src="image_236.png" /></figure>

01

02

03

04

05

06

07

08

09

10

11

12
# 6. Sea Ice modelling

<figure><img src="image_237.png" /></figure>

## 6.1. General introduction to sea ice models

6.1.1. Objective, applications and beneficiaries

6.1.2. Fundamental theoretical background

## 6.2. Sea Ice forecast and multi-year systems

6.2.1. Architecture singularities

6.2.2. Input data: available sources and data handling

6.2.3. Modelling component

6.2.3.1. Basic equations and modelling choices

6.2.3.2. Sea ice rheology

6.2.3.3. Community sea-ice models

6.2.3.4. Coupling of sea ice to atmosphere and ocean

6.2.3.5. Model setup

6.2.4. Ensemble Modelling

6.2.5. Data Assimilation systems

6.2.5.1. Ensemble-based methods

6.2.5.2. Variational methods

6.2.5.3. Challenges with coupled data assimilation

6.2.6. Validation strategies

6.2.7. Outputs

6.2.8. Examples of operational sea ice forecasting systems

## 6.3. References

<figure><img src="image_238.png" /></figure>
## 6.1.
# General introduction to sea ice models

### 6.1.1. Objective, applications and beneficiaries

The main objective of an operational sea ice forecasting system is to provide users with a reliable estimate of the state of the ice cover and its temporal evolution. To meet this purpose, the system needs to be coupled to, or use data from, ocean and atmosphere forecasting systems. Some form of data assimilation is also required to counteract errors due to the chaotic nature of the atmosphere-ocean-ice system. Users of sea ice forecasting systems are either stakeholders operating in the Arctic or downstream service providers who use the information as an input to their own services. With a changing climate and a warming Arctic, the number of stakeholders interested in operating in that region is growing.

The Arctic is getting warmer with temperatures rising at approximately twice the rate of the global average (Overland et al., 2016) but also more attractive for business as its natural resources are becoming available for exploitation and transport for the first time in our history. These include about 13% of the world's oil and gas resources as estimated by the United States Geological Survey (Gautier et al., 2009), gold and other metals, and 5.5% of the freshwater resources stored on Greenland (Kundzewicz et al., 2007). Changing environmental conditions are modifying ecosystems in diverse ways. In the Barents Sea, the cod are thriving thanks to warming conditions (Kjesbu et al., 2014). A migration behaviour of boreal

generalist fishes to cooler waters is also observed in the Bering Sea (Mueter and Litzow, 2008). These changes have implications for fisheries management and more generally for the Arctic ecosystem. Cruise tourism in the Arctic is also developing fast since operators can offer comfortable icebreaker cruises all the way to the North Pole.

The NSR along the Russian coast of the Arctic, which was heavily used by the Soviet Union until the 1990's, could again become an attractive alternative to reach East Asia from Western Europe. The route is indeed shorter than the one crossing Suez Passage (17000 km instead of 22000 km for a Rotterdam-Shanghai voyage) and would save fuel. However, in case of accidents, cargo and fuel would pose serious threats for the Arctic environment. Coastguards and navies of the Arctic nations must then be prepared for assisting vessels, performing search and rescue operations, and remediating oil spills in ice-infested waters, with frequently poor communication capabilities that may hinder access to new information.

The oil and gas exploration and production need sea ice forecasting both on local scales, to simulate individual ice floes on the theatre of their operations, and on large scales, to predict the time of the freeze up and break-up of the ice. It is expected that the exploration and production activities will be more active in relatively mild ice conditions than in severe ice conditions, which means that forecasts will have higher value

<figure><img src="image_239.png" /></figure>

<figure><img src="image_240.png" /><figcaption>**Figure 6.1.** Pack ice showing a pressure ridge on the left; Marginal Ice Zone with ice floes on the right. (Photos: E. Storheim, INTAROS/NERSC).</figcaption></figure>

<figure><img src="image_241.png" /></figure>
for the MIZ than for the ice pack. The MIZ, defined as the ice-covered region under the influence of surface waves from the open ocean, is particularly in need of forecasts to prevent risks such as ice floe's projections under the action of waves.

There are fewer stakeholder interests in the Southern Ocean, due to the reduced commercial activities in that region. However, ice-ocean predictions can provide information for tourism or scientific operations in the region, including access to Antarctic research stations and support for scientific research vessels. The complex rescue of a joint tourist-research vessel stuck within the Antarctic sea ice in December 2013 (A. Luck-Baker, BBC News, 21 January 2014, 1), requiring assistance from two icebreakers and a helicopter, highlighted the need for reliable predictions even in such a remote region. On longer timescales, changing sea ice conditions have implications for ice-dependent wildlife in the region, such as emperor penguins (e.g., Jenouvrier et al., 2012), which raises associated wildlife management concerns.

The shipping industry is primarily concerned with very detailed ice concentration, thickness and compression (and marginally snow depths, because deep snow can also impede the

progression of an icebreaker). On the other hand, in the aftermath of oil spills in ice-infested waters, search and rescue operations and forecasting are both dependent on ice motion and their diffusive properties that increase the search radius with time. The question of spatial and temporal resolution is especially critical for the latter case because of the strong scale-dependence of sea ice deformation rates (Rampal et al., 2008). In addition, the diffusion is higher in the chaotic MIZ than in the ice pack (Figure 6.1). The oil industry would ultimately need a detailed forecast of the position of each ice floe surrounding their operations for the day-to-day management of their activities, which can be only delivered by discrete-element models (Herman 2015, Rabatel et al., 2015). How to nest discrete-element models into the continuum sea ice models, considered in this chapter, remains an open question.

### 6.1.2. Fundamental theoretical background

The physical processes simulated by sea ice models are commonly split into two: vertical processes, related to thermodynamic growth and melt, and mechanical and dynamical processes giving rise to horizontal movement of ice (Figure 6.2).

<figure><img src="image_242.png" /><figcaption>**Figure 6.2.** A CICE Consortium graphic of sea-ice physics illustrates the complexity and breadth of variables at play (From 2).</figcaption></figure>

1. https://www.bbc.com/news/science-environment-25833307

2. https://www.lanl.gov/discover/science-briefs/2021/March/0322-cice.php

<figure><img src="image_243.png" /></figure>
The thermodynamic growth and melt of ice can be thought of as the result of the diffusion of heat between ocean and atmosphere, through the ice. Additional complications arise primarily due to the presence of salt or brine pockets in the ice, and the presence of snow. The brine pockets affect the heat conductivity and heat capacity of the ice, while both heat conductivity and heat capacity of the snow, as well as its density, are affected by the state and type of snow, as well as snow metamorphosis.

The basics of thermodynamic modelling of sea ice have been well established since the early 70s (Maykut and Untersteiner, 1971), with the notable improvement in theoretical understanding brought by the application of mushy-layer theory to sea ice (Feltham et al., 2006), and substantial work relating to the dynamics of brine drainage and the multi-phase nature of sea ice (Vancoppenolle et al., 2007; Notz and Worster, 2009; Griewank and Notz, 2013). In terms of model development though, progress has been made in improving numerical performance and in technical aspects, such as conservation of heat, energy, and enthalpy (e.g., Semtner, 1976; Bitz and Libscomb, 1999; Winton, 2000; Huwald et al., 2005). Recently, the more advanced multi-phase physics have also found its way into large-scale sea-ice models (Turner et al., 2013; Turner and Hunke, 2015).

The fundamentals of ice dynamics modelling are less firmly rooted in basic theoretical understanding. While most of the terms of the momentum equation are well understood and follow the basic formulation of the Navier-Stokes equation on a rotating sphere, the formulation of internal stresses is less certain. These describe the response of the ice to external forcing and are, as such, at the heart of sea ice dynamical modelling.

Sea ice is a solid material and, as such, can only move once fractured or broken. In most sea ice models this is taken into

account by assuming a rate-independent (von Mises) plasticity. This approach was originally proposed by Coon et al. (1974) but reshaped into a more computationally tractable form in the viscous-plastic model proposed by Hibler (1979), in which the ice is assumed to deform in a linear-viscous manner until it reaches a plastic threshold, representing the fracturing or breaking of the ice. The fracturing process is, as such, simulated explicitly at the grid scale.

However, the process of ice fracturing has been shown to be the result of the propagation of fracturing events from small spatial scales to large ones (Weiss and Marsan, 2004). This results in fractal characteristics of the deformation rates (e.g., Marsan et al., 2004; Rampal et al., 2008; Stern and Lindsay, 2009, Schulson and Hibler, 2017). It means that a sea ice model hoping to correctly capture the deformation of the ice must account for this propagation of fracturing events from small to large scales. As the propagation starts at very small spatial and temporal scales (Oikkonen et al., 2017), a geophysical scale model must account for this through a sub-grid scale parameterisation.

The role and importance of fracture dynamics is still a hotly debated subject within the sea ice modelling community. The fractal nature of sea ice deformation is generally accepted and the scaling of deformation rates is recognised as a potential tool and metric for model evaluation and improvement (Rampal et al., 2016; Spreen et al., 2017; Hutter et al., 2018; Rampal et al., 2019; Bouchat et al., 2021). At the same time, it is still unresolved the question of whether to explicitly simulate the fracturing process at a very high resolution (Hutter et al., 2019) or to use a sub-grid scale parameterisation of the fracturing process at a more modest resolution (Dansereau et al., 2016; Rampal et al., 2016).

## 6.2.

# Sea Ice forecast and multi-year systems

### 6.2.1. Architecture singularities

This section and the next one focus on the “forward integration” spot in the centre of Figure 4.1, designing the architecture of an OOFS.

Sea ice drift forecasts are affected by multiple sources of uncertainties. The surface winds are one of the most important external forces driving the motion of the sea ice in the central Arctic (Thorndike and Colony, 1982). More-

over, the uncertainties in the atmospheric reanalysis in the Arctic are higher than those at the mid-latitudes, and observations are insufficient to estimate the statistical characteristics (scale, amplitudes) of the errors. Rabatel et al. (2018) investigated the sensitivity of sea ice drift using neXtSIM-EB for the uncertainties of the surface winds. They concluded that, in regions of highly compact ice cover, the accuracy of surface wind forcing and sea ice rheology are both important for the probabilistic forecast skill of sea ice trajectories.

<figure><img src="image_244.png" /></figure>
<table><caption>**Table 6.1.** Overview of operating and approved satellites and sensors for the sea ice observations grouped into: ESA and Eumetsat missions (yellow), 3rd Party Missions (green) and new approved missions from ESA and NASA/CNES (blue). Spectrometers and infrared radiometers are only sensing in cloud free conditions. Note (*) that CFOSAT flies a combined altimeter and real aperture radar at five distinct incidence angles up to 10 degrees. Harmony (**) comprises two bi-static satellites that will fly in convoy with Sentinel-1.</caption><thead><tr><th>Passive-μwaves</th><th>Scatterometer</th><th>SAR</th><th>Altimeter</th><th>Spectrometer</th><th>Infrared Radiometer</th></tr></thead><tbody><tr><td rowspan="2">SMOS</td><td rowspan="2">Metop-B/C ASCAT</td><td rowspan="2">Sentinel-1A/B</td><td>CRYOSAT-2</td><td>Sentinel-3 A/B</td><td>Sentinel-3 A/B</td></tr><tr><td>Sentinel-3 A/B</td><td>Sentinel-2 A/B</td><td>Metop AVHRR</td></tr><tr><td>AMSR-2</td><td>CFOSAT</td><td>Radarsat2</td><td>Altika</td><td>Aqua MODIS</td><td>Aqua MODIS</td></tr><tr><td rowspan="2">SMAP</td><td rowspan="2">Oceansat2</td><td>Radarsat constellation</td><td>ICESat/ICESat 2</td><td></td><td></td></tr><tr><td>CFOSAT*</td><td>CFOSAT*</td><td></td><td></td></tr><tr><td rowspan="4">CIMR</td><td rowspan="4">Metop second generation</td><td>Sentinel-1 C/D</td><td>SWOT</td><td>Sentinel-3 C/D</td><td>Sentinel-3 C/D</td></tr><tr><td>Rose-L</td><td>Cristal</td><td>Sentinel-2 C/D</td><td></td></tr><tr><td>HARMONY**</td><td></td><td></td><td></td></tr><tr><td>BIOMASS</td><td></td><td></td><td></td></tr></tbody></table>

The ocean below the ice contains large quantities of heat and momentum, enough to melt the sea ice and to cause ice drift and deformations. Uncertainties in ocean temperature, vertical mixing, and currents are then very meaningful for the sea ice. The surface ocean salinity is important, as the melting point temperature depends on it. However, measuring ocean properties and particularly currents below the sea ice is challenging and uncertainties are rather high.

Uncertain initial conditions, particularly the sea ice thickness, persist a long time (Chevallier and Salas-Mélia, 2012). Blockley and Peterson (2018) showed that the sea ice conditions in spring persist typically a few months into the summer and are an important source of large-scale predictability. Errors in the position of the ice edge at the beginning of a forecast are usually persistent throughout the forecast run and ought to be post-processed for practical use.

Finally, sea ice models are dependent on their numerous model parameters, both in the sea ice dynamics and thermodynamics (Urrego-Blanco et al. 2016).

### 6.2.2. Input data: available sources and data handling

Initialized forecasts are critically dependent on the observations used for their initialization. To be useful for operational systems, observations are needed in near real-time for short-term forecasts and with limited time lag for seasonal and longer forecasts. There are unique challenges involved in polar observations because of its remoteness, harsh conditions, and long polar night. However, forecasting systems are making use of satellite observations for initialization, most routinely for sea ice concentration. Additionally, new products, such as sea ice thickness and drift, are becoming available and may ultimately improve the predictive capabilities.

Sea ice reconnaissance flights were mostly occasional until after the second world war, with the exception of the USSR which started systematic flights with Polar Aviation as early as 1929 to monitor the Northern Sea Route. The USA and Japan gradually increased the frequency of their flights at the turn of the 1950's and adopted the WMO sea ice charting standard proposed in 1952 (WMO, 1970). These flights are still used nowadays, mostly in Canada, but have elsewhere been superseded by satellite data.

<figure><img src="image_245.png" /></figure>
Regular and routine sea ice observations are today performed by a variety of satellite-based measurements provided by several space agencies (as grouped in the matrix in Table 6.1 and organised by satellite sensor classes). This has been accomplished thanks to a large number of major technical and scientific milestones and achievements over more than 40 years, as further addressed below. Note that spectrometers and radiometers are only sensing in cloud free conditions.

Started in 1978, the longest satellite record to cover the whole Arctic comes from polar orbiting passive microwave sensors onboard the satellites SMMR, SSM/I, AMSR-E and AMSR2 (Cavalieri and Parkinson, 2012) which provide the sea ice areal concentration. Their main advantage is that they can see through clouds but still a few issues remain, especially with the summer ice, because the sensor does not properly discriminate between open water and signatures from wet snow and melt ponds. This and other technical issues are accommodated differently in a multitude of algorithms that calculate sea ice concentrations from the raw passive microwave retrievals (Ivanova et al., 2014, 2015). This is an important matter for data assimilation as we will see in Section 6.2.5. The resolution of passive microwaves depends on the frequency band used, with the most precise low-frequency channels having the largest footprint (as large as 60 km). However, gridded sea ice concentration data can be found at resolutions between 6 and

25 km, which is consistent with current operational models of the whole Arctic but still coarse with respect to the needs of any operational users navigating in ice-infested waters. SAR and satellite data in the visible channels (VIIRS, AVHRR, MODIS, SPOT) provide much more detail at spatial resolutions finer than 1 km, which is what the ship captain would need, for example to detect and sail along a lead. However, both types of data suffer from poor coverage, SAR images because the acquisition frequency may be limited, and visible data because they are impaired by the frequent cloud coverage and by Arctic winter darkness.

For short-term forecasts, it is important to assess how the sea ice is moving. Various sea ice drift products are obtained from different satellites and can be split in two types: 1) the coarse resolution, full spatial coverage products using passive microwave radiometers and scatterometers (most accurate retrievals in winter because of the aforementioned limitations of passive microwave data during summer; see review by Sumata et al., 2014); and 2) the high resolution but reduced coverage SAR-based products (Kwok 2006). The SAR coverage has recently significantly improved by the launch of the ESA Sentinel-1 A/B missions offering full daily coverage in high latitudes (Korosov and Rampal. 2017). In comparison, drifting buoys on sea ice still provide the longest (more than 40 years) data record of the IABP but with limited spatial coverage.

<table><caption>**Table 6.2.** Overview of sea ice variables observed per group of satellite sensors listed in Table 6.1. Note (*) that spectrometers and radiometers are only sensing in cloud free conditions.</caption><thead><tr><th>Satellite Sensors</th><th>Extent and concentration</th><th>Sea Ice type</th><th>Sea Ice thickness</th><th>Snow depth</th><th>Sea Ice drift</th><th>Open leads</th><th>Melt ponds</th><th>Waves in ice</th></tr></thead><tbody><tr><td>Passive microwaves</td><td>X</td><td>First and multi-year ice</td><td>Less than 50 cm thick from L-band</td><td></td><td>X</td><td></td><td></td><td></td></tr><tr><td>Scatterometer</td><td>(X)</td><td>(X)</td><td></td><td></td><td>X</td><td></td><td></td><td></td></tr><tr><td>SAR</td><td>(X)</td><td>(X)</td><td></td><td></td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>Altimeter (radar, laser)</td><td></td><td></td><td>Thicker than ~ 50 cm</td><td>X</td><td></td><td>X</td><td></td><td></td></tr><tr><td>Spectrometer *</td><td>(X)</td><td>Ridges, rough and flat ice</td><td></td><td></td><td>(X)</td><td>X</td><td>X</td><td>(X)</td></tr><tr><td>Infrared Radiometer *</td><td>(X)</td><td></td><td></td><td></td><td></td><td>X</td><td>(X)</td><td></td></tr></tbody></table>

<figure><img src="image_246.png" /></figure>
Sea ice thickness observations from satellites have recently become routinely available. These use different principles to obtain either: sea ice freeboard of thick ice, for example from the satellite altimetry missions ICESat, ICESat2 (Kwok et al., 2007), CryoSat-2 (Laxon et al., 2013) and Sentinel-3; and the thickness of thin ice, derived from the SMOS (Tian-Kunze et al., 2014). These observations are quite complex and come with relatively high uncertainties (Zygmuntowska et al., 2014, Tian-Kunze et al., 2014). As discussed above, sea ice thickness is an important source of sea ice predictability on seasonal and longer timescales. Other aspects of the sea ice, such as snow cover, snow thickness and melt pond characteristics, may also be important for sea ice forecasts on seasonal and longer timescales. Remote sensing is used to characterise these aspects of the sea ice. For example, snow depth information is being provided through the NASA Operation Ice-Bridge airborne campaign (Kurtz et al., 2013), and melt pond fractions have been derived from satellite data in the visible channels (Rösel et al., 2012). Combined use of IceSAT-2, CryoSat-2 and Altika has also demonstrated promising capabilities to recover reliable snow depth estimates during winter (Guerreiro et al., 2016). Ice mass buoys are also providing in-situ measurements of snow depth and other sea ice characteristics (Richter-Menge et al., 2006; Perovich et al., 2008). However, only limited work has been done to quantify the possible influence of these types of observations for forecasting systems.

As already indicated in Table 6.1, the continuity of sea ice observations from satellites are indeed assured by the approved future satellite missions such as CIMR, Cristal, and ROSE-L Coopericus Sentinel Expansion missions under preparation by ESA. However, more dedicated field campaigns are still needed to assess the uncertainties of the satellite-based retrievals, as well as to harvest the multi-sensor synergies as can be noted

<figure><img src="image_247.png" /><figcaption>**Figure 6.3.** Illustration of a vertical temperature profile in a column consisting of sea ice of thickness hi and topped by snow of depth hs. The heat conduction equation is discretized in the vertical by Δz-thick levels. (adapted from Lisæter 2009).</figcaption></figure>

from Table 6.2. In turn, this would improve the quality and use of satellite data, and expectedly advance the forecast skill of sea ice on seasonal to interannual timescales.

### 6.2.3. Modelling component

#### 6.2.3.1. Basic equations and modelling choices

Most modern large-scale sea ice models are based on very similar foundations. The ice is generally modelled as a continuum using a Eulerian perspective, with the sea ice moving in a horizontal plane, subject to both external and internal forces. The dynamic evolution of the sea ice cover is described using two continuity equations and the momentum equation. The thermodynamic evolution is modelled within each column of the grid and is modelled as a heat diffusion process within the slab of sea ice. There are substantially varying degrees of complexity in the treatment of the thermodynamic processes, ranging from treating all the ice as being of a single thickness (Hibler, 1979) and treating the heat diffusion without resolving the temperature profile (Semtner, 1976), to using multiple thickness categories (Hibler, 1979, and numerous later variations) and treating the heat diffusion using mushy-layer dynamics (Feltham, et al., 2006).

The main equations for a simple dynamic model of sea ice with two categories (ice and open water) are the two continuity equations and the momentum equation. The continuity equation for mass is:

$$ \frac{\partial m}{\partial t} + \nabla \cdot (vm) = S_m \quad (6.1) $$

where m is the sea ice mass per unit area, $S_m$ is a thermodynamic source/sink term and v is velocity. In the case of a single sea ice category the continuity equation for the sea ice distribution takes the basic form:

$$ \frac{\partial A}{\partial t} + \nabla \cdot (vA) = S_A \quad (6.2) $$

with A is the sea ice concentration and $S_A$ is a source/sink term. In addition, the condition $A \le 1$ is imposed. This can be interpreted as a ridging condition since m can increase even if A does not (Hibler, 1979). Together these equations describe the advection of the sea ice in a given velocity field.

The momentum equation is generally written as (Connolley et al., 2004):

$$ m \frac{\partial v}{\partial t} = A(\tau_a + \tau_w) - m f \hat{k} \times v - mg \nabla H - \nabla \cdot \sigma \quad (6.3) $$

Here $\hat{k}$ is a unit vector normal to the surface, $\tau_a$ and $\tau_w$ are the air and water stresses, f is the Coriolis parameter, g is the gravitational acceleration, $\nabla H$ is the gradient of the sea surface height and $\sigma$ is the sea ice stress tensor. The acceleration term on the left-hand side may be set to zero, depending on

<figure><img src="image_248.png" /></figure>
<figure><img src="image_249.png" /><figcaption>**Figure 6.4.** Left: a schematic example of an ice thickness distribution, the thickness classes are in x-axis and the y-axis relates to the area concentrations. The continuous distribution is shown with a solid line, while the discretized version is shown in filled bars. Right: illustration of the subgrid-scale ice thickness distribution in a sea ice model (only two classes of 3 and 4 m thickness for the sake of illustration) (adapted from Lisæter, 2009).</figcaption></figure>

the model implementation. The last term on the right hand side $\nabla \cdot \sigma$, describes forces due to internal stress while the other terms are all external factors. Wind and water stresses are generally treated as quadratic drag (e.g., McPhee, 1975). In the absence of internal stress, the sea ice is in “free drift” and the model simplifies drastically. Free drift forecasts have therefore been used for a long time (Grumbine 1998) are still used operationally.

The thermodynamic equation is the heat diffusion equation:

$$ \rho c \frac{\partial T}{\partial t} = k \frac{\partial^2 T}{\partial t^2} \qquad (6.4) $$

where $\rho c$ is the heat capacity of sea ice or snow and $k$ is the heat conductivity. This equation can be solved in various ways (see Figure 6.3) (e.g., Maykut and Untersteiner, 1971; Semtner, 1976; Bitz and Libscomb, 1999; Winton, 2000; Huwald et al., 2005), discretized in the vertical (Figure 6.3). These take into account different physical properties and numerical solutions in solving the equation.

In addition to these grid-scale quantities, many models consider various sub-grid scale information and parameterisations. The most important of those is arguably the ice thickness distribution (ITD). This assumes that each grid cell of

the model contains not is of uniform thickness, but of varying thicknesses described by an ice thickness distribution $g$ (see Figure 6.3). This is in principle a continuous distribution of thicknesses, which is modified through dynamic and thermodynamic processes. The governing equation of evolution of the ice thickness distribution is (e.g., Thorndike et al., 1975):

$$ \frac{\partial g}{\partial t} = -\nabla \cdot (vg) - \frac{\partial}{\partial h}(fg) + \Psi \qquad (6.5) $$

where $f$ is the thermodynamic growth or melt rate, $h$ is the ice thickness, and $\Psi$ is a mechanical redistribution function.

In practice, sea ice models must use a discretized version of the ice thickness distribution, resulting in models with several distinct thickness categories (Figure 6.4 right, for a top view of a grid cell). The thickness distribution then becomes (Bitz et al., 2001):

$$ g(h) = g_0 \delta(h) + \sum_{i=1}^{M} g_i \delta(h - H_i) \qquad (6.6) $$

where $M$ is the number of thickness categories, $M_i$ is the thickness of category $i$, and $\delta(h)$ is the Dirac delta function. Various implementations exist, but the one from Bitz et al., (2001) with five thickness categories remains a popular choice.

<figure><img src="image_250.png" /></figure>
<figure><img src="image_251.png" /><figcaption>**Figure 6.5.** Two yield curves commonly associated with sea ice rheology: in red the elliptic yield curve used in the (E)VP models, and in blue a Mohr-Coulomb yield curve, for instance used in the brittle rheology of neXtSIM.</figcaption></figure>

In addition to these two basic components, a large number of sub-grid scale processes can and should be represented, depending on the use cases for each model. These include simulation of melt points (Flocco et al., 2010; Hunke et al., 2013), changes in atmospheric and oceanic drag due to sea ice roughness (Tsamados et al., 2014), and salt rejection from freezing sea ice (Vancoppenolle et al. 2009).

### 6.2.3.2. Sea ice rheology

The relationship between the internal stress and resulting deformation is referred to as rheology. Basically, all continuum, geophysical-scale sea ice models currently employ the VP rheology proposed by Hibler (1979) or some direct descendant of that work. The VP rheology treats the sea ice as a continuum and assumes it deforms in a viscous manner with a high viscosity until the internal stress reaches a plastic threshold, determined by a yield curve which usually has an elliptic shape (see Figure 6.5). Several important improvements have been made to the original VP rheology (e.g. Hunke and Dukowicz, 1997; Lemieux et al., 2010; Bouillon et al., 2013; Kimmritz et al., 2017), but the physical principles remain the same.

The VP rheology has enjoyed tremendous success and is used for time scales from days to centuries and spatial scales from tens of kilometres to basin scales. However, the VP rheology is not without faults when it comes to both the underlying assumptions (see in particular Coon et al., 2007) and the results it produces. In model inter-comparison studies, there is generally a very large spread - well beyond observed variability - in key prognostic variables such as sea ice thickness, concentration, and drift (Chevallier et al., 2017; Tandon et al., 2018). The sharp gradients in velocities, which are known as LKFs that are related to ridge and lead formation, are also poorly reproduced in any VP-based model running at a coarser resolution than about 2 km - a resolution that is an order of magnitude higher than the observational data (Sreenet al., 2017; Hutter et al., 2019).

Therefore, several authors, such as Tremblay and Mysak (1997), Wilchinsky and Feltham (2004), Schreyer et al. (2006), Girard-Ardhuin and Ezraty (2012), Dansereau et al. (2016), and Ólason et al. (2022), have suggested alternative approaches to the VP rheology. The EAP rheology of Wilchinsky and Feltham (2004) was implemented in the CICE model (Tsamados et al., 2013) and has been used in several studies, although it was not widely adopted yet (Bouchat et al., 2021; Hutter et al. 2021). The brittle rheologies of Girard et al. (2011), Dansereau et al. (2016), and Ólason et al. (2022) have all been implemented in the neXtSIM model (Bouillon and Rampal 2015; Rampal et al., 2019; Ólason et al., 2022) and used for forecasting and scientific research by the team involved in the model. The current neXtSIM version uses the BBM rheology of Ólason et al. (2022).

### 6.2.3.3. Community sea-ice models

Practically all sea ice models used in modern forecasting platforms are based on the principles described above. They use a Eulerian reference frame and use some version of the VP or the Elastic-Viscous-Plastic (EVP) rheologies. The thermodynamic growth and melt of ice are described through the diffusion of heat between ocean and atmosphere, through the ice. As such, they all follow the same general design philosophy. The main differences exist in the form of different choices of parameterisation and differences in data assimilation approaches.

Today, the CICE model (e.g. Hunke et al., 2021) is likely the most widely used sea ice model for operational forecasts. This model was developed at the Los Alamos National Laboratory and was originally designed to be part of the CCSM. Thanks to its clean and modular design, the model has been used in other multiple modelling systems, as a stand-alone model, part of sea ice-ocean models, and part of climate and earth-system models. The LIM and SI3 models (Rousset et al., 2015), which are part of the NEMO modelling system, are also very widely used, but only within the NEMO modelling system. Other sea-ice models include the SIS (Adcroft et al., 2019), which is part of the GFDL ocean modelling system MOM, the MITgcm sea-ice model (Losch et al., 2010) and the FESOM sea-ice model (FESIM, Danilov et al., 2015).

In contrast, only a few stand-alone sea ice models have used moving Lagrangian coordinates (Hopkins 2004), among which the neXtSIM model (Rampal et al., 2016; Rampal et al., 2019; Ólason et al., 2022) and the DEMSI model (Turner et al., 2022). neXtSIM-F is unique among forecasting models as it uses both a moving Lagrangian mesh and a newly developed brittle rheology, the Brittle Bingham-Maxwell (Williams et al., 2021; Ólason et al., 2022). This setup gives results that are clearly different from the classical systems, and arguably more realistic (Rampal et al., 2016; Rampal et al., 2019; Ólason et al., 2020; Ólason et al., 2022). The key improvement is a

<figure><img src="image_252.png" /></figure>
much more realistic representation of the deformation sta-
tistics of the ice cover, which gives more realistic leads and
ridges in the model. Sea ice drift simulated by neXtSIM is also
very realistic, and the pan-Arctic ice-thickness distribution is
also quite good (Williams et al., 2021).

6.2.3.4. Coupling of sea ice to atmosphere and ocean

Sea ice models are integral parts of Earth system models. The reason for this is that at high latitudes sea ice insulates the relatively warm ocean from the cold atmosphere. Over an unbroken sea ice cover, the atmosphere can therefore cool much more than it could if it was not insulated by the presence of sea ice. This has an impact on all ocean-atmosphere interactions in the polar regions, and therefore a global climate or Earth system model without a sea ice model can simply not function.

Sea ice interacts with the atmosphere through heat, moisture, and momentum exchanges. In summer incoming short-wave radiation melts the ice surface but would warm up the ocean surface in the absence of sea ice. In winter, heat conduction from the ocean and through the ice only results in a very modest amount of heat flux to the atmosphere. However, the dominant heat flux source is radiant cooling through long wave radiation from the surface. This happens because surface cooling through long wave radiation is much more efficient than the heat conduction through ice from the ocean, resulting in a surface that is colder than the lowest layers of the atmosphere. The result is a predominant temperature inversion and a stable atmospheric boundary layer. This reduces even further the latent and sensible heat fluxes from the surface. However, openings in the ice (leads and polynyas) expose the relatively warm ocean surface to the atmospheric boundary layer, which causes mixing and breaks down the stable boundary layer.

Momentum transfer between ice and atmosphere happens through wind stress at the surface of the ice. This is the main driver of ice movement and exerts a drag on the atmosphere, slowing down the wind. The amount of momentum transferred between ice and atmosphere is determined primarily by the stability of the atmospheric boundary layer (Gryanik and Lüpkes, 2017) and the roughness of the ice. While parameterisations and studies on the ice surface roughness have been proposed (Lüpkes et al., 2012, Castellani et al., 2014), consistent and basin-scale observations of the atmospheric drag coefficient over sea ice are currently unavailable (Petty et al., 2017). In a modelling context, our ability to predict ice surface roughness is severely limited, as most ice-atmosphere coupled models do not take surface roughness into account when calculating atmosphere-ice momentum exchanges.

Sea ice interacts with the ocean through heat, fresh-water,
and salt exchanges, as well as momentum exchanges. During

summer, the mixed layer may warm up due to shortwave
heating through openings in the ice. This makes the ice melt-
ing from below, causing release of both freshwater and salt
into the ocean. In winter, the atmosphere extracts heat from
the ocean through the ice, causing new ice to form at the
bottom of the existing ice pack. This causes a net heat and
freshwater flux out of the ocean. However, most of the salt
present in the ocean cannot enter the ice, since the ice is
much fresher than the ocean (ca. 15 vs 30 PSS for newly
formed ice in the Arctic). This results in a layer of very salty
water forming below the ice, which then sinks into the mixed
layer. The resulting salt plumes generally reach the bottom of
the halocline but may also be mixed into the mixed layer in
the presence of turbulence.

Momentum transfer between ice and ocean happens through interface stress between ocean and ice. The momentum coupling of ice and ocean is much stronger than that of ice and atmosphere, and the ice can be considered as the first layer in the ocean's Ekman spiral. Ice-ocean stress drives most geostrophic flows in ice covered areas, as well as some larger scale circulation patterns.

It is also worth mentioning the interaction between sea ice and ocean waves. Waves entering the ice pack may mechanically fracture it into smaller sea ice floes. This can widen the MIZ, which may also be viewed as the area where the ice is fractured by waves. Fracturing the ice into smaller floes increases the mobility of the ice cover, the momentum transfer between atmosphere, ocean, and ice, and this may cause enhanced melting of the ice through lateral melt. The ice in turn dampens the waves causing an attenuation of the wave amplitude, so that waves will only penetrate a limited distance into the ice pack, depending on the size of the waves and the compactness of the pack. Wave-ice interactions are of major importance in the Southern Ocean, but less so in the Arctic, where much less of the ice edge is exposed to open ocean.

Virtually all climate or earth-system models today include sea-ice models of the classical description above, i.e., a Eule-rian reference frame, VP or EVP rheology, and thermodynam-ics and column physics of varying complexity. They gener-ally include very simplistic formulations for the momentum transfer between atmosphere, ocean, and ice, and no ice-wave interactions. This is true for all the CMIP6 models. In fact, the sea ice models used in today's forecasting models were all designed for climate modelling, the only current exception is the above mentioned neXtSIM model. It is not clear how this lineage of the models affects the quality of their short-term predictions. It could be argued that a good large-scale sea ice model should be able to represent scales from ca. 1 km up to the basin scales and from hours to centu-ries. This is not the current case, but the discussion of why it is this way and how to address it is still in its infancy (Hunke et al., 2020; Blockley et al., 2020).

<figure><img src="image_253.png" /></figure>
<figure><img src="image_254.png" /><figcaption>**Figure 6.6.** Left: example of ice trajectories from an ensemble of 10 members of 7-days sea ice drift of synthetic floats in an area of the Barents Sea from the TOPAZ system using randomly perturbed winds. The mean sea ice thickness is indicated above. Right: illustration of the ensemble spread in end point positions increasing as a function of the forecast length. The uncertainty growth depends strongly on the region (from Bertino et al., 2015).</figcaption></figure>

### 6.2.3.5. Model setup

In nearly all forecasting platforms, the sea ice model is coupled to an ocean model. There are platforms that use fully coupled atmosphere-ocean-sea ice models and only a few platforms using a stand-alone sea ice model. The reasons for this are partly historical: most sea ice models are written as parts of ocean models. Ocean forecasting and re-analysis platforms have tended to include a sea ice model from the start, making a dedicated sea ice forecasting platform redundant. In addition, the coupling between sea ice and ocean is quite strong, so running a separate sea ice forecasting platform can bring its own set of challenges. On the other hand, a stand-alone sea ice forecasting platform can be run at a higher resolution and can be used as a technology preview, as in the case of the neXtSIM-F platform.

### 6.2.4. Ensemble Modelling

Probabilistic forecasts, which are widely used in weather forecasting (Molteni et al., 1996; Leutbecher and Palmer, 2008), are still in their infancy in sea ice forecasting. Probabilistic predictions rely on an ensemble of model simulations (e.g. a Monte Carlo simulation) used to describe the forecast uncertainty stemming from errors in the model parameters, initial and boundary conditions, and any external forcing. The resulting range of model outputs is used to retrieve statistical information, such as the ensemble mean and its spread (i.e. the standard deviation), which are thus used instead of the deterministic forecast to estimate the associated uncertainty (see Figure 6.6). The multiple simultaneous sources of errors usually make the forecast accuracy of the ensemble mean ex-

ceed that of the single deterministic prediction (Leith, 1974), although the spread often underestimates the actual forecast error when the sources of error are not all adequately accounted for (Buizza et al., 2005). Monte Carlo techniques are already common practice in different areas (e.g. Dobney et al., 2000; Hackett et al., 2006; Breivik and Allen, 2008; Melsom et al., 2012; Motra et al., 2016; Duraisamy and Iaccarino, 2017) and a common tool for sensitivity analysis.

### 6.2.5. Data Assimilation systems

As introduced in the previous section, a sea ice forecast needs to regularly assimilate operational observations, which, at present, are mostly satellite data. The most tempting way forward is to insert directly the satellite observed concentrations and thicknesses into the model. However, this is not as easy as it sounds in a complex sea ice code where a large number of model variables are dependent on each other. Hence, various data assimilation methods are used for sea ice models, similar to those used for ocean physical, biogeochemical models or weather models. The most common method is nudging, which is less disruptive than direct insertion: the data are introduced gradually over a given time scale (Lindsay and Zhang, 2006). The nudged model is then assumed to adjust itself progressively using its own equations. But how much can we rely on such adjustments?

When the ocean mixed layer is too warm to sustain sea ice but observations show the presence of sea ice, a data assimilation system updating only sea ice would add sea ice on top of the warm waters, but the huge heat capacity of the ocean would then melt the added sea ice almost immediately. The

<figure><img src="image_255.png" /></figure>
ocean mixed layer temperature and salinity must be adjusted accordingly. This suggests that when used in a coupled ice-ocean system, assimilation of sea ice observations ought to be coupled in the sense that it should update both the sea ice and the ocean properties consistently. In data assimilation jargon, this means that the sea ice observation should be projected down to the ocean column using a multivariate forecast error covariance matrix.

### 6.2.5.1. Ensemble-based methods

Dynamical model ensembles are a practical way to estimate the covariances mentioned above. In data assimilation terminology, the state vector must include all prognostic variables of the coupled model (ocean and sea ice variables) and the ensemble of model runs can be used to calculate empirically the cross-covariances between sea ice and ocean variables. Similarly, observations of the ocean are used to update sea ice variables, although this situation is less common. Using an EnKF (see section 5.5.2), Lisæter et al. (2003) demonstrated that the coupled assimilation of sea ice properties can modify the ocean surface temperatures in rather systematic ways (adding sea ice cools down the water), but not ocean salinities. However, according to sea ice halodynamics, the freezing of sea ice injects salty brines to the ocean mixed layer and the melting releases fresher water, but these simple relationships explain only a part of the sea ice-salinity cross-covariances and a relationship may arise in other situations without the intervention of sea ice thermodynamics: the wind may occasionally blow the sea ice on top of more saline water. Sakov et al. (2012) showed how the sea ice-salinity cross-covariance can change sign on either side of the ice edge in the Barents Sea: the sea ice-salinity correlation turns negative on the ocean side because the main process responsible for melting is the advection of warm and saline Atlantic water near the surface, thus the sea ice-salinity correlation is made through the intermediate of the surface temperature variable. The last finding does not hold in locations where the sea ice is isolated from the Atlantic water, but such isolation may not remain forever if the open water mixing reaches these warm waters (Rippeth et al., 2015). The assimilation of sea ice concentrations with the EnKF described in Lisæter et al. (2003) was included in the near-real-time TOPAZ forecasts in 2003.

### 6.2.5.2. Variational methods

An alternative to ensemble methods is the use of an adjoint model as in the 4D-variational (4D-Var) data assimilation method. The adjoint model and the tangent linear model calculate the sensitivity of observed variables to the control variables within the duration of the assimilation window. If tangent linear and adjoint models are available both for the ocean and the sea ice models, they can exchange information about the interface variables, like the heat, salt, and mo-

mentum fluxes. Since these correlations are usually monovariate at the beginning of the assimilation window, the length of the assimilation window should be as long as possible. The most recent experiments report successful applications of the 4D4D-VarVAR in an Arctic regional configuration for durations of one year or longer (Fenty and Heimbach 2013; Fenty et al., 2015); an adjoint model for the EVP sea ice rheology has been introduced later (Toyoda et al. 2019). The advantage of the 4D-Var method is that it returns one optimised model trajectory, which is very useful for oceano-graphic interpretation (Kauker et al., 2009) and quantitative network design (Kaminski et al., 2015) but, to our knowledge, 4D-Var is not used for operational sea ice-ocean forecasting.

A computationally simplified variant of 4D-Var is known as 3D-Var, in which the same increment is used to compute the model equivalent of the observation-minus-reference state differences at all times in the assimilation. Owing to the relatively low cost of the scheme compared with the full 4D-Var, 3D-Var is commonly used by operational forecasting centres around the world (Usui et al., 2006; Mogensen et al., 2012; Hebert et al., 2015; Tonani et al., 2015; Waters et al., 2015; Smith et al., 2016, see Table 6.3 below).

### 6.2.5.3. Challenges with coupled data assimilation

Coupled multivariate covariances do not necessarily cure all the troubles of assimilating sea ice observations. Another source of problems is the lack of respect of the traditional Gauss-linear assumptions underlying classical data assimilation methods. By definition, sea ice concentrations have bounded values between zero and one, while other sea ice tracer variables (thickness, snow depth) have positive values only. Ocean temperatures are not allowed below the freezing point. While it should be easy for a monovariate assimilation method, based on a heuristic covariance function, to preserve monotonicity and therefore the bounds of variables (Wackernagel, 2003), an ensemble-based covariance (or a tangent linear model) may generate values out of bounds. Honouring the bounds can be forced by different means, either by nonlinear transformations of the variables (a method called Gaussian anamorphosis in geostatistics; Bertino et al., 2003, Barth et al., 2015) or by including inequality constraints in the cost function (Lauvernet et al., 2009; Simon et al., 2012; Janjic et al., 2014). Altogether, the benefits of multivariate flow-dependent covariances still outbeat the inconvenience of values out of bounds.

There are continuous improvements to data assimilation methods in chaotic high-dimensional systems, such as coupled sea ice-ocean models. But new models and new observations always call for further developments in data assimilation. In particular, sea ice models expressed in Lagrangian grids with automatic remeshing are uncommon targets for data assimilation. Ensemble Kalman Filtering techniques rely

<figure><img src="image_256.png" /></figure>
<figure><img src="image_257.png" /><figcaption>**Figure 6.7.** Left: example of areas of excess ice (A+) and missing ice (A-) for a given TOPAZ forecast in the European Arctic; the validation data is the Norwegian Ice Charts. Right: associated Fraction Skill Scores showing that the forecast beats persistence at +5 and +9 days lead time irrespective of spatial scales for that specific forecast (both from 3rd of Copernicus Marine Service).</figcaption></figure>

on cross-covariances between observed and unobserved variables, which implies that the grid cells have to be uniquely identified across different members of the ensemble. This also becomes difficult when adaptive remeshing is turned on, unless the Lagrangian model output is interpolated on a fixed grid, at the risk of smoothing the very localised kinematic features (long cracks, ridges and leads) that they are meant to simulate (Aydoğdu et al., 2019). Lagrangian models do not offer any easy differentiation/automatic adjoint capabilities, thus preventing the use of variational techniques. It should also be noted that a coupling framework such as CESM is sufficiently flexible to allow several instances of a model component to be run (for example, the atmosphere) for each instance of another (for example, the ocean), thus allowing to use different data assimilation methods for the sea ice, ocean, land, and atmosphere. An important aspect in view of coupled data assimilation and ensemble forecasting is that the uncertainties are consistent across these compartments; the error statistics at the base of the atmosphere are consistent with those at the surface of the sea ice and similarly between the bottom of the sea ice and the ocean surface. This is possible to enforce if all components of the coupled system use an ensemble to represent the errors.

### 6.2.6. Validation strategies

Since a measure of RMS errors of sea ice concentrations depend on arbitrary choices made by the person doing the scoring (these errors diminish as more open ocean is included in the validation area), more targeted sea ice validation metrics express the skill as distance of the forecast from the observed ice edge (Dukhovskoy et al., 2015). In the Arctic, the skill of the 24-hour forecast of ice edge location is about 50 km for the TOPAZ (including both seasonal biases and RMS errors, updated at http://cmems.met.no/ARC-MFC/) and 40 km down to 30 km depending on the input data sources in the ACNFS (Hebert et al. 2015; Posey et al., 2015), although both methods may differ and be sensitive to special configurations of the ice edge. The area of discrepancy is accepted as an objective metric with the IIEE introduced by Goessling et al. (2016). The dependence of such metrics on spatial scales can be further included in the evaluation using the FSS (Melsom et al., 2019) and an extension of the IIEE metric proposed by Goessling and Jung (2018) for the evaluation of ensemble forecasts of the ice edge. Examples of IIEE and FSS are shown in Figure 6.7.

It should be noted that the metrics related to isolines (like the ice edge, classically defined at the 15% ice concentration isoline, or at other critical values such as 50% and 85%) are

<figure><img src="image_258.png" /></figure>

<footnote n="3">3. https://cmems.met.no/ARC-MFC/</footnote>
<figure><img src="image_259.png" /><figcaption>**Figure 6.8.** Typical Example of a two-day sea ice drift from satellite observations (OSI-SAF, left) and a model (TOPAZ4, right) (courtesy of A. Melsom, MET Norway).</figcaption></figure>

play also to other isolines, like the frontier between FYI and MYI, or to the MIZ/pack boundary. Contingency tables are also a valuable approach to the validation of sea ice concentrations (Smith et al., 2016), as well as the threat scores or the Heidke Skill Score.

The forecast skills for sea ice drift have received comparatively less attention but errors in sea ice drift are important, both for their contribution to the displacement of the sea ice edge and for their cumulative contribution to the sea ice thickness distribution. Long climate simulations hint for a seasonal dependence of the forecast skills, also noted in free drift simulations (Grumbine 1998). Biases of sea ice drift have been revealed in IPCC simulations (Tandon et al., 2018) related to the seasonal cycle. To our knowledge, there are no signs that these shortcomings are corrected in recent forecast models (Xie et al., 2017, for the TOPAZ system; Hebert et al., 2015 for the ACNFS), although a review of global reanalysis systems shows that some models simulate correctly the minimum sea ice drift in March (Chevallier et al., 2017). Hebert et al. (2015) also noted that the forecast of drifter positions beats persistence although the forecast of drift speed does not, indicating that the drift direction is better forecasted than the drift speed. How to remedy these shortcomings? Although adjusting the mean speed of sea ice can be easily achieved by tuning the drag coefficients, there is no simple tuning that can make the sea ice accelerate over years or shift its seasonal cycle.

The assimilation of sea ice drift data has been so far less successful than that of sea ice concentrations: Stark et al. (2008) showed a 50% reduction of errors in ice speed but no

benefit to ice concentrations and Sakov et al. (2012) indicated a low sensitivity of the sea ice drift to external perturbations in the wind forcing, which points to a shortcoming of the TOPAZ4 version of the EVP sea ice rheology. Qualitatively, the large-scale patterns of sea ice drift can be reproduced by such a model (see a typical situation in Figure 6.8) but the observed gradients between areas of low sea ice drift (North of Greenland) and strong sea ice drift (North of the Barents Sea) are smoothed by the model, which tends to simulate intermediate values of the sea ice drift speed. The forecast of 24-hours ice trajectories and locations exhibits an RMS error of 6.3 km in TOPAZ4 (Melsom et al., 2015), which does not seem to beat a simple free drift predictor (5 km in Grumbine, 1998). It should be noted that the validation is done against different data sources (sea ice drift from satellite SAR images versus IABP buoys) and at different periods (years 2012-2015 versus 80's and 90's decades). The SIDFEx (Goessling et al., 2020), carried out in the framework of the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) ice camp, has been the first to collect forecasts from international systems and has shown that a consensus forecast could be successfully used to order detailed satellite images of the ice camp in advance. Beyond the use of RMS errors, several alternative metrics for sea ice drift validation have been reviewed by Grumbine (2013). Validation metrics for an ensemble of trajectories from a probabilistic ice drift forecast have been proposed by Rabatel et al. (2018) and refined in Cheng et al. (2020) based on an analogy with Search and Rescue operations, in which the ensemble of trajectories define a search ellipse; the success of the forecast is the probability of containment of the target inside the ellipse.

<figure><img src="image_260.png" /></figure>
The forecast of sea ice thickness also suffers from excessive smoothness: thick sea ice is too thin and thin sea ice is too thick (Johnson et al., 2012) and errors reach easily one to two metres. There is a dynamical contribution to these errors with the too high sea ice drift speed North of Greenland ex-aggerating the transport of thick sea ice into the Beaufort Gyre. However, thermodynamic contributions cannot be ex-cluded either (in particular from snow and melt ponds). More generally, any error in the model initial conditions, atmo-spheric and ocean boundary conditions or its inherent param-eters will eventually accumulate in sea ice thickness biases, which means that different errors can cancel each other and yield a correct sea ice thickness for the wrong reasons. It is worth stressing the important role of snow on sea ice as an effective insulator, its presence can inhibit both the growth and melt of sea ice and thus reduce its seasonal cycle. Snow predictions in sea ice-ocean models are very dependent on the quality of precipitation from weather analyses and fore-casts which are difficult to validate and usually vary from one product to another (Lindsay et al., 2014).

### 6.2.7. Outputs

Information on formats and types of outputs by all kinds of OOFS can be found in Chapter 4. In this Section, we list the variables related to sea ice forecasts:

* sea ice concentration (SIC)

* sea ice thickness (SIT)

* sea ice drift velocity in x- and y- directions (SIUV)

* snow depth (SNOW)

* sea ice age

* sea ice albedo (SIALB)

* sea ice temperature

Sea ice forecasting systems generally comply with CF standards. The CF metadata conventions are a widely used standard for atmospheric, ocean, and climate data. Standard names are de-fined in a CF Standard Name Table (see 4). Standard vari-able names from the CMIP nomenclature can be found in Notz et al. (2016).

### 6.2.8. Examples of operational sea ice forecasting systems

Most present day short-term forecast systems (listed in Table 6.3) assimilate sea ice concentration and are therefore expect-ed to perform well at forecasting the ice edge. These systems include the Canadian RIOPS (Smith et al., 2021), the United States ACNFS/GOFS3.1 (Hebert et al., 2015), the Italian GOFS16 (Iovino et al., 2016) the Global and the Arctic Marine Forecasting System (TOPAZ, Sakov et al. 2012) by the Coperni-cus Marine Services. Stand-alone sea ice models, like neXtSIM-F (Williams et al., 2021), are also used for forecasting purposes and, given that their control vector excludes the ocean, they can be initialised more flexibly than coupled ice-ocean sys-tems. Baltic forecasting systems are omitted for brevity. Ocean data assimilated are also omitted from Table 6.3.

<table><caption>**Table 6.3.** List of present-day short-term global and Arctic forecast systems. Note that the output spatial resolution may differ from the native resolution.</caption><thead><tr><td>Area</td><td>System</td><td>Resolution (km)</td><td>Model</td><td>Assimilation (method and sea ice data)</td><td>Variables</td><td>Website</td></tr></thead><tbody><tr><td>Arctic</td><td>ArclOPS</td><td>18 km</td><td>MITgcm</td><td>LESTKF SIC, SIT</td><td>SIC, SID, SIT</td><td>http://www.oceanguide.org.cn/IceIndexHome/ThicknessIce</td></tr><tr><td>Arctic</td><td>NOAA (Bob Grumbine)</td><td>N/A</td><td>Free drift</td><td>N/A</td><td>SIUV</td><td>https://mag.ncep.noaa.gov/model-guidance-model-parameter.php?group=Model%20Guidance&amp;model=ICE-DRIFT&amp;area=POLAR&amp;ps=area#</td></tr><tr><td>Global</td><td>RTOFS</td><td>3.5 km</td><td>HYCOM-CICE5</td><td>3DVAR SIC</td><td>SIC, SIT, SIUV</td><td>https://polar.ncep.noaa.gov/global/</td></tr></tbody></table>

4. http://cfconventions.org/standard-names.html

<figure><img src="image_261.png" /></figure>
<table><thead><tr><th>Area</th><th>System</th><th>Resolution (km)</th><th>Model</th><th>Assimilation (method and sea ice data)</th><th>Variables</th><th>Website</th></tr></thead><tbody><tr><td>Arctic</td><td>TOPAZ4</td><td>12.5 km</td><td>HYCOM-CICE3</td><td>EnKF SIC, SIUV, SIT</td><td>SIC, SIT, SIUV, SNOW, SIALB</td><td>https://marine.copernicus.eu/</td></tr><tr><td>Arctic</td><td>neXtSIM-F</td><td>7.5 km*</td><td>neXtSIM</td><td>Nudging SIC</td><td>SIC, SIT, SIUV, SNOW</td><td>https://marine.copernicus.eu/</td></tr><tr><td>Global</td><td>MOi</td><td>3.5 km</td><td>NEMO-LIM2</td><td>SEEK SIC</td><td>SIC, SIT, SIUV</td><td>https://marine.copernicus.eu/</td></tr><tr><td>Global</td><td>GIOPS</td><td>12 km</td><td>NEMO-CICE4</td><td>3DVAR SIC</td><td>N/A</td><td>CONCEPTS - Science.gc.ca</td></tr><tr><td>Arctic</td><td>RIOPS</td><td>3.5 km</td><td>NEMO-CICE4</td><td>3DVAR SIC</td><td>N/A</td><td>https://science.gc.ca/eic/site/063.nsf/eng/h_97620.html</td></tr><tr><td>Global</td><td>GOF3.1</td><td>3.5 km</td><td>HYCOM-CICE4</td><td>3DVAR SIC</td><td>SIC, SIT, SIUV</td><td>https://www7320.nrlssc.navy.mil/GLBhycomcice1-12</td></tr><tr><td>Global</td><td>ECMWF</td><td>12 km</td><td>NEMO-LIM2</td><td>3DVAR SIC</td><td>SIC, SIT</td><td>https://www.ecmwf.int/en/forecasts/datasets/set-i</td></tr><tr><td>Arctic</td><td>DMI</td><td>10 km</td><td>HYCOM-CICE4</td><td>Nudging SIC</td><td>N/A</td><td>http://ocean.dmi.dk/models/hycom.uk.php</td></tr><tr><td>Global</td><td>Met Office coupled DA</td><td>12 km</td><td>NEMO-CICE5</td><td>3DVAR SIC</td><td>SIC, SIT, SIUV</td><td>https://marine.copernicus.eu/</td></tr><tr><td>Global</td><td>Met Office FOAM</td><td>3.5 km</td><td>NEMO-CICE5</td><td>3DVAR SIC</td><td>N/A</td><td></td></tr><tr><td>Arctic**</td><td>VENUS</td><td>2.5km</td><td>IcePOM</td><td>N/A</td><td>SIC, SIT</td><td>https://ads.nipr.ac.jp/venus.mirai/#/mirai</td></tr></tbody></table>

<figure><img src="image_262.png" /></figure>

<footnote>* Note that the resolution of a Lagrangian triangular mesh is not comparable to square grids, thus the output resolution is 3 km.</footnote>

<footnote>** VENUS is deployed on demand.</footnote>
## 6.3. References

Adcroft, A., Anderson, W., Balaji, V., Blanton, C., Bushuk, M., Dufour, C. O., et al. (2019). The GFDL global ocean and sea ice model OM4.0: Model description and simulation features. *Journal of Advances in Modelling Earth Systems*, 11, 3167-3211, https://doi.org/10.1029/2019MS001726

Aydoğdu, A., Carrassi, A., Guider, C.T., Jones, C. K. R. T., and Rampal, P. (2019). Data assimilation using adaptive, non-conservative, moving mesh models. *Nonlinear Processes in Geophysics*, 26(3), 175-193, https://doi.org/10.5194/npg-26-175-2019

Barth, A., Canter, M., Van Schaeybroeck, B., Vannitsem, S., Massonnet, F., Zunz, V., Mathiot, P., Alvera-Azcárate, A., Beckers, J.-M. (2015). Assimilation of sea surface temperature, ice concentration and ice drift in a model of the Southern Ocean. *Ocean Modelling*, 93, 22-39, https://doi.org/10.1016/j.ocemod.2015.07.011

Bertino, L., Bergh, J., and Xie, J. (2015). Evaluation of uncertainties by ensemble simulation, ART JIP Deliverable 3.3. NERSC, Bergen, Norway.

Bertino, L., Evensen, G., Wackernagel, H. (2003). Sequential Data Assimilation Techniques in Oceanography. *International Statistical Review* 71(2), 223-241, http://www.jstor.org/stable/1403885

Bitz, C. M. and Lipscomb, W. H. (1999). An energy-conserving thermodynamic model of sea ice. *Journal of Geophysical Research: Oceans*, 104(C7), 15669-15678, https://doi.org/10.1029/1999JC900100

Bitz, C. M., Holland, M. M., Weaver, A. J., and Eby, M. (2001). Simulating the ice-thickness distribution in a coupled climate model. *Journal of Geophysical Research: Oceans*, 106(C2), 2441-2463, https://doi.org/10.1029/1999JC000113

Blockley, E. W. and Peterson, K. A. (2018). Improving Met Office seasonal predictions of Arctic sea ice using assimilation of CryoSat-2 thickness. *The Cryosphere*, 12, 3419-3438, https://doi.org/10.5194/tc-12-3419-2018

Blockley, E., Vancoppenolle, M., Hunke, E., Bitz, C., Feltham, D., Lemieux, J., Losch, M., Maisonnave, E., Notz, D., Rampal, P., Tietsche, S., Tremblay, B., Turner, A., Massonnet, F., Ólason, E., Roberts, A., Aksenov, Y., Fichefet, T., Garric, G., Iovino, D., Madec, G., Rousset, C., Salas y Melia, D., Schroeder, D. (2020). The future of sea ice modeling: Where do we go from here? *Bulletin of the American Meteorological Society*, 101, E1302-E1309, https://doi.org/10.1175/BAMS-D-20-0073.1

Bouchat, A., Hutter, N.C., Chanut, J., Dupont, F., Dukhovskoy, D. S., Garric, G., Lee, Y. J., Lemieux, J.-F., Lique, C., Losch, M., Maslowski, W., Myers, P. G., Ólason, E., Rampal, P., Rasmussen, T.A.S., Talandier, C., Tremblay, B. (2021). Sea Ice Rheology Experiment (SIREX), Part I: Scaling and statistical properties of sea-ice deformation fields. *Journal of Geophysical Research: Oceans*, 127(4), e2021JC017667, https://doi.org/10.1029/2021JC017667

Bouillon, S., Fichefet, T., Legat, V., and Madec, G. (2013). The elastic-viscous-plastic method revisited. *Ocean Modelling*, 71, 2-12, https://doi.org/10.1016/j.ocemod.2013.05.013

Bouillon, S., and Rampal, P. (2015). On producing sea ice deformation data sets from SAR-derived sea ice motion. *The Cryosphere*, 9, 663-673, https://doi.org/10.5194/tc-9-663-2015

Breivik, Ø. and Allen, A. A. (2008). An operational search and rescue model for the Norwegian Sea and the North Sea. *Journal of Marine Systems*, 69(1-2), 99-113, https://doi.org/10.1016/j.jmarsys.2007.02.010

<figure><img src="image_263.png" /></figure>
Buizza, R., Houtekamer, P. L., Pellerin, G., Toth, Z., Zhu, Y., and Wei, M. (2005). A Comparison of the ECMWF, MSC, and NCEP Global Ensemble Prediction Systems. *Monthly Weather Review*, 133(5), 1076-1097, https://doi.org/10.1175/MWR2905.1

Castellani, G., Lüpkes, C., Gerdes, R., Hendricks, S. (2014). Variability of Arctic sea ice topography and its impact on the atmospheric surface drag. *Journal of Geophysical Research: Ocean*, 119(10), 6743-6762, https://doi.org/10.1002/2013JC009712

Cavalieri, D. J., Parkinson, C. L. (2012). Arctic Sea ice variability and trends, 1979-2010. *The Cryosphere*, 6, 881-889, https://doi.org/10.5194/tc-6-881-2012

Cheng, S., Aydoğdu, A., Rampal, P., Currass, A., Bertino, L. (2020). Probabilistic Forecasts of Sea Ice Trajectories in the Arctic: Impact of Uncertainties in Surface Wind and Ice Cohesion. *Oceans*, 1(4), 326-342, https://doi.org/10.3390/oceans1040022

Chevallier, M., and Salas-Melia, D. (2012). The role of sea ice thickness distribution in the Arctic sea ice potential predictability: A diagnostic approach with a coupled GCM. *Journal of Climate*, 25(8), 3025-3038, https://doi.org/10.1175/JCLI-D-11-00209.1

Chevallier, M., Smith, G. C., Dupont, F., Lemieux, J.-F., Forget, G., Fujii, Y., et al. (2017). Intercomparison of the Arctic sea ice cover in global ocean-sea ice reanalyses from the ORA-IP project. *Climate Dynamics*, 49, 1107-1136, https://doi.org/10.1007/s00382-016-2985-y

Connolley, W. M., Gregory, J. M., Hunke, E., McLaren, A. J. (2004). On the Consistent Scaling of Terms in the Sea-Ice Dynamics Equation. *Journal of Physical Oceanography*, 34(7), 1776-1780, https://doi.org/10.1175/1520-0485(2004)034<1776:OTCSOT>2.0.CO;2

Coon, M. D., Maykut, G. A., Pritchard, R. S., Rothrock, D. A., Thorndike, A. S. (1974). Modeling the pack ice as an elastic-plastic material. *AIDEX Bulletin*, 24, 1-106, Univ. of Wash., Seattle, 1974.

Coon, M., Kwok, R., Levy, G., Pruis, M., Schreyer, H., Sulsky, D. (2007). Arctic Ice Dynamics Joint Experiment (AIDEX) assumptions revisited and found inadequate. *Journal of Geophysical Research: Ocean*, 112, C11S90, https://doi.org/10.1029/2005JC003393

Danilov, S., Wang, Q., Timmermann, R., Iakovlev, N., Sidorenko, D., Kimmritz, M., Jung, T., and Schröter, J. (2015). Finite-Element Sea Ice Model (FESIM), version 2. Geoscientific Model Development, 8, 1747-1761, https://doi.org/10.5194/gmd-8-1747-2015

Dansereau, V., Weiss, J., Saramito, P., and Lattes, P. (2016). A Maxwell elasto-brittle rheology for sea ice modelling. *The Cryosphere*, 10, 1339-1359, https://doi.org/10.5194/tc-10-1339-2016

Dobney, A., Klinkenberg, H., Souren, F., and Van Borm, W. (2000). Uncertainty calculations for amount of chemical substance measurements performed by means of isotope dilution mass spectrometry as part of the PERM project. *Analytica chimica acta*, 420, 89-94.

Dukhovskoy, D. S., Ubnoske, J., Blanchard-Wrigglesworth, E., Hiester, H. R., and Proshutinsky, A. (2015). Skill metrics for evaluation and comparison of sea ice models. *Journal of Geophysical Research: Oceans*, 120(9), 5910-5931. https://doi.org/10.1002/2015JC010989

Duraisamy, K. A., and Iaccarino, G. (2017). Assessing turbulence sensitivity using stochastic Monte Carlo analysis. *Physical Review Fluids*, https://arxiv.org/abs/1704.05187

Feltham, D. L., Untersteiner, N., Wettlaufer, J. S., and Worster, M. G. (2006). Sea ice is a mushy layer. *Geo-physical Research Letters*, 33(14), L14501, ISSN 0094-8276, https://doi.org/10.1029/2006GL026290

<figure><img src="image_264.png" /></figure>
Fenty, I., and Heimbach, P. (2013). Coupled sea ice-ocean state estimation in the Labrador Sea and Baffin Bay. *Journal of Physical Oceanography*, 43(5), 884-904, https://doi.org/10.1175/JPO-D-12-065.1

Fenty, I., Menemenlis, D., and Zhang, H. (2015). Global coupled sea ice-ocean state estimation. *Climate Dynamics*, 49, 931-956, https://doi.org/10.1007/s00382-015-2796-6

Flocco, D., Feltham, D., Turner, A. (2010). Incorporation of a physically based melt pond scheme into the sea ice component of a climate model. *Journal of Geophysical Research: Oceans*, 115, C08012, https://doi.org/10.1029/2009JC005568

Gautier, D. L., Bird, K. J., Charpentier, R. R., Grantz, A., Houseknecht, D. W., Klett, T. R., ... Wandrey, C. J. (2009). Assessment of Undiscovered Oil and Gas in the Arctic. *Science*, 324, 1175-1179, https://doi.org/10.1126/science.1169467

Girard, L., Bouillon, S., Weiss, J., Amitrano, D., Fichefet, T., and Legat, V. (2011). A new modeling framework for sea-ice mechanics based on elasto-brittle rheology. *Annals of Glaciology*, 52, 123-132, doi:10.3189/172756411795931499

Girard-Ardhuin, F., and Ezraty, R. (2012). Enhanced Arctic Sea Ice Drift Estimation Merging Radiometer and Scatterometer Data. *IEEE Transactions on Geoscience and Remote Sensing*, 50(7), 2639-2648, doi:10.1109/TGRS.2012.2184124

Goessling H., Tietsche, S., Day, J.J., Hawkins, E., Jung, T. (2016). Predictability of the Arctic sea ice edge. *Geophysical Research Letters*, 43(4), 1642-1650, https://doi.org/10.1002/2015GL067232

Goessling, H. F., and Jung, T. (2018). A probabilistic verification score for contours: Methodology and application to Arctic ice-edge forecast. *Quarterly Journal of the Royal Meteorological Society*, 144, 735-743, https://doi.org/10.1002/qj.3242

Goessling, H. F. and the SIDFEx Team (2020). Making Use and Sense of 75,000 Forecasts of the Sea Ice Drift Forecast Experiment (SIDFEx), EGU General Assembly 2020, Online, 4-8 May 2020, EGU2020-18867, https://doi.org/10.5194/egusphere-egu2020-18867

Griewank, P. J., Notz, D. (2013). Insights into brine dynamics and sea ice desalination from a 1-D model study of gravity drainage. *Journal of Geophysical Research: Oceans*, 118, 3370-3386, doi:10.1002/jgrc.20247

Grumbine, R. W. (1998). Virtual Floe Ice Drift Forecast Model Intercomparison. *Weather and Forecasting*, 13, 886-890, https://doi.org/10.1175/1520-0434(1998)013<0886:VFIDFM>2.0.CO;2

Grumbine, R. W. (2013). Long Range Sea Ice Drift Model Verification. Camp. Springs, Maryland.

Guerreiro, K., Fleury, S., Zakharova, E., Rémy, F., and Kouraev, A. (2016). Potential for estimation of snow depth on Arctic sea ice from CryoSat-2 and SARAL/AltiKa missions. *Remote Sensing of Environment*, 186, 339-349, https://doi.org/10.1016/j.rse.2016.07.013

Gryanik, V. M., and Lüpkes, C. (2017). An Efficient Non-iterative Bulk Parametrization of Surface Fluxes for Stable Atmospheric Conditions Over Polar Sea-Ice. *Boundary-Layer Meteorology*, 166, 301-325, https://doi.org/10.1007/s10546-017-0302-x

Hackett, B., Breivik, Ø., Wettre, C. (2006). Forecasting the drift of objects and substances in the ocean. In E.P. Chassignet and J. Verron (ed.s) "Ocean Weather Forecasting, An Integrated View of Oceanography", 507-524, Springer, Dordrecht.

<figure><img src="image_265.png" /></figure>
Hebert, D. A., Allard, R. A., Metzger, E. J., Posey, P. G., Preller, R. H., Wallcraft, A. J., Phelps, M. W., and Smedstad, O. M. (2015). Short-term sea ice forecasting: An assessment of ice concentration and ice drift forecasts using the U.S. Navy's Arctic Cap Nowcast/Forecast System. *Journal of Geophysical Research: Oceans*, 120, 8327-8345, https://doi.org/10.1002/2015JC011283

Herman, A. (2015). Discrete-Element bonded particle Sea Ice model DESIgn, version 1.3a - model description and implementation. *Geoscientific Model Development*, 9, 1219-1241, https://doi.org/10.5194/gmd-9-1219-2016

Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. *Journal of Physical Oceanography*, 9, 817-846.

Hopkins, M. (2004). A discrete element Lagrangian sea ice model. *Engineering Computations*, 21, 409-421, https://doi.org/10.1108/02644400410519857

Hunke, E., and Dukowicz, J. (1997). An elastic-viscous-plastic model for sea ice dynamics. *Journal of Physical Oceanography*, 27(9), 1849-1867, https://doi.org/10.1175/1520-0485(1997)027<1849:AEVPMF>2.0.CO;2

Hunke, E., Hebert, D., Lecomte, O. (2013). Level-ice melt ponds in the Los Alamos sea ice model, CICE. *Ocean Modelling*, 71, 26-42, https://doi.org/10.1016/j.ocemod.2012.11.008

Hunke, E., Allard, R., Blain, P. et al. (2020). Should Sea-Ice Modeling Tools Designed for Climate Research Be Used for Short-Term Forecasting? *Current Climate Change Reports*, 6, 121-136, https://doi.org/10.1007/s40641-020-00162-y

Hunke, E., Allard, R., Bailey, D.A., Blain, P., Craig, A., Dupont, F., DuVivier, A., Grumbine, R., Hebert, D., Hol-land, M., Jeffery, N., Lemieux, J.F., Osinski, R., Rasmussen, T., Ribergaard, M., Roberts, A., Denise, W. (2021). CICE Version 6.3.0. http://doi.org/10.5281/zenodo.1205674

Hutter, N., Losch, M., and Menemenlis, D. (2018). Scaling properties of Arctic sea ice deformation in a high-resolution viscous-plastic sea ice model and in satellite observations. *Journal of Geophysical Research: Oceans*, 123, 672-687, https://doi.org/10.1002/2017JC013119

Hutter, N., Zampieri, L., Losch, M. (2019). Leads and ridges in Arctic sea ice from RGPS data and a new tracking algorithm. *The Cryosphere*, 13, 627-645, https://doi.org/10.5194/tc-13-627-2019

Huwald, H., Tremblay, L.-B., Blatter, H. (2005). A multilayer sigma-coordinate thermodynamic sea ice model: Validation against Surface Heat Budget of the Arctic Ocean (SHEBA)/Sea Ice Model Intercomparison Project Part 2 (SIMIP2) data. *Journal of Geophysical Research: Oceans*, 110, C05010, https://doi.org/10.1029/2004JC002328

Iovino, D., Masina, S., Storto, A., Cipollone, A., and Stepanov, V. N. (2016). A 1/16° eddying simulation of the global NEMO sea-ice-ocean system. *Geoscientific Model Development*, 9, 2665-2684, https://doi.org/10.5194/gmd-9-2665-2016

Ivanova, N., Johannessen, O. M., Pedersen, L. T., and Tonboe, R. T. (2014). Retrieval of Arctic Sea Ice Parameters by Satellite Passive Microwave Sensors: A Comparison of Eleven Sea Ice Concentration Algorithms, IEEE Transactions on Geoscience and Remote Sensing, 52, 7233-7246, doi:10.1109/TGRS.2014.2310136

Ivanova, N., Pedersen, L. T., and Tonboe, R. T. (2015). Inter-comparison and evaluation of sea ice algorithms: towards further identification of challenges and optimal approach using passive microwave observations. *The Cryosphere*, 9, 1797-1817, https://doi.org/10.5194/tc-9-1797-2015

Janjić, T., McLaughlin, D., Cohn, S. E., and Verlaan, M. (2014). Conservation of Mass and Preservation of Positivity with Ensemble-Type Kalman Filter Algorithms. *Monthly Weather Review*, 142(2), 755-773, https://doi.org/10.1175/MWR-D-13-00056.1

<figure><img src="image_266.png" /></figure>
Jenouvrier, S., Holland, M., Strœve, J., Barbraud, C., Weimerskirch, H., Serreze, M., Caswell, H. (2012). Effects of climate change on an emperor penguin population: analysis of coupled demographic and climate models. Available at https://core.ac.uk/download/pdf/9584173.pdf

Johnson, M., Proshutinsky, A., Aksenov, Y., Nguyen, A.T., Lindsay, R., Haas, C., Zhang, J., Diansky, N., Kwok, R., Maslowski, W., Häkkinen, S., Ashik, I., de Cuevas, B. (2012). Evaluation of Arctic sea ice thickness simulated by Arctic Ocean Model Intercomparison Project models. *Journal of Geophysical Research: Oceans*, 117(C8), https://doi.org/10.1029/2011JC007257

Kaminski, T., Kauker, F., Eicken, H., and Karcher, M. (2015). Exploring the utility of quantitative network design in evaluating Arctic sea ice thickness sampling strategies. *The Cryosphere*, 9, 1721-1733, https://doi.org/10.5194/tc-9-1721-2015

Kauker, F., Kaminski, T., Karcher, M., Giering, R., Gerdes, R., and Voßbeck, M. (2009). Adjoint analysis of the 2007 all time Arctic sea-ice minimum. *Geophysical Research Letters*, 36, L03707, doi:10.1029/2008GL036323

Kimmritz, M., Losch, M., Danilov, S. (2017). A comparison of viscous-plastic sea ice solvers with and without replacement pressure. *Ocean Modelling*, 115, 59-69, https://doi.org/10.1016/j.ocemod.2017.05.006

Kjesbu, O.S., Bogstad, B., Devine, J.A., Gjøsæter, H., Howell, D., Ingvaldsen, R.B., Nash, R.D. M., and Skjæraasen, J. E. (2014). Synergies between climate and management for Atlantic cod fisheries at high latitudes. *Proceedings of the National Academy of Sciences of the United States of America*, 111(9), 3478-3483, https://doi.org/10.1073/pnas.1316342111

Korosov, A. A., and Rampal, P. (2017). A Combination of Feature Tracking and Pattern Matching with Optimal Parametrization for Sea Ice Drift Retrieval from SAR Data. *Remote Sensing*, 9(258), https://doi.org/10.3390/rs9030258

Kundzewicz, Z. W., Mata, L. J., Arnell, N. W., Doll, P., Kabat, P., Jimenez, B., Miller, K., Oki, T., Zekai, S. and Shiklomanov, I. (2007) Freshwater resources and their management. In: Parry, M. L., Canziani, O. F., Palutikof, J. P., van der Linden, P. J. and Hanson, C. E. (eds.) "Climate Change 2007: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change". Cambridge University Press, 173-210, ISBN 9780521880091

Kurtz, N., Farrell, S. L., Studinger, M., Galin, N., Harbeck, J. P., Lindsay, R., Onana, V. D., Panzer, B., and Sonntag, J. G. (2013). Sea ice thickness, freeboard, and snow depth products from Operation IceBridge airborne data. *The Cryosphere*, 7, 1035-1056, https://doi.org/10.5194/tc-7-1035-2013

Kwok, R. (2006). Contrasts in sea ice deformation and production in the Arctic seasonal and perennial ice zones. *Journal of Geophysical Research: Oceans*, 112(C2), https://doi.org/10.1029/2005JC003172

Kwok, R., Cunningham, G. F., Zwally, H. J., and Yi, D. (2007). Ice, Cloud, and land Elevation Satellite (ICESat) over Arctic sea ice: Retrieval of freeboard. *Journal of Geophysical Research: Oceans*, 112(C12), https://doi.org/10.1029/2006JC003978

Lauvernet, C., Brankart, J.-M. M., Castruccio, F., Broquet, G., Brasseur, P., and Verron, J. (2009). A truncated Gaussian filter for data assimilation with inequality constraints: Application to the hydrostatic stability condition in ocean models. *Ocean Modelling*, 27(1-2), 1-17. https://doi.org/10.1016/j.oceemod.2008.10.007

Laxon S. W., Giles, K. A., Ridout, A. L., Wingham, D. J., Willatt, R., Cullen, R., Kwok, R., Schweiger, A., Zhang, J., Haas, C., Hendricks, S., Krishfield, R., Kurtz, N., Farrell S., and Davidson, M. (2013). CryoSat-2 estimates of Arctic sea ice thickness and volume. *Geophysical Research Letters*, 40, 732-737, https://doi.org/10.1002/grl.50193

<figure><img src="image_267.png" /></figure>
Leith, C. E. (1974). Theoretical skill of Monte Carlo forecasts. Monthly Weather Review, 102(6), 409-418, https://doi.org/10.1175/1520-0493(1974)102<0409:TSOMCF>2.0.CO;2

Lemieux, J.-F., Tremblay, B., Sedlácek, J., Tupper, P., Thomas, S., Huard, D., and Auclair, J.-P. (2010). Improving the numerical convergence of viscous-plastic sea ice models with the Jacobianfree Newton Krylov method. Journal of Computational Physics, 229(8), 2840-2852, https://doi.org/10.1016/j.jcp.2009.12.011

Leutbecher, M., and Palmer, T. N. (2008). Ensemble forecasting. Journal of Computational Physics, 227(7), 3515-3539, https://doi.org/10.1016/j.jcp.2007.02.014

Lindsay, R. W., and Zhang, J. (2006). Assimilation of Ice Concentration in an Ice-Ocean Model. Journal of Atmospheric and Oceanic Technology, 23(5), 742-749, https://doi.org/10.1175/JTECH1871.1

Lindsay, R., Wensnahan, M., Schweiger, A., and Zhang, J. (2014). Evaluation of seven different atmospheric reanalysis products in the Arctic. Journal of Climate, 27(7), 2588-2606, https://doi.org/10.1175/JCLI-D-13-00014.1

Lisæter, K. A., Rosanova, J., Evensen, G. (2003). Assimilation of ice concentration in a coupled ice-ocean model, using the Ensemble Kalman filter. Ocean Dynamics, 53, 368-388, https://doi.org/10.1007/s10236-003-0049-4

Lisæter, K. A. (2009). A multi-category sea-ice model. Technical report 304. NERSC, Thormøhlensgt. 47, N-5006 Bergen, Norway.

Losch, M., D. Menemenlis, J.-M. Campin, P. Heimbach, and C. Hill, 2010. On the formulation of sea-ice models. Part 1: Effects of different solver implementations and parameterizations. Ocean Modelling, 33(1-2), 129-144, https://doi.org/10.1016/j.ocemod.2009.12.008

Lüpkes, C., Gryanik, V. M., Hartmann, J., Andreas, E. L. (2012). A parameterization, based on sea ice morphology, of the neutral atmospheric drag coefficients for weather prediction and climate models. Journal of Geophysical Research: Atmospheres, 117(D13), https://doi.org/10.1029/2012JD017630

Marsan, D., Stern, H., Lindsay, R., Weiss, J. (2004). Scale Dependence and Localization of the Deformation of Arctic Sea Ice. Physical Review Letter, 93(17), https://link.aps.org/doi/10.1103/PhysRevLett.93.178501

Maykut, G. A., and Untersteiner, N. (1971). Some results from a time-dependent thermodynamic model of sea ice. Journal of Geophysical Research (1896-1977), 76, 1550-1575, https://doi.org/10.1029/JC076i006p01550

McPhee, M. G. (1975). Ice-ocean momentum transfer for the AIDJEX ice model. AIDJEX Bulletin, 29, 93-111.

Melsom, A., Counillon, F., LaCasce, J. H., Bertino, L. (2012). Forecasting search areas using ensemble ocean circulation modeling. Ocean Dynamics, 62(8), 1245-1257, https://doi.org/10.1007/s10236-012-0561-5

Melsom, A., Palerme, C., Müller, M. (2019). Validation metrics for ice edge position forecasts. Ocean Science, 15(3), 615-630, https://doi.org/10.5194/os-15-615-2019

Mogensen, U. B., Ishwaran. H., Gerds, T. A. (2012). Evaluating Random Forests for Survival Analysis using Prediction Error Curves. Journal of Statistical Software, 50(11), 1-23, doi: 10.18637/jss.v050.i11

Molteni, F., Buizza, R., Palmer, T. N., Petroliagis, T. (1996). The ECMWF Ensemble Prediction System: Methodology and validation. Quarterly Journal of the Royal Meteorological Society, 122, 73-119, https://doi.org/10.1002/qj.49712252905

Motra, H. B., Stutz, H., Wuttke, F. (2016). Quality assessment of soil bearing capacity factor models of shallow foundations. Soils and Foundations, 56(2), 265-276, https://doi.org/10.1016/j.sandf.2016.02.009

<figure><img src="image_268.png" /></figure>
Mueter, F. J., Litzow, M. A. (2008). Sea ice retreat alters the biogeography of the Bering Sea continental shelf. *Ecological Applications*, 18(2), 309-20, https://doi.org/10.1890/07-0564.1

Notz, D., and Worster, M. G. (2009). Desalination processes of sea ice revisited. *Journal of Geophysical Research: Oceans*, 114(C5), https://doi.org/10.1029/2008JC004885

Notz, D., Jahn, A., Holland, M., Hunke, E., Massonnet, F., Stroeve, J., ... Vancoppenolle, M. (2016). The CMIP6 Sea-Ice Model Intercomparison Project (SIMIP): Understanding sea ice through climate-model simulations. *Geoscientific Model Development*, 9(9), 3427-3446, https://doi.org/10.5194/gmd-9-3427-2016

Oikkonen, A., Haapala, J., Lensu, M., Karvonen, J., Itkin P. (2017). Small-scale sea ice deformation during N-ICE2015: From compact pack ice to marginal ice zone. *Journal of Geophysical Research: Oceans*, 122(6), 5105-5120, https://doi.org/10.1002/2016JC012387

Ólason, E., Rampal, P., and Dansereau, V. (2021). On the statistical properties of sea-ice lead fraction and heat fluxes in the Arctic. *The Cryosphere*, 15, 1053-1064, https://doi.org/10.5194/tc-15-1053-2021

Ólason, E., Boutin, G., Korosov, A., Rampal, P., Williams, T., Kimmritz, M., Dansereau, V. (2022). A new brittle rheology and numerical framework for large-scale sea-ice models. Submitted to Journal of Modelling the Earth System (JAMES), https://doi.org/10.1002/essoar.10507977.2

Overland, J.E, Dethloff, K., Francis, J.A., Hall, R.J., Hanna, E., Kim, S.-J., Screen, J.A., Shepherd, T.G., Vihma, T. (2016). Nonlinear response of mid-latitude weather to the changing Arctic. *Nature Climate Change*, 6, 992-999, https://doi.org/10.1038/nclimate3121

Perovich, D. K., Richter-Menge, J. A., Jones, K. F., Light, B. (2008). Sunlight, water, and ice: Extreme Arctic sea ice melt during the summer of 2007. *Geophysical Research Letters*, 35(11), https://doi.org/10.1029/2008GL034007

Petty, A. A., Markus, T., Kurtz, N. (2017). Improving our understanding of Antarctic sea ice with NASA's Operation IceBridge and the upcoming ICESat-2 mission. *US CLIVAR Variations*, 15(3), doi:10.5065/D6833QQP

Posey, P. G., Metzger, E.J., Wallcraft, A.J., Hebert, D.A., Allard, R.A., Smedstad, O.M., Phelps, M.W., Fetterer, F., Stewart, J.S., Meier, W.N., and Helfrich, S.R. (2015). Improving Arctic sea ice edge forecasts by assimilating high horizontal resolution sea ice concentration data into the US Navy's ice forecast systems. *The Cryosphere*, 9, 1735-1745, https://doi.org/10.5194/tc-9-1735-2015

Rabatel, M., Labbé, S., and Weiss, J. (2015). Dynamics of an Assembly of Rigid Ice Floes. *J. Journal of Geophysical Research: Oceans*, 120(9), 5887-5909, https://doi.org/10.1002/2015JC010909

Rabatel, M., Rampal, P., Carrassi, A., Bertino, L., and Jones, C.K.R T. (2018). Impact of rheology on probabilistic forecasts of sea ice trajectories: Application for search and rescue operations in the Arctic. *The Cryosphere*, 12, 935-953, https://doi.org/10.5194/tc-12-935-2018

Rampal, P., Weiss, J., Marsan, D., Lindsay, R., and Stern, H. (2008). Scaling properties of sea ice deformation from buoy dispersion analysis. *Journal of Geophysical Research: Oceans*, 113(C3), C03002, https://doi.org/10.1029/2007JC004143

Rampal, P., Bouillon, S., Ólason, E., and Morlighem, M. (2016). NeXtSIM: A new Lagrangian sea ice model. *The Cryosphere*, 10, 1055-1073, https://doi.org/10.5194/tc-10-1055-2016

Rampal, P., Dansereau, V., Olason, E., Bouillon, S., Williams, T., Korosov, A., and Samaké, A. (2019). On the multi-fractal scaling properties of sea ice deformation. *The Cryosphere*, 13, 2457-2474, https://doi.org/10.5194/tc-13-2457-2019

<figure><img src="image_269.png" /></figure>
Richter-Menge, J. A., Perovich, D. K., Elder, B. C., Claffey, K., Rigor, I., Ortmeyer, M. (2006). Ice mass-balance buoys: a tool for measuring and attributing changes in the thickness of the Arctic sea-ice cover. *Annals of Glaciology*, 44(1), 205-210, doi:10.3189/172756406781811727

Rippeth, T.P. et al. (2015). Tide-mediated warming of Arctic halocline by Atlantic heat fluxes over rough topography. *Nature Geoscience*, 8(3), 191-194, https://doi.org/10.1038/ngeo2350

Rousset, C., Vancoppenolle, M., Madec, G., Fichefet, T., Flavoni, S., Barthélemy, A., Benshila, R., Chanut, J., Levy, C., Masson, S., and Vivier, F. (2015). The Louvain-La-Neuve sea ice model LIM3.6: global and regional capabilities. *Geoscientific Model Development*, 8, 2991-3005, https://doi.org/10.5194/gmd-8-2991-2015

Rösel, A., Kaleschke, L., and Birnbaum, G. (2012). Melt ponds on Arctic sea ice determined from MODIS satellite data using an artificial neural network. *The Cryosphere*, 6, 431-446, https://doi.org/10.5194/tc-6-431-2012

Sakov, P., Counillon, F., Bertino, L., Lisæter, K. A., Oke, P., and Korablev, A. (2012). TOPAZ4: an ocean-sea ice data assimilation system for the North Atlantic and Arctic. *Ocean Science*, 8, 633-656, https://doi.org/10.5194/os-8-633-2012

Schreyer, H. L., Sulsky, D. L., Munday, L. B., Coon, M. D., Kwok, R. (2006). Elastic-decohesive constitutive model for sea ice. *Journal of Geophysical Research: Oceans*, 111(C11), https://doi.org/10.1029/2005JC003334

Schulson, E. M., Hibler, W. D. (2017). The fracture of ice on scales large and small: Arctic leads and wing cracks. *Journal of Glaciology*, 27(127), 319-322, https://doi.org/10.3189/S0022143000005748

Semtner, A. J. (1976). A Model for the Thermodynamic Growth of Sea Ice in Numerical Investigations of Climate. *Journal of Physical Oceanography*, 6(3), 379-389. https://doi.org/10.1175/1520-0485(1976)006<0379:AMFTTG>2.0.CO;2

Simon, E., Samuelsen, A., Bertino, L., and Dumont, D. (2012). Estimation of positive sum-to-one constrained zooplankton grazing preferences with the DEnKF: a twin experiment. *Ocean Science*, 8(4), 587-602, https://doi.org/10.5194/os-8-587-2012

Smith, G. C., Roy, F., Reszka, M., Surcel Colan, D., He, Z., Deacu, D., Belanger, J.-M., Skachko, S., Liu, Y., Dupont, F., Lemieux, J.-F., Beaudoin, C., Tranchant, B., Drévillon, M., Garric, G., Testut, C.-E., Lellouche, J.-M., Pellerin, P., Ritchie, H., Lu, Y., Davidson, F., Buehner, M., Caya, A. and Lajoie, M. (2016). Sea ice forecast verification in the Canadian Global Ice Ocean Prediction System. *Quarterly Journal of the Royal Meteorological Society*, 142(695), 659-671, https://doi.org/10.1002/qj.2555

Smith, G. C., Liu, Y., Benkiran, M., Chikhar, K., Surcel Colan, D., Gauthier, A.-A., Testut, C.-E., Dupont, F., Lei, J., Roy, F., Lemieux, J.-F., and Davidson, F. (2021). The Regional Ice Ocean Prediction System v2: a pan-Canadian ocean analysis system using an online tidal harmonic analysis. *Geoscientific Model Development*, 14, 1445-1467, https://doi.org/10.5194/gmd-14-1445-2021

Sreen, G., Kwok, R., Menemenlis, D., and Nguyen, A. T. (2017). Sea-ice deformation in a coupled ocean-sea-ice model and in satellite remote sensing data. *The Cryosphere*, 11, 1553-1573, https://doi.org/10.5194/tc-11-1553-2017

Stark, J. D., Ridley, J., Martin, M., and Hines, A. (2008), Sea ice concentration and motion assimilation in a sea ice-ocean model. *Journal of Geophysical Research: Oceans*, 113(C5), https://doi.org/10.1029/2007JC004224

Stern, H., Lindsay, R. (2009). Spatial scaling of Arctic sea ice deformation. *Journal of Geophysical Research: Oceans*, 114(C10), https://doi.org/10.1029/2009JC005380

Sumata, H., Lavergne, T., Girard-Ardhuin, F., Kimura, N., Tschudi, M. A., Kauker, F., Karcher, M., Gerdes, R. (2014). An intercomparison of Arctic ice drift products to deduce uncertainty estimates. *Journal of Geophysical Research: Oceans*, 119(8), 4887-4921, https://doi.org/10.1002/2013JC009724

<figure><img src="image_270.png" /></figure>
Tandon, N. F., Kushner, P. J., Docquier, D., Wettstein, J. J., and Li, C. (2018). Reassessing sea ice drift and its relationship to long-term Arctic sea ice loss in coupled climate models. *Journal of Geophysical Research: Oceans*, 123, 4338-4359. https://doi.org/10.1029/2017JC013697

Thorndike, A. S., Rothrock, D. A., Maykut, G. A., Colony, R. (1975). The thickness distribution of sea ice. *Journal of Geophysical Research: Oceans*, 80(33), 4501-4513, https://doi.org/10.1029/JC080i033p04501

Thorndike, A. S., Colony, R. (1982). Sea ice motion in response to geostrophic winds. *Journal of Geophysical Research: Oceans*, 87(C8), 5845-5852, https://doi.org/10.1029/JC087iC08p05845

Tian-Kunze, X., Kaleschke, L., Maaß, N., Mäkynen, M., Serra, N., Drusch, M., and Krumpen, T. (2014). SMOS-derived thin sea ice thickness: algorithm baseline, product specifications and initial verification. *The Cryosphere*, 8, 997-1018, https://doi.org/10.5194/tc-8-997-2014

Tonani M., et al. (2015). Status and future of global and regional ocean prediction systems. *Journal of Operational Oceanography*, 8(sup2), s201-s220, http://doi.org/10.1080/1755876X.2015.1049892

Toyoda, T., Hirose, N., Urakawa, L. S., Tsujino, H., Nakano, H., Usui, N., Fujii, Y., Sakamoto, K., and Yamanaka, G. (2019). Effects of Inclusion of Adjoint Sea Ice Rheology on Backward Sensitivity Evolution Examined Using an Adjoint Ocean-Sea Ice Model. *Monthly Weather Review*, 147(6), 2145-2162, https://doi.org/10.1175/MWR-D-18-0198.1

Tremblay, L.-B., and Mysak, L. A. (1997). Modeling Sea Ice as a Granular Material, Including the Dilatancy Effect. *Journal of Physical Oceanography*, 27(11), 2342-2360. https://doi.org/10.1175/1520-0485(1997)027<2342:M-SIAAG>2.0.CO;2

Tsamados, M., Feltham, D. L., and Wilchinsky, A. V. (2013). Impact of a new anisotropic rheology on simulations of Arctic sea ice. *Journal of Geophysical Research: Oceans*, 118, 91-107, https://doi.org/10.1029/2012JC007990

Tsamados, M., Feltham, D. L., Schroeder, D., Flocco, D., Farrell, S. L., Kurtz, N., Bacon, S. (2014). Impact of Variable Atmospheric and Oceanic Form Drag on Simulations of Arctic Sea Ice. *Journal of Physical Oceanography*, 44, 1329-1353, https://doi.org/10.1175/JPO-D-13-0215.1

Turner, A. K., and Hunke, E. C. (2015). Impacts of a mushy-layer thermodynamic approach in global sea-ice simulations using the CICE sea-ice model. *Journal of Geophysical Research: Oceans*, 120, 1253-1275, https://doi.org/10.1002/2014JC010358

Turner, A. K., Hunke, E. C., and Bitz, C. M. (2013). Two modes of sea-ice gravity drainage: A parameterization for large-scale modeling. *Journal of Geophysical Research: Oceans*, 118, 2279-2294, https://doi.org/10.1002/jgrc.20171

Turner, A. K., Peterson, K. J., and Bolintineanu, D. (2022). Geometric remapping of particle distributions in the Discrete Element Model for Sea Ice (DEMSI v0.0). *Geoscientific Model Development*, 15, 1953-1970, https://doi.org/10.5194/gmd-15-1953-2022

Urrego-Blanco, J. R., Urban, N. M., Hunke, E. C., Turner, A. K., and Jeffery, N. (2016). Uncertainty quantification and global sensitivity analysis of the Los Alamos sea ice model. *Journal of Geophysical Research: Oceans*, 121, 2709-2732, https://doi.org/10.1002/2015JC011558

Usui, N., Ishizaki, S., Fujii, Y., Tsujino, H., Yasuda, T., Kamachi, M. (2006). Meteorological Research Institute multivariate ocean variational estimate (MOVE) system: Some early results. *Advances in Space Research*, 37/4), 806-822, https://doi.org/10.1016/j.asr.2005.09.022

<figure><img src="image_271.png" /></figure>
Vancoppenolle, M., Bitz, C. M., and Fichefet, T. (2007). Summer landfast sea ice desalination at Point Bar- row, Alaska: Modeling and observations. *Journal of Geophysical Research: Oceans*, 112(C4), https://doi. org/10.1029/2006JC003493

Vancoppenolle, M., Fichefet, T., Goosse, H., Bouillon, S., Madec, G., and Maqueda, M. A. M. (2009). Simu- lating the mass balance and salinity of Arctic and Antarctic sea ice. 1. Model description and validation. *Ocean Modelling*, 27(1-2), 33-53, https://doi.org/10.1016/j.ocemod.2008.10.005

Wackernagel, H. (2003). Multivariate Geostatistics. An Introduction with Applications. Springer Science & Business Media, 12 feb 2003 - 388 pp.

Waters, J., Lea, D.J., Martin, M.J., Mirouze, I., Weaver, A., and While, J. (2015). Implementing a variational data assimilation system in an operational 1/4 degree global ocean model. *Quarterly Journal of the Roy- al Meteorological Society*, 141(687), 333-349, https://doi.org/10.1002/qj.2388

Weiss, J., and Marsan, D. (2004). Scale properties of sea ice deformation and fracturing. *Comptes Rendus Physique*, 5(7), 735-751, https://doi.org/10.1016/j.crhy.2004.09.005

Wilchinsky, A. V., and Feltham, D. L. (2004). A continuum anisotropic model of sea-ice dynamics. *Proceed- ings of the Royal Society A*, 460(2047), https://doi.org/10.1098/rspa.2004.1282

Williams, T., Korosov, A., Rampal, P., Ólason, E. (2021). Presentation and evaluation of the Arctic sea ice forecasting system neXtSIM-F. *The Cryosphere*, 15, 3207-3227, https://doi.org/10.5194/tc-15-3207-2021

Winton A. (2000). A Reformulated Three-Layer Sea Ice Model. *Journal of Atmospheric and Oceanic Tech- nology*, 17(4), 525-531, https://doi.org/10.1175/1520-0426(2000)017<0525:ARTLSI>2.0.CO;2

WMO (1970). Sea-ice nomenclature, terminology, codes and illustrated glossary. WMO Rep. WMO/OMM/ BMO 259, 121 pp.

Xie, J., Bertino, L., Counillon, F., Lisæter, K. A., and Sakov, P. (2017). Quality assessment of the TOPAZ4 reanalysis in the Arctic over the period 1991-2013. *Ocean Science*, 13(1), 123-144, http://doi.org/10.5194/ os-13-123-2017

Zygmuntowska, M., Rampal, P, Ivanova, N., and Smedsrud, L. H. (2014). Uncertainties in Arctic sea ice thickness and volume: new estimates and implications for trends. *The Cryosphere*, 8, 705-720, https:// doi.org/10.5194/tc-8-705-2014

<figure><img src="image_272.png" /></figure>
# 7. Storm surge modelling

<figure><img src="image_273.png" /></figure>

<figure><img src="image_274.png" /></figure>

CHAPTER COORDINATOR

Fujiang Yu

CHAPTER AUTHORS (in alphabetical order)

David Byrne, Jianxi Dong, Begoña Pérez Gómez, and Shichao Liu

<figure><img src="image_275.png" /></figure>
# 7. Storm surge modelling

## 7.1. General introduction to storm surge

### 7.1.1. Overview of storm surge disaster

7.1.1.1. Disasters and forecasting

7.1.1.2. The impact of climate change on storm surge

### 7.1.2. Basic description of storm surge phenomena

#### 7.1.3. Physics of storm surge

7.1.3.1. Meteorological forcing

7.1.3.2. The influence of topography and bathymetry

## 7.2. Storm surge modelling

### 7.2.1. Architecture components and singularities

### 7.2.2. Input data: available sources and data handling

7.2.2.1. Bathymetry and geometry

7.2.2.2. Tidal boundaries

7.2.2.3. Meteorological inputs

### 7.2.3. Modelling component

7.2.3.1. Governing equations

7.2.3.2. 2D barotropic and 3D baroclinic models for storm surge

7.2.3.3. Wetting and drying scheme

7.2.3.4. Grid types

7.2.3.5. Discretization method

7.2.3.6. Existing models for storm surge modelling

### 7.2.4. Data assimilation systems

7.2.4.1. Sources of error in storm surge models

7.2.4.2. Assimilated data sources for storm surge modelling

7.2.4.3. Application of data assimilation to real time forecasting systems

7.2.4.4. Examples from real operational systems

### 7.2.5. Ensemble modelling

### 7.2.6. Validation strategies

### 7.2.7. Outputs

7.2.7.1. Time series outputs

7.2.7.2. Maximum elevation field

7.2.7.3. Ensemble forecast field

7.2.7.4. Animation output

### 7.2.8. Existing operational storm surge forecasting systems

## 7.3. References

<figure><img src="image_276.png" /></figure>

<figure><img src="image_277.png" /></figure>
## 7.1. General introduction to storm surge

Many natural phenomena can cause the sea to rise and fall, such as wind, air pressure, celestial gravity, earthquakes, etc. The sea level changes caused by different phenomena have different periods. For example, wind waves have a period of several seconds, tsunami waves of few minutes to tens of minutes, and the period of storm surge and astronomical tide is about several hours to several days (Figure 7.1).

Among them, the storm surge brings huge economic losses and risks to coastal countries every year (Murty, 1988). In order to reduce the impact of storm surge disasters on coastal residents, understanding and forecasting storm surge have always been an important objective for marine forecasters.

This chapter will introduce the main overview and elements of storm surge modelling, to guide technical personnel to engage in related work and give full play to the role of storm surge numerical models in various fields.

<figure><img src="image_278.png" /><figcaption>**Figure 7.2.** The impact of the storm surge caused by the super typhoon Haiyan on the Philippines, the coastal villages of Tacloban were destroyed (Credits: Photography Marcel Crozet, ILO, 11-2013).</figcaption></figure>

### 7.1.1. Overview of storm surge disaster

#### 7.1.1.1. Disasters and forecasting

Storm surge refers to the phenomenon of abnormal water level rise in a coastal or inland body caused by strong atmospheric disturbances, such as tropical cyclones (typhoons, hurricanes), extratropical cyclones, strong winds from cold fronts, and sudden change in atmospheric pressure.

As a complex coastal dynamic process of major coastal marine disasters, storm surge has received much attention by major affected countries all over the world. Storm surge disasters are mainly caused by the abnormal water level rise and by flooding. The disaster causing factors include not only the storm surge, but also coupling with the effect of astronomical tide and nearshore waves. Storm surge disasters (Figure 7.2, Figure 7.3, and Figure 7.4) not only include the damage to ports,

<figure><img src="image_279.png" /><figcaption>**Figure 7.1.** Frequencies and periods of the vertical motions of the ocean surface (adapted from Pérez et al., 2013).</figcaption></figure>

<figure><img src="image_280.png" /></figure>
<figure><img src="image_281.png" /><figcaption>**Figure 7.3.** People walk among debris next to a ship washed ashore in the aftermath of super typhoon Haiyan in Tacloban, Philippines, 11 November 2013. (Credits: ILO, 11-2013).</figcaption></figure>

wharfs, dykes, but also include the disasters caused by flooding houses, farmland, and aquaculture facilities.

Areas with severe storm surges are shown on a global map (Figure 7.5). The Gulf Coast of North America and the eastern coast of the United States are affected by storm surges generated by Atlantic hurricanes. In Europe, the North Sea coast is often affected by extratropical cyclones, which bring storm surge disaster. The coast of the Bay of Bengal in the Indian Ocean is threatened by storm surges caused by typhoons in the Indian Ocean. On the western Pacific Ocean coast, China, Japan, and the Philippines are frequently affected by storm surges caused by typhoons, and the north coast of China is also affected by extratropical cyclones.

In addition to the areas severely affected by the storm surge mentioned above, other countries or regions may also be affected by the storm surge. Areas with low elevation may face the threat of storm surge inundation, and the approaching channel may not meet the navigation requirements due to the drop in water level. For example, 20% of the land in the Netherlands is below mean sea level, and large areas of flooding may be caused without a very severe storm surge. For this reason, they built the famous Storm Surge Barriers

(Mooyaart and Jonkman, 2017). In Spain, surges of 60 cm contribute significantly to inundation processes.

<figure><img src="image_282.png" /><figcaption>**Figure 7.4.** The impact of the storm surge caused by the super typhoon Haiyan on the Philippines, the coastal villages of Tacloban were inundated with water (Credits: Photography Marcel Crozet, ILO, 11-2013).</figcaption></figure>

<figure><img src="image_283.png" /></figure>
<figure><img src="image_284.png" /><figcaption>**Figure 7.5.** Areas severely affected by storm surge.</figcaption></figure>

Storm surge forecasting is an important means of reducing disasters and losses, and a very necessary link in disaster prevention and mitigation. The methods of storm surge forecasting can be divided into two categories: empirical statistical forecasting and numerical forecasting. With the rapid development of computer technology, numerical models play an increasingly important role in storm surge forecasting. The establishment of a storm surge numerical model will provide strong support for storm surge forecasting. In addition to providing help for disaster prevention and mitigation, the numerical model of storm surge can also be used in offshore engineering design and marine disaster risk assessment of coastal cities.

In recent years, with the rapid economic development of coastal cities and the urgent needs of disaster prevention and mitigation, more and more ocean forecasting centres have started to establish operational storm surge models to provide relevant services for the above activities and purposes (more information in Section 7.2.8).

### 7.1.1.2. The impact of climate change on storm surge

Coastal cities are directly affected by global warming, sea temperature continues to increase, sea level fluctuates and rises, and natural disasters such as storm surges and huge waves show an increasing trend. Statistics show that there is a significant increase in global super typhoons (or category 4 and 5 hurricanes). In the 1970s, the number of super

typhoons accounted for 20% of the total tropical cyclones, while it rose to 35% in the 1990s. Among them, the most evident increase was in the North Pacific, Indian Ocean, and Southwest Indian Ocean, while the increase was the least in the North Atlantic (Webster et al., 2005). Therefore, storm surge disasters caused by typhoons showed an increasing trend, as well as the risk of storm surge disasters in coastal cities. The tide observation data also shows this characteristic. After the storm surge of typhoon Hato (2017) and typhoon Mangkhut (2018) affected coastal cities such as Zhuhai and Shenzhen in China, the return period of coastal tide levels changed significantly. The Hengmen Station, located on the west bank of the Pearl River Estuary (China), has shown that the tide level return period has been reduced from 200 years to 50 years, as well as the original 100 years tide level return period been reduced to 50 years at the Sanzao Station, and the Chiwan Station on the east bank of the Pearl River Estuary.

Sea level rise directly leads to the expansion of storm surge inundation area, increases the mean sea level, and various characteristic tide levels. The increased water depth and enhanced nearshore waves further strengthen the impact of storm surges.

### 7.1.2. Basic description of storm surge phenomena

Storm surges have periods of several hours to several days, and are usually superimposed on tides, wind waves and

<figure><img src="image_285.png" /></figure>
<figure><img src="image_286.png" /><figcaption>**Figure 7.6.** Storm surge components and drivers.</figcaption></figure>

swells (with a period of several seconds). Combination of these three factors causes extreme rise of coastal water level and often leads to huge storm surge disasters. However, sometimes the opposite situation can also be encountered: the wind blowing away from the direction of the open coast for a long time causes the water level to drop sharply along the shore and shoals exposed. In this case, the normal navigation is seriously affected, as well as anchoring of ships, especially large oil tankers.

The spatial range of storm surges is generally between tens and thousands of kilometres, and the time scale or period is about several to hundreds hours, which is between a tsunami and the astronomical tide. Since the area affected by storm surges moves with the movement of meteorological forcing, sometimes a storm surge process can affect a coastal area of 1000-2000km, and the impact time can be up to several days. In addition, the period of water level change by a storm surge ranges from some hours to several days, excluding seiches, tsunamis and wind waves.

According to its standard definition, a storm surge is caused by atmospheric disturbance, specifically abnormal alterations in water surface due to strong winds and atmospheric pressure changes. Storm surge can also occur in inland bodies, such as the Great Lakes in the US. In recent years, studies have shown that the nearshore waves breaking can also cause rise of the water level, in the range of tens of centimetres to metres, called wave setup. With the perspective of changes occurred in modern times, the definition of storm surges should be revised as the following: "storm surge refers to strong atmospheric disturbances, such as tropical cyclones (typhoons and hurricanes), extratropical cyclones, strong wind due to cold fronts, and sudden changes in atmospheric pressure inducing

abnormal water level rise combined with nearshore wave set-up" (Yu et al., 2020). See representation of storm surge components and drivers in Figure 7.6.

Meteorological tsunami, or meteotsunami, is caused by strong winds and sudden changes in atmospheric pressure and its period is equivalent to a tsunami. In the Tsunami Glossary (Tsunami) by the IOC's ITIC, meteotsunami is defined as following: "Tsunami-like phenomena generated by meteorological or atmospheric disturbances. These waves can be produced by atmospheric gravity waves, pressure jumps, frontal passages, squalls, gales, typhoons, hurricanes and other atmospheric sources. Meteotsunamis have the same temporal and spatial scales as tsunami waves and can similarly devastate coastal areas, especially in bays and inlets with strong amplification and well defined resonant properties (e.g. Ciutadella Inlet, Baleric Islands; Nagasaki Bay, Japan; Longkou Harbour, China; Vela Luka, Stari Grad and Mali Ston Bays, Croatia)."

The water level recorded at coastal or estuarine tide stations usually include a combination of changes caused by astronomical tides, storm surges, tsunamis, and other long waves. Generally, tide gauges filter out sea surface fluctuations caused by short-period waves in the order of seconds. The separation of storm surge phases is obtained by linear subtracting the harmonic analysis forecast astronomical tide from the hourly data (Figure 7.7 and Figure 7.8).

<figure><img src="image_287.png" /></figure>

<footnote n="1">1. http://itic.ioc-unesco.org/index.php?option=com_content&view=category&id=2339&Itemid=2433</footnote>
<figure><img src="image_288.png" /><figcaption>**Figure 7.7.** Observed water level, astronomical tide, and storm surge (water level subtracting astronomical tide); data from Zhanjiang tide station (China).</figcaption></figure>

<figure><img src="image_289.png" /><figcaption>**Figure 7.8.** Observed water level, astronomical tide, and storm surge (water level subtracting astronomical tide); data from Nandu tide station (China).</figcaption></figure>

<figure><img src="image_290.png" /></figure>
### 7.1.3. Physics of storm surge

#### 7.1.3.1. Meteorological forcing

Meteorological forcing is the main driver for storm surges. When a storm passes over the open sea, the low centre pressure of the storm will cause the water level to rise. The height of the surge is related to the barometric pressure drop of the storm, i.e. 1 mbar decrease corresponds approximately to 1 cm increase in sea level (Schalkwijk, 1948; Myers, 1954; Pore, 1964). The raised sea surface will propagate with the movement. At the same time, a kind of free long wave, induced by raised sea surface, could spread outward from the storm centre. This process will typically take place near the coast when interactions with bathymetry changes become relevant. If the pressure disturbance is moving at a speed comparable to the shallow water wave speed, the water level disturbance may be greatly amplified by resonance (Harris, 1957).

Compared with the long wave effect, the wind shear stress is the dominating forcing of storm surges in shallow water of nearshore and estuaries (Miller, 1958; Pore, 1964). With the wind blowing continuously, water accumulates at the coastal line causing the water level to rise. This phenomenon is referred to as "wind set-up" and its magnitude is inversely proportional to water depth. The wind set-up is particularly evident in semi-enclosed seas, such as Bohai Bay in China.

#### 7.1.3.2. The influence of topography and bathymetry

Storm surge is not only influenced by astronomical tide and waves, but also by topography and bathymetry. Due to the shoreline block, storm surge propagates from ocean to nearshore. The surge is generated by water accumulation at the shoreline. The magnitude of the surge is controlled by the

shape of the shoreline. In case of onshore direction, semi-enclosed bays or estuaries contribute to intensify storm surge than straight shoreline. That is because the shape of the semi-enclosed bay and estuary hinder water flow out. The water accumulates in the semi-enclosed bay or estuary continuously, resulting in a greater storm surge.

Another factor that can impact storm surge is the variation of bathymetry from the continental shelf to estuaries and coasts. Generally, the water depth of estuaries and coasts is shallower than the continental shelf, and the propagation speed of the storm surge wave is approximately proportional to the square root of the water depth. Therefore, the speed of the wave propagation at estuaries and coasts is slower than at the continental shelf. The storm surge waves converge at estuaries and coasts, causing the water level to increase.

On the other hand, in the process of storm surge wave propagation, the water depth at the crest is greater than at the preceding trough, and the movement of the crest is faster. So, the more waves move inland, the smaller the interval between the crests. This is more pronounced where the continental shelf is longer, for example in the North Sea, and hence larger storm surges will be caused due to the long continental shelf extension.

The propagation speed of storm surge waves is controlled by the water depth: it moves faster at high tide than at low tide. The wind effect is inversely proportional to the total water depth, and the same wind speed will produce a greater surge at low tide than at high tide. Combining the two effects above, surge in an estuary tends to be greater on the rising stage of the tide (Doodson, 1929; Doodson, 1956; Rossiter, 1961).

Extremely accurate topography and bathymetry, especially for shallow water areas, is key to storm surge modelling.

## 7.2.

# Storm surge modelling

### 7.2.1. Architecture components and singularities

Storm surge models are generally based on the two-dimensional shallow water equation to compute the water level and velocity. According to different modelling purposes, the storm surge model can be divided into: i) storm surge model without tide; ii) storm surge model including astronomical tide; and iii) storm surge flooding model considering inundation.

In a storm surge model without tide, only the effect of the meteorological forcing needs to be considered. Obtained result is only the water level rising and falling under the effect of atmospheric pressure and wind. The role of this model is generally to provide forecasters with a reference for the magnitude of storm surge when making forecasts, and it cannot truly reflect the fluctuations of water level.

Based on the former model, the tidal boundary conditions and tidal potential are introduced in the total water level storm surge model, in which the nonlinear interaction between
<figure><img src="image_291.png" /><figcaption>**Figure 7.9.** Storm surge modelling flow chart.</figcaption></figure>

storm surge and astronomical tide can be fully considered in the model. This model can be used for releasing storm surge numerical forecast products and providing reference for coastal response to storm surge disasters. It can also provide the water level changes in the target area under the influence of extreme weather conditions for coastal engineering (such as harbours, wharfs, seawalls, offshore wind farms, etc.).

The storm surge flooding model considering inundation is more complicated than the previous two models. The interaction of storm surge and astronomical tide, and the interaction of storm surge and wave are also considered in the model. The inundation range and depth of the coastal area can be obtained by a storm surge simulation. Figure 7.9 shows the main storm surge modelling components used by a forecasting system and that will be detailed in the next sections.

## 7.2.2. Input data: available sources and data handling

### 7.2.2.1. Bathymetry and geometry

Reliable and accurate shoreline and bathymetry data are the basis for modelling of storm surge, tidal, and storm surge flooding models. In the process of model setup, the computational grid should be determined according to the demand, and then the data of shoreline and bathymetry should be collected according to the location and scope.

Sources of bathymetric data can be found on Section 4.2.4. These data need to be used with caution when establishing high-resolution storm surge models. As the bathymetry of coastal or estuary changes rapidly over time, these upstream data may not be able to be updated in time. Therefore, the correctness of the data needs to be verified before it is used.
When the published upstream data do not meet the requirements, it is preferable to use the latest sea chart or high-resolution DEM data. In addition, the datum of the data needs to be confirmed as, due to different sources of data, the datum could be different. In order to truly reflect the effects of bathymetry on storm surges, data from different sources need to be unified on the same datum.

### 7.2.2.2. Tidal boundaries

In the storm surge and storm surge flooding modelling, tidal waves are generally used as open boundary conditions. Open boundary conditions can be velocity, water level, or harmonic constant.

Tidal height at any location and time can be written as a function of harmonic constituents according to the following general relationship:

$$H(t) = H_0 + \sum_{n=1}^{N} f_n H_n \cos[a_n t + (V_0 + u)_n - g_n] \quad (7.1)$$

where

$H(t)$ = height of the tide at any time $t$

$H_0$ = mean water level above some defined datum such as mean sea level

$H_n$ = mean amplitude of tidal constituent $n$

$f_n$ = factor for adjusting mean amplitude $H_n$ values for specific times

$a_n$ = speed of constituent $n$ in degrees/unit time

$t$ = time measured from some initial epoch or time, i.e., $t=0$ at $t_0$

$(V_0 + u)_n$ = value of the equilibrium argument for constituent $n$ at Greenwich and when $t=0$

$g_n$ = epoch of constituent $n$, i.e., phase shift from tide-producing force to high tide from $t_0$

In Eq. 7.1, for a tidal component, the $f_n$ and the $(V_0 + u)_n$ will change with time, and the $H_n$ and the $g_n$ will change with the geographical location. Therefore, according to the start time of simulation, the $f_n$ and $(V_0 + u)_n$ of the proposed tidal component should be set in the model, and the $H_n$ and the $g_n$ boundary conditions should be given according to the position of the boundary grid node.

The location of the tidal open boundary is critical. First of all, it is necessary to ensure that the grid can completely cover the definition area. Secondly, it is better to set the open

boundary near the tidal station, because the tidal constituent data near the tidal station is more accurate. Finally, it is preferable not to set the open boundary at the amphidromic point nearby, because too small amplitude will bring uncertainty to the simulation.

The data for tidal boundaries can be obtained by downloading publicly available data on tidal harmonic constants covering most of the oceans. TPXO (Egbert et al., 1994; Egbert and Erofeeva, 2002) is a widely used global tidal data. It is a series of fully global models of ocean tides, which best fits, in a least squares sense, the Laplace Tidal Equations and altimetry data. The TPXO models include eight primary (M2, S2, N2, K2, K1, O1, P1, Q1), two long period (Mf, Mm) and 3 non-linear (M4, MS4, MN4) harmonic constituents (plus 2N2 and S1 for TPXO9 only). More detailed information can be found at $\mathbb{S}^2$. In addition, also the NAO.99b tide model (Matsumoto et al., 2000 and 2001), which is developed by assimilating TOPEX/POSEIDON Altimeter Data into Hydrodynamical Model, can provide global tide data. This model provides 16 short-period harmonic constituents (M2, S2, K1, O1, N1, N2, P1, K2, Q1, M1, J1, O01, 2N2, Mu2, Nu2, L2, L2, and T2), 7 long-period harmonic constituents (Mf, Mm, Msf, Msm, Mtm, Sa, Ssa) data with of 0.5 degrees, and provides 16 short-period harmonic constants around Japan with a resolution of 5 minutes. More detailed information can be found at $\mathbb{S}^3$.

Once the model and tidal boundaries have been established, they need to be tuned and calibrated before using. The model can be run with tidal boundaries and tidal potential. The tidal results are more sensitive to changes in the bottom friction coefficient and the bathymetry. By calibrating the tidal simulation, a reasonable bottom friction coefficient can be set for the model. At the same time, the difference in the tidal results caused by the bathymetry from dissimilar sources helps to find more suitable bathymetry data for the model.

### 7.2.2.3. Meteorological inputs

In storm surge or storm surge flooding modelling, the input from meteorological forcing mainly includes surface wind shear stress and atmospheric pressure at sea surface level. In deep water, the sea level rises are mainly caused by the atmospheric pressure gradient, i.e. the water level rises approximately 1 cm at every reduction in pressure of 1 mbar. In shallow water, estuary and nearshore, wind shear stress is the dominant force of the storm surge, and the sea level rise is proportional to the square of wind speed. Therefore, accurate meteorological inputs, especially sea surface wind, is essential for storm surge modelling. The accuracy of storm surge results depends largely on the quality of meteorological data.

<figure><img src="image_292.png" /></figure>

<footnote n="2">2. https://www.tpxo.net/global</footnote>

<footnote n="3">3. https://www.miz.nao.ac.jp/staffs/nao99/index_En.html</footnote>
Depending on the storm surge modelling purposes, the types of meteorological forcing input are different. At present, there are two main sources of wind field for storm surge modelling: atmospheric model and empirical formula. Atmospheric models can provide global or regional meteorological fields; the main elements required for storm surge calculation are sea level pressure and 10 metres wind. The horizontal resolution of these data is between ten and tens of kilometres, and the forecast period can reach up to 240 hours. This kind of meteorological field is mainly used in the calculation of storm surge caused by extratropical cyclones. Compared with extratropical cyclones, tropical cyclones are smaller in scale but stronger in intensity, and atmospheric models cannot resolve the structure of tropical cyclones well. Therefore, the meteorological field from atmospheric models is not applicable to the typhoon storm surge calculation. Empirical formulas for tropical cyclone pressure and wind are often applied to create meteorological forcing fields for tropical storm surge models. Since the wind speed and pressure structure of tropical cyclones are close to axisymmetric, their distribution can be reasonably represented with an empirical formula for the radial distribution of wind or pressure.

The commonly used empirical formulas for pressure distribution mainly include the following:

Takahashi (1939):

$$ \frac{P(r) - P_0}{P_{\infty} - P_0} = 1 - \frac{1}{1 + \frac{r}{R}} \qquad (7.2) $$

Myers (1954):

$$ \frac{P(r) - P_0}{P_{\infty} - P_0} = 1 - \frac{1}{\sqrt{1 + 2\left(\frac{r}{R}\right)^2}} \qquad (7.3) $$

Myers (1954):

$$ \frac{P(r) - P_0}{P_{\infty} - P_0} = e^{-\frac{r}{R}} \qquad (7.4) $$

Jelesnianski (1965):

$$ \frac{P(r) - P_0}{P_{\infty} - P_0} = \frac{1}{4} \left(\frac{r}{R}\right)^3, \quad 0 \le r < R \qquad (7.5) $$

$$ \frac{P(r) - P_0}{P_{\infty} - P_0} = 1 - \frac{3R}{4r}, \quad R \le r < \infty \qquad (7.6) $$

Bjerknes (1921):

$$ \frac{P(r) - P_0}{P_{\infty} - P_0} = 1 - \left[1 + \left(\frac{r}{R}\right)^2\right]^{-1} \qquad (7.7) $$

Holland (1980):

$$ \frac{P(r) - P_0}{P_{\infty} - P_0} = e^{-\frac{A}{R^B}} \qquad (7.7) $$

where:

$P_\infty$ = the ambient pressure;

$P_0$ = the pressure at the tropical cyclone centre;

$P(r)$ = the pressure at $r$ from tropical cyclone centre;

A and B = empirical parameters that control the tropical cyclone size.

The tropical cyclone wind field is formed by the superposition of two vector fields (Ueno, 1981). The first vector field is a wind field symmetrical to the centre of the cyclone. The wind vector passes through the isobar and points to the left with a 20° deflection angle.

The wind speed is proportional to the gradient wind, which can be expressed by the following formula:

$$ V_v = V_r \frac{2Rr}{R^2 + r^2} \qquad (7.8) $$

$V_r$ = the maximum wind speed;

R = the radius of the maximum wind.

The second vector field caused by the movement of cyclone is superimposed on the wind system for the stationary symmetric cyclone, and that the wind velocity $\vec{V}_f$ is:

$$ \vec{V}_f = \vec{V}_t \exp \left( -\frac{\pi |r - R|}{4R} \right) \qquad (7.9) $$

$\vec{V}_t$ = velocity of typhoon;

R = the radius of the maximum wind.

Consequently, the wind velocity $\vec{W}$ is:

$$ \vec{W} = \vec{V}_v + \vec{V}_f \qquad (7.10) $$

The typhoon centre pressure, maximum wind speed, moving speed and other parameters can be obtained from websites of national meteorological agencies.

### 7.2.3. Modelling component

#### 7.2.3.1. Governing equations

The governing equations of numerical simulation were determined by the flow dynamic theory.

The three-dimensional flow equations are as follows:

<figure><img src="image_293.png" /></figure>
* Continuity equation:

$$ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0 \qquad (7.11) $$

* Momentum equation:

$$ \frac{\partial u}{\partial t} + \frac{\partial(uu)}{\partial x} + \frac{\partial(uv)}{\partial y} + \frac{\partial(uw)}{\partial z} - f_v = -\frac{1}{\rho} \frac{\partial p}{\partial x} + \frac{1}{\rho} \left( \frac{\partial \tau_{xx}}{\partial x} + \frac{\partial \tau_{xy}}{\partial y} + \frac{\partial \tau_{xz}}{\partial z} \right) \qquad (7.12a) $$

$$ \frac{\partial v}{\partial t} + \frac{\partial(vu)}{\partial x} + \frac{\partial(vv)}{\partial y} + \frac{\partial(vw)}{\partial z} - f_u = -\frac{1}{\rho} \frac{\partial p}{\partial y} + \frac{1}{\rho} \left( \frac{\partial \tau_{xy}}{\partial x} + \frac{\partial \tau_{yy}}{\partial y} + \frac{\partial \tau_{yz}}{\partial z} \right) \qquad (7.12b) $$

$$ \frac{\partial p}{\partial z} = -\rho g \qquad (7.12c) $$

Based on hydrostatic approximations and incompressible assumption of fluid, the depth integrated two-dimensional storm surge governing equations can be written as:

* Continuity equation:

$$ \frac{\partial \xi}{\partial t} + \frac{\partial}{\partial x}[(\xi + h)u] + \frac{\partial}{\partial y}[(\xi + h)v] = 0 \qquad (7.13) $$

* Momentum equation:

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} - f_v + g \frac{\partial \xi}{\partial x} \pm \frac{\tau_{x,S} - \tau_{x,b}}{pH} + \frac{1}{\rho} \frac{\partial P}{\partial x} = \epsilon \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \qquad (7.14a) $$

$$ \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} - f_u + g \frac{\partial \xi}{\partial y} \pm \frac{\tau_{y,S} - \tau_{y,b}}{pH} + \frac{1}{\rho} \frac{\partial P}{\partial y} = \epsilon \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \qquad (7.14b) $$

where:

$\xi$= free surface elevation relative to the geoid;

$h$ = water depth;

$f$ = Coriolis coefficient;

$\rho$ = density of water;

$g$ = gravitational acceleration;

($\tau_{x,S}$, $\tau_{x,b}$) = free-surface shear stress in x and y direction;

$$ (\tau_{x,b}, \tau_{y,b}) = C_d \rho_a |W|(W_x, W_y) $$

$W$ = wind speed at 10 metres above sea surface;

$W_x$, $W_y$ = wind speed components in x and y direction;

$C_d$ = wind Drag coefficient which is relevant to wind speed;

$\tau_{x,b}$, $\tau_{y,b}$ = bottom shear stress in x and y direction;

$$ (\tau_{x,b}, \tau_{y,b}) = \rho g \sqrt{u^2 + v^2} \frac{C^2}{C} (u, v) $$

$C$ = Chezy coefficient, $C = \frac{h^{1/6}}{n}$ n is roughness coefficient;

$u$, $v$ = depth-averaged horizontal velocity components in x and y direction;

$$ u = \xi + h \int_{-h}^{\xi} u dz $$

$$ v = \xi + h \int_{-h}^{\xi} v dz $$

$P$ = atmospheric pressure at the free surface;

$\epsilon$ = depth-averaged horizontal eddy viscosity.

### 7.2.3.2.2D barotropic and 3D baroclinic models for storm surge

Hydrodynamic models are generally divided into vertically averaged 2D barotropic, 3D barotropic, or 3D baroclinic models. These types of models can be used for storm surge modelling. Minato (1998) first studied the effect of a 3D model on storm surge results and simulated the water level change caused by Typhoon 7010 in Tosa Bay, Japan. The results showed that the difference between the 3D model and the 2D model is about 2 to 10%. Weisberg and Zheng (2008) found that under the condition of setting the same bottom friction coefficient, the 3D model simulates higher extreme water level than the 2D model. Subsequently, simulating the storm surge of Hurricane Ike in the Gulf of Mexico, Zheng et al. (2013) pointed out that the 2D and the 3D models have differences in the trend and peak values of water level, but the calibration of the bottom friction is more important. Ye et al. (2020) studied the influence of baroclinic models on storm surge simulation results. Sensitivity tests show that the impact of the baroclinic model on storm surge is not significant, but it has a greater impact on current.

<figure><img src="image_294.png" /></figure>
Since the vertical velocity distribution structure of the 2D model is different from that of the 3D model, the vertical average velocity is greater than that near the bottom layer, so that the bottom shear stress of the 2D model will be greater. Satisfactory results can be obtained for both types of models by calibrating the bottom friction.

Regarding operational storm surge modelling, computational efficiency is a factor that must be considered. 3D models generally divide the water body into multiple layers vertically, so they require more computing time than 2D models. If the purpose of storm surge modelling is to obtain water level, rather than currents, the 2D models are the best choice for balancing calculation efficiency and accuracy. Most of the existing operational storm surge forecasting systems are 2D models.

#### 7.2.3.3. Wetting and drying scheme

During the simulation of storm surge inundation, the grid points close to the coastline in the model will be either wet or dry due to the fluctuation of water level. Therefore, the model needs to determine the dry and wet state of a grid point according to the state of the surrounding grid points.

Assuming that the state of a certain grid point is wet, if the calculated water level makes the water volume less than zero, then this grid point will become a dry grid point and will not participate in the next calculation. In practical calculations, to prevent the negative value of the water volume at a grid point, which makes the momentum equation and continuity equation meaningless, a threshold value greater than zero is usually selected. When the water volume is less than the threshold value, the state of the grid point is defined as dry.

Assuming that the state of a certain grid point is dry, the first step is to check how many of the surrounding grid points are wet grid points. If more than one is wet grid point, then the water level is averaged over these wet grid points. If the averaged water level is greater than the threshold water level, the state of this grid point may become wet, otherwise it still remains as a dry one. In the second step, it is necessary to further check the transport over the cross-section area between the grid point and the surrounding wet grid points. If these transport cross-section areas are all positive, then the dry grid point becomes a wet grid point and participates in the next calculation; otherwise, it is still a dry grid point.

#### 7.2.3.4. Grid types

The grids used in most sea level models are mainly divided into two categories: structured and unstructured grids. The structured grid nodes are arranged in an orderly manner, and the connection relationship between adjacent nodes is fixed. In contrast, the unstructured grid nodes are arranged in an

unordered manner and the adjacent nodes have no fixed connection relationship. Differently from structured grids, unstructured grids with triangular elements allow to adjust the resolution flexibly to depict complex shapes of coastline and estuary (Figure 7.10).

<figure><img src="image_295.png" /><figcaption>**Figure 7.10.** Structured grid and unstructured grid at coastal area.</figcaption></figure>

Structured grids mainly have two forms: rectangular grids and curvilinear grids. Relatively speaking, the curvilinear grid can adapt better to the complex shape of the coastline, and it can also realise the change of grid resolution that is more advantageous for storm surge simulation in estuary area (Figure 7.11).

<figure><img src="image_296.png" /><figcaption>**Figure 7.11.** An example of curvilinear grid.</figcaption></figure>

Triangular grids were the main forms of unstructured grids until a new form of unstructured grids, the SCVTs (Figure 7.12), appeared in ocean modelling a decade ago (Ringler et al., 2013). Like the triangular grids, the SCVTs can adapt well to complex coastline, and facilitate a smooth transition from coarse resolution grid cell to high resolution grid cell. In addition, they also solve the computational instability problem caused by small acute angles in the triangular grid. At present, the SCVTs model has been used in the China Sea (as shown in Figure 7.13).

<figure><img src="image_297.png" /></figure>
### 7.2.3.5. Discretization method

The main discretization methods include finite difference method, finite element method, and finite volume method. Discretization methods also correspond to the types of grids used. The finite difference method is generally used for structured grids, while the finite element method and finite volume method are generally used for unstructured grids.

The FDM is one of the simplest and oldest methods to solve the storm surge problems, and it is still widely used. This method divides the solution domain into differential grids, replacing the continuous solution domain with a finite number of grid nodes. By using the Taylor series expansion, the derivative of the governing equation is discretized by the difference quotient of the function value on the grid node, so as to establish the algebraic equations with the value on the grid node as the unknown quantity. This method is an approximate numerical solution that directly transforms a differential problem into an algebraic problem. The mathematical concept is intuitive and simple to express. It is an earlier and relatively mature numerical method.

<figure><img src="image_298.png" /><figcaption>**Figure 7.12.** An example of Voronoi tessellation schematic diagram: centroid in red, Voronoi circle in green, edges in grey.</figcaption></figure>

<figure><img src="image_299.png" /><figcaption>**Figure 7.13.** A storm surge model grid based on SCVTs applied to the China Sea.</figcaption></figure>

<figure><img src="image_300.png" /></figure>
<figure><img src="image_301.png" /><figcaption>**Figure 7.14.** Geometric interpretation of difference expression.</figcaption></figure>

($\frac{\partial u}{\partial x}$)$_i \approx \frac{u_{i+1}-u_i}{\Delta x}$ forward difference

($\frac{\partial u}{\partial x}$)$_i \approx \frac{u_i-u_{i-1}}{\Delta x}$ backward difference

($\frac{\partial u}{\partial x}$)$_i \approx \frac{u_{i+1}-u_{i-1}}{2\Delta x}$ central difference

The basic difference expression mainly has three forms (Figure 7.14): forward difference, backward difference, and central difference. The first two formats are first-order derivatives, while the last format is second-order derivatives. Different computational schemes can be obtained through the combination of several different formats of time and space.

According to the precision of the scheme, we can distinguish among first-order, second-order, and high-order accuracy schemes. Depending on the influence of the time factor, the difference scheme can also be divided into explicit scheme, implicit scheme, explicit and implicit alternate scheme.

The FEM is based on variational principle and weighted residual method. The basic solution idea is to divide the computational domain into a finite number of non-overlapping elements. In each element, some appropriate nodes are selected as interpolation points of the solution function. The variable in the differential equation is changed into a linear expression composed of the node value of each variable or its derivative and the selected interpolation function. The differential equation is solved discretely by means of variational principle or weighted residual method.

According to the difference of the weight function and the interpolation function, the finite element method is divided into several computational schemes. For the choice of weight function, there are collocation methods, moment method, least square method, and Galerkin method. According to the shape of the computing cell grid, there are triangular grid, quadrilateral grid, and polygonal grid. Triangular grids are commonly used in storm surge modelling. The accuracy of the interpolation function is divided into linear interpolation functions and high-order interpolation functions. Different combinations also constitute different finite element calculation schemes.

The FVM is also called the control volume method. The basic idea is to divide the computational domain into a series of non-overlapping control volumes, and make a control volume around each grid point. Integrate the differential equa-

tions to be solved for each control volume to obtain a set of discrete equations. The unknown is the value of the dependent variable at the grid point. In order to obtain the integral of the control volume, it is necessary to assume the changing law of the value between grid points, i.e. the distribution profile (continuous or segmented) of the assumed value. From the selection method of the integral region, the finite volume method belongs to the subregion method in the weighted residual method. From the approximate method of the unknown solution, the finite volume method is a discrete method using local approximation. The physical meaning of the discrete equation is the conservation principle of the dependent variable in a finite controlled volume, just as the differential equation expresses the conservation principle of the dependent variable in an infinitely small controlled volume. The discrete equation obtained by the finite volume method requires the integral conservation of the dependent variable to be satisfied for any set of control volumes, and naturally also for the entire computation area. This is the attractive advantage of the finite volume method. Some discrete methods, such as the finite difference method, only satisfy the integral conservation if the grid is extremely fine. The finite volume method shows accurate integral conservation even in the case of coarse grids. As far as the discrete method is concerned, the finite volume method can be regarded as an intermediate between the finite element method and the finite difference method.

### 7.2.3.6. Existing models for storm surge modelling

Numerical simulation of storm surges began in the 1950s. After decades of development, it has emerged that many numerical models can be used to storm surge simulations. Commercial models include MIKE21 and TuFlow, while examples of free models are: ADCIRC, Delft3D-FLOW, POM, FVCOM, ROMS, and SCHISM. Free numerical models generally provide the source code of the model, so that the model can be modified as needed when establishing a forecasting system. The models listed in Table 7.1 can be used to establish a complete operational storm surge forecasting system.

<figure><img src="image_302.png" /></figure>
<table><caption>**Table 7.1.** Geometric interpretation of difference expression.</caption><thead><tr><td>Model</td><td>Grid topology</td><td>Numerical methods</td><td>Nesting capabilities</td><td>Website</td></tr></thead><tbody><tr><td>MIKE21</td><td>Structured curvilinear grid and unstructured grid</td><td>Alternating direction implicit method for structured grid. Cell-centered finite volume method for unstructured grid</td><td>Nesting is not possible</td><td>https://www.mikepoweredbydhi.com/products/mike-21-3</td></tr><tr><td>TuFlow</td><td>Structured grid and unstructured grid</td><td>2nd order semi-implicit matrix solver for structured grid. Finite-Volume for unstructured grid</td><td>Sub-Grid</td><td>https://www.tuflow.com/</td></tr><tr><td>ADCIRC</td><td>Unstructured grid</td><td>Finite element method in space and finite difference method in time</td><td>Not available</td><td>https://adcirc.org/</td></tr><tr><td>Delft3D-FLOW</td><td>Structured curvilinear grid</td><td>Alternating direction implicit method</td><td>Nested boundary conditions</td><td>https://oss.deltares.nl/web/delft3d/</td></tr><tr><td>POM</td><td>Structured curvilinear grid</td><td>Finite difference scheme</td><td>Not available</td><td>http://www.ccpo.edu/POMWEB/index.html</td></tr><tr><td>FVCOM</td><td>Unstructured grid</td><td>Finite volume method</td><td>Nesting at the boundaries</td><td>http://codfish.smast.umass.edu/fvcom/</td></tr><tr><td>ROMS</td><td>Structured curvilinear grid</td><td>Second-order finite differences</td><td>One-way nesting</td><td>https://www.myroms.org/</td></tr><tr><td>SCHISM</td><td>Unstructured mixed triangular/ quadrangular grid</td><td>Semi-implicit Galerkin finite element method</td><td>One-way nesting</td><td>http://ccrm.vims.edu/schismweb/</td></tr></tbody></table>

### 7.2.4. Data assimilation systems

Data assimilation techniques are used to combine model and observed data to obtain the best estimate of the state of a system (see Chapter 4.3 for more details). Statistical techniques are often employed to find a solution which, ideally, minimises some error metric. For storm surges, this is done to obtain fields of sea surface height that can help us to better understand past events or to improve the quality of forecasts. An overview of the application of data assimilation to storm surge modelling and forecasting is provided in this section. Henceforth, references to errors mean some metric distance between the model and observations.

#### 7.2.4.1. Sources of error in storm surge models

In order to reduce errors in storm surge models, especially for forecasting in which the accuracy may influence real-time decision making, it is important to understand the sources of error. Some of the main sources are given below.

**Quality of input datasets**

This includes atmospheric surface forcing and tidal forcing at the boundaries (see Section 7.2.2). Storm surges are largely forced phenomena; therefore, the accuracy of forcing is key and errors in the related datasets may be transferred into the

<figure><img src="image_303.png" /></figure>
storm surge component of the modelled sea surface height. Errors in input datasets may arise from similar sources as the storm surge model, including model and instrument errors.

For example, errors in tidal amplitudes and phases at the boundaries will propagate with the tidal waves into the domain. As discussed in Section 7.2.1, interactions with the tides can influence both timing and height of a storm surge, so it is important to have an accurate tidal component in the model. The accuracy of the surface atmospheric forcing is important as well, especially the components of wind and surface pressure. Due to the forced nature of storm surges, this is one of the largest sources of error in storm surge forecasts (Horsburgh et al., 2011).

## Tuning of model parameters

It is common practice to adjust various model parameters to obtain a better solution. For example, it could include the tuning of bottom/surface friction coefficients. It is very unlikely to find a perfect parameter set, and the iterative processes often used can lead to non-optimal solutions. See Section 7.2 for more information about parameters used in storm surge models.

## Representativity errors

Representativity errors arise from models' ability to represent variables and processes such as resolution (Daley, 1991). For example, a coarse model may not be able to resolve finer scale features, which is present in the observations. For storm surges, this is particularly important nearby complex coastlines and estuaries. Similarly, coarser models may mean smoother bathymetry in these areas, which can significantly affect the modelled surge.

A model may not simulate all processes required to accurately model a storm surge. For example, if tides are not included in the model, only the atmospherically forced component of sea surface height is being generated, and contributions from tide-surge interactions will be missing. Other examples include the lack of tidal processes such as self-attraction and loading, or not including the inverse barometer effect in the model.

### 7.2.4.2. Assimilated data sources for storm surge modelling

An attempt to reduce the impact of the errors described in the previous section can be made using data assimilation. For storm surge modelling, assimilation of observations may occur directly into the model or indirectly via input datasets.

Datasets used as atmospheric forcing often contain assimilated observations. The generation of the storm surge

is highly dependent on the model's interaction with these datasets and it is vital that they are accurate. For example, the forecasted atmospheric fields used at the UK Met Office use assimilation of atmospheric observations. There are also many reanalysis datasets available, such as the ECMWF ERA5 dataset (Hersbach et al., 2020) that assimilates observations from multiple sources to generate atmospheric data. While these examples are suitable for extratropical storm surges, they may not sufficiently resolve intense tropical cyclones, meaning that parametric methods may be a better option (see Section 7.2.2.3). There are also assimilative alternatives, such as the MTCSWA datasets (Knaff et al., 2011), which blend together parametric fields and observations. These have been shown to have some benefit for forecasting (Byrne et al., 2017). The same is true for the datasets used to derive tidal signals at the model boundaries. Examples of such datasets include TPXO (Egbert et al., 2002) and FES (Lyard et al., 2021), which incorporate data from satellite altimetry and tide gauges. See Section 7.2.2.2 for more information on these tidal datasets.

Sea surface height may also be assimilated directly into the modelled sea surface. There are two sources used, both with advantages and disadvantages: tide gauges and satellite altimetry. Tide gauges (and other fixed instruments such as bottom pressure recorders) offer information that is frequent and consistent in time, making them useful for capturing ocean processes of all frequencies (including storm surges). However, they are generally spatially sparse. On the other hand, altimetry data offer good spatial coverage but are less consistent in time, as a satellite only returning to the same location once every number of days. This makes the data useful for longer periods of periodic ocean processes.

Tide gauge data are currently assimilated for storm surge forecasting at some institutions (see Section 7.2.4.4). There are representativity challenges that must be considered when using these data. Most importantly, modelled sea surface height variables and observed variables must represent the same physical quantity. For example, do both datasets contain the same components of sea level such as tides and inverse barometer effect? The datum on which the data are based must also be considered. The sea level anomaly can be used to overcome these problems if a long enough record is available.

### 7.2.4.3. Application of data assimilation to real time forecasting systems

For real time forecasting, data assimilation is used to generate an improved initial condition for a forecast model run. The forward propagation of errors can be reduced by creating a more realistic initial condition. This is important to improve the lead times over which good forecasts may be given. The use of data assimilation for storm surge forecasting has been shown to offer improvements over short lead

<figure><img src="image_304.png" /></figure>
<figure><img src="image_305.png" /><figcaption>**Figure 7.15.** Illustration of two distinct stages of sequential data assimilation for forecasting storm surges.</figcaption></figure>

times (Heemink, 1986; Verlaan et al., 2005; Madsen et al., 2015; Zijl et al., 2015; Byrne, 2021). However, the duration of these improvements may be limited to a few hours of forecast, due to the forced nature of storm surges.

This improved initial condition can be generated by running the model for some historical period up until today, including forcing with assimilated observations and potentially direct assimilation of sea surface height. This model run can effectively be seen as a continuous simulation with assimilative steps at some predefined frequency, for example every 6 hours. When a forecast is desired, the most recent state from this model can be taken and used as the initial condition. A forecast simulation is then done using no assimilation, as no observations are available in the future. This means that the atmospheric forcing used is also a forecast. Figure 7.15 illustrates this process.

There are several methods that have been used with success in storm surge modelling, including Optimal Interpolation (Gandin, 1966; Lorenc, 1981 and 1986; Daley, 1991), variational assimilation (Lorenc, 1986), Kalman filters (Kalman, 1960) and Ensemble Kalman Filters, (Evensen, 2004). In all methods, a key step is the estimation of spatial error covariances in both the model and the observations. This can be parametrically, as shown in the example in Figure 7.16, or by deriving covariances from an ensemble of model states. An example of the latter is described in Section 7.2.4.4.

Data assimilation has the potential to add significantly to the computation and time resources required, especially for en-

semble systems. As noted, for real-time forecasting systems is vital that a balance is made between accuracy and speed, i.e. useful forecasts need to be delivered in a timely manner (Horsburgh et al., 2011).

<figure><img src="image_306.png" /><figcaption>**Figure 7.16.** An example of correlation length scale estimation for assimilation of sea surface height into a barotropic storm surge model of the North Sea (Byrne et al. 2021). Such a length scale could be used to assimilate tide gauge observations.</figcaption></figure>

<figure><img src="image_307.png" /></figure>
### 7.2.4.4. Examples from real operational systems

The examples below are correct at the time of writing.

#### UKMO

The UK Met Office provides storm surge forecasting for the United Kingdom. Its 2D operational model does not currently assimilate any data into the model sea surface height. The atmospheric forcing used does include assimilated observations and comes from UKMO or the European Centre for Medium-Range Weather Forecasts (ECMWF) models, depending on the system. They also have more general operational 3D models that assimilate temperature, salinity and sea level anomaly, but the last is only done in deep water.

#### Rijkswaterstaat

Rijkswaterstaat provides storm surge forecasts for The Netherlands. Its system assimilates information from tide gauges around the Northwest European Shelf, especially in the North Sea (Verlaan et al., 2005; Zijl et al., 2015). They use a steady-state Kalman Filter (SSKF), which uses a stationary Kalman gain derived from an ensemble of states, such as might be used for the ensemble Kalman filter (EnKF). SSKF offers more computational efficiency than EnKF, and potentially a better representation of the error covariance than the standard Kalman Filter. Additional localization steps are also applied to the assimilation, to limit the distance from observations over which information is assimilated.

### 7.2.5. Ensemble modelling

Like any other forecasts, sea level predictions have an associated uncertainty. The threat to life and property of extreme sea level events makes estimation of this uncertainty and the generation of a range of possible water levels (probabilistic forecast) particularly important for risk managers and decision-makers. The error of a single forecast time series can be assessed by comparison with in-situ tide gauges at specific locations and grid points. However, uncertainty of the forecast and its dependence on the forcing, model characteristics, and set up is usually unknown.

Uncertainty of sea level forecasts depends on several factors and may contain errors in both the tide and the non-tidal residuals (storm surge) components. During a storm, the storm surge is mainly driven by the weather conditions at the sea surface. This is considered to be the dominant source of uncertainty in sea level forecasts and may change significantly depending on the meteorological conditions. For this reason, ensemble storm surge forecasts based on weather ensemble prediction systems (EPS) are the most common approach to generate probabilistic forecasts (forecast plus a confidence interval).

The weather is a chaotic system highly sensitive to the initial state (Lorenz, 1965) that can only be deterministically predicted for about 10 days. Therefore, the standard procedure for dealing with forecasts uncertainty, i.e. the combination of different model solutions or ensemble modelling, was initially applied to meteorological forecasts (Leith, 1974; Hamill et al., 2000). Conceptual background of ensembles is chaos theory; they are a valuable tool to deal with equations in which several nonlinear processes and interacting variables are present. This is the case of meteorological models but also of ocean models and, particularly, of tide and surge models. Hence, their application is today strongly recommended in oceanography.

An EPS is based on the combination of a set of forecasts with different controlled changes in the initial conditions, the model physics or the open boundary conditions (Palmer and Williams, 2010; Flowerdew et al., 2010). All these modifications are designed to represent the uncertainties in the knowledge of the weather state. For example, different initial conditions allow to include those perturbations that grow most rapidly in time, in a context where slight changes to the initial conditions may lead to significantly different forecasts (Buizza and Palmer, 1995). Slight modifications of the set of equations (including different values in the parameterization constants representing different processes) also provide estimation of model uncertainties contributing to the forecast error.

Deviation of wind and sea level pressure fields from their actual evolution will determine a corresponding deviation on the predicted sea level. Therefore, their uncertainty derived from weather EPS will cause an uncertainty in the evolution of the sea level, linked to the meteorological forcing and affecting mainly the storm surge component. If different weather predictions are used to drive different sea level simulations, the probability distribution function of the forecast sea level values allows estimating the uncertainty of the sea level forecast and the probability of exceeding a given sea level threshold.

The ECMWF EPS has been operational since 1992 (Molteni et al., 1996). It was first applied to storm surge operational forecasts in the North Sea by Flowerdew et al. (2010 and 2013), who provided skilled probabilistic forecasts of sea level and showed that ensemble spread was a reliable indicator of uncertainty during large surge events.

Storm surge ensemble predictions have been used to forecast sea level in Venice by Mel and Lionello (2014a). They used a 50 members ensemble to simulate 10 events showing that EPS slightly increases the accuracy of the prediction with respect to the deterministic forecast, and that the probability distribution of maximum sea level produced by the EPS is acceptably realistic. They also showed that the storm surge peaks correspond to maxima of uncertainty and that the uncertainty of such maxima increases linearly with the forecast range. The same procedure was used by Mel and Lionello (2014b) for the simulation

<figure><img src="image_308.png" /></figure>
of the operational forecast practice for a three-month peri-
od (fall 2010). It revealed that uncertainty for short and long
lead times of the forecast is mainly caused by the uncertain-
ty of the initial condition and of the meteorological forcing,
respectively. The probability forecast based on this ensemble
technique has a clear skill in predicting the actual probability
distribution of sea level. A computationally cheap alternative,
called ensemble dressing method, has been proposed by Mel
and Lionello (2016). It replaces the explicit computation of un-
certainty by ensemble forecast with an empirical estimate. In-
stead of performing multiple forecasts, the procedure "dress-
es" the forecast of sea level with an error distribution form,
which includes, on one hand, a dependence of the uncertainty
on surge level and lead time and, on the other hand, of the
uncertainty of the meteorological forcing. This computational-
ly cheap alternative also provides acceptably realistic results.

Apart from the meteorological input, other sources of error
on sea level forecasts can be attributed to the ocean model
characteristics and/or to the setup of the system: bathymet-
try, spatial resolution, model domain, tidal forcing, temporal
resolution of the meteorological input, barotropic or baro-
clinic models, ocean open boundary conditions, etc. Cur-
rently, sea level variations on timescales of hours/days are
operationally forecasted through different barotropic and
baroclinic models, sometimes over the same area. Therefore,
another option is the combination of existing operation-
al models with different characteristics, forcings and even
physics (multi-model forecast).

A multi-model storm surge forecast was first implemented by
Deltares (an independent Dutch institute) in 2008, combin-
ing existing operational storm surge forecasts from differ-
ent countries in the North Sea. The system included the use
of the Bayesian Model Average (BMA) statistical technique

for validation of the different members and generation of
a combined improved prediction, with a confidence inter-
val (Beckers et al., 2008). In the same year, this methodol-
ogy was tested for the Spanish coast by Puertos del Estado
(Spain) (Pérez et al., 2012), combining the output of Nivmar
(Álvarez-Fanjul et al., 2001), an existing storm surge fore-
casting system, with circulation (baroclinic) models already
operating in the region. Nowadays, at Puertos del Estado is
operational a multi-model surge forecast named ENSURF
that combines Nivmar with two Copernicus Marine Service
regional operational models, IBI-MFC (Sotillo et al., 2015) and
MedFS (Clementi et al., 2021).

The BMA technique requires near-real time access to tide
gauge data and automatic quality control of this data (as
required by the Nivmar system as well), and specific data
tailoring of model outputs. It is applied to the surge or
non-tidal residual component of sea level because this can
be approximated by a normal distribution (which is not the
case for total sea level including tides, especially for strong
semidiurnal regimes). So, observations from tide gauges and
model data for those models providing total sea level must
be previously decided. This could be considered a limitation
but, in practice, it is the best way of optimising the final to-
tal sea level forecast by using the tidal component obtained
from historical tide gauge observations at each site. ENSURF
is also a valuable operational validation tool that allows a
detailed assessment of the skills of different models to fore-
cast coastal sea levels. A first deterministic forecast is pro-
vided by the old Nivmar solution early in the morning every
day and, when later the Copernicus Marine Service forecasts
become available, they are integrated with the tide gauge
data and, by means of the BMA technique, a probabilistic
forecast band is generated for each harbour (Pérez-González
et al., 2017, Pérez-Gómez et al., 2019) (Figure 7.17).

<figure><img src="image_309.png" /><figcaption>**Figure 7.17.** Example of sea level probabilistic forecast generated by the multi-model ENSURF for the Barcelona harbour, validated against Barcelona tide gauge (hourly data). Top panel: total sea level; bottom panel: surge component. Blue: tide gauge data; orange: tide prediction; black: BMA forecast; grey: BMA confidence interval.</figcaption></figure>

<figure><img src="image_310.png" /></figure>
Generally, use of the probabilistic methodology improves the forecast and gives significant added value to existing operational systems, as there is no single model that outperforms at all tidal stations and synoptic conditions. However, further work must be done with the BMA technique to predict the storm peaks which, in some weather conditions, are better captured by single systems.

A multi-model ensemble forecasting system has been recently developed for the Adriatic Sea combining 10 models predicting sea level height (either storm surge or total water level) and 9 predicting waves characteristics (Ferrarin et al., 2020). Other examples of this technique can be found in New York (Di Liberto et al., 2011) and the North Sea (Siek and Solomatine, 2011).

### 7.2.6. Validation strategies

Storm surge models have been traditionally validated with time series of coastal sea level measured by tide gauges. These data allow assessing the skills of the model to reproduce observed water heights at specific points along the coast. Note the advantages this application presents with respect to sea level data from satellite altimetry, less reliable along the coastal strip and with a lower temporal resolution. Fortunately, there are hundreds of tide gauges around the world that become a very valuable and reliable dataset for storm surge validation (Muis et al., 2016 and 2020, Fernández-Montblanc et al., 2020). In some cases, these stations provide ancillary meteorological data, such as wind and atmospheric pressure, which can also be used to validate the model meteorological forcing. In addition, tide gauges transmit data in near-real time that can be integrated in an operational validation of the forecasting system (Álvarez-Fanjul et al., 2001).

In most cases, the forecast will provide the tide and storm surge signals (hourly to daily timescales), usually dominant at the tide gauge records, but will not be able to reproduce higher-frequency sea level oscillations such as infragravity waves, seiches or meteotsunamis, with periods of the order of a few minutes. It is important to define which observed "sea level" product will best fit the validation purpose, according to the physical processes included in the system. The most adequate standard product for existing operational storm surge forecasting systems are filtered hourly values from tide gauges. As new models include additional processes (e.g. fully coupled models including wave effects), higher resolution bathymetries, and forcing fields, the use of lower temporal sampling data will become more important and the validation process more challenging.

Normally, the model output at the grid point closer to the tide gauge is selected. However, the validation results will

depend on the resolution and quality of the bathymetry data, as well as on the location of the tide gauge: if it is in an open site or inside a harbour or bay with important local effects, it may not be resolved by the model.

A careful validation of both tide and surge components should be performed to verify not only the final total sea level forecast, but also the quality of the tidal signal in the model, which can be an important source of error especially on shallow waters with high tidal range. For these reasons, model and tide gauge data must be de-tided, applying a harmonic analysis to both time series, as well as computing the tide prediction for the analysed period, and the surge or non-tidal residual after tide subtraction. The performance assessment can then be made in terms of comparison of the harmonic constants (amplitude and phase) from model and observations, and in terms of model data time series comparison of tide, surge and total sea level.

It is important to mention that sea levels measured by the tide gauge will be related to a local, regional or national datum. The model forecast is theoretically referred to mean sea level, though this mean sea level may be affected by the model setup implementation, domain, etc. Therefore, the mean should be subtracted from both time series at each location before comparison.

Metrics for time series validation can be found at Section 4.5.1. Most of these metrics describe the overall performance of the model for a time period of several days (the storm duration), months or years. However, they do not reflect the predictive skill for extreme surge events. This can be better evaluated, for example, in terms of the differences in the highest percentiles (e.g.: 95th, 99th percentiles) or the maximum observed and modelled value. For validation of long time series (multi-decadal hindcasts), it is possible to use annual maxima, annual percentiles, and extreme sea levels for specific return periods, obtained through extreme sea level analysis (Muis et al., 2020).

Taylor diagrams can be used to graphically indicate the performance of different competing models or solutions, providing information of the Pearson correlation coefficient, the standard deviation and the RMSE at each tide gauge.

Progressively, the storm surge models will consider inundation, and additional validation of the flooding extent will be required. This is less straightforward and requires other types of data, such as locations of flooded points, marks left by the water or reports about the flood chronology (Le Roy et al., 2015). This information is commonly available after the event for a delayed mode validation; e.g. for validation of inundation, Loftis et al. (2017) used crowdsourced GPS data and maps of flooded areas obtained by drones.

<figure><img src="image_311.png" /></figure>
### 7.2.7. Outputs

The main outputs of storm surge models are: time series output, maximum elevation field output (extreme values at every time step for water surface elevation), ensemble forecast elevation field output, animation output.

#### 7.2.7.1. Time series outputs

The time series output is usually plotted in a two-dimensional rectangular coordinate system, the abscissa is time and the ordinate is water level. The time series output of storm surge models are the water level changes at a certain location. In order to facilitate the comparison between the

results of simulations and the observation data, multiple result curves can be plotted in the same coordinate system.

Generally, the results of storm surge models (without tide) can be directly used for plotting time series diagrams, but sometimes attention should be paid to the change of the total water level at a certain point. Therefore, there is a way to output the total water level, that is the results of the storm surge model without astronomical tide directly superimposed on the astronomical tide from harmonic analysis, obtaining in this way the time series results of the total water level. Figures 7.18 and 7.19 show examples of time series model storm surge result and storm surge superimpose on harmonic analysis tide.

<figure><img src="image_312.png" /><figcaption>**Figure 7.18.** Time series model storm surge result (blue line) and observed storm surge (red circle).</figcaption></figure>

<figure><img src="image_313.png" /><figcaption>**Figure 7.19.** Time series model storm surge result superimposed on predicted tide (blue line) and observed water level (red circle).</figcaption></figure>

<figure><img src="image_314.png" /></figure>
### 7.2.7.2. Maximum elevation field

Among the output methods of storm surge models, there is an output form called maximum elevation field. This kind of field output is not the water level field at a certain time, but extracts the highest water level value of each grid point in the calculation process to form the maximum elevation field. The maximum water level field can be used to grasp the distribution of the maximum water level during a Typhoon and to identify the more severely affected areas along the coast. Figure 7.20 shows the maximum storm surge field during 2019 Typhoon Mitag (1918).

<figure><img src="image_315.png" /><figcaption>**Figure 7.20.** The maximum storm surge field during 2019 Typhoon Mitag (1918).</figcaption></figure>

### 7.2.7.3. Ensemble forecast field

The storm surge ensemble forecast usually uses the respective meteorological forcing fields of the ensemble members to calculate the storm surge separately. The output of ensemble forecast fields mainly includes the following forms: (a) ensemble mean field; (b) probability field; and (c) postage stamp maps.

#### a. Ensemble mean storm surge field

In order to obtain a definite forecast result, it is necessary to synthesise the respective results of the ensemble members. It is generally used to assign different weights to the results of each member and to superimpose the results of all members, i.e. the weighted average method (Wang et al., 2010). The superimposed result is output in the form of the elevation field, and the ensemble forecast field is obtained as result.

The track map of Typhoon Mitag can be seen in Figure 7.21. The storm surge results of the subjective typhoon forecast track, fast track, slow track, left track, and right track are used to synthesise the ensemble forecast water level field applying the weighted average method (Figure 7.22). In this example, the weight of the storm surge result of the subjective forecast typhoon track is 60%, while the weight of the storm surge result of the other tracks are all 10%.

#### b. Probability storm surge field

The typhoon ensemble forecasting tracks for storm surge numerical simulation can describe the surge field

<figure><img src="image_316.png" /><figcaption>**Figure 7.21.** Track map of Typhoon Mitag. Red line: middle track; black line: fast track; cyan line: slow track; magenta line: left track; blue line: right track.</figcaption></figure>

<figure><img src="image_317.png" /></figure>
<figure><img src="image_318.png" /><figcaption>**Figure 7.22.** Synthesis of ensemble forecast water level field.</figcaption></figure>

under different typhoon track scenarios. By equitably assigning weights to different track results, the probability field distributions under different extreme values of storm surges can be clearly displayed, and the intensity probability of coastal storm surges can be more intuitively presented (Liu et al., 2020). Figures from 7.23 to 7.25 show the probability distribution of forecast storm surge over 0.5m, 1.0m and 2.0m of Typhoon Mitag.

### c. Postage stamp maps

A postage stamp map is a set of small storm surge maps drawn by the results of the individual members (Figure 7.26). Forecasters can learn about the possible situation of each ensemble member through the postage stamp map, thereby estimating the magnitude of the maximum storm surge and the range of impact. (WMO, 2012)

<figure><img src="image_319.png" /><figcaption>**Figure 7.24.** Distribution of probability forecasting of storm surge over 1.0m.</figcaption></figure>

<figure><img src="image_320.png" /><figcaption>**Figure 7.23.** Distribution of probability forecasting of storm surge over 0.5m.</figcaption></figure>

<figure><img src="image_321.png" /><figcaption>**Figure 7.25.** Distribution of probability forecasting of storm surge over 2.0m.</figcaption></figure>

<figure><img src="image_322.png" /></figure>
<figure><img src="image_323.png" /><figcaption>**Figure 7.26.** Postage stamp map of storm surge field of Typhoon Mangkhut.</figcaption></figure>

#### 7.2.7.4. Animation output

The output of the water level field is very helpful for grasping the distribution of the storm surge process over a whole region. The output of the storm surge water level field is the elevation value of all grid points at a certain time. The elevation field figure at each moment is taken as a frame of the animation, and all the frames are connected to form the elevation field animation (Figure 7.27). The elevation field animation can intuitively reflect the changes in the water level of the entire area during the impact of the storm surge.

<figure><img src="image_324.png" /><figcaption>**Figure 7.27.** Synthesis of elevation field animation.</figcaption></figure>

### 7.2.8. Existing operational storm surge forecasting systems

After decades of development of storm surge numerical models, many countries have established their own operational storm surge forecasting models. For example, in the United States a storm surge forecasting system is operating through the SLOSH model, in which the wind field is established based on the cyclone path, maximum wind speed radius, storm centre, and environmental pressure difference; it provides operational forecast products and storm surge inundation guidance products (Jarvinen and Lawrence, 1985). China established ver3.0 of the PMOST forecasting system, which is based on a depth-averaged two-dimensional shallow water equation in the vector invariant form, and uses a SCVTs grid. It can enhance the resolution in key areas and fit the coastline. The system is able to couple astronomical tides and simulate flooding processes. With the GPU acceleration technology, the efficiency of storm surge simulation along the coast of China can reach 60 sec/day. The system can also perform ensemble forecasts based on multiple typhoon events and storm surge probability forecasting. The Indian Institute of Technology storm surge model was developed in the 1980s and applied to storm surge forecasts in the Indian Ocean and the Arabian Sea. It uses rectangular Cartesian coordinates and separates land and water during calculations. It has been applied throughout the north Indian Ocean (Dube et al., 1984, Dube et al., 1985).

<figure><img src="image_325.png" /></figure>
The two-dimensional storm surge model developed by the Japan Meteorological Agency also uses rectangular coordinates (Hasegawa et al., 2015). In the numerical calculation, the water and land are separated with the flexible mesh, i.e. fine grid is used in shallow water and coarse grid is used in deep water. The system can also provide ensemble forecasting products. In the mid-1980s, the Netherlands developed a numerical fluid dynamics model called the DCSM, which uses a depth integrated shallow water equation. The driving force of the model is provided by a high-resolution regional meteorological model. In the early 1990s, the Kalman filtering method was used in DCSM to assimilate the water level (Verlaan et al., 2005; de Vries, 2009). The UKMO developed the storm surge forecasting model CS3X, which is a tide-storm

surge model. In this operational storm surge forecasting system, the interaction of tide and storm surge is considered. In recent years, a storm surge ensemble forecasting has been developed in this system (Flowerdew et al., 2013). A storm surge model covering the French overseas territories has been operated since the 1990s by Meteo-France (Daniel et al., 2009); it was established based on the spherical nonlinear shallow water equation. In order to solve the problem of the shore boundary, the C-grid difference format was adopted, with meteorological forcing provided by the Holland model (see Section 7.2.2.3). Table 7.2 provides a list and features of storm surge forecasting systems currently operating in various countries.

<table><caption>**Table 7.2.** List and features of operational storm surge forecast models.</caption><thead><tr><td>Model</td><td>Area</td><td>Type</td><td>Grid</td><td>Country</td></tr></thead><tbody><tr><td>HAMSOM, Nivmar</td><td>Mediterranean Sea and Iberian Peninsula</td><td>Vertically integrated barotropic</td><td>10 minutes</td><td>Spain</td></tr><tr><td>Mike 21 pre-operational 3-D 2-D finite element MOG2D</td><td>North Sea, Baltic Sea</td><td>2-D hydrodynamic</td><td>Finite difference 9 nmi, 3 nmi, 1 nmi, 1/3 nmi</td><td>Denmark</td></tr><tr><td>Coupled ice-ocean NPAC</td><td>Grand Banks, Newfoundland, Labrador<br>NE Pacific, 120°W-160°W, 40°N-62°N</td><td>3-D circulation based on the Princeton Ocean Model</td><td>Approximately 20 km x 20 km<br>Finite difference curvilinear C-grid 1/8 degree</td><td>Canada</td></tr><tr><td>Mike 21 pre-operational 3-D 2-D finite element MOG2D</td><td>North Sea, Baltic Sea</td><td>2-D hydrodynamic</td><td>Finite difference 9 nmi, 3 nmi, 1 nmi, 1/3 nmi</td><td>Denmark</td></tr><tr><td>JMA Storm Surge</td><td>23.5°N-46.5°N, 122.5°E-146.5°E</td><td>2-D linearized shallow water</td><td>Staggered Arakawa C-grid. 1 minute latitude/longitude</td><td>Japan</td></tr><tr><td>KMA Storm Surge</td><td>20°N-50°N, 115°E-150°E</td><td>2-D barotropic surge and tidal current based on the Princeton Ocean Model</td><td>Approximately 8 km x 8 km, finite difference curvilinear C-grid 1/12 degree</td><td>Republic of Korea</td></tr></tbody></table>

<figure><img src="image_326.png" /></figure>
<table><thead><tr><td>Model</td><td>Area</td><td>Type</td><td>Grid</td><td>Country</td></tr></thead><tbody><tr><td>NIVELMAR</td><td>Portuguese mainland coastal</td><td>Shallow water</td><td>1 minute latitude x 1 minute longitude</td><td>Portugal</td></tr><tr><td>SMARA storm surge</td><td>Shelf sea 32°S-55°S, 51°W-70°W. Rio de la Plata</td><td>2-D depth-averaged</td><td>Geographical Arakawa C-grid, 1/3 degree latitude x 1/3 degree longitude<br>1/20 degree latitude x 1/20 degree longitude</td><td>Argentina</td></tr><tr><td>BSH circulation (BSHcmod)<br>BSH surge (BSHsmod)</td><td>North-east Atlantic, North Sea, Baltic</td><td>3-D hydrostatic circulation<br>2-D barotropic surge</td><td>Regional spherical, North Sea, Baltic 6 nmi, German Bight Western Baltic, 1 nmi, surge North Sea, 6 nmi, north-east Atlantic 24 nmi</td><td>Germany</td></tr><tr><td>Caspian Storm Surge</td><td>Caspian Sea 36°N-48.5°N, 45°E-58°E<br>North Caspian Sea 44.2°N-48°N, 46.5°E-55.1°E</td><td>2-D hydrodynamic, based on MIKE 21 (DHI Water & Environment)</td><td>10 km x 2 km</td><td>Kazakhstan</td></tr><tr><td>HIROMB/NOAA</td><td>North-east Atlantic, Baltic</td><td>3-D baroclinic</td><td>C-grid, 24 nmi</td><td>Sweden</td></tr><tr><td>WAQUA-in-Simona/DCSM98</td><td>Continental shelf 48°N-62°N, 12°E-13°E</td><td>2-D shallow water, ADI method, Kalman filter data assimilation</td><td>1/8 degree longitude x 1/12 degree latitude</td><td>Netherlands</td></tr><tr><td>Derived from MOTHY oil spill drifts model</td><td>Near-Europe Atlantic (Bay of Biscay, Channel and North Sea)<br>8.5°E-10°E, 43°N-59°N<br>West Mediterranean basin (from the Strait of Gibraltar to Sicily)<br>Restricted area in overseas departments and territories</td><td>Shallow-water equations</td><td>Arakawa C-grid 5' of latitude x 5' of longitude Finer meshes</td><td>France</td></tr><tr><td>SLOSH</td><td>Sea area south of Hong Kong within 130 km</td><td>Finite difference</td><td>Polar, 1 km near to 7 km, South China Sea</td><td>Hong Kong, China</td></tr><tr><td>Short-term sea-level and current forecast</td><td>Caspian Sea and near-shore low-lying zones</td><td>3-D hydrodynamic baroclinic</td><td>3 nmi horizontal, 19 levels</td><td>Russian Federation</td></tr></tbody></table>

<figure><img src="image_327.png" /></figure>
<table><thead><tr><td>Model</td><td>Area</td><td>Type</td><td>Grid</td><td>Country</td></tr></thead><tbody><tr><td>IIT-Delhi, IIT-Chennai, NIOT-Chennai</td><td>East and west coasts of India and high-resolution areas</td><td>Non-linear, finite element, explicit finite element</td><td>For inundation model average spacing of 12.8 km offshore direction and 18.42 km along shore</td><td>India</td></tr><tr><td>CS3 tide surge</td><td>North-west European shelf waters</td><td>Finite difference, vertically averaged</td><td>C-grid 12 km, nested finer resolution</td><td>United Kingdom</td></tr><tr><td>SLOSH</td><td>Atlantic and gulf coasts</td><td>2-D depth integrated</td><td>625 meters</td><td>United States</td></tr><tr><td>PMOST 3.0</td><td>China sea<br>10°N-45°N, 105°E-140°E</td><td>2-D depth-averaged barotropic</td><td>SCVTs unstructured grid, 10km at boundary and 500m along shoreline</td><td>People's Republic of China</td></tr></tbody></table>

<figure><img src="image_328.png" /></figure>

## 7.3. References

* Álvarez-Fanjul, E., Perez Gomez, B., Sanchez Arevalo, I.R. (2001). Nivmar: A storm surge forecasting system for Spanish Waters. *Scientia Marina*, 65(S1), https://doi.org/10.3989/scimar.2001.65s1145

* Beckers J. V. L., Sprokkereef, E., and Roscoe, K. L. (2008). Use of Bayesian model averaging to determine uncertainties in river discharge and water level forecasts, Proceeding of “4th International Symposium on Flood Defence: Managing Flood Risk, Reliability and Vulnerability”, Toronto, Ontario, Canada, 6–8 May, 2008.

* Bjerknes, V. (1921). On the dynamics of the circular vortex with applications to the atmosphere and atmospheric vortex and wave motions. *Quarterly Journal of the Royal Meteorological Society*, 48(204), 375-376. https://doi.org/10.1002/qj.49704820414

* Buizza, R., and Palmer, T.N. (1995). The Singular-Vector Structure of the Atmospheric Global Circulation. *Journal of the Atmospheric Sciences*, 52(9), 1434-1456, https://doi.org/10.1175/1520-0469(1995)052<1434:TSVSOT>2.0.CO;2

* Byrne, N. J. (2017). Deterministic models of Southern Hemisphere circulation variability. Ph.D. thesis, University of Reading, 105 pp.

* Byrne, D., Horsburgh, K., and Williams, J. (2021). Variational data assimilation of sea surface height into a regional storm surge model: Benefits and limitations. *Journal of Operational Oceanography*, 1-14, https://doi.org/10.1080/1755876X.2021.1884405

* Clementi, E., Aydogdu, A., Goglio, A. C., Pistoia, J., Escudier, R., Drudi, M., Grandi, A., Mariani, A., Lyubartsev, V., Lecci, R., Cretí, S., Coppini, G., Masina, S., and Pinardi, N. (2021). Mediterranean Sea Physical Analysis and Forecast (CMEMS MED-Currents, EAS6 system) (Version 1) set. Copernicus Monitoring Environment Marine Service (CMEMS), https://doi.org/10.25423/CMCC/medsea_analysisforecast PHY_006_013_EAS6

<figure><img src="image_329.png" /></figure>
Daley, R. (1991). Atmospheric Data Analysis. Cambridge University Press: Cambridge, UK

Daniel, P., Haie, B., and Aubail, X. (2009). Operational Forecasting of Tropical Cyclones Storm Surges at Meteo-France. *Marine Geodesy*, 32(2), 233-242. https://doi.org/10.1080/01490410902869649

de Vries, H. (2009). Probability Forecasts for Water Levels at the Coast of The Netherlands. Marine Geodesy, 32(2), 100-107, https://doi.org/10.1080/01490410902869185

Di Liberto, T., and Colle, B.A. (2011). Verification of a Multimodel Storm Surge Ensemble around New York City and Long Island for the Cool Season. *Weather and Forecasting*, 26, 922-939.

Doodson, A. T. (1929). Report on Thames Floods. *Geophysical Memoirs*, Great Britain Meteorological Office, 5(47), 1-26.

Doodson, A. T. (1956). Tides and storm surges in a long uniform gulf. *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences*, 237(1210), 325-343.

Dube, S. K., Sinha, P. C., Rao, A. D., and Rao, G. S. (1985). Numerical modelling of storm surges in the Arabian Sea. *Applied Mathematical Modelling*, 9(4), 289-294, https://doi.org/10.1016/0307-904X(85)90067-8

Dube, S. K., Sinha, P. C., and Roy, G. D. (1984). Numerical simulation of storm surges induced by tropical storms impinging on the Bangladesh coast. *Coastal Engineering Proceedings*, 1(19), https://doi.org/10.9753/icce.v19.13

Egbert, G. D., Bennett, A. F., and Foreman, M. G. G. (1994). TOPEX/POSEIDON tides estimated using a global inverse model. *Journal of Geophysical Research: Oceans*, 99(C12), 24821-24852, https://doi.org/10.1029/94JC01894

Egbert, G. D., and Erofeeva, S. Y. (2002). Efficient Inverse Modeling of Barotropic Ocean Tides. *Journal of Atmospheric and Oceanic Technology*, 19(2), 183-204, https://doi.org/10.1175/1520-0426(2002)019<0183:Eimobo>2.0.Co;2

Evensem, G. (2004). Sampling strategies and square root analysis schemes for the EnKF. *Ocean Dynamics*, 54, 539-560, https://doi.org/10.1007/s10236-004-0099-2

Fernández-Montblanc, T., Vousdoukas, M.I., Mentaschi, L., Ciavola, P. (2020). A Pan-European high resolution storm surge hindcast. *Environment International*, 135, 105367, https://doi.org/10.1016/j.envint.2019.105367

Ferrarin, C., Valentini, A., Vodopivec, M., Karic, D., Massaro, G., Bajo, M., De Pascalis, F., Fadini, A., Ghezzo, M., Menegon, S., Bressan, L., Unguendoli, S., Fettich, A., Jerman, J., Licer, M., Fustar, L., Papa, A., Carraro, E. (2020). Integrated sea storm management strategy: the 29 October 2018 event in the Adriatic Sea. *Natural Hazards and Earth System Sciences*, 20, 73-93, https://doi.org/10.5194/nhess-20-73-2020

Flowerdew, J., Horsburgh, K., Wilson, C., and Mylne, K. (2010). Development and evaluation of an ensemble forecasting system for coastal storm surges. *Quarterly Journal of the Royal Meteorological Society*, 136(651), 1444-1456, https://doi.org/10.1002/qj.648

Flowerdew, J., Mylne, K., Jones, C., and Titley, H. (2013). Extending the forecast range of the UK storm surge ensemble. *Quarterly Journal of the Royal Meteorological Society*, 139(670), 184-197, https://doi.org/10.1002/qj.1950

Fujita, T. (1952). Pressure distribution within typhoon. *Geophysical Magazine*, 23, 437-451.

<figure><img src="image_330.png" /></figure>
Gandin, L.S. (1963). Objective analysis of meteorological fields. Translated from Russian.

Hamill, T. M., Snyder, C., and Morss, R.E. (2000). A comparison of probabilistic forecasts from bred, singular vector, and perturbed observation ensembles. *Monthly Weather Review*, 128(6), 1835-1851, https://doi.org/10.1175/1520-0493(2000)128<1835:ACOPFF>2.0.CO;2

Harris, D.L. (1957). The Effect of a Moving Pressure Disturbance on the Water Level in a Lake. In: "Interaction of Sea and Atmosphere". Meteorological Monographs, vol 2. American Meteorological Society, Boston, MA. https://doi.org/10.1007/978-1-940033-15-0_4

Hasegawa, H., Kohno, N., and Itoh, M. (2015). Development of Storm Surge Model in Japan Meteorological Agency. Available at http://www.waveworkshop.org/14thWaves/Papers/JCOMM_2015_j4.pdf

Heemink, A.W. (1986). Storm surge prediction using Kalman filtering. Dissertation TU Twente, Hydraulic Engineering Reports. Available at http://resolver.tudelft.nl/uuid:8a7c65f5-f8f5-4033-b796-adaef2eb39a7

Hersbach, H. Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R.J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., Thépaut Hersbach, J.N. (2020). The ERA5 global reanalysis. *Quarterly Journal of the Royal Meteorology Society*, 146(730), 1999-2049, https://doi.org/10.1002/qj.3803

Holland, G. J. (1980). An Analytic Model of the Wind and Pressure Profiles in Hurricanes. Monthly Weather Review, 108(8), 1212-1218. https://doi.org/10.1175/1520-0493(1980)108<1212:Aamotw>2.0.Co;2

Horsburgh, K., and de Vries, H. (2011). Guide to Storm Surge Forecasting. WMO-No. 1076. ISBN 978-92-63-11076-3

Jarvinen, B. R., and Lawrence, M. B. (1985). An evaluation of the SLOSH storm-surge model. *Bulletin of the American Meteorological Society*, 66(11), 1408-1411.

Jelesnianski, C. P. (1965). A numerical calculation of storm tides induced by a tropical storm impinging on a continental shelf. *Monthly Weather Review*, 93(6), 343-358. https://doi.org/10.1175/1520-0493(1993)093<0343:Ancos>2.3.Co;2

Kalman, R. E. (1960). A New Approach to Linear Filtering and Prediction Problems. *Journal of Basic Engineering*, 82(1), 35-45. https://doi.org/10.1115/1.3662552

Knaff, J.A., De Maria, M., and Molenar, D.A. (2011). An automated, objective, multiple satellite platform tropical cyclone surface wind analysis. *Journal of Applied Meteorology and Climatology*, 50(10), pp. 2149-2166, https://doi.org/10.1175/2011JAMC2673.1

Le Roy, S., Pedreros, R., André, C., Paris, F., Lecacheux, S., Marche, F., and Vinchon, C. (2015). Coastal flooding of urban areas by overtopping: dynamic modelling application to the Johanna storm (2008) in Gâvres (France). *Natural Hazards and Earth System Sciences*, 15, 2497-2510, https://doi.org/10.5194/nhess-15-2497-2015

Leith, C.E. (1974). Theoretical Skill of Monte Carlo Forecasts. Monthly Weather Review, 102 (6), 409-418, https://doi.org/10.1175/1520-0493(1974)102<0409:TSOMCF>2.0.CO;2

Liu, Q., Jiang, J., Yu, F., Zhang, C., Dong, J., Song, X., and Wang, Y. (2020). Typhoon storm surge ensemble forecast based on GPU technique. Acta Oceanologica Sinica, 39(5), 77-86, https://doi.org/10.1007/s13131-020-1570-8

<figure><img src="image_331.png" /></figure>