the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An automatic mesh generator for coupled 1D–2D hydrodynamic models
Younghun Kang
Ethan J. Kubatko
Twodimensional (2D), depthaveraged shallow water equation (SWE) models are routinely used to simulate flooding in coastal areas – areas that often include vast networks of channels and floodcontrol topographic features and/or structures, such as barrier islands and levees. Adequately resolving these features within the confines of a 2D model can be computationally expensive, which has led to coupling 2D simulation tools to less expensive onedimensional (1D) models. Under certain 1D–2D coupling approaches, this introduces internal constraints that must be considered in the generation of the 2D computational mesh used. In this paper, we further develop an existing automatic unstructured mesh generation tool for SWE models, ADMESH+, to sequentially (i) identify 1D constraints from the raw input data used in the mesh generation process, namely the digital elevation model (DEM) and land–water delineation data; (ii) distribute grid points along these internal constraints, according to feature curvature and userprescribed minimum grid spacing; and (iii) integrate these internal constraints into the 2D mesh size function and mesh generation processes. The developed techniques, which include a novel approach for determining the socalled medial axis of a polygon, are described in detail and demonstrated on three test cases, including two inland watersheds with vast networks of channels and a complex estuarine system on the Texas, USA, coast.
 Article
(24935 KB)  Fulltext XML
 BibTeX
 EndNote
