Leveraging mesh modularization to lower the computational cost of localized updates to regional 2D hydrodynamic model outputs

Abstract

Hydrodynamic model outputs are used in urban flood risk modelling, flood alert systems, and Monte Carlo hazard assessment. This study tackles an under-explored challenge wherein regular updates to the spatial characteristics of the watershed – due to factors such as changing land use – alter the watershed’s response to rainfall forcing, thus rendering existing model outputs obsolete. Because state-of-the-art hydrodynamic models are computationally expensive, frequently re-running simulations can be costly. Modularization addresses this problem by requiring re-computation only for a limited domain affected by the land use changes. This article introduces a novel approach by modularizing the 2D domain into independent sub-domains before (‘discrete’) and after (‘abstract’) the numerical computations. Using the Hydrologic Engineering Center River Analysis System (HEC-RAS) 2D model of a large urban watershed in Houston as an illustrative and generalizable testbed, we show that both the discrete and abstract modularization closely approximates the results from re-running the entire model. The computational cost of modularization scales linearly with model size for memory requirements as storing the solution on the interior boundaries (discrete) or throughout the domain (abstract) are necessary. This trade-off of memory for computation may facilitate advances in surrogate modelling or Monte Carlo flood risk assessment.

KEYWORDS:

• HEC-RAS
• urban hydrodynamics
• flood hazard

1. Introduction

Fast, detailed, and accurate, but scalable, hydrodynamic models are usually infeasible without High Performance Computing (HPC), yet the information they provide in resolving a floodplain response both temporally and spatially could influence flood hazard assessment, stormwater management, and emergency planning/response in urban systems. As the scale (both spatial and temporal), level of detail (i.e. resolution), and complexity of processes (e.g. presence of hydraulic structures) represented by the hydrodynamic model increase, so does the computational cost of performing a simulation. If complex large-scale models were viable, building a regional integrated flood information system that provides real-time flood warning (e.g. sub-hourly alerts on inundation hotspots and road network accessibility) and flood risk and flood hazard prediction (e.g. probabilistic floodplain mapping and damage estimates) might become feasible.

Recently, various studies have implemented Machine Learning (ML) methods in hydrodynamic applications to address this scale issue by leveraging the availability of big datasets to approximate a solution with significantly reduced computing time. Examples include using support vector machines (Bermúdez et al., 2019; Liu & Pender, 2015), random forests (Mobley et al., 2021; Wang et al., 2015), gradient boosting decision trees (Wu et al., 2020), extreme gradient boosting or XGBoost (Sanders et al., 2022), as well as deep learning methods such as convolutional neural networks (Hosseiny, 2021; Kabir et al., 2020) and Long Short-Term Memory (LSTM) (Nevo et al., 2022; Yang et al., 2022; Zou et al., 2023). While the use of ML and deep learning methods in flood risk, flood hazard, and flood prediction have grown increasingly popular owing to their promising results and endless potential applications, they could be dynamically unstable (Mosavi et al., 2018) and might require frequent retraining due to changes in boundary conditions, which is an issue common in urban hydrology.

To minimize the computational burden of constant model updates/retraining, we propose a mesh modularization framework that divides the modelling domain/mesh into multiple predefined subsets or modules. Since these subsets act as surrogates, one would only need to update or retrain specific modules that have encountered change instead of re-running the entire modelling domain. This article details an original proof of concept with discrete, pre-determined modules in a two-dimensional (2D) hydraulic model. Although we used specific hydraulic software – the Hydrologic Engineering Center River Analysis System (HEC-RAS) (Brunner, 2021) – the framework introduced in this study would be applicable to other hydrodynamic models and solvers. Additionally, the article also presents and discusses a pilot study using abstract modules generated after computation, with boundaries intended for eventual use in a real-time flood alert system for Harris County, Texas.

2. Background

In the last decade, 2D hydrodynamic solvers have gone from academic codes like TELEMAC-2D (World Bank, 2022) that solve governing Partial Differential Equations (PDEs) to widespread adoption in Hydrologic and Hydraulic (H/H) applications with many fully graphical interfaced models such as TUFLOW (BMT, 2022), InfoWorks-ICM (Innovyze, 2022), and HEC-RAS. These programs can compute detailed and accurate spatio-temporal floodplains and can show how these computed water surface elevations affect critical infrastructure. However, these tools are still computationally expensive, which poses a challenge to end-users. While numerous methods have been proposed to address this computational burden, we summarize them into three classes of method: Simplification, Supercomputing, and ML.

