Global Storm Tide Modeling with ADCIRC v55: Unstructured Mesh Design and Performance

This paper details and tests numerical improvements to ADCIRC, a widely used finite element method shallow water equation solver, to more accurately and efficiently model global storm tides with seamless local mesh refinement in storm landfall locations. The sensitivity to global unstructured mesh design was investigated using automatically generated triangular meshes with a global minimum element size (MinEle) that ranged from 1.5 km to 6 km. We demonstrate that refining resolution based on topographic seabed gradients and employing a MinEle less than 3 km is important for the global 5 accuracy of the simulated astronomical tide. Our recommended global mesh design (MinEle = 1.5 km) based on these results was locally refined down to two separate MinEle (500 m and 150 m) at the coastal landfall locations of two intense storms (Hurricane Katrina and Super Typhon Haiyan) to demonstrate the model’s capability for coastal storm tide simulations and to test the sensitivity to local mesh refinement. Simulated maximum storm tide elevations closely follow the lower envelope of observed high water marks (HWMs) measured near the coast. In general, peak storm tide elevations along the open coast are 10 decreased and the timing of the peak occurs later with local coastal mesh refinement. However, this mesh refinement only has a significant positive impact on HWM errors in straits and inlets narrower than the MinEle, and in bays and lakes separated from the ocean by these passages. Lastly, we demonstrate that the computational performance of the new numerical treatment is one-to-two orders of magnitude faster than studies using previous ADCIRC versions because gravity-wave based stability constraints are removed allowing for larger computational time steps. 15

more, barotropic global storm tide models can be used as components of Earth System Models to analyze risks posed by the long-term response of extreme sea level and coastal flooding to climate change in far greater detail than currently possible (Bouwer, 2018;Vousdoukas et al., 2018).
A key practical advantage of ocean models discretized using FEM as compared to FVM is that they are usually less sensitive to mesh quality (e.g., element skewness). Specifically, ocean models discretized using FMV often use staggered C-grid arrangements (e.g., Delft-FM) that have strict grid orthogonality requirements for numerical accuracy (Danilov, 2013;Fringer 30 et al., 2019). The orthogonal requirement makes mesh generation over wide areas with fractal shoreline boundaries an arduous task and is difficult to automate, although progress has been made (Herzfeld et al., 2020;Hoch et al., 2020). Despite the difficulties, the FVM Delft-FM (Flexible Mesh) based Global Tide and Surge model (GTSM) (Verlaan et al., 2015), has been meticulously developed and widely used to generate reanalysis datasets, describe historical trends, and make projections of extreme sea levels (Muis et al., 2016;Vousdoukas et al., 2018;Muis et al., 2019;Dullaart et al., 2019). The minimum (coastal) 35 resolution of GTSM has been historically limited to ∼5 km but recently upgraded to 2.5 km (1.25 km in Europe) (Dullaart et al., 2019).
In the absence of constraints on orthogonality or element skewness, automatically generating unstructured triangular meshes on the spherical Earth that accurately conform to the coastline and cover a wide-range of spatial sales [O(10 m)-O(10 km)] is completely realisable (Legrand et al., 2000;Gorman et al., 2006;Lambrechts et al., 2008;Roberts et al., 2019a). In one study, 40 automatically generated ocean basin-scale meshes with variable element sizes (50 m to 10 km) were used to conduct dozens of numerically stable FEM simulation experiments without mesh hand-edits or numerical limiters (Roberts et al., 2019b). The ability of FEM models to handle rapid transitions in mesh element sizes combined with the ease of mesh generation with modern technologies (e.g., Roberts et al., 2019a) enables the application of seamless local refinement directly into the global mesh where required, potentially on-the-fly based on atmospheric and ocean conditions that indicate a risk of coastal flooding. The ADCIRC storm tide model used in this study is a FEM solver that has been extensively used for detailed hurricane inundation studies at local or regional scales (e.g., Westerink et al., 2008;Bunya et al., 2010;Hope et al., 2013), and as an operational storm tide forecast model run by the U.S. National Oceanic and Atmospheric Administration (NOAA) (Funakoshi et al., 2011;Vinogradov et al., 2017). ADCIRC solves the SWEs that are composed of primitive continuity and non-conservative 60 depth-averaged momentum equations under astronomical and atmospherical forcing. After neglecting radial velocity terms, we formulate these equations in spherical coordinates as follows (Kolar et al., 1994), where ζ is the free surface, U and V are the depth-averaged zonal and meridional velocities respectively, H is the total water depth, t is time, λ is longitude, φ is latitude, R is the radius of the Earth, and ρ 0 is the reference density of water. Additional terms are defined as follows, 70 Ψ = p s ρ 0 − gη : external pressure and astronomical forcing f = 2Ω sin φ + tan φ R U : Coriolis + component of advection expanded in spherical coordinates : quadratic surface stress due to winds τ b = ρ 0 C f U 2 + V 2 : quadratic bottom stress due to friction 3 https://doi.org/10.5194/gmd-2020-123 Preprint. Discussion started: 28 July 2020 c Author(s) 2020. CC BY 4.0 License.
where p s is the surface air pressure, g is the gravitational acceleration, η is the summation of the equilibrium tidal potential and self-attraction and loading (SAL) tide (Ray, 1998), Ω is the angular speed of the Earth, R is the radius of the spherical Earth, C D is the quadratic wind drag coefficient, ρ s is the density of air at the ocean surface, U w and V w are the zonal and meridional 10-m wind velocities respectively, and C f is the quadratic bottom friction coefficient. C is the internal wave drag tensor that accounts for the energy conversion from barotropic to baroclinic modes through internal tide generation in the deep 85 ocean (Garrett and Kunze, 2007). Here, the local generation formulation is used (cf. Pringle et al., 2018a, b), in which C it is a global tuning coefficient, N b and N m are the seabed and depth-averaged buoyancy frequencies respectively, ω is set to the angular frequency of the M 2 tide, and ∇h λ , ∇h φ are the zonal and meridional topographic gradients respectively. Last, τ denotes the lateral stress tensor with ν t denoting the lateral mixing coefficient. The components τ λφ and τ φλ can be chosen to be either symmetric or non-symmetric as desired (Dresback et al., 2005). For this study we choose the symmetric option.

