In its original incarnation CEST (Zhang et al., 2008, 2012) handled advection via a straightforward, fully explicit, Eulerian finite difference scheme. This approach often leads to a prohibitive restriction on the size of the time–step that the original CEST model can employ. This is because, for numerical stability, the time–step used for the entire model must be chosen such that advection satisfies the well-known CFL condition (Courant et al., 1967). Here, in the spirit of Casulli's original approach (Casulli, 1990), we employ a semi–Lagrangian (SL) methodology for the velocity advection. Importantly, however, we extend the approach to high–order accuracy in both space and time by employing second–order Runge Kutta time integration and monotonic cubic spline interpolation in space. In theory, this type of advection is unconditionally stable. When computing the velocity advection, we work purely on the computational (image) grid; working on the curvilinear (physical) grid adds complexity as scales vary arbitrarily between cells and grid curvature is not necessarily constant. In order to effect the interpolation, and limiting, the two–dimensional processes are broken down into a sequence of one–dimensional processes along each co–ordinate axis. This is possible because the computational (image) grid is a regular Cartesian grid. Thus, without loss of generality, when discussing the base point interpolation we need only consider the 1D case. For the spatial interpolation we use Hermite cubic interpolation made monotone by use of the limiter proposed by Nair et al. (1999); from hereonin we shall refer to this as the NCS99 limiter. The NCS99 limiter is applied in each spatial dimensional in turn and works as follows: first find the local maximum and minimum surrounding the particle path base point x_b

Next, reset the interpolated value (obtained using Hermite cubic interpolation) at the base point f|_{x =xb} = f_b such that

Note that this amounts to a simple clamping operation and can lead to the suppression of certain genuine physical waves. For this reason the NCS99 limiter includes additional checks for monotonicity that must be satisfied before the limiter is applied. First NCS99 introduce the global minimum and global maximum of f denoted by f_min and f_max, respectively. Thus, if all of the following four inequalities hold the limiter should not be applied

This stipulates that the signal contains only one extremum in a five–mesh–length interval (therefore suppressing the Gibbs phenomenon). We have found this limiter to be particularly robust when compared with alternative formulations such as that proposed in Fedkiw et al. (2001). Figure 3 shows a comparison of SL advection, using a variety of base point interpolation schemes, for the advection of a 1D function with compound waves in a uniform velocity field. The 1D function comprises a combined Gaussian, triangle and square wave. Results are plotted after 100 time–steps for a grid comprising 500 points; clearly Hermite cubic interpolation with the NCS99 limiter gives the best performance in this complex case. Whilst this high–order accurate advection scheme is monotone, it is unsuitable for flows that contain or develop discontinuities as it is not conservative. Moreover, CEST is only suitable for smooth flows as the governing equations are themselves not in a (divergence) form that permits discontinuities. Thus, unless explicit shock fitting is utilized (Moretti, 2002), flow discontinuities cannot be expected to propagate at the correct strength or speed. We mention that, when using SL advection, in wet areas close to dry land, care should be taken to ensure that the particle path is not allowed to project too far back into an area that is completely dry. If this is allowed to happen then the velocity advection can cause a false zeroing of the velocity field in such cells. This issue can be avoided through the use of locally controlled time–stepping.

The physical diffusion terms are treated using a simple forward–in–time centered–in–space (FTCS) scheme (Press et al., 1992). We note that the use of a simple explicit scheme for diffusion introduces its own stability requirements; however, these are far less stringent than those associated with the advective terms. The Coriolis and wind stress terms are updated using first–order explicit Euler time integration and the bottom friction term is treated implicitly in the manner detailed in Kelly et al. (2015). After the velocities have been updated using the method of fractional steps the free surface can be updated. The operator that updates both u and v by SL advection and explicit diffusion, evaluation of curvature terms, pressure drop, wind forcing and finally an implicit bed friction update is denoted by F. Thus, we have

Evolution of the free surface requires solution of Equation (5). This is achieved using an implicit scheme which results in the following 5–diagonal linear system of equations:

where

and

The right hand side of Equation (13) is given by