One class of method uses models that reduce computational cost by representing a dominant subset of physical processes. An example of this approach is LISFLOOD-FP (Bates, 2022), which uses a 1D solver on a 2D structured mesh similar to the diffusive wave approximation of the Shallow Water Equations (SWEs) (vector form shown in Equation 1). This approximation decreases the computational cost and provides good results when those assumptions can be applied (Alonso et al., 2008). However, these simplified structured mesh solvers cannot account for control structures like reservoirs or other major retention owing to the irregular geometries of real applications, leading to mesh leakage and other root issues. In addition, structured solvers are often poor at preventing mesh leakage for coarser grids as the model scales up (Shustikova et al., 2019). Of particular relevance to the problem of urban flood risk management, representing the effects of control structures and producing estimates at the scale of buildings and infrastructure present challenges to flood risk modelling. As an example, releases from the Addicks and Barker reservoirs during Hurricane Harvey 2017 significantly contributed to the flooding experienced in the Houston region (Garcia et al., 2020). Although these simplified tools are useful when applicable, excluding relevant physical processes can limit their suitability for many problems.

The second major class of approach for PDE solvers is the implementation of supercomputing, which is used in hurricane forecasting and other complex fields. This could be a viable option in theory but might not be fiscally or logistically feasible for many municipalities in practice, where most floodplain modelling resides. Supercomputing works well but does not address the fundamental trade-off of computational expenses against model detail and scope.

The third class of approach is to train surrogate models or emulators. Surrogate models use machine learning to generate a computationally cheap approximation of a more complex model by using output from the complex model as training data. However, surrogate methods require the problem being learned to be static (L’Heureux et al., 2017). Consequently, every time there is an edit to the modelled domain, the entire training process needs to be redone to learn a new approximation. Since urban landscapes are inherently dynamic systems that change frequently, a static system assumption is ill-suited for certain applications (e.g. urban flood modelling). Adding land use and other geometric edits as another degree of freedom for the the surrgoate to learn adds a multiplier for the data requirements. These requirements already present a nearly infeasible task, which may see improvements as computer science methods and data storage continue to improve but are still far from implementation.

3. Methods

The introduction of modularization starts with an applied definition, and then a more formalized definition of the method. Next will be the sequence of tests performed, starting with the discrete case, or pre-determined modules in the domain setup followed by a full abstraction of the method, or modules created after the modelling has been completed. Both tests comprise three model runs with the same input rainfall data, one with the original terrain (the existing condition) and two with the edited terrain (the edited condition). One of the edited conditions is identical to the 2D domain of the existing condition, and the other is a subset of that 2D domain using model data from the existing condition as external boundary conditions, this setup ending with the two edited condition runs showing the same results as a comparative run being the goal.

3.1. Mesh modularization

Modularization draws from multi-domain methods intended to solve varied mesh refinement detail in large simulations (Tezduyar & Osawa, 2001) and domain decomposition methods for SWEs intended to parallelize PDEs (Navon & Cai, 1993). Conceptually, modularization involves defining the contiguous portion of the domain containing the drainage outlet where an edit to the model would affect the output water depth (i.e. outside the portion of the domain selected, or module, the change in the solution due to the edit is far less than the error of the base model and in practical use is unchanged). Below in Figure 1, a watershed bounded by a 2D domain $\mathrm{\Omega }$ is discretized on a grid shown as squares. Figure 1(A) shows the existing condition with the location of a future edit being highlighted in red, and the channel being highlighted in blue.

Figure 1. Graphical overview of modularization: (A) the red square is the location of an edit in the domain $\Omega$ and the blue represents the main channel; (B) the red edit will change the water surface elevation in the green region and will change the channel volumes along the purple; (C) the boundary of the changed regions in B are denoted by ${\Omega }^{\ast }$. Using this boundary as an upstream boundary condition, a new solution in the pink region can be solved in place of $\Omega$ for the red edit. The forming of the new domain without prior modelling discretization will be defined later as abstract modularization. If the model had been built having $\Omega$ and ${\Omega }^{\ast }$ components from the start, the modularization would be discrete.

3.1.1. Practical explanation of modularization