90
To properly compute the governing equations on the spherical Earth in the FEM framework used by ADCIRC, we have upgraded the model formulation and code as detailed in Sect. S1 of the supplementary document. This involves rotating the Earth so that the pole singularity is removed (Sect. S1.3) before applying a rectilinear mapping projection to transform the governing equations into a Cartesian form with spherical-based corrections to the spatial derivatives (Sect. S1.1). Here, the continuity equation is multiplied by a factor dependent on the choice of cylindrical projection used (e.g., Mercator) to produce 95 a conservative form that leads to discrete mass conservation and stability (cf. Hervouet, 2007;Castro et al., 2018) (Sect. S1.2).
The stability of the ADCIRC solution scheme (Sect. S1.4) was analyzed in one-dimensional linear form (Sect. S1.5) to provide guidelines for the choice of numerical parameters that can be chosen to remove the gravity-wave based (CFL) constraint. The validity of this analysis in the 2D nonlinear form has been demonstrated through the numerical simulations presented in this paper. With a semi-implicit time integration scheme, computational time steps (up to 120s) permitted are much larger than the 100 CFL constraint and as a result facilitate computationally efficient global simulations. Hereafter, the updated code in this study refers to the new release, ADCIRC v55. Solutions using the uncorrected model formulation are referred to by the previous version, ADCIRC v54.