Finally, after the free surface has been updated, the final velocity update is performed using the pressure gradient. Use is made of the staggered grid to obtain second–order accuracy for the spatial gradient of the free surface in this step.

CEST employs a straightforward wetting–drying algorithm that is based on an accumulated water volume (Zhang et al., 2008). Free surface elevation and water depth at both the cell center and its four boundaries are all used to calculate the accumulated water volume. At the beginning of a model timestep, cells are assigned as being either wet or dry based on whether the water depth at the cell center is above or below a threshold depth H_TOL. During wetting, if the free surface elevation at the center of a wet cell is higher than that at an adjacent dry cell, and the water depth at the shared boundary between these two cells H_k (obtained by linear interpolation) is greater than a predefined threshold, the water is allowed to flow from the wet cell into the dry cell and accumulate there. The flux of water crosses a maximum of four shared boundaries between a dry cell, and any wet neighbors. The water interchange velocities (u_k, where k = 1...4 represents the four cell boundaries) are approximated by solving a simplified 1D momentum equation which ignores the contribution of advection, Coriolis force, air pressure drop and wind shear giving

with x_k being the direction of u_k and ζ_k is a linear reconstruction of the free surface elevation at the kth cell interface. The accumulated water volume ΔQ in dry cells is computed as

where ΔQ(i,j)n is the accumulated volume from the previous time step and Δ_k denotes the cell length (which is Δξ or Δη depending on the value of k). If the water depth, obtained from the accumulated volume, in a dry cell exceeds H_TOL the cell is reflagged as being wet. During drying a cell is set to be dry if the water depth at the cell center falls below H_TOL. Note that, if the water depth at a cell boundary is less than H_TOL, water will stop flowing across this boundary even before the cell itself is completely dry. Also, if it is the case that the linearly reconstructed water depths at all four boundaries of a cell are less than H_TOL then the cell is set to be dry. This simple algorithm conserves water mass, but not momentum. The approach has proven extremely robust in a huge number of storm surge simulations carried out at the IHRC over the last decade.

Parallelization of the CEST model is achieved using the Message Passing Interface (MPI) library, see: https://www.open-mpi.org/. To achieve parallelization a horizontal strip–type domain decomposition on the computational (image) grid is employed. The computational domain is split into a number of horizontal strips and strips are allocated to each process on a single strip per process basis. Two layers of halo regions (ghost cells) is employed to transfer information between processes. It should be noted that no load balancing, or process optimization, is currently implemented.

Figure 4 shows a schematic representation of the domain decomposition employed by CEST. Assuming a domain that runs from y_min to y_max and comprises im × jm cells, only J index is split to save on computational time and simplify the domain decomposition. In other words, the computational domain is split vertically among all available processes np, including the root process p=0, according to Algorithm 1. The horizontal direction, I index for each process keeps same.

Whilst the parallelization of the F operator is straightforwards, solution of the implicit equation for the free surface (Equation 13) is not. To facilitate the paralellization of the code, and avoid the need for multicoloring, the pre–conditioned conjugate gradient (PCCG) method employed to solve the continuity equation for the free surface (Equation 13) in CEST is replaced with simple Jacobi iteration. This allows for a straightforward parallelization of the implicit part of Casulli's algorithm. Whilst this approach leads to slower convergence this is more than offset by the ability to employ multiple processes (and the associated decrease in the linear system size that each process is required to solve).

Figure 5 shows the model domain and the measurement locations for Hurricane Irma (2017), including tide gauges from NOAA, High Water Marks (HWMs) from US Geological Survey (USGS), and river gauges from USGS. We employ a very high–resolution curvilinear grid that comprises approximately 640,000 computational cells. The cell size at the open ocean area is approximately 1,470 × 1,470 m, which gradually reduces to 200 × 200 m at the shoreline areas and further smaller inland. It should be noted that, for display purposes, the grid lines are not shown in Figure 5. All model boundaries that lie in water are specified as open; details of the boundary conditions can be found in the section below.