As a more concrete example of modularization, if a large new development offset their runoff with a retention pond (red square in Figure 1(A)), the only region of the floodplain, or Water Surface ELevation (WSEL), that would differ from this change is the immediate region surrounding the terrain edit (green in Figure 1(B)) and the volumetric changes to the channels downstream of it (purple in Figure 1(B)). Anything draining into the region affected by the new development region or into the channels downstream of the development region would not be influenced by it, and as a result would not have any effect on the computed floodplain solution assuming the same external forcings. Defining this as an internal boundary where inside the new developed pond will change the floodplain, and outside of it there will be no influence (green and purple in Figure 1(B)), a contiguous subset is produced of the full domain with a known solution along the boundary except for the already dependent boundary condition of the outlet for the watershed (pink in Figure 1(C)).

This subset of the full domain is called a module (${\mathrm{\Omega }}^{\ast }$ in Figure 1(C)) as it can be edited independently, replacing the existing portion of the full model that it occupied (the red dashed outlined regions in Figure 1(C)). The naming was then extended to the entire method of splitting the physics-based domain, called modularizing the domain. In normal engineering applications, the storage of the solution information along these modular boundaries (or all cell faces in the abstract case) would be excessive, but in the same scenario of ML methods trading off computation for stored information with the training dataset, this storage minimizes the necessary computations to update the full model as it only computes the subset of the domain where edits to the output will occur.

3.1.2. Formally defining a module by its boundary

As a more formal definition of this module, we will define it as the boundary, or ${\mathrm{\Omega }}^{\ast }$, in Figure 1(C). Given that the governing PDE, in this case the SWE, is well posed and solved under reasonable physical conditions (Duruk Mutlubas et al., 2020; Ghader & Nordström, 2014), any edit (defined as a localized change to topography or land use) has a continuous relationship to the change in the solution, or water depth, on a discretized 2D domain. Below is a walkthrough of this statement that implies the existence of the previously defined internal boundary line in the domain where there is an insignificant solution difference from the edit to form a module.

Beginning with the vector form of the momentum equation from the HEC-RAS hydraulic reference manual (Brunner, 2021), (2–150) with (2–153) and (2–154) are included as Equation (1) below: $\overrightarrow{V}$ is the velocity vector, $f_c$ is the Coriolis parameter, is the virtual unit vector, $g$ is the gravitational acceleration, $z_s$ is the water surface elevation, $h$ is the water depth, $v_t$ is the eddy viscosity tensor, $n$ is Manning’s roughness coefficient, $R$ is the hydraulic radius, and ${\tau }_s$ is the wind surface stress vector. (1) $\begin{array}{rl}& \frac{\mathrm{\partial }\overrightarrow{V}}{\mathrm{\partial }t}+(\overrightarrow{V}\cdot \mathrm{\nabla })\overrightarrow{V}+f_{c}k\times \overrightarrow{V}\\ & =-g\mathrm{\nabla }{z}_s+\frac{1}{h}\mathrm{\nabla }\cdot (v_{t}h\mathrm{\nabla }\overrightarrow{V})-\frac{n^{2}g||\overrightarrow{V}||\overrightarrow{V}}{R^{4/3}}+\frac{{\tau }_s}{\rho h}\end{array}$(1) The PDE in Equation (1) can be rearranged to isolate the portion influenced by the topographic and roughness changes to the right-hand side (RHS), putting the other terms on the left-hand side (LHS) in Equations (2): (2) $\begin{array}{rl}& \frac{\mathrm{\partial }\overrightarrow{V}}{\mathrm{\partial }t}+(\overrightarrow{V}\cdot \mathrm{\nabla })\overrightarrow{V}+f_{c}k\times \overrightarrow{V}-\frac{1}{h}\mathrm{\nabla }\cdot (v_{t}h\mathrm{\nabla }\overrightarrow{V})\\ & -\frac{{\tau }_s}{\rho h}=-g\mathrm{\nabla }{z}_s-\frac{n^{2}g||\overrightarrow{V}||\overrightarrow{V}}{R^{4/3}}\end{array}$(2) Next, each side is renamed as a function of its dependent variables in Equations (3) and (4): (3) $\begin{array}{rl}S(\overrightarrow{V},h)& =\frac{\mathrm{\partial }\overrightarrow{V}}{\mathrm{\partial }t}+(\overrightarrow{V}\cdot \mathrm{\nabla })\overrightarrow{V}+f_{c}k\times \overrightarrow{V}\\ & -\frac{1}{h}\mathrm{\nabla }\cdot (v_{t}h\mathrm{\nabla }\overrightarrow{V})-\frac{{\tau }_s}{\rho h}\end{array}$(3) (4) $\begin{array}{rl}b(\overrightarrow{V})& =-g\mathrm{\nabla }{z}_s-\frac{n^{2}g|\overrightarrow{V}|\overrightarrow{V}}{R^{4/3}}\end{array}$(4) with the following three points: the RHS and LHS remain equal, the total energy of the system outside of the edit remains constant, and the total mass of the full system is unchanged. It can be implied that a change to $b$ has a continuous relationship to a change in the solution of $h$ for a given velocity. This implication also suggests that the norms of these changes are continuously related, as shown in Equation (5): (5) $\begin{array}{r}S(\overrightarrow{V},h)=b(\overrightarrow{V})\Rightarrow ||\mathrm{\Delta }b{||}^2\mapsto ||\mathrm{\Delta }S(h){||}^2\end{array}$(5) This continuity further implies that there exists some maximum distance L in the 2D domain $\Omega$ of x, y coordinates, where the local change to the solution $\epsilon$ is less than the error of the model $\text{e}$ for each $\mathrm{\Delta }{b}_{x,y}$, or the x, y location of a topographic or roughness edit, as shown in Equations (6a) and (6b): (6a) $\begin{array}{rl}\Rightarrow \mathrm{\exists }L;L_{\mathrm{\Delta }{b}_{x,y}}& =max(||\mathrm{\Delta }{b}_{x,y}-(\mathrm{\Delta }S(h),\epsilon ){||}_{L^2(\Omega )})\end{array}$(6a) (6b) $\begin{array}{rl}\text{s}\text{.t}\text{. }\epsilon & \le |\text{e}|;e,\epsilon \in R;x,y\in \Omega \end{array}$(6b) In total, this implies there exists a line in the 2D domain where the solution along that line is unchanged with respect to the model error. Therefore, if all the cell faces or discrete module boundaries lie beyond this line L away from the domain edit, the solution along this line can be used as an external boundary condition for the portion of the domain lying inside it that contains the edit. The resulting module, or ${\Omega }^{\ast },$ produces the same model solution for the PDE as the full domain when it is run with the same edit.