Unstructured Triangular Mesh Generation on the Earth
The global unstructured meshes in this study are generated automatically using scripts with Version 3.0.0 of OceanMesh2D 105 (Roberts et al., 2019a;Pringle and Roberts, 2020). No post-processing hand-edits of any mesh were necessary to facilitate numerically stable simulations. Meshes are built in the stereographic projection centered at the North Pole to maintain angle conformity on the sphere and have the elements wrap around the Earth seamlessly including an element placed over the North Nominal element-to-element gradation limit on resolution ⇒ |∇ER| < α MaxEle Nominal maximum element size bound ER ≤ α ds: shortest distance to the shoreline, TM 2 : period of the M2 tidal wave, h: still-water depth, h * : low-pass filtered h, in which F lp (L) is the low-pass filter with cutoff length, L.
Pole (Lambrechts et al., 2008). Interested readers can execute "Example_7_Global.m" contained within the OceanMesh2D package to generate their own global mesh in a similar fashion.

110
Mesh design is handled through mesh size (resolution distribution) functions that are defined on a regular structured grid, usually that of the topo-bathymetric digital elevation model (DEM). In this study we use functions based on distance-toshoreline, bathymetric depths, and topographic gradients (see Table 1 for definitions). The final mesh size function is found by taking the minimum of all individual functions, and applying nominal minimum and maximum mesh resolution bounds and a element-to-element gradation limiter to bound the transition rate (Roberts et al., 2019a). The effects of the individual mesh size 115 functions and bounds on barotropic tides in a regional model have been previously detailed in Roberts et al. (2019b), which we use to guide our experiments exploring mesh design.
Additionally, this study makes use of the OceanMesh2D "plus" function which seamlessly merges two arbitrary meshes together keeping the finer resolution in the overlapping region. We use this function to apply local mesh refinement to a global mesh in storm-affected coastal regions to better resolve semi-enclosed bays and lakes, inlets, backbays, channels, and other 120 small-scale shoreline geometries.

Experimental Design
The experimental design we pursue is composed of two distinct steps, both with the purpose to maximize model efficiency while maintaining a threshold of accuracy. First, we begin with a mesh design that we assume is a highly-refined discretization of the Earth in a global sense, and systematically relax the mesh size parameters to reduce the total number of mesh vertices 125 while trying to minimize any negative impacts on global model accuracy. Second, we take the resulting recommended global mesh design from the previous step, and apply local mesh refinement to increase the coastal resolution in the storm landfalling region and potentially improve local model accuracy. for global ocean tides (Greenberg et al., 2007), and an element-to-element gradation limit/distance expansion rate (G and D) as high as 0.35 is tolerable as long as the TLS mesh size function is applied (Roberts et al., 2019b). As a note, the FL parameter is used in the construction of the TLS mesh size function to filter out small-scale topographic features (c.f., Table 1) that are potentially unimportant to the local barotropic physics, thus disregarding those features for the application of higher resolution.
This study is the first to test the effect of the FL parameter in detail. In each mesh design perturbation, only one of the three To assess the effect on the global model accuracy as the mesh designs are coarsened, we compare simulated astronomical tidal solutions to the data-assimilated TPXO9-Atlas. We focus on astronomical tides in this step because they can be reduced to a series of harmonic constituents of well-defined frequencies to make systematic global comparisons (Roberts et al., 2019b;  Table 2).

Step 2: Local Mesh Refinement
In this step, the recommended global mesh design from Step 1 (Sect. 2.3.1) is used but with additional patches of high resolution (local mesh refinement) near the landfall location of two storm events. We choose to focus on a particularly significant historical To assess the accuracy and inter-compare mesh designs, we primarily use the output of the maximum simulated storm tide elevation from each event. For validation we compare to surveys of high-water marks (HWMs) in the landfall regions that are 170 located close to the coast (within 1.5 km of the nearest 150-m locally refined mesh vertex) and attributed primarily to surge for both Katrina (URS Group Inc, 2006a, b) and Haiyan (Mori et al., 2014). Note that for Katrina we add a value of 0.23 m to the simulated storm tide elevations to account for a steric offset and the conversion to NAVD88 vertical datum from local mean sea level (Bunya et al., 2010)

