On the automatic and a priori design of unstructured mesh resolution

15 This study investigates the design of unstructured mesh resolution and its impact on the modeling 16 of barotropic tides along the United States East Coast and Gulf Coast (ECGC). A discrete repre17 sentation of a computational ocean domain (mesh design) is necessary due to finite computational 18 resources and an incomplete knowledge of the physical system (e.g., shoreline and seabed topogra19 phy). The selection of mesh resolution impacts both the numerical truncation error and the approx20 imation of the system’s physical domain. To increase confidence in the design of high-resolution 21 coastal ocean meshes and to quantify the e cacy of current mesh design practices, an automated 22 mesh generation approach is applied to objectively control resolution placement based on a priori 23 information such as shoreline geometry and seabed topographic features. The simulated harmonic 24 tidal elevations for each mesh design are compared to that of a reference solution, computed on a 25 10.8 million vertex mesh of the ECGC region with a minimum shoreline resolution of 50-m. Our 26 key findings indicate that existing mesh designs that use uniform resolution along the shoreline and 27 slowly varying resolution sizes on the continental shelf ine ciently discretize the computational 28 domain. Instead, a targeted approach that places fine resolution in narrow geometric features, along 29 steep topographic gradients, and along pronounced submerged estuarine channels, while aggres30 sively relaxing resolution elsewhere, leads to a mesh with an order of magnitude fewer vertices 31 than the reference solution with comparable accuracy (within 3% harmonic elevation amplitudes in 32 99% of the domain). 33

adjustments to frictional and dissipative parameterizations in order to agree with measured data. 140 However, when the mesh underresolves shoreline and seabed features, the system's response may 141 become distorted leading to an inability to correctly produce solutions across the entire domain 142 and energy spectrum. An example of this would be tuning the model to agree with observations of 143 dominant semi-diurnal elevation tidal constituent regionally but this may not lead to a good agree-144 ment globally nor for the other tidal constituents. Instead, by gathering knowledge on how tidal 145 solution depends on mesh resolution in realistic coastal modeling problems, we can enable e cient 146 and uniformly more accurate mesh designs that can then facilitate more dynamically correct cali-147 brations of friction parameterizations. 148 Our premise is that the numerical modeling of the circulation and flow of water is largely 149 driven and controlled by the representation of the physical system and the representation of the 150 physical system is integrally related to the mesh sizing functions. Thus, the sizing functions need 151 to be carefully considered for ensuring high fidelity coastal ocean hydrodynamic simulations that 152 have a relatively low associated computational cost. This is particularly relevant for operational/real-153 time forecast systems in order to be practically computationally feasible. Many of the previously 154 used a priori mesh size heuristics (e.g., topographic-length scale, and distance-to-shoreline) have 155 proven useful in practice for producing accurate solutions for tides and storm surges. Thus, we 156 have devised an approach that combines and builds on such mesh size heuristics to variably resolve excluded from calculations for numerical stability purposes; however, all terms were included in 239 the calculations. Further wetting/drying was enabled although a minimum depth is enforced on the 240 shoreline of 1 m below sea level to ensure flow through narrow channels on the scale of the min-241 imum resolution. A constant quadratic bottom friction was used with the standard coe cient of 242 0.0025. Horizontal dissipation was parameterized through a constant lateral eddy viscosity term of 243 50 m 2 s 1 . The GWCE mass matrix is solved using an explicit time discretization with mass lump-244 ing instead of the consistent implicit method. This choice was not found to a↵ect the simulation 245 results at the 2 second simulation timesteps we are using here with the Courant-limited explicit 246 timestepping scheme. Therefore, the explicit method was preferred due to improved computation- The construction of regional coastal ocean meshes for hydrodynamic simulations in mod-250 els such as ADCIRC is an involved process with many degrees of variation. In order to analyze 251 how mesh resolution may a↵ect numerical simulations, it is vital to have an automated and repro-252 ducible workflow to systematically control aspects of the mesh design. By reproducible we mean 253 that given the exact same inputs and options, the vertex locations of a new instance of the mesh 254 will be approximately the same having vertex/elemental densities within a fraction of the target 255 density function and leading to negligible di↵erences between simulation results repeated on vari-256 ous instances of the mesh. The approximate similarity of meshes is evidenced in results throughout 257 the manuscript: nearly similar mesh designs exhibit the smallest relative di↵erences between their 258 solutions. 259 Some approaches and tools have been developed recently to make these workflows feasible 260 [Engwirda, 2017;Gorman et al., 2008;Candy and Pietrzak, 2018;Roberts et al., 2018].  et al., 2018] with support for mesh gener-265 ation using map projections to ensure that meshes on the sphere conform to Earth's curvature and 266 obey user-defined resolution requests which are specified in meters. Any map projection that is 267 featured in the m_map mapping package [Pawlowicz, 2018] can be selected.