3.2. Discrete modularization test

The initial modularization testing was done on a highly urbanized, upstream portion of the White Oak Bayou watershed in Houston (Figure 2). This section was chosen owing to the limited amount of natural infiltration and dense quantity of flood control structures. The scale of the modelling domain was minimized to both limit computation time and limit the number of modelling elements that were still significant with 25 individual 2D modelling domains or discrete modules and 97 individual connections between them. The discrete modules can be seen outlined in Figure 2.

Figure 2. Overview of testing and case study modelling in Texas, USA. The yellow region shows the discrete testing region zoomed in on the upper right with the black lines delineating the sub-domains of 2D models. A detail of the terrain edits for a local pond fill is shown in the lower right for location T1. The blue outline is the location of the abstract testing region for the case study.

Some of these modules were part of a segmented 2D channel section implemented for downstream volume conveyance changes. Each module was interconnected with every adjacent module with HEC-RAS ‘SA/2D Area Connections’ (or on-terrain weirs highlighted purple in Figure 3). Each 2D module was discretized using an internal mesh generation program in HEC-RAS with 61 m (200 ft) cell spacing and used mesh refinement along the channels with 30.5 m (100 ft) cell spacing (Figure 3). The delineations for the cuts defining the modules were determined from the Harris County 1 m LiDAR Digital Elevation Map (DEM) (HGAC, 2018) of the region using an ArcMap™ program called ArcHydro (now an integrated tool in ArcGIS^® Pro) (ESRI, 2013) effectively defining sub-basins of the modelling domain.

Figure 3. Stage and flow comparisons of the discrete modularization tests performed. The rows display the stage and flow of the outlet for the tested region. The columns of existing conditions and edited condition differences highlight the module edit run in the context of the change the geometry caused.