180
The full-resolution Global Self-consistent Hierarchical High-resolution Shorelines (GSHHS) dataset (Wessel and Smith, 1996) is used to define the shoreline boundary when making the mesh. The most recent global bathymetry DEM, GEBCO_2019 8 https://doi.org/10.5194/gmd-2020-123 Preprint. Discussion started: 28 July 2020 c Author(s) 2020. CC BY 4.0 License. (GEBCO Compilation Group, 2019), which has an equatorial resolution of ∼500 m, is used to prescribe the bathymetry for the model (Sect. S2.2). Underneath Antarctica ice shelves, the RTopo-2 DEM (Schaffer et al., 2016), which has an equatorial resolution of ∼1 km, is used to prescribe ocean depths taking into account the ice shelf thickness. The SAL tide is specified  August 25 to August 31 (Fig. 4b). Super Typhoon Haiyan meteorology is described by the best-track Holland parametric vortex model (Holland, 1980) from November 4 to November 10 ( Fig. 4c).
When only astronomical tides are simulated, we force the model with only astronomical forcing (η) for the five leading 195 astronomical tidal constituents (M 2 , S 2 , N 2 , K 1 , O 1 ), and analyze for the corresponding constituent amplitude and phases using a 28-day harmonic analysis. These five constituents are chosen so that we can use a relatively short 28-day harmonic analysis period (Ngodock et al., 2016), which would otherwise need to be extended to around 180 days if other constituents are included because of the closeness in their frequencies (e.g., K 1 and P 1 ) (Pringle et al., 2018a). When simulating storm tides, both atmospheric (τ w and p s ) and astronomical forcings are invoked, this time using the following ten tidal potential-generating 200 constituents to obtain a more complete tidal signal: M 2 , S 2 , N 2 , K 2 , O 1 , P 1 , K 1 , Q 1 , MF, MM. In both cases the model is spun- When simulating on the Ref mesh using ADCIRC v55, significant improvements to the prediction of astronomical tidal constituents were measured compared to ADCIRC v54. In particular, M 2 amphidromes in the high latitude regions are largely corrected such that any disparities between TPXO9-Atlas and our updated model solutions are qualitatively hard to discern 210 from a global perspective (Fig. 5). Moreover, M 2 RMSE t is less than 2.5 cm over most of the ocean with the largest remaining deep ocean hotspot in the North Atlantic (Fig. 6). Indeed, deep ocean M 2 RMSE t = 2.87 cm (Table 3), which is smaller than those previously computed for non-assimilated barotropic tidal models (Stammer et al., 2014;Schindelegger et al., 2018), and within the range of errors computed for solutions obtained by embedding a state ensemble Kalman Filter (perturbed data assimilation) into a forward ocean circulation model (Ngodock et al., 2016). Furthermore, the 5-constituent total tidal error, 215 RMSE t|tot , is less than 4 cm in the deep ocean. In shallow regions, the M 2 RMSE t is 13.9 cm, which is slightly smaller than obtained in Schindelegger et al. (2018), and RMSE t|tot is 17.2 cm. Note that the area-weighted median value of RMSE t|tot (c.f., Appendix A) in shallow water is just 6.63 cm.  Table 2 for mesh design details.