Hydrodynamic models are routinely used to simulate, analyze, and assess the effects of physical phenomena that result in coastal flooding, such as tsunamis and hurricane storm surges. Typically, the twodimensional (2D), depthaveraged shallow water equations, equipped with a suitable wetting and drying algorithm, are used to model inundation within the coastal floodplain in the main (see, for example, Luettich and Westerink, 1999; Bunya et al., 2009; Dawson et al., 2011). However, these coastal regions often include vast networks of smallscale channels that, while playing a significant part in the conveyance of flood waters propagating into and through the floodplain, are often left underresolved in practice, due to the computational expense of adequately resolving them within the confines of a strictly 2D modeling approach. Moreover, in many situations, flow in open channels (and storm sewer systems in urban areas) can be adequately described by simpler onedimensional (1D), sectionaveraged flow equations (as demonstrated in, for example, Pramanik et al., 2010; Gichamo et al., 2012; Timbadiya et al., 2014; Bhuyian et al., 2015; Price, 2018). This situation has led to the development of a number of coupled 1D–2D modeling approaches over the past decade (see references below) that aim to more accurately and efficiently simulate flooding events, and coastal hydrodynamics in general, than traditional 2D modeling approaches.
These types of coupled models have been applied in a number of different hydrodynamic scenarios but have been most widely used to simulate river flooding problems in various settings (e.g., Liu et al., 2015; D'Alpaos and Defina, 2007; Kuiry et al., 2010; Martini et al., 2004; Marin and Monnier, 2009; Gejadze and Monnier, 2007; Timbadiya et al., 2015; Stelling and Verwey, 2006; Li et al., 2021; MoralesHernández et al., 2016). Some recent investigations have focused on inundation in urban areas, where storm sewer systems and channels, which are approximated in 1D, interact with a 2D flooding model (Vojinovic and Tutulic, 2009; Adeogun et al., 2012, 2015; Delelegn et al., 2011; Seyoum et al., 2012). Other coupled models have been developed and applied for coupled riverine–estuarine flows near coastal areas (Bakhtyar et al., 2020; Lin et al., 2006); river–lake flows (Chen et al., 2012; Pham Van et al., 2016); supercritical flow in crossroads (Ghostine et al., 2015); overland–open channel flows (West et al., 2017); and river closure projects (Lin et al., 2020). More examples of coupled models, including commercial models, are summarized elsewhere (Néelz and Pender, 2009; Teng et al., 2017; Woodhead et al., 2007).
The coupled models are often categorized based on the types of interactions that occur between the 1D–2D flows, but here we categorize them into three types, based on the way the 1D and 2D domains are connected (see Fig. 1). The first type is boundaryconnected domains (e.g., Chen et al., 2012; Liu et al., 2015; Bakhtyar et al., 2020; Ghostine et al., 2015; Pham Van et al., 2016). These types of domains are widely used for river–lake or river–estuary systems (Chen et al., 2012; Bakhtyar et al., 2020; Pham Van et al., 2016; Ghostine et al., 2015), where the longitudinal (or frontal) flows from 1D domains enter the 2D domains as a boundary condition. However, in some cases, the interaction can be made by lateral flow through breach between the 1D–2D domains (Liu et al., 2015). The second type is internally connected domains (e.g., West et al., 2017; Kuiry et al., 2010; D'Alpaos and Defina, 2007; Martini et al., 2004; Marin and Monnier, 2009; Gejadze and Monnier, 2007; Stelling and Verwey, 2006; Timbadiya et al., 2015; Vojinovic and Tutulic, 2009; Bunya et al., 2023). These types of domains are widely used for river–floodplain systems. The interaction is made along the whole 1D domain, where the discharge from 2D domains can enter the 1D domains or the 1D channel flows exceeding bank level can enter the 2D domain. The third type is vertically connected domains (e.g., Vojinovic and Tutulic, 2009; Fan et al., 2017; Adeogun et al., 2012, 2015; Delelegn et al., 2011). These types of domains are mostly used for urban inundation with storm sewer systems. The interaction is made at points where 1D–2D domains are connected vertically, where the surcharged overflow can enter the 2D domains.
Depending on the type of 1D–2D domain connections, the mesh generation can be straightforward or extremely complicated. The simplest case is boundaryconnected domains. In this case, the 2D computational meshes can be generated with automatic mesh generators that have been developed for 2D hydrodynamic models (see, for example, Persson and Strang, 2004; Conroy et al., 2012; Koko, 2015; Roberts et al., 2019; Engwirda, 2017; Hagen et al., 2002; Bilgili et al., 2006; Candy and Pietrzak, 2018; Avdis et al., 2018; Remacle and Lambrechts, 2018; Gorman et al., 2006). It would then remain to generate 1D meshes and to connect them to the 2D mesh at the boundary, which is straightforward once the 1D domain is generated. For the vertically connected domains, staggered 1D–2D computational meshes are widely used, i.e., 1D line elements are not constrained to be colocated with 2D mesh edges, given that the connections between the 1D and 2D domains are limited to points. For example, the meshes used in Delelegn et al. (2011); Adeogun et al. (2012, 2015) are generated with standard twodimensional mesh generators, such as Gaja3D (Rath, 2007) and Triangle (Shewchuk, 1996), without any consideration of the locations of links.
The most complicated case is the internally connected domains. Unlike vertically connected domains, collocated meshes, which here means that 1D line elements are aligned with 2D mesh edges, are desirable. The difficulty of mesh generation for this type comes from the fact that the resolutions of the 1D–2D domains are obviously closely intertwined with each other; however, the desired mesh resolutions for each domain may be quite different. Additionally, the 1D line elements must serve as internal constraints in the mesh generation process, which must be carefully identified and preprocessed to avoid overconstraining the 2D elements, leading to triangles of poor quality.
In this paper, we present an automatic mesh generator for internally connected 1D–2D hydrodynamic models that is an extension of an existing mesh generator for shallow water equation (SWE) models, ADMESH+ (Conroy et al., 2012), which built upon the ideas and methodology of Persson's DistMesh program – a simple, opensource mesh generator implemented in MATLAB (Persson and Strang, 2004). The ADMESH+ mesh generation process can be briefly outlined as follows. First, a mesh or element size function, h, is constructed so that is used to prescribe element sizes h(x) throughout a given 2D domain. These element sizes are based on a number of geometric factors, such as shoreline and/or boundary curvature and bathymetric/topographic gradients, as well as userdefined inputs, such as target maximum and minimum element sizes and mesh grading specifications (i.e., the ratio of neighboring elements should not exceed some specified factor). Given the element size function, a Delaunay triangulation of an initial set of mesh nodes with a density proportional to $\mathrm{1}/h(x{)}^{\mathrm{2}}$ is then generated and the nodes of this initial mesh are repositioned by solving for force equilibrium (iteratively) around each node, making use of a spring mechanics analogy (see Conroy et al., 2012, for complete details).
The primary improvements incorporated into ADMESH+ in this work are twofold. First, an automatic identification of 1D domains is developed. This requires, as input to the mesh generation process, a digital elevation model (DEM) and a socalled land and water mask that identifies the (initially) dry and wet (respectively) portions of the domain. Using these input datasets, separate methods of identification of narrow channels (i.e., those below a userspecified 2D minimum element size) that define 1D model domains over land and water are developed. These methods involve tightly integrating ADMESH+ with a MATLABbased topographic toolbox, TopoToolbox (Schwanghart and Kuhn, 2010), and identifying the medial axis of the water portion of the domain for which a novel approach is developed. Second, from the extracted 1D model domains, target 1D mesh node distributions are computed along smooth spline approximations of the 1D line segments according to channel curvature; i.e., mesh node density is increased in areas of high curvature and relaxed in straight segments, as well as the underlying 2D mesh size function computed during the standard mesh generation process. Similar to the 2D force equilibrium approach mentioned above, the actual 1D mesh node distribution is then determined from the target mesh size through the use of a 1D spring mechanics approach. The final internal constraint obtained is then used within the 2D mesh generation process to obtain meshes suitable for coupled (specifically, internally connected) 1D–2D hydrodynamic models.
The rest of this paper is organized as follows. In the next section, the framework of the proposed methodology and the primary input datasets of the mesh generation are discussed. Details of the algorithms developed to identify internal constraints are then described in Sect. 3, with illustrative examples for simple and complex geometries. The 1D and 2D mesh generation process with the identified internal constraints is then described in Sect. 4. Finally, the application of the methodologies developed are demonstrated in Sect. 5, and the paper is concluded by providing a brief summary of the work and some possible future directions are identified.
Consider a hydrodynamic model domain Ω_{2D}⊂ℝ^{2} defined by a simple (i.e., not selfintersecting) polygon, possibly with holes. A socalled mesh size function $h:{\mathrm{\Omega}}_{\text{2D}}\to \mathbb{R}$ that assigns a target element size $\mathrm{\Delta}=h(x,y)$ to each point $(x,y)\in {\mathrm{\Omega}}_{\text{2D}}$ plays a fundamental role in the construction of a triangulation 𝒯_{h} of Ω_{2D}. Here, the mesh size function h is represented as a bilinear interpolant constructed on a rectilinear background grid that consists of a set of points X, defined such that Conv(Ω_{2D})⊂X, where Conv(Ω_{2D}) denotes the convex hull of Ω_{2D}. In the existing ADMESH+ framework, several factors can be considered in the construction of the mesh size function, including userspecified minimum and maximum element sizes, boundary curvature (see Conroy et al., 2012, for a complete description).
In our previous work (Conroy et al., 2012), the polygonal domain Ω_{2D} and the mesh size function h (together with its corresponding background grid) have served as the two primary inputs to the mesh generation process of ADMESH+. In this work, the incorporation and construction of a third primary input is described that is fundamental to the generation of suitable meshes for coupled 1D–2D hydrodynamic models, namely a set of internal constraints that consists of a set of line segments interior to Ω_{2D}, along which the edges in the 2D triangulation must be constrained. Three types of internal constraints are considered. The first type is line segments that represent centerlines of narrow channels; i.e., those channels that cannot be accurately resolved under the userspecified minimum element size. These line segments make up the aforementioned 1D hydrodynamic domain. The second type is line segments that align with subgridscale topographic features/structures, such as narrow barrier islands, levees, and weirs, along which certain internal boundary conditions are enforced in the 2D hydrodynamic model (see, for example, Dawson et al., 2011). Finally, the third type of internal constraint is line segments that represent the boundary between the land and water subdomains, i.e., the shorelines. Note that this type of internal constraint is neither part of the 1D domain nor an internal boundary but is desired in order to provide a clear distinction between the land and water subdomains by the generated 2D mesh.
We note that the described internal constraints can be provided directly from relevant data sources; for example, channel centerlines are available from the U.S. Geological Survey (USGS) National Hydrography Dataset (NHD) (U.S. Geological Survey, 2016) but are more generally obtained based on two functions defined over Ω_{2D}. One is a function $f:{\mathrm{\Omega}}_{\text{2D}}\to \mathbb{R}$ that assigns the (bare) Earth surface elevation $z=f(x,y)$ to a point $(x,y)\in {\mathrm{\Omega}}_{\text{2D}}$. This is typically provided by a combination of topographic and bathymetric (gridded) DEMs. A second is a socalled indicator function ${\mathbb{1}}_{{\mathrm{\Omega}}_{\text{w}}}:{\mathrm{\Omega}}_{\text{2D}}\to \left\{\mathrm{0},\mathrm{1}\right\}$ of a subset Ω_{w}⊂Ω_{2D}, defined as
which indicates whether a point $(x,y)\in {\mathrm{\Omega}}_{\text{2D}}$ is (initially) wet (1) or dry (0), i.e., a socalled land–water mask. We refer to Ω_{w} as the water subdomain of Ω_{2D} (see Fig. 2). This is typically provided by a dataset of closed polygons whose interiors indicate the water subdomain (see, for example, Wessel and Smith, 1996). Below, we describe our methodology for identifying internal constraints from the datasets that inform these two functions.
In this section, we present our methodology for extracting the three types of internal constraints described above. The methodologies vary for the land and water portions of the domain. The key of the extraction of open channels from the land portion of the domain is the drainage area computed using the gradient of the input DEM(s). In the water portion of the domain, central to the extraction of open channels, is the width measurement of relevant features of the domain. Note that this methodology is further applied to identify internal boundaries (the second type of internal constraint) from the land portion of the domain. The third type of internal constraint – namely the boundaries between the land and water domains – is obtained as a byproduct from this methodology. As two distinct methodologies are applied, they are described in separate subsections below.
3.1 Open channels from the land subdomain
In the land subdomain, only the first type of internal constraint mentioned above, namely channel centerlines, are identified. From input DEM(s), internal constraints representing dry channels (flow paths) in the land subdomain are detected by integrating the existing ADMESH+ code with TopoToolbox – a widely used, MATLABbased topographic toolbox (Schwanghart and Kuhn, 2010).
The channel detection algorithm in TopoToolbox is a flow accumulation algorithm based on the gradient of the DEM. At each grid point of the DEM, flow direction is computed for eight connected neighbors to create a global flow direction matrix M, which is defined by
where z_{i} is elevation at cell i, and d_{ij} is the distance between cells i and j. Starting with a uniform unit water depth on each grid point of the DEM, expressed as a vector w, the global flow direction matrix is multiplied with the water depth vector to compute water depth at the next time step
where w^{k} is water depth vector at kth iteration. This operation is repeated until water completely leaves the system. Then, drainage area at each grid point is defined as the summation of water depth over all iterations:
where a is drainage area vector. Finally, the channels are identified by a threshold of minimum drainage area. Note that while the algorithm is developed based on an iteration method, in TopoToolbox, the drainage area is computed directly by using geometric procession:
Note that the quality of extracted open channels using the flow accumulation algorithm depends on the quality of the DEM. For example, if there is a sink in the midst of an open channel, then the numerical flow computed using the given algorithm cannot propagate at the sink, which results in an erroneous channel network. In order to improve the quality of openchannel extraction, sinks are filled as preprocess using the fillsinks
function provided in TopoToolbox.
As mentioned above, while other datasets of openchannel centerlines (for example, USGS NHD) can be used, in our experience, open channels extracted from TopoToolbox show better alignment with input DEMs.
3.2 Open channels from the water subdomain and internal boundaries from the land subdomain
The water subdomain exhibits a wide range of scales from large open waterbodies to narrow channels and islands, the latter of which are represented as holes in the water mask. For the purpose of identifying narrow channels and small islands in the water subdomain, we describe our methodology for performing a widthbased decomposition of the water subdomain that makes use of a userdefined minimum width δ_{w}; i.e., features of width <δ_{w} are identified as narrow.
3.2.1 Widthbased mask decomposition
We describe the stepbystep procedure for determining the widthbased mask decomposition with the aid of the simple example shown in Fig. 3. For masks with more complex boundaries and/or holes, additional processes are required, which are described in Sect. 3.2.2.

Step 1. For a given polygon (Fig. 3a), find the medial axis (Fig. 3b).
The medial axis MA(P), of any closed polygon P, is defined as the set of interior points that have equal distances to two or more points on the boundary of P. Our computation of the medial axis is based on the vector distance transform (VDT) originally proposed by Mullikin (1992). Given a polygon P with the boundary ∂P, the value of the VDT, V(x,P), at point $\mathit{x}=(x,y)\in P$ is defined as
$$\begin{array}{}\text{(7)}& \mathit{V}(\mathit{x},P)=\mathbf{\Phi}(\mathit{x},\partial P)\mathit{x},\end{array}$$where
$$\begin{array}{}\text{(8)}& \mathbf{\Phi}(\mathit{x},\partial P)=\underset{\mathit{y}\in \partial P}{\mathrm{arg}\phantom{\rule{0.125em}{0ex}}inf}\phantom{\rule{0.25em}{0ex}}\left\right\mathit{x}\mathit{y}\left\right.\end{array}$$Using the VDT, the medial axis of P can be obtained as (see Appendix A)
$$\begin{array}{}\text{(9)}& \mathrm{MA}\left(P\right)=\mathit{\{}\mathit{x}\in P:\mathrm{\nabla}\cdot \mathit{V}(\mathit{x},P)>\mathrm{0}\mathit{\}}.\end{array}$$Note that the VDT defined by Eq. (7) is not well defined on the medial axis points, as multiple arguments of the infimum are produced from Eq. (8). In numerical implementation, we arbitrarily choose one of the arguments of infimum of Eq. (8).
Finally, this algorithm produces medial axis points, which is a subset of background grid points along the medial axis. However, the following methodology is implemented based on medial axis branches, which consists of line segments. For details of the construction of the medial axis branches from medial axis points, see Appendix A.

Step 2. Prune the medial axis (Fig. 3b).
After obtaining the medial axis, pruning near corners is required to improve the quality of the width function, which is described later. First, we construct a hierarchy of medial axis branches; i.e., branches with free ends are order 1, branches connected to order 1 are order 2, etc. Next, for each medial axis point p on order 1 branches, we define width and angle of medial axis as
$$\begin{array}{}\text{(10)}& \mathrm{\ell}=\underset{i=\mathrm{1},\mathrm{\dots},\mathrm{4}}{max}\left\right(\mathit{p}+\mathit{v})({\mathit{p}}_{i}+{\mathit{v}}_{i})\left\right\end{array}$$and
$$\begin{array}{}\text{(11)}& \mathit{\theta}=\underset{i=\mathrm{1},\mathrm{\dots},\mathrm{4}}{max}{\mathrm{cos}}^{\mathrm{1}}\left({\displaystyle \frac{\mathit{v}\cdot {\mathit{v}}_{i}}{\left\right\mathit{v}\left\right\phantom{\rule{0.25em}{0ex}}\left\right{\mathit{v}}_{i}\left\right}}\right),\end{array}$$where ${\mathit{p}}_{i}\phantom{\rule{0.125em}{0ex}}(i=\mathrm{1},\mathrm{\dots},\mathrm{4})$ are four neighboring background grid points of p, and v and v_{i} are VDTs corresponding to p and p_{i}; i.e., $\mathit{v}=\mathit{V}(\mathit{p},P)$ and ${\mathit{v}}_{i}=\mathit{V}({\mathit{p}}_{i},P)$ (see Fig. 4). As mentioned in the previous step, the VDTs on the medial axis are arbitrarily chosen in numerical implementation. This may cause differences for lengths of the medial axis, but these differences are insignificant (see Fig. 4). Note that inverse cosine function computed with
acos
function in MATLAB gives ${\mathit{\theta}}_{i}\in [\mathrm{0},\mathit{\pi}]$. Finally, the medial axis points near corners are identified and pruned, based on the thresholds$$\begin{array}{}\text{(12)}& \mathit{\theta}<{\mathit{\delta}}_{\mathit{\theta}}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}\mathrm{\ell}<{\mathit{\delta}}_{\mathrm{\ell}},\end{array}$$where δ_{θ} and δ_{ℓ} are specified threshold values used for pruning. Based on our experiments, the values δ_{θ}=0.9π and δ_{ℓ}=2δ_{w} provide reasonable results.

Step 3. Compute distance functions and define the width function (Fig. 3c, d, and e).
Now we define two distance functions; one measures the closest distance from a point x to the boundary $\mathrm{d}(\mathit{x},\partial P)$, and the other measures the closest distance from a point to the x medial axis d(x,MA(P)). The former can be measured using the VDT:
$$\begin{array}{}\text{(13)}& \mathrm{d}(\mathit{x},\partial P)=\left\right\mathit{V}(\mathit{x},P)\left\right.\end{array}$$Then, the width function is defined as twice the sum of the two distance functions (see Fig. 5); i.e.,
$$\begin{array}{}\text{(14)}& {f}_{\mathrm{w}}(\mathit{x},P):=\mathrm{2}\left(\mathrm{d}\right(\mathit{x},\partial P)+\mathrm{d}(\mathit{x},\mathrm{MA}\left(P\right)\left)\right).\end{array}$$Note that, without pruning the medial axis, the width function results in f_{w}(x)≈0 in the vicinity of every corner. This is because both boundary lines and the (unpruned) medial axis exist at every corner of the polygon boundary, and thus $\mathrm{d}(\mathit{x},\partial P)\approx \mathrm{0}$ and $\mathrm{d}(\mathit{x},\mathrm{MA}(P\left)\right)\approx \mathrm{0}$. This results in the area around every corner being identified as narrow, even though they are corners of large waterbodies. This is the rational for pruning the medial axis around corners, which was presented in Step 2.

Step 4. Decompose the polygon with a userdefined minimum width δ_{w} (Fig. 3f).
With the width function Eq. (14), the given polygon P can be decomposed as
$$\begin{array}{}\text{(15)}& {P}_{\mathrm{1}}:=\mathit{\{}\mathit{x}\in P:{f}_{\mathrm{w}}(\mathit{x},P)<{\mathit{\delta}}_{\mathrm{w}}\mathit{\}}\end{array}$$and
$$\begin{array}{}\text{(16)}& {P}_{\mathrm{2}}:=\mathit{\{}\mathit{x}\in P:{f}_{\mathrm{w}}(\mathit{x},P)\ge {\mathit{\delta}}_{\mathrm{w}}\mathit{\}}.\end{array}$$We call the masks P_{1} and P_{2} level 1 and level 2 masks, respectively. Note that the level 1 and level 2 masks can be used for 1D and 2D domains for simple cases, but additional processes are required for complex geometries. The level 1 and level 2 masks are not identical to the 1D and 2D domain in that case, and thus we use the term “level” to avoid confusion.
3.2.2 Additional processes for complex geometries
The challenge in real applications comes from complex boundaries and the presence of small narrow islands within large waterbodies. The complexities result in the noise of the width function and do not provide a clear distinction between level 1 and level 2 masks, and additional processes are applied to resolve the noise. In order to catch the narrow islands, the mask decomposition needs to be applied to the land mask in a similar way to that of the water mask decomposition. Applying the mask decomposition for both land and water masks also requires an additional process, which is described in this section. Hereinafter, the land and water masks are denoted by M_{L} and M_{W}, respectively. These processes are described below in the context of an example (see Figs. 6 and 7).

Step 1. Remove small islands (Fig. 6a and b).
Small islands, by which we mean islands with areas that are much smaller than the square of minimum element size, are removed at the first step for two reasons. First, it is expected that small islands will not have significant effects on the hydrodynamic models. Second, mask decomposition with small islands tends to result in poor internal constraints. Mask decomposition converts small islands to internal constraints with the length of their major axis, and therefore, small islands tend to create internal constraints shorter than the minimum element size. Note that short internal constraints can be identified and removed after decomposing masks as well. In this example, islands with area less than 1000 m^{2} are removed, where the minimum element size is set as 45 m (see Fig. 6a and b).

Step 2. Apply widthbased decomposition with filling to land and water masks (Fig. 6e and f).
While the mask decomposition gives a clear distinction between narrow and wide regions in the simple example (see Fig. 3), more complex cases can result in noise in the width function, which can be removed through a socalled filling method.
The filling method is based on the idea of the maximal disk (see Fig. 8a), which is defined as
$$\begin{array}{}\text{(17)}& D(\mathit{x},P)=\mathit{\{}\mathit{y}\in P:\mathit{x}\mathit{y}\le \mathrm{d}(\mathit{x},\partial P\left)\mathit{\right\}}.\end{array}$$And the level 1 and level 2 masks are updated with the maximal disks centered in level 2 masks
$$\begin{array}{}\text{(18)}& {M}_{\stackrel{\mathrm{\u0303}}{\text{i2}}}:={M}_{\text{i2}}\cup {\mathit{D}}_{i}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{M}_{\stackrel{\mathrm{\u0303}}{\text{i1}}}:={M}_{\text{i1}}\backslash {\mathit{D}}_{i},\end{array}$$where
$$\begin{array}{}\text{(19)}& {\mathit{D}}_{i}=\mathit{\left\{}D\right({\mathit{x}}_{j},{M}_{\text{i}})\text{for all}{\mathit{x}}_{j}\in \mathrm{MA}({M}_{\text{i}})\cap {M}_{\text{i2}}\mathit{\}},\end{array}$$and M_{i}=M_{L} and M_{W}. Note that the maximal disks are sought on level 2 masks only. Also, there are some regions of ${M}_{\stackrel{\mathrm{\u0303}}{\text{i1}}}$ that do not include any medial axis (see Fig. 8b). These regions are redundant, as they cannot be represented by the medial axis. Therefore, we transfer such regions to level 2 masks and update the masks (see Fig. 8c).
$$\begin{array}{}\text{(20)}& {M}_{\stackrel{\mathrm{\u203e}}{\text{i2}}}:={M}_{\stackrel{\mathrm{\u0303}}{\text{i1}}}\cup {\mathit{R}}_{i}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{M}_{\stackrel{\mathrm{\u203e}}{\text{i1}}}:={M}_{\stackrel{\mathrm{\u0303}}{\text{i1}}}\backslash {\mathit{R}}_{i},\end{array}$$where
$$\begin{array}{}\text{(21)}& {\mathit{R}}_{i}=\mathit{\{}P\in {M}_{\stackrel{\mathrm{\u0303}}{\text{i1}}}:P\cap \mathrm{MA}({M}_{\text{i}})=\varnothing \mathit{\}}.\end{array}$$Note that with these procedures, the level 2 masks have a smoother boundary and provide a better representation of the large waterbodies; see Fig. 6e and f, which are updated masks of panels Fig. 6c and d, respectively, after filling.

Step 3. Find regions of a level 1 land mask that are surrounded by a level 2 water mask (Fig. 7g).
The purpose of this step is to identify regions of a level 1 land mask that are retained as internal constraints in the 2D hydrodynamic model. These narrow land regions are typically modeled as internal boundaries over which simple subgridscale flow parameterizations are performed, using, for example, a simple weirbased formula (see, for example, Dawson et al., 2011).
First, the thinness of land regions is determined by the socalled isoperimetric ratio (IPR), which is defined as ${\text{perimeter}}^{\mathrm{2}}/\text{area}$. Note that the IPR is a dimensionless number, which is higher for thinner regions. Here, we set a threshold of 30. Second, in order to identify if each region is surrounded by level 2 water mask, we first set a buffer for each region. The buffer size is set as half of the length of the minor axis of the ellipse that has the same normalized second central moments as the region (see Fig. 9a and b). The minor axis length is computed using the
regionprops
function in MATLAB. We then check the area of the level 2 masks within the buffer. We define the land region to be surrounded by level 2 water mask if the area of level 2 water mask is greater than twice the area of level 1 land mask within the buffer (see Fig. 9c and d, which show selected regions of level 1 land mask shown in Fig. 9a and b, respectively). Note that level 1 water masks will be internal constraints (channel centerlines). This means that if a region of the level 1 land mask is surrounded by a level 1 water mask, then there are too many internal constraints too close to each other. Therefore, we only retain regions of the level 1 land mask that are surrounded by a level 2 water mask, denoted by M_{L1W2}, which can be expressed as$$\begin{array}{}\text{(22)}& \begin{array}{rl}{M}_{\text{L1W2}}& =\mathit{\{}P\in {M}_{\stackrel{\mathrm{\u203e}}{\text{L1}}}:\text{IPR}(P)>\mathrm{30}\\ & \text{and}{\displaystyle \frac{\left\left{M}_{\stackrel{\mathrm{\u203e}}{\text{W2}}}\cap B\left(P\right)\right\right}{\left\left{M}_{\stackrel{\mathrm{\u203e}}{\text{L2}}}\cap B\left(P\right)\right\right}}\mathrm{2}\mathit{\}},\end{array}\end{array}$$where IPR(P) is the isoperimetric ratio, B(P) is the buffer of region $P\in {M}_{\stackrel{\mathrm{\u203e}}{\text{L1}}}$, and ⋅ denotes area.

Step 4. Transfer regions of the level 1 land mask identified in Step 3 to the water mask (Fig. 7h).
As described in Step 3, the narrow regions identified will now serve as internal constraints (specifically, internal boundaries as described above) in the mesh generation process and no longer need to be in land mask. Therefore, these regions are transferred into the water mask and updated land and water masks are defined as
$$\begin{array}{}\text{(23)}& {\displaystyle}{\stackrel{\mathrm{\u203e}}{M}}_{\text{L}}& {\displaystyle}:={M}_{\text{L}}\backslash {M}_{\text{L1W2}}\text{(24)}& {\displaystyle}{\stackrel{\mathrm{\u203e}}{M}}_{\text{W}}& {\displaystyle}:={M}_{\text{W}}\cup {M}_{\text{L1W2}}.\end{array}$$Note that the centerlines of narrow land regions are used as internal constraints (internal boundaries) for mesh generation:
$$\begin{array}{}\text{(25)}& {\mathit{s}}_{\mathrm{1}}=\mathrm{MA}\left({M}_{\text{L}}\right)\cap {M}_{\text{L1W2}}.\end{array}$$ 
Step 5. Apply the widthbased mask decomposition to the updated water mask (Fig. 7i).
The widthbased decomposition of the water mask that now includes the narrow land regions described above will be different from the widthbased decomposition of the original input water mask. Therefore, the widthbased decomposition must be applied again. Note that, in this step, the decomposition for the land mask is not required. The centerlines of updated level 1 water mask ${\stackrel{\mathrm{\u203e}}{M}}_{\text{W1}}$ form the domain of a 1D hydrodynamic model; i.e.,
$$\begin{array}{}\text{(26)}& {\mathrm{\Omega}}_{\text{1D}}:=\mathrm{MA}\left({\stackrel{\mathrm{\u203e}}{M}}_{\text{W}}\right)\cap {\stackrel{\mathrm{\u203e}}{M}}_{\text{W1}}.\end{array}$$Here, we have the second type of internal constraint (1D domain) for mesh generation
$$\begin{array}{}\text{(27)}& {\mathit{s}}_{\mathrm{2}}={\mathrm{\Omega}}_{\text{1D}}.\end{array}$$The domain of the 2D model is the entire domain, except the domain of the 1D model. Note that it includes the updated level 1 water mask, as well as the updated level 2 water mask and land mask; i.e.,
$$\begin{array}{}\text{(28)}& {\mathrm{\Omega}}_{\text{2D}}:={\stackrel{\mathrm{\u203e}}{M}}_{\text{L}}\cup {\stackrel{\mathrm{\u203e}}{M}}_{\text{W2}}\cup ({\stackrel{\mathrm{\u203e}}{M}}_{\text{W1}}\backslash {\mathrm{\Omega}}_{\text{1D}}).\end{array}$$It is preferred that the mesh elements are aligned along boundaries between waterbodies and land. This can be ensured by passing the boundary of ${\stackrel{\mathrm{\u203e}}{M}}_{\text{W2}}$ as internal constraints, even though it is neither an open channel nor an internal boundary, which is the third type of internal constraint for mesh generation:
$$\begin{array}{}\text{(29)}& {\mathit{s}}_{\mathrm{3}}=\partial {\stackrel{\mathrm{\u203e}}{M}}_{\text{W2}}.\end{array}$$ 
Step 6. Construct mainstreams of 1D domain and internal boundaries.
The 1D domains and internal boundaries contain the centerlines of narrow regions of the water and land mask. These centerlines are obtained from the medial axis of the mask, which has been pruned and ordered into hierarchic medial axis branches. This branchwise ordering creates several short medial axis segments (see, for example, Fig. 10), which is undesirable for computing the internal constraint curvature that is used to help determine the size of the 1D elements, as described in the next section. Therefore, a procedure to construct channel mainstreams is applied as follows. First, at each joint, the pair of segments, or branches, that have minimum (absolute) curvature are merged together to form a new segment (the mainstream). If there are more than three branches at a joint, another pair of branches with minimum curvature, excluding the mainstream, is merged and set as a substream. The substreams are collected until there is no pair at a joint.
To conclude this section, the algorithm described generally provides a good identification of narrow regions and the three types of internal constraints. However, this identification is not perfect and can result in some narrow regions being falsely identified, or misclassified, as level 2 regions. There are two possible reasons for this misclassification. The first one is related to the quality of the medial axis calculation. Recall that the medial axis is obtained by computing the divergence of the VDT and is subsequently pruned based on specified tolerances. This approach can result in an inaccurate estimation of the width function, especially within small regions. The second reason relates to the use of the dimensionless width parameter. Specifically, the identification of level 1 land regions surrounded by level 2 water masks (as described in Step 3) depends on the IPR threshold and the ratio of the surrounding level 2 water and land mask. This can result in some level 1 land regions being omitted, which we want to see represented by their centerlines. Note that a higher quality of the medial axis can be achieved by using other algorithms (see, for example, Lee, 1982) or using an adaptive background grid such as octree (Yerry and Shephard, 1983). In a simpler way, a higher resolution of the (uniform Cartesian) background grid can be utilized, which will cause memory and computational inefficiencies. In order to further enhance the quality of classification, a complex algorithm of criteria choice may be required, such as machine learning.
Additionally, there can be some internal constraints, even if they are correctly identified, that result in elements that are too small or of poor quality. For example, there can be internal constraints that are too close to each other. Note that, by choosing level 1 land regions surrounded by a level 2 water mask (in Step 3), it is unlikely that internal boundaries are too close to open channels. However, there are some internal boundaries/open channels which are too close to the third internal constraint, namely boundaries between waterbodies and land.
While these problems are related to identification of narrow and wide regions and internal constraints, it is easier to resolve them during the mesh generation process itself. Therefore, treatment of the problems mentioned here will be described in the next section (specifically, see Sect. 4.4).
The ADMESH+ is built using the force equilibrium algorithm, originally proposed by Persson and Strang (Persson and Strang, 2004). The basic idea of the force equilibrium algorithm is to model each element edge as a spring. It assigns a repelling force to each element edge; for example, a force of the form $f=k({l}_{h}{l}_{i})$, where k is a spring constant, l_{i} is the length of the element edge in the current triangulation 𝒯_{i}, and l_{h} is the desired length determined by the target mesh size h. Total repelling forces for each node are computed and applied to reposition nodes. The algorithm iterates this process and seeks to equilibrium state. For complete details, see Persson and Strang (2004).
In this section, we introduce a methodology to generate 2D finite element meshes given the identified internal constraints. The goal of mesh generation with internal constraints is to use efficient mesh resolution along the internal constraints so that it preserves the geography of the study areas with reduced computational demand. This is achieved by assigning a mesh size inversely proportional to the curvature of internal constraints. Also, additional processes are applied to ensure robustness of the force equilibrium algorithm with internal constraints.
4.1 Smoothing internal constraints
The internal constraints, which are given by user input or extracted from a DEM, are often based on a structured grid. This can produce, for example, a choppy set of channel centerlines (see Fig. 11) that do not provide a good basis for computing channel curvature, which is used in the process of determining mesh node placement and element size.
Therefore, in order to provide a smooth curve from which we can compute curvature, a cubic spline smoothing is applied based on the csaps
function in MATLAB. The csaps
function returns a smooth spline interpolation f_{p} to the N data points (r_{i},y(r_{i})), $i=\mathrm{1},\mathrm{\dots},N$ that minimizes
where p is a smoothing parameter, D^{2}f_{p} is the second derivative of f_{p}, and w_{i} and λ denote error measure weights and a weight function, respectively (see the MATLAB csaps
help documentation for more details). Note that internal constraints are twodimensional curves $\mathit{s}=({x}_{i},{y}_{i})$. Thus, for each internal constraint, we define a parametric curve as (x(r_{i}),y(r_{i})), where r_{i}=i for $i=\mathrm{1},\mathrm{\dots},N$, and find smoothed curves ${\stackrel{\mathrm{\u0303}}{x}}_{p}$ and ${\stackrel{\mathrm{\u0303}}{y}}_{p}$ individually. Note that the parameter r_{i} is a set of arbitrary values, and the csaps
function depends not only on the smoothing parameter p but also on the parameter r_{i}. In order to get standardized smoothing, we find a smoothing parameter with root mean square error (RMSE) closest to userdefined RMSE:
Then a smoothed curve is defined as
Note that the userdefined RMSE (RMSE_{desired}) should be carefully selected. If it is too high, then the smoothed curve is close to a straight line. If it is too low, then it does not give enough smoothing. In our numerical experiments, a smoothing with RMSE_{desired} between 1 and 10 m is generally appropriate in representing the overall curvature of the 1D constraints (see Fig. 11). It is noteworthy that RMSE_{desired} may be determined based on the resolution of the DEM.
4.2 Initial target mesh size and 2D gradient limiting
Given a smoothed internal constraint segment ${\stackrel{\mathrm{\u0303}}{\mathit{s}}}_{i}$, initial target mesh sizes along the curve are computed by
where K is the number of elements per radian (a userdefined parameter) and $\mathit{\kappa}\left({\stackrel{\mathrm{\u0303}}{\mathit{s}}}_{i}\right(r\left)\right)$ is the curvature of smoothed internal constraints ${\stackrel{\mathrm{\u0303}}{\mathit{s}}}_{i}\left(r\right)$.
It is often desired to ensure that the 2D element sizes grade properly in the final mesh. Common approaches to achieving this are marching methods (e.g., Persson, 2006; Roberts et al., 2019) and gradient limiting (e.g., Conroy et al., 2012; Persson, 2006). In this paper, we adopt the gradient limiting approach. Briefly, we find a steadystate solution to the socalled gradient limiting equation:
where g is related to a userdefined parameter that controls the ratio of neighboring element sizes (see Conroy et al., 2012 for details). The gradient limiting equation is solved with the initial condition
The gradient limiting is applied with the initial mesh size, h_{0}(x), on internal constraints being defined by Eq. (33). Note that a onedimensional gradient limiting equation could be solved for each internal constraint. However, this may result in inappropriate target mesh sizes if internal constraints are close to each other (see Fig. 12).
4.3 Generation of 1D meshes on internal constraints
In order to generate collocated 1D elements and 2D edges, 1D meshes are first generated on internal constraints and then used as fixed points in the 2D mesh generation. The target element sizes of the 1D meshes on each internal constraint are defined by projecting the gradientlimited mesh size h_{2D}, which is the solution of Eq. (34) with ${\stackrel{\mathrm{\u0303}}{h}}_{\text{2D}}$ as the initial condition. That is,
Then, applying 1D force equilibrium with the target size on each internal constraint provides 1D nodes ${r}_{i},i=\mathrm{1},\mathrm{\dots},N$. Now, fixed points of the 2D mesh generation are defined by
Note that the positioning of the fixed points is based on the original internal constraints s(r) instead of smoothed curves $\stackrel{\mathrm{\u0303}}{\mathit{s}}\left(r\right)$ because smoothing can result in deviations from the original set of points defining the internal constraints, as described in Sect. 4.1. Also, it is necessary to keep junctions of the curves so that the physical connections are not missed. This can simply be ensured by using junction points as fixed points of the 1D force equilibrium.
4.4 Postprocesses for 1D mesh
As noted in Sect. 3.2.2, postprocesses for the generated 1D meshes are applied to improve identification of internal constraint types and to improve 2D mesh quality.
Note that regions of the level 2 land mask M_{L2} and the updated level 2 water mask ${\stackrel{\mathrm{\u203e}}{M}}_{\text{W2}}$ represent wide regions that will be represented with 2D elements. Such regions are required to completely include at least one 2D element. However, for the regions falsely identified as level 2, their perimeter is not long enough to have three or more edges. This is likely to happen for regions with small areas that are round in shape. The IPR is lower for these types of areas, which results in those regions not being selected with the IPR filter when level 1 land regions surrounded by level 2 water mask are identified (see Eq. 22). Again, the boundaries of the level 2 regions are the third type of internal constraint. When 1D meshes are generated along the internal constraints corresponding to the falsely identified regions, then there is only one element per region, and such 1D elements are transferred to the first or second type of internal constraint.
As also noted in Sect. 3.2.2, there can be some of the first and second type of internal constraints, open channels and internal boundaries, that are too close to the third type of internal constraints, boundaries between waterbodies and land. This will result in 2D elements that are too small between such internal constraints. Therefore, if the all 1D mesh nodes on an internal constraint are closer than ${h}_{\mathrm{min}}/\mathrm{2}$ to any other 1D mesh on the third type of internal constraint – the boundaries between the land and water domains – then the 1D meshes on the other two types of internal constraint are removed.
As we fixed junction points to preserve physical connections, there can be 1D mesh node clusters, which are sets of 1D mesh nodes located within ${h}_{\mathrm{min}}/\mathrm{4}$ each other. Since this will result in 2D elements that are too small, the 1D mesh nodes are merged into the centroid of the cluster.
Finally, note that the mesh size of 1D elements (the distance between the 1D nodes in Eq. 37) is not identical to the target mesh size provided by Eq. (36). The mesh generated from the force equilibrium algorithm has mesh sizes relative to target mesh sizes (see Persson and Strang, 2004, for details). Due to the nature of force equilibrium algorithm, there can be 1D elements for which the length is shorter than ${h}_{\mathrm{min}}/\mathrm{2}$. In order to obtain highquality 2D elements, such 1D elements are merged into their neighboring 1D elements.
4.5 Generation of 2D meshes with fixed points
The 2D force equilibrium is applied after 1D meshes are generated by adopting 1D mesh nodes as fixed points. Again, note that the mesh size of the 1D mesh (the distance between the 1D nodes in Eq. 37) is not identical to the target mesh size Eq. (36). Due to the discrepancy between the target mesh size (of the 2D mesh) and the actual mesh size (of the 1D mesh), there are some nonconverging nodes near internal constraints. Note that, in general, the displacement of nodes decreases during the force equilibrium iterations and nodes converge to their final node locations. However, the nonconverging nodes keep moving back and forth near internal constraints, and an additional treatment is applied to resolve the nonconverging node situation.
This treatment consists of the following steps (see Fig. 13): (1) find nodes near the internal constraints that are not fixed points but have distances to the internal constraints that are shorter than ${h}_{\mathrm{min}}/\mathrm{2}$. (2) Compute the length of the closest 1D element to this node. (3) If the length of closest element is greater than 2h_{min}, then add the node to a new fixed point in the middle of the 1D element (Fig. 13a1 and a2). Otherwise, remove the node (Fig. 13b1 and b2).
This density control is applied while the 2D force equilibrium is being applied. However, note that there might be a number of nodes near the internal constraints which are located within ${h}_{\mathrm{min}}/\mathrm{2}$ at early stages of the force equilibrium process. Therefore, the density control is applied for later stages of the force equilibrium, which starts from 0.8×maximum iteration number.
The mesh generation algorithm presented in this paper is applied to three test cases. The first two test cases are for inland watersheds without water subdomains. The third test case is applied for a coastal basin to highlight the performance of identification of 1D domains in water subdomains.
There are several measures used to assess the quality of a mesh (see Field, 2000). The measure used in this paper is twice the ratio of the inradius r and circumradius R of each triangular element, i.e.,
where a, b, and c are the edge lengths of the triangular element. Note that this measure gives q=1 for an equilateral triangle and q=0 for a completely degenerate triangle (see Fig. 14).
Since mathematical error bounds for numerical methods are influenced by the smallest angle in the mesh, it is desired that a mesh consists of nearly equilateral triangles (Field, 2000). Meshes created with ADMESH+ typically have a mean element quality measure of q≥0.90 and a minimum element quality measure of q>0.30, where q=0.30 corresponds to an element with a minimum angle of 20° (Conroy et al., 2012).
5.1 The Middle Bosque River watershed
The Middle Bosque River watershed (MBRW), located in central Texas, has been the subject of numerous computational hydrological studies (see, for example, Bailey et al., 2021; Park et al., 2019; Tefera and Ray, 2023). Given this interest and the complex network of channels that must be represented for accurate model studies (see Fig. 15), the MBRW presents an ideal test case for the developed 1D–2D mesh generation process described.
The MBRW covers an area of approximately 516 km^{2} within the much larger Brazos River basin (≈119 174 km^{2}) – the secondlargest river basin by area within Texas. The boundary of the MBRW is obtained from the USGS Watershed Boundary Dataset (WBD) (U.S. Geological Survey, 2014), which provides the input (polygonal) domain Ω_{2D} for the mesh generation process. In addition to this input, a DEM covering the MBRW is available from the USGS 3D Elevation Program (3DEP) (U.S. Geological Survey, 2017), with the highest resolution available for the whole watershed being $\mathrm{1}/\mathrm{3}$ arcsec (approximately 10 m). While channel centerlines are available for the MBRW from the aforementioned USGS National Hydrography Dataset (NHD), in this test case, the channel centerlines, which constitute the 1D domain Ω_{1D}, are extracted within our mesh generator using TopoToolbox, as described in Sect. 3.1. The MBRW domain boundary, DEM, and extracted channels are shown in Fig. 15.
Given the inputs described above, the mesh is automatically generated using the procedure outlined with the following userdefined parameters: minimum and maximum element sizes are set to 30 and 500 m, respectively; the number of elements per radian K in Eq. (33) is set to 20; the grading limit g in Eq. (34) is set to 0.15; and the smoothing RMSE in Eq. (31) is set to 10 m. The resulting mesh is shown in Figs. 16, 17, and 18, where several qualities of the mesh can be visually noted. First, the dashed blue lines of panels (a1) and (b1) of Fig. 16 show closeups of the smooth spline approximation of the channel centerlines that have been extracted from the input DEM. The accompanying panels (a2) and (b2) of Fig. 16 show the node distribution of the 1D mesh that is generated along these channels, where it can be noted that smaller element sizes are present in highly curved areas and where elements sizes are relaxed in straighter channel segments (this is also visible, perhaps more so, in the zoomins of Fig. 18). Given these 1D channel elements, the 2D mesh is then generated and postprocessed as described; see Figs. 17 and 18, where it can be noted that the 2D elements of the generated mesh are constrained along the channel centerlines, grading out to larger 2D element sizes away from the channels, all while maintaining high quality. Specifically, the generated 2D mesh has a mean element quality of q=0.97, with approximately 99 % of the 183 610 elements of the mesh having a quality of q>0.83 and only 18 elements (corresponding to 0.01 % of the elements) having a quality of q<0.50, with the minimum element quality being q=0.33. While the minimum and maximum element sizes are set to 30 and 500 m, the minimum and maximum element sizes of the mesh are 18 and 572 m. Note that, again, the mesh generated from the force equilibrium algorithm has mesh sizes relative to the target mesh size. The mesh was generated in 13.16 min.
5.2 The Walnut Gulch Experimental Watershed
The Walnut Gulch Experimental Watershed (WGEW), established in southeastern Arizona in the 1950s, operates as an outdoor laboratory for studying hydrologic and erosion processes. Over the years, an extensive database of precipitation, runoff, and sediment records has been collected (Renard et al., 2008; Goodrich et al., 2008; Stone et al., 2008), making it, like the previous test case, the subject of numerous studies (see, for example, Meng et al., 2008; Goodrich et al., 2012; Yu and Duan, 2017) and an ideal test case for the developed mesh generator.
The WGEW covers approximately 149 km^{2} in Cochise County in southeastern Arizona. As with the previous test case, the boundary of the domain is obtained from USGS WBD (U.S. Geological Survey, 2014), which provides the input (polygonal) domain Ω_{2D} for the mesh generation process. Additionally, for the WGEW, both finescale (1 m) and coarsescale ($\mathrm{1}/\mathrm{3}$ arcsec) DEMs are available through the USGS 3DEP. The NHD dataset for channel centerlines is also available, but as with the previous test case, TopoToolbox is used to extract open channels. The top panel of Fig. 19 shows the domain boundary, the DEM, and the extracted channel networks of the WGEW.
Given these inputs, the mesh is generated with the following userdefined parameters: minimum and maximum mesh sizes of 30 and 500 m, respectively; the number of elements per radian K is 20; the grading limit is 0.15; and the smoothing RMSE is set to 10 m. The mesh that is generated is shown in Figs. 19 and 20. Like the previous test case, the 2D elements of the generated mesh are constrained along the channel centerlines, grading out to larger 2D element sizes away from the channels, while maintaining high quality throughout the mesh. Again, as with the previous test case, the mesh has a mean element quality of q=0.97. Furthermore, out of 115 459 elements, only 1155 elements (corresponding to 1 percentile) have quality lower than q=0.83 and only 12 elements (corresponding to 0.01 percentile) have quality below q=0.52, with the minimum element quality being q=0.29. The minimum and maximum element sizes of the mesh are 20 and 283 m, and the mesh was generated in 9.61 min.
5.3 Neches River tidal watershed
The Neches River flows southeast for approximately 670 km, entering into Sabine Lake and then into the Gulf of Mexico near Port Neches (see Texas Parks and Wildlife Department, 1974). This case study is focused on the Neches River tidal segment, which stretches approximately 45 km from the Salt Water Barrier to Sabine Lake and whose drainage area is approximately 546 km^{2} (Schramm and Jha, 2020). This area is routinely included in the computational domains used for storm surge simulations in Texas and southwestern Louisiana (see, for example, Dawson et al., 2011; Bunya et al., 2010) and includes rivers and streams of widely varying scales. For example, the Neches River, which is the main channel of the study area, has a width of approximately 300 m, which has small tributaries with widths on the order of 10 to 30 m. The complex geometry in the study area is not limited to the channels. Additionally, it includes a number of islands, whose areas range from 10 m^{2} to 300 km^{2}.
For simplicity, a rectangular study area (see black line in Fig. 21) is chosen as the domain boundary, where the identification of internal constraints associated with the land–water mask decomposition will be applied. The water mask of the study area is obtained from shoreline data provided by the National Oceanic and Atmospheric Administration (NOAA) Continually Updated Shoreline Product (CUSP) (see National Oceanic and Atmospheric Administration (NOAA), 2011). The NOAA CUSP provides a set of line segments as polylines in a shapefile format (see Fig. 21). Successive line segments, which are connected to each other, are merged to construct the boundary of the water mask (see Fig. 22). The land mask is then obtained by subtracting the water mask from the rectangular study area. As with the previous two test cases, a DEM for the domain is obtained through the USGS 3DEP.
First, the water mask is preprocessed following the steps described in Sect. 3.2.2, with small islands with an area smaller than 2000 m^{2} being filtered out and by using δ_{w}=100 m. The mesh is then generated with the following userdefined parameters: minimum mesh size of 100 m, maximum of 1000 m, the number of elements per radian K is 20, the grading limit is 0.15, and the smoothing RMSE is 5 m.
The 2D mesh generated from the water mask (shoreline data) has three types of internal constraints, namely open channels (1D domain), internal boundaries, and boundaries between waterbodies and land (see Figs. 23 and 24). The algorithm automatically identified 333 open channels, 180 internal boundaries, and 68 boundaries between waterbodies and land. This identification allows the preservation of most of the channel networks (in particular, see panels a and b of Fig. 24) and smallscale islands (in particular, see panel c of Fig. 24) in the waterbodies without using extremely small elements and provides a sharp delineation between land and water. It should be noted that there are a few narrow channels that were not identified as 1D domains. This occurs for open channels with free ends, which correspond to an order 1 medial axis branch, as a result of medial axis pruning (see Step 2 in Sect. 3.2.1). Likewise, there are a few smallscale islands that are not identified as internal boundaries. However, overall, the algorithm does an exceptional job of automatically identifying the internal constraints based on the specified width parameter, while maintaining elements of high quality. The 2D mesh has a mean element quality of q=0.93, with only 624 (out of 62 403) elements having element quality lower than q=0.65 (corresponding to 1 percentile) and only 6 elements (corresponding to 0.01 percentile) having element quality lower than q=0.25. The minimum element quality in this case is q=0.16, which is lower than the previous two test cases. This is a result of internal constraints being close to one another and is a tradeoff for preserving geographic features (see Fig. 25). While the minimum and maximum element sizes are set to 100 and 1000 m, the minimum and maximum element sizes of the mesh are 25 and 1384 m. The minimum element size is significantly smaller than target minimum element size, which is also caused by internal constraints located close to one another.
An automatic mesh generation algorithm with internal constraints, especially for coupled 1D–2D hydrodynamic models, is presented in this paper. The main objectives of the proposed algorithm are to automatically identify internal constraints (mainly channel centerlines) in the domain and to generate collocated meshes along the internal constraints with efficient sizing. The identification of internal constraints is developed for both land and water subdomains. TopoToolbox is used to extract channel centerlines from land subdomains, and an additional smoothing is applied to estimate appropriate curvature of the lines. The extraction of internal constraints from water subdomains is based on a width function and a userdefined threshold of thin channels. This enables the identification of three types of internal constraints if a water mask is given. Several additional processes are developed for complex water subdomains, including the representation of thin islands as internal boundaries. The meshes generated with the proposed algorithm have a precise alignment along the given internal constraints with the efficient sizing of highquality 2D elements. This is obtained by assigning proper target mesh sizes to the 1D–2D force equilibrium algorithms and applying postprocessing of 1D elements and density control.
While the test cases presented in this paper have, in general, elements of high quality, there are still a few elements of poor quality. This can occur when internal constraints are located too close to each other; for example, thin channels that are very close to each other or thin islands that are located very close to a shoreline. Note that these cases can be resolved by ignoring the object or allowing very small element sizes if the element quality is of higher priority. However, the proposed algorithm places a higher priority on keeping geographical features with relatively low resolution. Alternatively, poor quality of elements can be resolved with postprocess software. For example, MeshGUI (Blain et al., 2008) offers the operations of Smooth Angles
to reduce minimal angle of elements and Reduce Connectivity
to resolve high valence. However, such postprocessing software should be applied carefully. as the poorquality elements are likely created near the internal constraints, which are not considered in the postprocessing software.
Future work may include the following two objectives. First, an efficient background grid such as an octree or unstructured grid can be used to improve the computational efficiency, especially for the identification of internal constraints in water subdomains. A key factor for the identification of thin regions is the computation of the width function, which requires that the background grid is fine enough to span thin regions. It is expected that the use of an efficient background grid such as an octree or unstructured grid can improve the computational efficiency. Second, an automatized algorithm to retrieve crosssectional profiles from channels would be beneficial. Note that channel crosssectional representations, which are typically specified as triangles, rectangles, or trapezoids, are required for most of the coupled 1D–2D hydrodynamic models. While the cross sections of channels in land subdomains can be detected from the DEM, there is some ambiguity for the width of the channels. On the other hand, the width of channels in the water subdomain is relatively clear, as it can be identified with the water mask. However, the crosssectional information would need to be provided by supplemental bathymetric survey data, as standard DEMs do not contain bathymetric elevations.
In this section, we provide details of the medial axis computation briefly described in Sect. 3.2.1 and Eq. (9). In particular, we demonstrate the fact that the divergence of the VDT has positive values only on the medial axis. This is based on Voronoi polygons and their properties, as described in Lee (1982). A given polygon P can be partitioned into a set of Voronoi polygons, with boundaries referred to as Voronoi edges (see Fig. A2, for example). Voronoi polygons can be categorized into two types. One consists of Voronoi polygons whose boundaries include a segment of the external boundary ∂P. Another consists of Voronoi polygons whose boundaries do not include any segments of the external boundary. We refer to these two types as lateral and wedge types, respectively (the white and gray polygons, respectively, in Fig. A2).
In the case of lateraltype Voronoi polygons, the terminal points of the VDT are on the corresponding external boundary segment (see Lemma 1 in Lee, 1982, and Fig. A2). Note that by the definition of VDT, the VDT is perpendicular with the corresponding external boundary segment. This VDT can be represented, once the corresponding external boundary segment is projected to y=0 (see Fig. A1), as follows:
In the wedgetype Voronoi polygons, the terminal points of the VDT are vertices of a given polygon (see Lemma 1 in Lee, 1982, and Fig. A2). This VDT can be represented, once the corresponding vertex is projected to the origin (see Fig. A1), as follows:
Therefore, we have
Now, ∇⋅V is computed on the Voronoi edges. From Corollary 3 in Lee (1982), the medial axis is a subset of the Voronoi edges. Let us call the Voronoi edges that are not part of the medial axis the extra Voronoi edges. Note that, according to Corollary 3 in Lee (1982), the extra Voronoi edges are a subset of the boundaries between lateral and wedgetype Voronoi polygons. With the projection of external boundaries and vertices described above and Eqs. (A1) and (A2), it can be shown that VDT is continuous but nondifferentiable on extra Voronoi edges. Therefore, the divergence of VDT cannot be analytically defined on extra Voronoi edges, but here we show that the numerical divergence of VDT is between −1 and −2 on extra Voronoi edges. For simplicity, let us assume that the extra Voronoi edge is projected to x=0, which is the case shown in Fig. A1. Since the VDT is continuous and differentiable on each lateral and wedgetype Voronoi polygon, forward and backward difference schemes on $\mathit{x}=(\mathrm{0},y)$ give $\mathrm{\nabla}\cdot \mathit{V}(\mathit{x},P)=\mathrm{2}$ and −1, respectively. A central difference scheme on $\mathit{x}=(\mathrm{0},y)$ gives
Therefore, the numerical divergence from forward, backward, and central difference schemes is between −1 and −2 on extra Voronoi edges. Finally, on the medial axis, the VDT is discontinuous, and the divergence cannot be computed explicitly. However, given the fact that the VDT diverges from the medial axis to external boundaries, the numerical divergence must be positive. The ∇⋅V is numerically computed with divergence
function in MATLAB and shown in Fig. A3.
However, computation of ∇⋅V on a background grid results in positive values for multiple grid points near the medial axis (see Fig. A4). As the medial axis branches play a key role in the proposed methodology, medial axis branches are constructed from the cluster of medial axis points. The morphological operator bwmorph
in MATLAB is utilized for the construction of medial axis. The cluster of medial axis points are thinned first, and end points and branch points are identified. Then, starting from an end point or branch point, all traversable points based on 8connectivity are identified until it reaches to another end point or branch point. After repeating this process, medial axis branches are constructed.
One advantage of this method is that it requires low additional computational cost. In our mesh generation algorithm, the VDT is computed as part of computing the distance map (note that $\mathrm{d}(\mathit{x},\partial P)=\left\right\mathit{V}(\mathit{x},P)\left\right$), which is as essential requirement. Therefore, medial axis can be obtained with a simple additional step, the computation of divergence of VDT, which incurs a small computational cost. The medial axis can also be found with criteria based on the gradient of the distance map $\mathrm{d}(\mathit{x},\partial P)$ (see Koko, 2015; Roberts et al., 2019), i.e.,
where ϵ is a userspecified parameter (typically taken to be 0.9); however, based on our experiments, the medial axis computed from Eq. (9) tends to be more accurate than the medial axis computed from Eq. (A7).
The current version of ADMESH+, the mesh generator presented in this study, is available from the Zenodo archive at https://doi.org/10.5281/zenodo.10242565 (Kang et al., 2023). The ADMESH+ code is under active development, and the latest version will be made available upon request to the corresponding author.
EJK made an initial suggestion and guided this work. YK and EJK equally contributed to the development of the framework. YK developed the model code and performed it for the given test cases.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This research has been supported by U.S. National Science Foundation (grant no. ICER1854991).
This paper was edited by Simone Marras and reviewed by Shintaro Bunya and Chris Massey.
Adeogun, A., Pathirana, A., and Daramola, M.: 1D2D Hydrodynamic Model Coupling for Inundation Analysis of Sewer Overflow, J. Eng. Appl. Sci., 7, 356–362, https://doi.org/10.3923/jeasci.2012.356.362, 2012. a, b, c
Adeogun, A. G., Daramola, M. O., and Pathirana, A.: Coupled 1D2D Hydrodynamic Inundation Model for Sewer Overflow: Influence of Modeling Parameters, Water Science, 29, 146–155, https://doi.org/10.1016/j.wsj.2015.12.001, 2015. a, b, c
Avdis, A., Candy, A. S., Hill, J., Kramer, S. C., and Piggott, M. D.: Efficient Unstructured Mesh Generation for Marine Renewable Energy Applications, Renewable Energy, 116, 842–856, https://doi.org/10.1016/j.renene.2017.09.058, 2018. a
Bailey, R. T., Tasdighi, A., Park, S., TavakoliKivi, S., Abitew, T., Jeong, J., Green, C. H., and Worqlul, A. W.: APEXMODFLOW: A New integrated model to simulate hydrological processes in watershed systems, Environ. Modell. Softw., 143, 105093, https://doi.org/10.1016/j.envsoft.2021.105093, 2021. a
Bakhtyar, R., Maitaria, K., Velissariou, P., Trimble, B., Mashriqui, H., Moghimi, S., Abdolali, A., der Westhuysen, A. J. V., Ma, Z., Clark, E. P., and Flowers, T.: A New 1D/2D Coupled Modeling Approach for a RiverineEstuarine System Under Storm Events: Application to Delaware River Basin, J. Geophys. Res.Oceans, 125, e2019JC015822, https://doi.org/10.1029/2019JC015822, 2020. a, b, c
Bhuyian, M. N. M., Kalyanapu, A. J., and Nardi, F.: Approach to Digital Elevation Model Correction by Improving Channel Conveyance, J. Hydrol. Eng., 20, 04014062, https://doi.org/10.1061/(ASCE)HE.19435584.0001020, 2015. a
Bilgili, A., Smith, K. W., and Lynch, D. R.: BatTri: A TwoDimensional BathymetryBased Unstructured Triangular Grid Generator for Finite Element Circulation Modeling, Comput. Geosci., 32, 632–642, https://doi.org/10.1016/j.cageo.2005.09.007, 2006. a
Blain, C., Linzell, R., and Massey, T.: MeshGUI: A Mesh Generation and Editing Toolset for the ADCIRC Model, Tech. Rep. NRL/MR/7322–089083, Naval Research Laboratory, https://apps.dtic.mil/sti/tr/pdf/ADA478174.pdf (last access: 28 Novmeber 2023), 2008. a
Bunya, S., Kubatko, E. J., Westerink, J. J., and Dawson, C.: A Wetting and Drying Treatment for the Runge–Kutta Discontinuous Galerkin Solution to the Shallow Water Equations, Comput. Method. Appl. M., 198, 1548–1562, https://doi.org/10.1016/j.cma.2009.01.008, 2009. a
Bunya, S., Dietrich, J. C., Westerink, J. J., Ebersole, B. A., Smith, J. M., Atkinson, J. H., Jensen, R., Resio, D. T., Luettich, R. A., Dawson, C., Cardone, V. J., Cox, A. T., Powell, M. D., Westerink, H. J., and Roberts, H. J.: A HighResolution Coupled Riverine Flow, Tide, Wind, Wind Wave, and Storm Surge Model for Southern Louisiana and Mississippi. Part I: Model Development and Validation, Mon. Weather Rev., 138, 345–377, https://doi.org/10.1175/2009MWR2906.1, 2010. a
Bunya, S., Luettich, R. A., and Blanton, B. O.: Techniques to Embed Channels in Finite Element Shallow Water Equation Models, Adv. Eng. Softw., 185, 103516, https://doi.org/10.1016/j.advengsoft.2023.103516, 2023. a
Candy, A. S. and Pietrzak, J. D.: Shingle 2.0: generalising selfconsistent and automated domain discretisation for multiscale geophysical models, Geosci. Model Dev., 11, 213–234, https://doi.org/10.5194/gmd112132018, 2018. a
Chen, Y., Wang, Z., Liu, Z., and Zhu, D.: 1D–2D Coupled Numerical Model for ShallowWater Flows, J. Hydraul. Eng., 138, 122–132, https://doi.org/10.1061/(ASCE)HY.19437900.0000481, 2012. a, b, c
Conroy, C. J., Kubatko, E. J., and West, D. W.: ADMESH: An Advanced, Automatic Unstructured Mesh Generator for Shallow Water Models, Ocean Dynam., 62, 1503–1517, https://doi.org/10.1007/s1023601205740, 2012. a, b, c, d, e, f, g, h
D'Alpaos, L. and Defina, A.: Mathematical Modeling of Tidal Hydrodynamics in Shallow Lagoons: A Review of Open Issues and Applications to the Venice Lagoon, Comput. Geosci., 33, 476–496, https://doi.org/10.1016/j.cageo.2006.07.009, 2007. a, b
Dawson, C., Kubatko, E. J., Westerink, J. J., Trahan, C., Mirabito, C., Michoski, C., and Panda, N.: Discontinuous Galerkin Methods for Modeling Hurricane Storm Surge, Adv. Water Resour., 34, 1165–1176, https://doi.org/10.1016/j.advwatres.2010.11.004, 2011. a, b, c, d
Delelegn, S. W., Pathirana, A., Gersonius, B., Adeogun, A. G., and Vairavamoorthy, K.: MultiObjective Optimisation of Cost–Benefit of Urban Flood Management Using a 1D2D Coupled Model, Water Sci. Technol., 63, 1053–1059, https://doi.org/10.2166/wst.2011.290, 2011. a, b, c
Engwirda, D.: JIGSAWGEO (1.0): locally orthogonal staggered unstructured grid generation for general circulation modelling on the sphere, Geosci. Model Dev., 10, 2117–2140, https://doi.org/10.5194/gmd1021172017, 2017. a
Fan, Y., Ao, T., Yu, H., Huang, G., and Li, X.: A Coupled 1D2D Hydrodynamic Model for Urban Flood Inundation, Adv. Meteorol., 2017, e2819308, https://doi.org/10.1155/2017/2819308, 2017. a
Field, D. A.: Qualitative Measures for Initial Meshes, Int. J. Numer. Meth. Eng., 47, 887–906, https://doi.org/10.1002/(SICI)10970207(20000210)47:4<887::AIDNME804>3.0.CO;2H, 2000. a, b
Gejadze, I. Y. and Monnier, J.: On a 2D “Zoom” for the 1D Shallow Water Model: Coupling and Data Assimilation, Comput. Method. Appl. M., 196, 4628–4643, https://doi.org/10.1016/j.cma.2007.05.026, 2007. a, b
Ghostine, R., Hoteit, I., Vazquez, J., Terfous, A., Ghenaim, A., and Mose, R.: Comparison between a Coupled 1D2D Model and a Fully 2D Model for Supercritical Flow Simulation in Crossroads, J. Hydraul. Res., 53, 274–281, https://doi.org/10.1080/00221686.2014.974081, 2015. a, b, c
Gichamo, T. Z., Popescu, I., Jonoski, A., and Solomatine, D.: River CrossSection Extraction from the ASTER Global DEM for Flood Modeling, Environ. Modell. Softw., 31, 37–46, https://doi.org/10.1016/j.envsoft.2011.12.003, 2012. a
Goodrich, D. C., Keefer, T. O., Unkrich, C. L., Nichols, M. H., Osborn, H. B., Stone, J. J., and Smith, J. R.: LongTerm Precipitation Database, Walnut Gulch Experimental Watershed, Arizona, United States, Water Resour. Res., 44, W05S04, https://doi.org/10.1029/2006WR005782, 2008. a
Goodrich, D. C., Burns, I. S., Unkrich, C. L., Semmens, D. J., Guertin, D. P., Hernandez, M., Yatheendradas, S., Kennedy, J. R., and Levick, L. R.: KINEROS2/AGWA: Model Use, Calibration, and Validation, T. ASABE, 55, 1561–1574, https://doi.org/10.13031/2013.42264, 2012. a
Gorman, G. J., Piggott, M. D., Pain, C. C., de Oliveira, C. R. E., Umpleby, A. P., and Goddard, A. J. H.: Optimisation Based Bathymetry Approximation through Constrained Unstructured Mesh Adaptivity, Ocean Model., 12, 436–452, https://doi.org/10.1016/j.ocemod.2005.09.004, 2006. a
Hagen, S. C., Horstmann, O., and Bennett, R. J.: An Unstructured Mesh Generation Algorithm for Shallow Water Modeling, Int. J. Comput. Fluid D., 16, 83–91, https://doi.org/10.1080/10618560290017176, 2002. a
Kang, Y., Kubatko, E. J., Conroy, C. J., and West, D. W.: YounghunKang/ADMESH: V3.0.1, Zenodo [code and data set], https://doi.org/10.5281/zenodo.10242565, 2023. a
Koko, J.: A Matlab Mesh Generator for the TwoDimensional Finite Element Method, Appl. Math. Comput., 250, 650–664, https://doi.org/10.1016/j.amc.2014.11.009, 2015. a, b
Kuiry, S. N., Sen, D., and Bates, P. D.: Coupled 1D–Quasi2D Flood Inundation Model with Unstructured Grids, J. Hydraul. Eng., 136, 493–506, https://doi.org/10.1061/(ASCE)HY.19437900.0000211, 2010. a, b
Lee, D. T.: Medial Axis Transformation of a Planar Shape, IEEE T. Pattern Anal., PAMI4, 363–369, https://doi.org/10.1109/TPAMI.1982.4767267, 1982. a, b, c, d, e, f
Li, D., Chen, S., Zhen, Z., Bu, S., and Li, Y.: Coupling a 1Dlocal Inertia 2D Hydraulic Model for Flood Dispatching Simulation in a Floodplain under Joint Control of Multiple Gates, Nat. Hazards, 109, 1801–1820, https://doi.org/10.1007/s1106902104899z, 2021. a
Lin, B., Wicks, J. M., Falconer, R. A., and Adams, K.: Integrating 1D and 2D Hydrodynamic Models for Flood Simulation, P. I. Civil Eng.Wat. M., 159, 19–25, https://doi.org/10.1680/wama.2006.159.1.19, 2006. a
Lin, J., Jin, S., Ai, C., and Ding, W.: Experimental and Numerical Investigation of River Closure Project, Water, 12, 241, https://doi.org/10.3390/w12010241, 2020. a
Liu, Q., Qin, Y., Zhang, Y., and Li, Z.: A Coupled 1D–2D Hydrodynamic Model for Flood Simulation in Flood Detention Basin, Nat. Hazards, 75, 1303–1325, https://doi.org/10.1007/s1106901413733, 2015. a, b, c
Luettich, R. A. and Westerink, J. J.: Elemental Wetting and Drying in the ADCIRC Hydrodynamic Model: Upgrades and Documentation for ADCIRC Version 34.XX, Contract Report, US Army Corps of Engineers, https://adcirc.org/wpcontent/uploads/sites/2255/2018/11/1999_Luettich01.pdf (last access: 24 March 2022), 1999. a
Marin, J. and Monnier, J.: Superposition of Local Zoom Models and Simultaneous Calibration for 1D–2D Shallow Water Flows, Math. Comput. Simulat., 80, 547–560, https://doi.org/10.1016/j.matcom.2009.09.001, 2009. a, b
Martini, P., Carniello, L., and Avanzi, C.: Two dimensional modelling of flood flows and suspended sedimenttransport: the case of the Brenta River, Veneto (Italy), Nat. Hazards Earth Syst. Sci., 4, 165–181, https://doi.org/10.5194/nhess41652004, 2004. a, b
Meng, H., Green, T. R., Salas, J. D., and Ahuja, L. R.: Development and Testing of a TerrainBased Hydrologic Model for Spatial Hortonian Infiltration and Runoff/On, Environ. Modell. Softw., 23, 794–812, https://doi.org/10.1016/j.envsoft.2007.09.014, 2008. a
MoralesHernández, M., Petaccia, G., Brufau, P., and GarcíaNavarro, P.: Conservative 1D–2D Coupled Numerical Strategies Applied to River Flooding: The Tiber (Rome), Appl. Math. Model., 40, 2087–2105, https://doi.org/10.1016/j.apm.2015.08.016, 2016. a
Mullikin, J. C.: The Vector Distance Transform in Two and Three Dimensions, CVGIPGraph. Model. Im., 54, 526–535, https://doi.org/10.1016/10499652(92)900726, 1992. a
National Oceanic and Atmospheric Administration (NOAA): Continually Updated Shoreline Product, NOAA [data set], https://www.ngs.noaa.gov/CUSP/ (last access: 23 July 2020), 2011. a
Néelz, S. and Pender, G.: Desktop Review of 2D Hydraulic Modelling Packages, Environment Agency, Bristol, ISBN 9781849110792, https://assets.publishing.service.gov.uk/media/6033a9888fa8f543294411a8/_SC080035_Desktop_review_of_2D_hydraulic_packages_Phase_1_Report.pdf (last access: 23 September 2021), 2009. a
Park, S., Nielsen, A., Bailey, R. T., Trolle, D., and Bieger, K.: A QGISbased graphical user interface for application and evaluation of SWATMODFLOW models, Environ. Modell. Softw., 111, 493–497, https://doi.org/10.1016/j.envsoft.2018.10.017, 2019. a
Persson, P.O.: Mesh Size Functions for Implicit Geometries and PDEbased Gradient Limiting, Eng. Comput., 22, 95–109, 2006. a, b
Persson, P.O. and Strang, G.: A Simple Mesh Generator in MATLAB, SIAM Rev., 46, 329–345, https://doi.org/10.1137/S0036144503429121, 2004. a, b, c, d, e
Pham Van, C., de Brye, B., Deleersnijder, E., Hoitink, A. J. F., Sassi, M., Spinewine, B., Hidayat, H., and SoaresFrazão, S.: Simulations of the Flow in the Mahakam River–Lake–Delta System, Indonesia, Environ. Fluid Mech., 16, 603–633, https://doi.org/10.1007/s1065201694454, 2016. a, b, c
Pramanik, N., Panda, R. K., and Sen, D.: One Dimensional Hydrodynamic Modeling of River Flow Using DEM Extracted River Crosssections, Water Resour. Manag., 24, 835–852, https://doi.org/10.1007/s1126900994746, 2010. a
Price, R. K.: Toward Flood Routing in Natural Rivers, J. Hydraul. Eng., 144, 04017070, https://doi.org/10.1061/(ASCE)HY.19437900.0001414, 2018. a
Rath, S.: Model Discretisation in 2D Hydroinformatics Based on High Resolution Remote Sensing Data and the Feasibility of Automated Model Parameterisation, PhD thesis, Hamburg University of Technology, 2007. a
Remacle, J.F. and Lambrechts, J.: Fast and Robust Mesh Generation on the Sphere–Application to Coastal Domains, Comput. Aided Design, 103, 14–23, https://doi.org/10.1016/j.cad.2018.03.002, 2018. a
Renard, K. G., Nichols, M. H., Woolhiser, D. A., and Osborn, H. B.: A Brief Background on the U.S. Department of Agriculture Agricultural Research Service Walnut Gulch Experimental Watershed, Water Resour. Res., 44, W05S02, https://doi.org/10.1029/2006WR005691, 2008. a
Roberts, K. J., Pringle, W. J., and Westerink, J. J.: OceanMesh2D 1.0: MATLABbased software for twodimensional unstructured mesh generation in coastal ocean modeling, Geosci. Model Dev., 12, 1847–1868, https://doi.org/10.5194/gmd1218472019, 2019. a, b, c
Schramm, M. and Jha, A.: Technical Support Document for Four Total Maximum Daily Loads for Indicator Bacteria in Neches River Tidal, Tech. Rep., Texas Water Resource Institute, Texas A&M University, https://www.tceq.texas.gov/downloads/waterquality/tmdl/nechesrivertidalrecreational118/118as471nechestidalbacteriatsd2020july.pdf (last access: 14 February 2024), 2022. a
Schwanghart, W. and Kuhn, N. J.: TopoToolbox: A Set of Matlab Functions for Topographic Analysis, Environ. Modell. Softw., 25, 770–781, https://doi.org/10.1016/j.envsoft.2009.12.002, 2010. a, b
Seyoum, S. D., Vojinovic, Z., Price, R. K., and Weesakul, S.: Coupled 1D and Noninertia 2D Flood Inundation Model for Simulation of Urban Flooding, J. Hydraul. Eng., 138, 23–34, https://doi.org/10.1061/(ASCE)HY.19437900.0000485, 2012. a
Shewchuk, J. R.: Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator, in: Applied Computational Geometry Towards Geometric Engineering, edited by: Lin, M. C. and Manocha, D., Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, 203–222, ISBN 9783540706809, https://doi.org/10.1007/BFb0014497, 1996. a
Stelling, G. S. and Verwey, A.: Numerical Flood Simulation, in: Encyclopedia of Hydrological Sciences, chap. 16, John Wiley & Sons, Ltd, ISBN 9780470848944, https://doi.org/10.1002/0470848944.hsa025a, 2006. a, b
Stone, J. J., Nichols, M. H., Goodrich, D. C., and Buono, J.: LongTerm Runoff Database, Walnut Gulch Experimental Watershed, Arizona, United States, Water Resour. Res., 44, W05S05, https://doi.org/10.1029/2006WR005733, 2008. a
Tefera, G. W. and Ray, R. L.: Hydrology and hydrological extremes under climate change scenarios in the Bosque watershed, NorthCentral Texas, USA, Environ. Sci. Pollut. R., https://doi.org/10.1007/s11356023274771, 2023. a
Teng, J., Jakeman, A. J., Vaze, J., Croke, B. F. W., Dutta, D., and Kim, S.: Flood Inundation Modelling: A Review of Methods, Recent Advances and Uncertainty Analysis, Environ. Modell. Softw., 90, 201–216, https://doi.org/10.1016/j.envsoft.2017.01.006, 2017. a
Texas Parks and Wildlife Department: An Analysis of Texas Waterways: A Report on the Physical Characteristics of Rivers, Streams, and Bayous in Texas, Tech. Rep., The Texas Agricultural Extension Service, Texas A&M University System, https://tpwd.texas.gov/publications/pwdpubs/pwd_rp_t3200_1047/ (last access: 24 October 2021), 1974. a
Timbadiya, P. V., Patel, P. L., and Porey, P. D.: OneDimensional Hydrodynamic Modelling of Flooding and Stage Hydrographs in the Lower Tapi River in India, Current Science, 106, 708–716, 2014. a
Timbadiya, P. V., Patel, P. L., and Porey, P. D.: A 1D–2D Coupled Hydrodynamic Model for River Flood Prediction in a Coastal Urban Floodplain, J. Hydrol. Eng., 20, 05014017, https://doi.org/10.1061/(ASCE)HE.19435584.0001029, 2015. a, b
U.S. Geological Survey: USGS National Watershed Boundary Dataset (WBD) Downloadable Data Collection, USGS [data set], https://www.usgs.gov/nationalhydrography/watershedboundarydataset (last access: 19 June 2023), 2014. a, b
U.S. Geological Survey: USGS National Hydrography Dataset (NHD) Downloadable Data Collection – National Geospatial Data Asset (NGDA) National Hydrography Dataset (NHD), USGS [data set], https://www.usgs.gov/nationalhydrography/nationalhydrographydataset (last access: 30 March 2022), 2016. a
U.S. Geological Survey: Digital Elevation Models (DEMs) – USGS National Map 3DEP Downloadable Data Collection: U.S. Geological Survey, USGS [data set], https://www.usgs.gov/3delevationprogram (last access: 19 June 2023), 2017. a
Vojinovic, Z. and Tutulic, D.: On the Use of 1D and Coupled 1D–2D Modelling Approaches for Assessment of Flood Damage in Urban Areas, Urban Water J., 6, 183–199, https://doi.org/10.1080/15730620802566877, 2009. a, b, c
Wessel, P. and Smith, W. H. F.: A Global, SelfConsistent, Hierarchical, HighResolution Shoreline Database, J. Geophys. Res.Sol. Ea., 101, 8741–8743, https://doi.org/10.1029/96JB00104, 1996. a
West, D. W., Kubatko, E. J., Conroy, C. J., Yaufman, M., and Wood, D.: A Multidimensional Discontinuous Galerkin Modeling Framework for Overland Flow and Channel Routing, Adv. Water Resour., 102, 142–160, https://doi.org/10.1016/j.advwatres.2017.02.008, 2017. a, b
Woodhead, S., Asselman, N., Zech, Y., SoaresFrazao, S., Bates, P., and Kortenhaus, A.: Evaluation of Inundation Models, FLOODsite Project Report T080701, http://www.floodsite.net/html/partner_area/project_docs/T08_08_01_inundation_modelling_ExecSum_v2_4_p01.pdf (last access: 23 September 2021), 2007. a
Yerry, M. A. and Shephard, M. S.: A Modified Quadtree Approach To Finite Element Mesh Generation, IEEE Comput. Graph., 3, 39–46, https://doi.org/10.1109/MCG.1983.262997, 1983. a
Yu, C. and Duan, J.: Simulation of Surface Runoff Using Hydrodynamic Model, J. Hydrol. Eng., 22, 04017006, https://doi.org/10.1061/(ASCE)HE.19435584.0001497, 2017. a
 Abstract
 Introduction
 Overview of coupled 1D–2D hydrodynamic domains
 Identification of internal constraints
 Force equilibrium with internal constraints
 Results
 Conclusions
 Appendix A: Computation of medial axis with vector distance transform
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Overview of coupled 1D–2D hydrodynamic domains
 Identification of internal constraints
 Force equilibrium with internal constraints
 Results
 Conclusions
 Appendix A: Computation of medial axis with vector distance transform
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Financial support
 Review statement
 References