ArcHydro determines the direction of drainage across the domain with a slope raster it generates from the DEM and delineates drainage boundaries from the user defined watershed outlet and an eight directional slope raster. A channel shapefile is similarly determined from a threshold value calculation from the same slope raster and accompanied flow accumulation raster made with ArcHydro. The generated channel lines were then buffered 122 m (400 ft) to ensure it was well outside the floodplain and manually cut perpendicular to the flow direction to divide into smaller modules. These manually made channel modules were then placed over the ArcHydro delineations of sub-basins superseding the existing boundaries to produce the complete set of modules after some manual dissolving of insignificantly sized sub-basins resulting in what is shown in the call out of Test in Figure 2.

Other geometry inputs for the discrete testing model included the land use data, which was the NLCD 2016 data from the Multi-Resolution Land Characteristics Consortium (MRLC, 2020), and the associated roughness values, which were derived from Kalyanapu et al. (2010). The roughness values were adjusted to 75% of those reported study values owing to previous calibration work done (Garcia et al., 2020) on a model in the neighbouring Addicks watershed. For the single rainfall input for this phase in testing, an approximately 100-year two-hour event of 17.78 cm (7.0 in.) was used as the rain-on-grid input for the model and was derived from the NOAA Atlas 14 data statistics (NOAA, 2018) of the location. The rainfall was temporally split into 6.35 cm (2.5 inch) and 11.43 cm (4.5 inch) volumes over the two one-hour periods to minimize the computation time of these tests by simulating a flash flood condition. The only boundary condition aside from the rainfall is the normal depth condition at the watershed outlet using the average slope of the region making the model a rain-on-grid simulation.

To test the modularization methodology, a comparison of the outputs from a contiguous set of modules with full domain outputs was conducted. To compare the edited condition runs effectively, significant changes to the output hydrograph needed to be created. In this study, significant terrain modifications were made through filling selected detention ponds to a constant DEM elevation across each scenario. The edited condition for these few tests consisted of various detention pond modifications and were each named after the pond(s) being edited. The ponds selected for testing were intended to sample variations effectively that could occur from both the size and location of the detention changes made in the domain. The three ponds selected were compared in a total of four evaluation runs, one condition for each pond being filled individually: T1, T2, and T3, and one of all three ponds being filled together, T123 (Figure 2). Each edited condition was compared to the existing condition run to ensure that significant changes were made to the outflow hydrographs at the watershed outlet and to validate their use for the comparative analysis.

3.3. Abstract modularization case study

Along with the implementation of a distributed net rainfall solver for the model, HEC-RAS 6.0.0 onward allowed for the abstraction of the module development such that a single domain could be used, and a single module could be generated after the fact from the internal stored data. The inputs to the HEC-RAS model were the same but with updated land use (MRLC, 2020) and the addition of soils data for HEC-RAS to use Green-and-Amp to solve for net rainfall in each cell. These input soils were from the USDA Natural Resources Conservation Service Soils database (USDA, 2022) and the parameters were based on USDA soils documentation as well (USDA, 2019). The rainfall rasters were scrubbed from the MTArchive hosted by Iowa State (Herzmann, 2022) and were processed into gridded [.dss] formats do be read easily by HEC-RAS using the Vortex API from the USACE HEC Github (Lehman, 2022). The storm chosen was Tax Day 2016 (HCFCD, 2022), which caused significant flooding in parts of Harris County including Brays Bayou where the edited condition will focus. As with the discrete case, the rain-on-grid simulation had a normal depth boundary condition at the outlet, but additionally some normal depth conditions on the perimeter as the terrain-based watershed was not computed. Water can freely exit the model on a majority of the boundary with the average terrain slope as constructed.

Since all of the data at each cell face, centre, and more can be pulled from the computation run data, the modules do not need to be defined explicitly in the model development but rather to be constructed after the domain edit using the stored data. This method has some requirements for necessary meta-data from the domain since the pulled internal data will not directly map to a newly constructed module. This includes which cells and faces are in the new domain, which faces are now on the domain boundary, what is the order of these boundary faces, which directions are the computed norms of these faces on the boundary, and what faces are being combined to a single external condition (since a corner no longer has external references). All this meta-data needs to be computed or pulled to allow an arbitrary cut to form a new module that can be updated independently.

Exactly how this was accomplished for HEC-RAS is not elaborated here as it would be very model/solver specific. Overall, the automation of this process will produce a set of boundary conditions each of which corresponds to a single cell face with the exception of corners and which has the flow passing through those faces for each computation timestep (or as near as feasible to each timestep). The dependent boundary conditions, such as the watershed outlet, will remain the same, and no cell centre will move or change volume from the original model to the module.