Solution Variability with Global Mesh Design Parameters
The distribution (ECDF curves) of RMSE t|tot degrades as the mesh size function parameters are relaxed (Fig. 7). Most of this 220 degradation occurs in the body of the distribution rather than the tails. This characteristic implies that the K test statistic is a good metric of disparity between mesh designs (Fig. 8) because it measures the greatest vertical distance between ECDF curves which has a greater chance of being larger in the body. K is positive for all mesh designs indicating that the Ref mesh is indeed statistically the best performing mesh.
Increasing MinEle has a clear but gradual degenerative effect on the solution as it is increased from 1.5 km to 3 km. The 225 disparities in the ECDF curves noticeably grow as MinEle is increased to 6 km; the value of K increases from 0.057 for 3 km to 0.175 for 6 km. In comparison, the ECDF curves and the value of K is changed comparatively little as the TLS parameter is decreased from 20 to 5 (K is just 0.032 for TLS = 5). In fact, the solutions are close to identical for TLS values of 10 and 20. However, as the TLS function is turned off (TLS = 0), the solution is severely degraded and the value of K is the greatest for any mesh design (= 0.257). The effect of using the FL function and relaxing the parameter on the ECDF curves is fairly 230 gradual overall, however the magnitude of this change is not trivial; K is increased to 0.070 for FL = 1/5.
In summary, the results demonstrate that decreasing the TLS parameter from 20 through to 5 substantially decreases the number of vertices while it has a relatively small effect on the tidal solution compared to the other experiments. On the other hand, increasing the FL parameter has a comparatively small impact on vertex count reduction, while increasing MinEle has a relatively large impact on the solution. The final choice of mesh design is dependent on one's tolerance for error, but in 235 general it is preferable to choose a mesh that is plotted close to bottom-left corner of the graph in Fig. 8. Following this logic, the TLS-B mesh design (MinEle = 1.5 km, TLS = 5, FL = 0) appears to be the most efficient one tested (see Fig. 6 for the spatial distribution of M 2 RMSE t on this mesh). The results also suggest that combinations of MinEle larger than 1.5 km (∼2.25-3 km) and a TLS parameter smaller than 20 (∼5-10) could be employed to potentially lead to an efficient mesh design.  Table 2).
Nevertheless, using a small MinEle is in and of itself useful to provide extra coastal resolution, which may be more important 240 as we consider local storm tide accuracy. Thus, TLS-B is chosen as the base mesh design in Sect. 3.2.

Validation on the 150-m Locally-Refined Meshes
The maximum storm tide elevation due to Hurricane Katrina approaches 8 m in the Hancock and Harrison Counties of Mississippi ( Fig. 9a), comparable to previously conducted high-fidelity simulation results (Dietrich et al., 2009)

245
Haiyan, the maximum storm tide elevation exceeds 6 m near Tacloban due to local amplification in the Leyte Gulf (Fig. 9b), similar to previous simulation results ( Mori et al., 2014). Qualitatively good agreement with the plotted HWMs is demonstrated for both storms with a few exceptions.
As pointed out by Mori et al. (2014), since inundation is not simulated and the effects of wave setup are ignored in these simulations we expect the maximum storm tide height to match the lower envelope of observed HWMs due to the amplification 250 by topography. By numbering the HWMs starting at the most southwest point and following the shoreline clockwise, it is indeed illustrated that the simulation results tend to mach the lower envelope of observed HWMs quite closely (Fig. 10)  Points are numbered by starting from the most southwest point and following the shoreline clockwise around.
Time series of Hurricane Katrina at NOAA tide gauges show that the timing of the peak storm tide elevation and the amplitude and phase of the tide signal prior to landfall are well represented by the model (Fig. 11). The modeled peak is underestimated by ∼0.4 m at gauge 8735180, and at gauge 8743281 where the largest peak storm tide occurred, tide gauge recording was interrupted as the storm was making landfall, but the simulation closely follows the observations up until this 260 point. A roughly constant discrepancy of ∼0.4-0.5 m between the simulation and observation develops at the gauges following the last high tide prior to storm landfall, likely explaining the underestimate at gauge 8735180 (by this logic the peak may be overestimated at 8761724). Comparing to time series of simulations with wind wave-coupling (Roberts and Cobell, 2017), it appears that this discrepancy can be attributed to wave setup effects that are ignored here. Time series of Super Typhoon Haiyan also show that the amplitude and phasing of the tide at the selected locations are fairly well represented as compared to 265 the TPXO9-Atlas (Fig. 11). Storm tide heights at Tacloban are dominated by the short and intense surge event but the duration of surge is likely underestimated because the parametric vortex model lacks background winds. Due to the timing of the storm landfall during the lower high tide, peak storm tides exceeded the higher high tide by just 0.5 m at Guiuan. In fact for Guiuan and Canuay Island the minimum storm tide levels were more severe than the peak levels. 0.2 m larger on MinEle = 500 m mesh than the MinEle = 150 m mesh (Fig. 12). Interestingly, neither the MinEle = 1.5 km mesh or the MinEle = 500 m mesh includes Lake Pontchartrain, while the MinEle = 150 mesh does (Fig. 2). This is because 275 the Rigolets strait connecting Lake Pontchartrain to Lake Borgne is approximately 500 m at its narrowest, thus the MinEle = 500 m mesh is at the cutoff point for meshing the strait and providing hydraulic connectivity between the two lakes. Yet, the maximum elevation difference from Lake Borgne across to Mobile Bay between MinEle = 1.5 km and MinEle = 500 m is still much greater than between MinEle = 500 m and MinEle = 150 m (refer time series at gauges 8735180 and 8743281 in Fig. 11 in addition to Fig. 12). Nevertheless, since the MinEle = 150 m is the only mesh to resolve Lake Pontchartrain, the 280 HWMs surrounding the lake (point numbers 7-23) can only be reasonably estimated by simulations on this mesh, resulting in the smallest HWM error statistics (Fig. 10a).