268
A number of meshes are automatically generated in Lambert conformal conic projection 269 space using the multiscale meshing approach [Roberts et al., 2018], whereby multiple boxes are 270 used to cover the region roughly indicated by the green and red colored zones in Figure 1(a)-(b).

271
Inside these boxes, the minimum resolution L min is specified to between 50 m and 250 m, depend-272 ing on the experiment (see Section 2.4). A larger box covering the whole study region is used to 273 mesh the rest of the domain with a minimum resolution of 1 km that is placed uniformly along 274 the shoreline. The result is one seamless unstructured mesh, in which the software automatically 275 smooths mesh resolution sizes between regions.

276
Topo-bathymetric data, available on a structured grid (DEM), is interpolated onto the mesh 277 vertices using the grid-scale averaging approach that is built into the mesh generation software 278 [Roberts et al., 2018]. Grid-scale averaging is used to minimize aliasing of the seabed topography 279 on the mesh vertices that would other arise from curve-fitting interpolation schemes (e.g., linear 280 interpolation). The minimization of sub-grid scale topo-bathymetric features in the interpolated 281 seabed topography is important in order to study the e↵ect of mesh resolution on the solution.  (Table 1). Within each experiment three meshes (categorized as 'fine', 'medium', and 286 'coarse' resolution) are generated by varying a single mesh size function parameter while holding 287 all other parameters constant. Note that the variation of mesh sizing parameters is a multi-faceted 288 problem and all the parameters interact (e.g., one parameter's value can mask e↵ects of another).

289
For example, a relatively higher feature size may cause finer resolution in deep o↵shore features 290 that can be largely influential on the simulation of tides, as later shown. All meshes require a min-291 imum mesh size and an element-to-element mesh size gradation rate (henceforth referred to as gra-292 dation), which are set to 50-m uniformly along the shoreline and 15%, respectively, unless other-293 wise stated. The maximum mesh size is set to 10 km for all meshes.

294
The e↵ect of the mesh size functions on the resulting triangulation's that are used in the vari-295 ous experiments (Table 1) (Table 1) was generated to act as a proxy for the 323 'true' solution against which our meshes in the experiments are compared. In this mesh, a set 324 of depth-based maximum element size constraints were used and a mesh size gradation of 15%. contains N = 10,746,955 vertices and represents a mesh design that we classify as 'overly-discretized' 331 in the sense that as this study will later demonstrate, it is possible to substantially reduce the vertex 332 count while maintaining solution accuracy. It is important to note that for all mesh configurations 333 the REF mesh is indeed finer than the other mesh designs except for the S20 mesh design in a lo-334 cal region on the Western side of the GOM.

335
Each mesh was used to perform a 122-day tidal simulation to assess the e↵ects on the astro-336 nomical tides due to variations in mesh design. In these simulations, ADCIRC is forced through 337 the tidal equilibrium potential and SAL terms throughout the domain and at the open ocean bound-338 aries with four major semi-diurnal (M 2 , N 2 , S 2 , K 2 ) and four major diurnal tidal constituents (K 1 , 339 O 1 , P 1 , Q 1 ). Open boundary elevations are obtained from TPXO9.1 (http://volkov.oce.orst.edu/ 340 tides/global.html) tidal solutions; SAL terms are obtained from FES2014 tidal loading solutions 341 (ftp://ftp.legos.obs-mip.fr/pub/FES2012-project/data/LSA/FES2014/). In the assessment of the re-342 sults of these simulations, a focus is placed primarily on the variation in the major semi-diurnal 10 35 L min : minimum mesh size L ma x : maximum mesh size L i : mesh size at i defined by the circumradius of each triangle X i : coordinate of i on grid to compute edgelengths i and j: adjacent elements g: gradation d s : shortest distance to the shoreline d m : shortest distance to the medial axis b: topo-bathymetric depth (positive below sea level) r: gradient operator tide (M 2 ) since this is the predominant tidal constituent along the ECGC. The major diurnal tide 344 (K 1 ) is also included where relevant.

345
The tidal elevations are decomposed into harmonic constituents using a least-squares method

RE =
where A is the harmonic elevation amplitude of the tidal constituent in the experiment (ID) and 351 the REF meshes. A focus is placed on the M 2 and K 1 elevation amplitudes as these represent the 352 predominant semi-diurnal and diurnal constituents in the ECGC domain (Section 2.1).

353
The calculation of the RE is proceeded in this manner to keep data extrapolation to a mini-354 mum so that the same shoreline geometric complexity as depicted in each mesh is present in both  Figure 1). To be consistent throughout, CAFE curves only consider 363 errors (A I D A RE F ) that exceed 1 mm or feature a RE greater than 0.1%. On these CAFE plots, greater (less) than the positive (negative) value on the x-axis. A solution that has "converged" in-366 dicates that 99% of the comparison region has a ±5% RE. This definition of convergence may be 367 arbitrary but it represents a statistic that can enable a consistent comparison between solutions and 368 a more strigent accuracy standard than currently set by the U.S. COASTAL Act1.

369
Last, in Section 3.5 we summarize the experiments through the standard deviation of the where E is the complex root-mean-square error of a tidal constituent for one cycle and account 379 for the amplitude and phase errors, A and ✓ are the amplitudes and phase lags of the tidal con- The representation of the shoreline determines the simulated accuracy in modeling the phys-394 ical interaction between forcing agents (e.g., tides, winds, and waves) with shoreline geometri-395 cal features (e.g., coves, headlands, back-bays, and lagoons). From a modeling standpoint, the 396 shoreline's representation must be simplified to satisfy computational resources by removing fine-  It is important to note that the variation in the minimum element size along the shoreline will [Massey, 2015;Shewchuk, 2002]. Further, the analysis by Hagen et al. [2000] for one dimensional 493 domains demonstrates that a high gradation value (g ⇡ 0.5) leads to the introduction of odd order 494 truncation error term, which lowers the order of the method to first-order accurate and/or degrades 495 the local/global accuracy.

496
Fundamentally, the gradation rate will impact many aspects of the mesh design at once. A 497 higher valued mesh size gradation will degrade the approximate representation of the seabed to-498 pography by creating comparatively coarser mesh sizes away from the targeted zones of fine res-499 olution. As was described in Table 1, the meshes that vary the gradation rate utilize a minimum 500 resolution of 50-m along the shoreline (L50). Note that the mesh generator is bounding the gra- where V is the total mesh volume for the mesh denoted by ID and is calculated as the sum of all it is apparent that there is a diminishing reduction in the total vertex count of the mesh with in-516 creased gradation. For the purposes of this study, we were not able to explore meshes with gra-517 dation greater than 35% due in Experiment 3 (Gx) due to the introduction of triangles with very 518 skewed aspect ratios and obtuse and acute angles that created numerical accuracy issues.

519
The increase in mesh size gradation from 15% to 35% leads to a highly amplified error pat-528 tern in the NA region for both M 2 and K 1 constituents as well as along the MAB for M 2 (Fig-529 ure 8). In the NA subdomain (Gulf of Maine), the M 2 RE is increased from 2-5% for G15 to 10-530 21% for G35 (colors are saturated in Figure 8a), in which the maximum RE is focused on the 531 Georges Bank. In contrast to the response in the M 2 's RE, the K 1 's RE is nearly uniformly de-532 graded from -3% for G15 to -6% for G35 in the NA subdomain. The M 2 RE in the MAB, SAB, 533 and eastern GOM tends to weakly deamplify by approximately 1% to 5% along the continental 534 shelf zones. In contrast to the shoreline approximation experiment, a relatively large deamplifica-535 tion of the M 2 RE occurs in both the Chesapeake Bay and Delaware Bay as the gradation is en-536 larged to 35% (Figure 8a,b). The M 2 RE reaches as high as 15% in this region for the G35 experi-537 ment (colors are saturated in Figure 8a).

538
As the mesh size gradation grows, the tidal elevation amplitudes start to diverge substantially 542 from the REF solution (Figure 9). In 99% of the comparison zone, the G15 mesh has an M 2 error 543 between -1.3% and +3.0% RE whereas G35 has between -5.0% and +15% RE. Furthermore, the 544 G35 mesh design exhibits between -5.5% and +6.5% K 1 RE in the 99% comparison zone, which 545 is compared to -3.0% and +0% K 1 RE for G15. Unlike the shoreline experiment where all meshes 546 converged (based on the ±5% threshold definition of convergence), only the G15 mesh converges 547 for the M 2 constituent, and the G15 and G25 meshes converges for K 1 . However, it's important 548 to note that the tendency of the solution is convergent as the RE reduces when the mesh sizes are 549 made finer with lower gradation bounds.   propose an alternative strategy to resolve documented later on. Insets are shown to reflect areas that are described in the text. The e↵ects of the estuarine channel mesh size function have been investigated in Experi-642 ment 5 (Cx) using a mesh size gradation of 35% (G35). A higher mesh size gradation motivates 643 the resolution targeting approach because mesh element sizes are relaxed quickly away from the 644 channel thalwegs where finer resolution is applied, thus obtaining a mesh with overall fewer ver-645 tices than without the targeting approach. Furhtermore, a lower mesh size gradation (e.g., 15%) 646 would lead to finer resolution in the center of the estuary where the thalweg may be located and 647 may already adequately resolve the channels cross-sectional profile. The mesh vertex count in the 648 finest Cx mesh (C1.0) is increased by more than two-fold from the G35 mesh to approximately 3.1 649 million vertices (Figure 14d-f), still approximately 60% of the G15 mesh vertex count.

650
The refinement of the estuarine channel network primarily impacts the M 2 elevation ampli- (c.f., Figure 1) in the MAB region is reduced to the point that the error changes sign for C1.0 (+1-662 2% RE) (Figure 14a-c). Similarly, the CAFE curves also demonstrate a consistent reduction in M 2 and K 1 RE in the  (Table 1). Panels (c) and (d) are the same but for the standard deviation of the dimensional errors ( AE ). Note the di↵erences smaller than the significance threshold defined in this paper are shown and that the colorbars are not the same between panels (a) and (b). Richardson extrapolation [Roache, 1994;Blain et al., 1998]. In order to use this approach to esti-714 mate numerical truncation error, it was first verified that the leading order error terms indeed con-715 trolled the numerical convergence (i.e., asymptotic regime), spatial errors were found to be much 716 greater than the time discretization errors, and the ADCIRC solver in the current configuration 717 was a second order accurate method in space and time [Luettich and Westerink, 2004]. The order 718 of convergence was verified to be 2nd order accurate by refining L250R1 once more producing 719 L250R2.

720
The Richardson extrapolation base error (REBE) following [Roache, 1994] is calculated to 721 estimate numerical error with the following formulas: 722 REBE[coarse mesh] = ✏r n (r n 1) wheref L250 andf L250R1 are the solutions computed on the original and refined meshes andf RE F 723 is the solution computed on the reference mesh. X L250 and X L250R1 denote the spatially varying 724 mesh sizes throughout the computational domain.

725
The REBE (herein the numerical error) for the L250 and L250R1 M 2 amplitude elevation is 728 presented in Figure 17c,d and compared against the total error that was calculated from the REF 729 solution using Equation 1 (i.e., RE) like was performed in the rest of the paper (Figure 17a,b).

730
There is a similarity in the numerical and total error estimates particularly in the NA subdomain 731 where the magnitude of both errors are 3-5% for the L250 mesh and diminish to 1-2% for the 732 L250R1 mesh. However, the estimate of the greatest magnitude numerical error is co-located with 733 the periphery of the Georges Bank near sharp seabed topographic gradients, while the total error is 734 spread across the entire Georges Bank. In general, a weaker reduction in the total error is observed 735 compared to the numerical error. In particular, the total error is not reduced over the Georges Bank 736 or along most of the SAB and MAB coastline (Figure 17a-b). However, the numerical error is reduced almost everywhere to below the significance threshold. For instance, the refinement of L250 738 to L250R1 reduces the numerical error estimate in the Chesapeake Bay estuary in the MAB re- Figure 17. An estimate of the numerical error calculated via Richardson extrapolation following [Roache, 1994] obtained by refining the L250 mesh using a four-to-one refinement strategy to preserve the approximate problem. The previously described mesh size functions (Table 1)  The idea behind this sequence of mesh combinations (COMBOx) is to proceed from a more 773 simple design and move towards a more complex design to test the additive e↵ects, i.e., start with 774 uniform shoreline resolution (COMBO1); use variable shoreline resolution (COMBO2); add addi-775 tional resolution along estuarine channels (COMBO3). COMBO1 begins with a high gradation rate was shown to reduce tidal error metrics as compared to both the reference solution and measured 897 data as inland waterway conveyances are improved and frictional resistance is reduced. We remark 898 that other mesh size heuristics, such as the slope mesh size function and using finer resolution 899 along the shoreline with a low gradation rate can implicitly, but ine ciently, capture these sub-900 marine channel features. Thus, the application of the estuarine channel mesh size function allows 901 the usage of a higher mesh size gradation so as to focus resolution only on the submarine channels 902 allowing us to more e ciently discretize the estuarine environment. 903 We tested the performance of mesh design strategies that involved using a steep mesh size derwater vertices, which is similar in number to our L250 mesh. In contrast, the COMBO3 mesh, 941 which spans the same ECGC study region, utilizes up to five times finer resolution nearshore (50 942 m compared to 250 m) and up to ten times finer resolution along the continental slope (1 km com-943 pared to 10 km), with only 1.6 times the total number of underwater vertices than HSOFS. 944 We highlight that an important first step in the coastal model development procedure is to 945 construct a mesh that minimizes the physical domain approximation error before model tuning oc-946 curs vis-a-vis varying bottom friction, other dissipative coe cients, viscous models, and manually 947 altering ocean depths and shoreline form. As was evident in this paper, by improving the accu-948 racy of the approximate problem (i.e., the representation of the shoreline and seabed topography of the ECGC in which the bottom friction coe cient are typically modified [Szpilka et al., 2016].