The defined single cut in the original domain this automation will produce (the blue line in Figure 2) is the same as described earlier and defines the boundary where everything inside is affected by a terrain edit, and everything outside is not. The addition of abstracting this cut to any selection of continuous cell faces rather than discretely defining it in the construction of the model allows any arbitrary contiguous subset of cells containing the outlet to be used. This not only refines the improvements modularization can yield, but also truly minimizes the future computations and maximizes the improvements the methodology can bring.

Once the module computer questions are answered in the automation and the abstract module is generated, it needs to be tested in the same manner as in the discrete case. The test setup for this case study adopted a similar approach to the original test, however it was done in a more realistic manner. There are two golf courses in the middle of Brays Bayou that have been discussed with a view to becoming detention ponds for regional mitigation strategies. For this study, artificial ponds were cut into the LiDAR over the area of the existing golf courses to the same depths as the adjacent channel. These new detention ponds were used as the terrain edit on a hindcast of the Tax Day 2016 event meteorology. The existing condition Tax Day hindcast was used as the input dataset for the generated module with the new detention ponds for the edited condition and the full model was run again the with the edited condition. This case study was intended to show that the abstract modularization can work in real-world conditions and be used in a future pilot study for the model.

4. Results

As stated before, the comparison of the two edited conditions shows the viability of modularization for either the discrete case or the abstract case. This section will begin with the initial testing done with the discrete case and multiple scenarios (T1, T2, T3, and T123) with different terrain edits and then show the case study of the abstract case for the Tax Day hindcast.

4.1. Discrete modularization test

To start the analysis, the difference hydrographs (existing minus the full model edited conditions) are highlighted in the top panels of Figure 3. These graphs show the washed response to the detention pond filling for each tested condition and will highlight that the conditions did significantly edit the hydrographs. The top panels of Figure 3 clearly show that the most downstream pond, T1, shows the most significant individual difference to watershed outlet response, and this is confirmed by the combined condition of T123 showing only a marginal difference.

To visualize the performance of the module run in isolation against the full model being re-run (full edit) for each condition, the bottom panels of Figure 3 show the outlet hydrograph differences of the edited conditions (full edit minus module edit). As these graphs in Figure 3 show, there is no longer a correlation to the terrain changes of the tested conditions and the differences in outlet flow. The consistent temporal shape of the bottom panels in contrast with the top panels shows that the differences are associated with the model setup and domain specific errors for this exact computation. Since the nominal differences do not impact the temporal shape of the graphs, the on-terrain weir approach seems to be a main driver of the errors. Further support for this theory can be produced when calculating the averages of the data displayed in Figure 3 as shown in Table 1.

Table 1. Tested condition average differences and average percentage differences.

	T1		T2		T3		T123
Same data as for Figure 3.	Stage	Flow	Stage	Flow	Stage	Flow	Stage	Flow
Average difference (m) or (m^3/s)	0.04	6.63	0.06	9.92	0.06	10.08	0.03	5.03
Average percentage difference	0.04%	0.08%	0.07%	0.16%	0.08%	0.17%	0.03%	0.07%
NSE	0.997	0.984	0.994	0.960	0.994	0.967	0.998	0.989

Note. The difference and percentage difference stage and flow hydrographs for the tested geometry edits were averaged. T1 was the closest to the drainage outlet increasing in distance sequentially, and T123 includes all three geometry edits.

The percentage difference data shown in Table 1 was computed by dividing the differences by the full edit condition value for each timestep producing the percentage with the convention that the full edit condition is true and then averaging that time series. As the results of Table 1 show, the average nominal difference is quite small, but most importantly the average percentage difference is at or less than 0.77% in all cases and is extremely consistent across the runs showing the deltas to be a modelling setup error rather than stemming from the modularization. In addition, the near zero percentage average difference rate on the stage values further supports the application for inundation mapping, and by extension real-time flood warning. These results show a strong case in favour of this methodology, even in this crude discrete implementation, being able to be quite successful when building a scalable ML training dataset with these potential errors being kept in mind from the start.

4.2. Harris county case study for abstract modules

As stated in the previous sections, the flow averaging of the physically implemented discrete module boundaries was a major driver of errors in the methodology. To provide a case study addressing this, the abstract implementation was performed, and the results are displayed below.