Solution Variability with Local Mesh Refinement
In the case of the Haiyan, the predominant maximum storm tide elevation difference (∼0.4-0.6 m) is located at the back of Leyte Gulf near Tacloban between the MinEle = 1.5 km and MinEle = 500 m meshes (Figs. 11 and 13). The decrease in elevations in the MinEle = 1.5 km mesh might be explained by the omission of the San Juanico Strait (∼800 m wide at its 285 narrowest, Figure 3). In contrast, the elevations near Tacloban and in the strait increase by up to 0.6 m as the MinEle = 500 m mesh is refined to MinEle = 150 m. There is also a reduction in the elevations (by up to 0.2 m) just offshore of the Leyte Gulf along the shelf break for the coarser meshes, which is better represented in the higher resolution renditions. Overall, the choice of local MinEle has a small effect on the representation of HWMs for Haiyan (Fig. 10b). The only major noticeable impact is for the HWMs near Tacloban in the San Juanico Strait (point numbers 25-63) on the MinEle = 1.5 km mesh. Since this mesh 290 does not resolve the strait, the simulated estimate of the HWMs here are all taken from the same or nearby mesh vertices at the back of the gulf.
Lastly, we mention two additional general observations. First, for both storms the far-field effects of local mesh refinement were found to be negligible. Second, storm tide elevation time series show that the not only does the peak elevation tend to decrease as mesh refinement is made, but the timing of the peak also tends to occur later (Fig. 11) -most clearly illustrated for 295 gauge 8735180.

Computational Performance
For all astronomical tide simulations performed in the global mesh design experiments (Sect. 3.1) the time step, ∆t, was set to 120 s, which is equivalent to a Courant (CFL) number of 5-22 on the global mesh designs used in this study (see Sect. S1.5 for model stability details). The resulting computational wallclock times for the astronomical tide simulations ranged from 5 s 300 to 30 s per simulation day depending on the total vertex count of the mesh (Fig. 14). All simulations were computed on 480 Haswell processors with a Mellanox FDR Infiniband network connection.
For the storm tide simulations performed in the local mesh refinement experiments (Sect. 3.2) the time step had to be reduced for finer meshes due to model instability related to nonlinear terms that are treated explicitly which become more significant as element sizes nearshore (in shallow depths) are reduced. For Katrina, the time step for the MinEle = 150 m mesh had to 305 be reduced to 40 s. For Haiyan, the time step had to be reduced for both the MinEle = 500 m (∆t = 90 s) and the MinEle = 150 m (∆t = 25 s) mesh. The resulting computational wallclock times were increased for the smaller time steps (Fig. 14).
However, the simulations on the MinEle = 150 m meshes were proportionally around two times faster per ∆t compared to the coarser meshes. This is attributed to the heavy I/O related to reading meteorological data during the simulation that limits computational speed-up as time steps are increased beyond a certain value. In fact, for the same ∆t, wallclock times for the 310 storm tide simulations were 2-3 times longer than the astronomical tide runs which require I/O only at the very beginning and end of the simulation. Furthermore, the Haiyan simulations were slightly slower than Katrina simulations because of the higher resolution of the reanalysis meteorology (c.f., Sect. 2.4).
In addition, the computational performance of our Katrina simulations are compared to previous ADCIRC model ones (Bunya et al., 2010;Dietrich et al., 2011), which employed ∆t = 1 s leading to wallclock times one-to-two orders of magnitude 315 greater than those in this study (Fig. 14). We remark that the Bunya et al. (2010);Dietrich et al. (2011) studies used finer mesh sizes (MinEle = 50 m) and included the floodplain leading to more wetting-drying, which can be a source of numerical instability. Our results nonetheless indicate the potential for substantial computational speed-up on suitably designed meshes.