The bathymetric and topographic data are required for calculating the water depths and elevations of the grid cells in a model basin. The topographic data used in this study mainly come from USGS, and the bathymetric data come from NOAA. Water depths for grid cells at the open ocean were calculated based on the ETOPO1 global relief dataset from NOAA, which has a resolution of 1 arc minute ( 1.8 km). Water depths for grid cells in coastal areas were interpolated from the U.S. coastal relief dataset from NOAA with a resolution of 3 arc second ( 90 m) (https://www.ngdc.noaa.gov/mgg/bathymetry/relief.html). The USGS 90, 30, 10, and 3 m digital elevation models (DEM) were used to calculate the elevation of grid cells on the land (http://viewer.nationalmap.gov/viewer/).

At the open boundaries the model is forced using a nine–component tide comprising the M₂, S₂, K₁, O₁, Q₁, K₂, N₂, M₄, and M₆ components. These constituents were obtained from the ADvanced CIRCulation model (ADCIRC) Tidal Databases East Coast 2015 database of tidal constituents (Szpilka et al., 2016). An inverse pressure adjustment is made to the water surface specification at the tidal boundaries. The inverse pressure approach partially accounts for the meteorological forcing at the boundary by imposing the inverted barometer effect (Blain et al., 1994). Thus, the free surface at the open boundaries is modified by the amount Δη(x,y,t)=-ps′(x,y,tt)ρwg where ps′ is an atmospheric pressure change and ρ_w is the sea water density.

A spatially varying value of Manning's n is employed. Figure 6 shows the spatial distribution of the Manning's coefficient used in this paper. This coefficient map was generated based on the national land cover dataset (NLCD) created by USGS in 2001, using the approach proposed in Zhang et al. (2012), where details can be found. Note that for the open ocean area, a constant coefficient n = 0.02 was employed.

Figure 7 shows a comparison for the data at the NOAA tide gauges (Figure 5). The parallel CEST model in general captures the major trend of the tide and the storm surge induced by Irma (2017). At Naples, Fort Myers, and Mckay Bay Entrance, the surface elevations were relatively well predicted by CEST, while at Virginia Key, St Petersburg, and Old Port Tampa, CEST tends to underestimate the surface elevation variation.

Figure 8 presents the comparison for the HWM data. The CEST predictions have a root mean square error (RMSE) of approximately 0.69 m against the observations. This error is contributed significantly by the underprediction of the HWMs at Florida Keys and the South Florida mangrove zones (where the USGS gauges were marked, see Figure 5). A reason for this underprediction may be that these areas are close to the domain boundary, hence there is a limited fetch for wind to push water in.

Figure 9 presents the computed maximum storm surge height across the entire computational domain. It can be seen that the most severe storm surge inundation occurs at the right hand side of the track near the landfall location, where the area is known as the South Florida mangrove zone. The computed maximum storm surge height is over 3 m, but the overall inundation is kept within the mangrove zone due to the resistance of mangrove trees (Zhang et al., 2012). Figure 10 further shows the computed maximum storm surge profile along the four profiles depicted in Figure 9b. The maximum storm surge height gradually reduces as it moves inland. The inundation extents are roughly around 10 km at Profiles 1 and 2, and approximately half of that at Profiles 3 and 4. It should be mentioned however that storm surge could move further upland in the rivers as seen in Figure 9b. The overall maximum surge pattern is comparable with ADCIRC model results (Kowaleski et al., 2020).

The parallel performance of CEST on the 4-day simulation of Hurricane Irma (2017) was examined using the parallel efficiency (E_p) and speedup (S_p) coefficients. The coefficients are defined following Chen et al. (2018):

where T₁ and T_p are the total CPU time when using 1 and p processors. The simulation presented in this paper were run on a 32-core i7 CPU (3.7GHz) workstation. T₁ is approximately 150 min. The parallel efficiency and the speedup coefficients are plotted in Figure 11. As can be seen, the parallel CEST model achieves linear and even super-linear speedup when the number of processors used are small (≤ 4). When the number of processors increases the speedup increases but also becomes flattened and the parallel efficiency drops linearly. Despite that 4 processors appear to be the optimal number for the current simulation in terms of parallel efficiency, as many as 10 processors can be used to reduce the total CPU time of simulation as a priority. In the current case, T₁₀ is approximately 22 min.