Figure 4 shows the depth hydrographs of the lower, middle, and upstream portions of the tested module labelled points C1, C2, and C3. The flood depth difference map is coloured such that all values between −0.15 and 0.15 m are blue, and the few outliers are either green or red. In these hydrographs the existing condition run for the Tax Day hindcast can be seen in light grey labelled as existing and the two edited conditions share the same naming conventions as the discrete test. In all of these graphs, it is evident that there are significant differences from the existing condition and the full and module edited conditions. This shows that the conversion of the golf courses in mid-stream Brays Bayou to detention ponds decreased peak flows during the Tax Day event, which is consistent with common practice in urban hydrology. An important point that can be observed in these figures is that the module edit and full edit are still not perfectly identical, although they may seem so from the image. This is due to the lack of energy conservation across the cut boundary. Since the external boundary condition can only be set as a flow hydrograph with an Energy Gradient (EG) slope to distribute (Brunner, 2021), only the flow, and by extension mass, from the shallow water equations is conserved across the boundary. The velocities of these flow time series are not identical to the original dataset since only the flow rate and a static EG slope were input. This can be observed with the mild amount of time delay and peak drop shown in Figure 4. This is a problem that can be corrected further through a decomposition of the face velocities and flow rates to calculate the average EG slope for each cell face and adding that into the automated processes. However, the differences observed between the full model and this abstract module are far less than the errors of the original model, so it was decided to show the results and address those final differences later if needed for the full surrogate development.

Figure 4. Spatial maximum and hydrograph outputs of the flood stage.

To show a distributed comparison of the module to the full domain, a table was generated using the WSEL time series at each cell in the edit condition module so as to compare it to the same cell in the edit condition for the full domain, which can be seen in Table 2.

Table 2. WSEL meta-analysis of individual cell errors.

	Tax Day 2016			Harvey 2017
RMSE (m)	RMSE	NSE	KGE	RMSE	NSE	KGE
Mean	0.100	0.952	0.973	0.159	0.967	0.978
Standard deviation	0.230	0.609	0.093	0.310	0.195	0.068
Sample variance	0.053	0.370	0.009	0.096	0.038	0.005
Minimum	0.000	−16.686	−0.360	0.001	−2.371	−0.146
Maximum	3.429	1.000	1.000	3.926	1.000	1.000

Note. The stage hydrograph at each computation cell was compared for RMSE, NSE, and KGE. The table displays the resulting statistics of that dataset.

After observing some significant outliers for all the commutated cells, some of those outliers can be seen in Figure 4.

These results show that modularization of the 2D domain can in practice reproduce the results of a solver with only a subset of the domain, meaning that a training dataset produced from a mass conservative solver can be updated for edits with minimal computation outside the locations where the solution will change. As a result, an application of this modularization is to facilitate the deployment of ML-based surrogate methods. Modularization can support this aim by reducing the computational cost of updating hydrodynamic model outputs, which ML models use as training data, as regional land use patterns shift. This can facilitate Monte Carlo risk assessments and advanced real-time flood warning and information systems.

5. Summary

Modularization is a way of leveraging the fact that data storage is far cheaper than computation. Using stored internal solve data on the faces of a finite volume solver can allow us to construct a boundary condition from this data. A set of these boundary conditions can be used to solve a small subset of the domain where edits are occurring (with some conditions), called a module. As the name implies, the module solution can be used in place of re-computing the full domain and will produce results with errors smaller than the original model making it viable for practical applications like flood warning systems or regional flood mapping. Owing to the minimization of the update cost (as only the edited portions are re-computed), dynamic systems like urban landscapes become worth the computational investment of developing large regional models of un-simplified physics-based solvers.

The methods shown in this article are intended to push flood warning technology forward now and enable connected fields to develop uses like ML based H/H surrogate models. Some fields, such as hazard research and risk assessment, already have a need for fast cheap fully distributed hydrodynamic outputs to enable more data intensive statistical modelling like Monte Carlo methods. Modularization could allow research to leap forward with feasible large-scale probabilistic modelling for areas like water resources management, emergency planning, or structural optimization. On top of these benefits, implementing these with ML methods for H/H can foster shorter time horizons for interconnections like real-time traffic alerts (Panakkal et al., 2019) and real-time decision support services for critical infrastructure. Pushing forward ML based flood warning systems that are implementable and scalable will have major short- and long-term impacts across many fields, and the modularization of the hydrodynamic model will make it feasible today.

Disclosure statement

No potential conflict of interest was reported by the authors.