Discussion
Results from the global mesh design experiments echo those found previously for a regional mesh of the Western North 320 Atlantic Ocean (Roberts et al., 2019b). Namely, that when employing a relatively large element-to-element gradation limit (G) computational processor. All simulations were computed on 480 computational processors (i.e., the variation in vertices per processor comes from the variation in the total number of vertices for each mesh design). The computational performance in this study using ADCIRC v55 is contrasted with previous ADCIRC model Katrina storm tide simulations (Bunya et al., 2010;Dietrich et al., 2011) that required very small time steps (∆t = 1 s).
of 35% used here, the TLS function is critical to ensure that water elevation solutions remain accurate by adequately resolving topographic features of importance. The relatively high value of G ensures that the mesh size rapidly enlarges in areas where no topographic gradient features exist to help avoid excessive mesh vertex counts. An interesting point of difference that we note in this study is the interplay between MinEle and TLS. MinEle not only controls the minimum coastal resolution but it 325 also controls the minimum resolution that can be applied along topographic gradient features in the TLS function. Here, the MinEle ranged between 1.5 km and 6 km which affected the strength of the TLS function in these mesh designs. Nevertheless, the MinEle-C design (MinEle = 6 km) was still superior to the TLS-C one (TLS = 0). There may be additional sensitivity to the TLS function at MinEle less than 1.5 km although it is likely quite small. In Roberts et al. (2019b), MinEle = 1 km was used along the outer shelf and shelf break and this appeared to be sufficient.

330
The results also show that the use of the barotropic Rossby radius-based low-pass filter (FL) in the TLS function is able to reduce mesh vertices without substantially degrading the solution. However, our results suggest that it is more efficient to simply reduce the TLS value to 5 compared to using TLS = 20 with FL. Combinations of TLS values between 5-10 with FL (=1/80-1/20) could also be applied depending on the tolerance to the solution accuracy. Using the FL parameter in more highly resolved local domains may also provide benefits to avoid over-resolving topographic slopes in shallow depths. How-335 ever, it must also be recognized that length scales of importance transition from Rossby radius scaling in the open ocean to those of shoreline features nearshore (LeBlond, 1991). Therefore, a more intelligent low-pass filter may be necessary for such applications.
The astronomical tide solution differences between global mesh designs were shown to be predominantly in the body of the area-weighted ECDF curves while the tails were almost identical. This implies that simulated tides over most of the area 340 of the ocean are affected by mesh resolution. However, in regions where the tidal range and error is large, which inevitably occurs on shallow shelves, all of the mesh designs have similarly large errors. This perhaps explains why a 1/12 • tidal model (Schindelegger et al., 2018) could obtain a similar M 2 RMSE t in shallow waters to our more highly resolved reference mesh (14.6 cm compared to 13.9 cm). There is likely an inherent uncertainty in the bathymetry and dissipation which prevents further decrease to the shallow water RMSE t which can be dominated by large errors (in shallow water the area-weighted 345 median RMSE t|tot is much smaller than RMSE t Last, it is widely recognized that sensitivities to local high resolution bathymetry datasets, and spatially varying bottom friction and surface ice friction are important (Lefevre et al., 2000;Le Bars et al., 2010;Zaron, 2017;Pringle et al., 2018a;Zaron, 2019), likely more so than the mesh resolution effects that we concentrate on here. We aim to develop a unified framework for calibrating bottom friction globally with improved local high resolution bathymetric datasets in future work. Moreover, we did not simulate inundation in this study, but it will be a crucial future step so that we can forecast and assess coastal flood risks 365 for all of Earth's coasts. Based on the results for global mesh design we recommend to aim for a minimum element size less than 3 km and to use the TLS function to resolve topographic gradient features with a TLS value of 5-10. Paired with the OceanMesh2D software, the ability to seamlessly apply local refinement allows the user to provide fine coastal resolution in the region of interest (e.g., the storm landfall region) without large increases to the total mesh vertex count (increase of 0.5-3% in this study). We found 375 that in general, peak storm tide elevations along the open coast are decreased and the timing of the peak occurs later with local coastal mesh refinement. When validated against observed high water marks measured near the coast, coastal mesh refinement only has a significant positive impact on errors in narrow straits and inlets, and in bays and lakes separated from the ocean by these passages.
The new ADCIRC v55 code capable of accurate global storm tide modeling with fine coastal resolution is computationally 380 efficient. For global meshes with minimum resolution as fine as 1.5 km, the computational wallclock time ranged from 5 s to 30 s per simulation day on 480 computational processors for astronomical tide simulations. Some improvements that we made to the numerical stability of the algorithm facilitated the application of relatively large 120 s time steps to achieve this efficiency. However, we found that the locally refined meshes often required smaller time steps (25 s-90 s). Nevertheless, these are still much larger than time steps used in previous studies with older versions of ADCIRC, resulting in one-to-two orders of 385 magnitude shorter computational times.
Code and data availability. The official release version of ADCIRC is available from the project website: http://adcirc.org/ under the terms stipulated there and is free for research or educational purposes. The exact version of the model used to produce the results used in this paper is archived on Zenodo (Pringle, 2020) along with the simulated tidal harmonic constituents on the various mesh designs and the simulated storm tide results with the corresponding model setup (Pringle, 2020).

390
The current version of OceanMesh2D is available from the project website: https://github.com/CHLNDDEV/OceanMesh2D under the GPL-3.0 license. The exact version of the model (V3.0.0) used to generate the meshes in this study is archived on Zenodo (Pringle and Roberts, 2020).
An application of the presented ADCIRC v55 model (on the TLS-B mesh design) providing 5-day forecasts of global storm surge is currently running in real-time and maximum surge elevations are available to view at https://wpringle.github.io/GLOCOFFS/. The forecast 395 is automatically updated every 6-hours based on GFS weather forecast schedules. Simulation wallclock time is approximately 10 min on 72 computational processors. Images are automatically archived by GitHub.
where A is the tidal amplitude and θ is the tidal phase lag, and the subscripts 'o' and 'm' refer to the observed (in this case TPXO9-Atlas) and modeled values respectively. To get the RMSE t of the total 5-constituent signal, RMSE t|tot , we use the approximation (Wang et al., 2012): where k indicates the arbitrary constituent number from the five leading tidal constituents (M 2 , S 2 , N 2 , K 1 , O 1 ).

405
The area-weighted ECDF curves of RMSE t|tot are computed through MATLAB's "ecdf" function with the 'frequency' input set equal to the elemental area [m 2 ]. With this construction, the area-weighted median value of RMSE t|tot is defined as the point of intersection of the ECDF curve with the 50% probability. Comparisons between ECDF curves can be made by taking the largest vertical distance between them to obtain K. We compute K one-sided in that it will be positive if the base (Ref) design is statistically superior and negative otherwise.

410
To obtain a single number to compare the solutions globally or in certain depths or regions, it is common to use the areaweighted root-mean square error of the RMSE t (Arbic et al., 2004), where dA indicates an area integral. RMSE t can be computed for a single constituent or for the 5-constituent total signal (= RMSE t|tot ).

415
Author contributions. WP prepared the manuscript, designed and implemented the coding upgrades into ADCIRC v55, designed and performed the experiments, and conducted the stability analysis. DW mathematically formulated and initially implemented most of the ADCIRC coding upgrades. KR and WP equally contributed to the coding and development of the OceanMesh2D software integral to this study, and KR provided critical feedback to improve the manuscript design. JW provided the research and computing resources and funding necessary to conduct this study.