A balanced watershed decomposition method for rain-on-grid simulations in HEC-RAS

Since the mid-nineties (Yang et al. 1998; Rodriguez et al. 2008), nonlinear partial differential equations (PDEs) and additional physical parameterizations have been considered to model more and more complex hydrological processes by using an increasingly accurate representation of topography. One-dimensional river models were replaced by two-dimensional (2D) hydrodynamic models operating in the riverbed. Recently, rain-on-grid software accounts for precipitation over the watershed as a source term. Best practices of hydrological models are concerned with the management of water resources, the forecasting of floods, the assessment of hazards (Falter et al. 2016; Fan et al. 2016), the simulation of extreme event scenarios (Sorribas et al. 2016), and landscape modifications (Bayer & Collischonn 2013), for instance.

Rain-on-grid software implements the solution of the nonlinear PDEs by means of a finite element or a finite volume method, for instance, and some structured or unstructured mesh that describes the watershed topography. Computing time and memory usage may be important as they depend on the catchment surface, the mesh size, and the stiffness of the nonlinear discrete equations to be solved in an iterative manner. A review (David & Schmalz 2021) of software, models, and parameterizations based on 20 very different rain-on-grid studies reports simulations comprising mesh elements ranging from 5,000 to 3,000,000 and run on appropriate (personal or parallel) computers. For the sake of generality, simulations are carried out by means of the ‘Hydrologic Engineering Center's River Analysis System’ (HEC-RAS) in the paper.

Two computing strategies may be implemented to balance the computational load over a multi-core computer. On the one hand, a sparse parallel solver may be implemented as in the HEC-RAS software (Brunner 2016) to split the linear system into small blocks in order to distribute them on different cores. This approach is particularly well suited to a shared memory multi-core processor. On the other hand, a domain decomposition method (DDM) may be set up on a partition of the catchment. GIS tools are used to split the watershed into a number of hydrological sub-units (Apostolopoulos & Georgakakos 1997; Vivoni et al. 2011). The PDE problem is decomposed into a number of boundary value sub-problems accordingly. A DDM then organizes the solution by gathering the solutions computed for the sub-units and by setting relevant boundary conditions at the interfaces. This approach is efficient on a shared memory multi-core processor and on a multi-processor cluster, assuming a good balance of the computational effort.

These computing strategies can be combined. A watershed partition is needed to implement a DDM. Then, boundary values corresponding to each of the sub-domains can be solved using a sparse parallel solver. The efficiency of this computing approach depends on the performance of its components. In particular, the watershed partitioning method should take into account a digital elevation model (DEM) and user-prescribed points for the modeling of the topography, gauging information for validation purposes, and DDM requirements for computational performances. State-of-the-art partition methods from Apostolopoulos & Georgakakos (1997); Vivoni et al. (2011); Pontes et al. (2017) propose to split the stream at confluences, then with respect to the reach length, ignoring any information about the area drained by a hydrological sub-unit. In the present paper, we use the flow accumulation information to propose a balanced watershed partition and decomposition method adapted to large rain-on-grid simulations.

The paper is organized as follows. A direct rainfall model as well as terrain and methodological constraints relevant to watershed partition methods are described in the Methodological framework section. We evaluate two state-of-the-art techniques and their limits. The section Area-balanced domain decomposition method introduces the area-balanced partition method and discusses a domain decomposition method adapted to rain-on-grid modeling. Our implementation builds on the HEC-RAS software and the HEC-RAS Controller tool. Rain-on-grid simulations carried out for two large catchments are reported and discussed in the section Results to demonstrate the abilities and performances of the proposed methods. The last section provides Conclusions.

The hydrological modeling of large-scale watersheds through the discretization of nonlinear PDEs usually involves a large number of variables and a long computational time.

Domain decomposition methods were proposed to face both memory and time shortages (Leandro et al. 2014; Kumar & Duffy 2015). Such a partition may also account for hydraulic structures that locally impact the flow in an important manner, gauging stations along the stream, or soil roughness heterogeneity, for instance.

This section introduces general equations for rain-on-grid modeling and evaluates memory requirements. State-of-the-art watershed partitioning methods are recalled to emphasize their capabilities and limitations, as well as the potential added value of such processes.

Rain-on-grid modeling is implemented by means of nonlinear PDEs. Without limiting the general scope, we focus on the two models implemented in the version 5.0 of HEC-RAS (Brunner 2010).

These are used to model:

• a flow Q_in entering in the domain at the inflow boundary ⁠:

  

  (5)

  where n is the unit normal vector to the boundary, s is the side boundaries, and the integral is taken over the boundary surface on which the boundary condition applies,
• an outflow condition at the outlet implemented as a water surface gradient ⁠:

  

  (6)
• a normal flux condition on the side boundary Ω_side:

  

  (7)

The nonlinear equations are discretized using an iterative scheme for the temporal integration and a finite element volume for the spatial component. This accounts for the watershed topography, which is provided as a DEM and a shapefile that delineates the area of interest. Other data are boundary values, Manning's roughness coefficient accounted in S_f, and hydrological data R.

It is worth mentioning that, using HEC-RAS, linear systems are solved by means of a sparse parallel solver that can afford medium-range discretization involving a few hundred thousand elements on a shared memory multi-core processor. In other words, due to memory shortage, the HEC-RAS sparse parallel solver abilities may be insufficient to perform rain-on-grid simulation on a unique multi-core processor when the discretization involves a very large-scale watershed and/or very fine mesh. Indeed, as discussed in David & Schmalz (2021), the application of the HEC-RAS rain-on-grid for the simulation of large-scale watersheds remains limited by the number of mesh elements in the computational domain. For example, the simulation of the Saar basin of 1,747 km² using a 25 m mesh generated about 3 million mesh elements, which exceeds the capacity of the available memory of the Intel Xeon 64 bits @3.80 GHz 6-core processor with 16 Gb of RAM memory we used.

In a nutshell, HEC-RAS 2D is mainly used on Windows platforms and the set-up of a multidomain simulation remains challenging. In this paper, a proof of concept is proposed to evaluate the benefits of the domain decomposition method. It is worth mentioning that HEC-RAS neither allows us to delineate watersheds nor performs a partition at confluences in an automatic manner.

Provided with a DEM and an outlet, GIS software tools such as ArcGis, QGIS, or TopoToolbox (Schwanghart & Kuhn 2010) allow for the delineation of the corresponding watershed and stream network through the calculation of a flow direction and a flow accumulation. As an illustration, the flow accumulation raster related to the case study DEM is plotted in Figure 1. It records the number of cells and, thus, the area drained by each of the DEM cells, outlining the stream network. The latter comprises a number of connection points called confluences and a number of links called river reaches.

Based on the confluence points, the basic watershed partition offered by GIS tools allows us to delineate sub-basins. When the stream network presents nearby confluences, this partition method may create very small sub-basins. On the other hand, long reaches may result in large sub-basins.

As proposed in Pontes et al. (2017) and implemented in IPH-Hydro Tools (Siqueira et al. 2016), the confluence method may be improved by dividing the longest river reaches to delineate smaller sub-basins of larger ones. However, this does not remedy the problem of very small sub-basins.

The former two methods may be viewed as zero-dimensional and one-dimensional since they are based on confluence points and reach lengths, respectively. This stream information cannot guarantee an area-balanced partition of the watershed, because the available flow accumulation information is not taken into account. This will be illustrated for a case study in the section Numerical results.

DEMs are only an approximation of the topography and may underestimate the significant impact of bridges and other engineered structures on the water flow.

Watershed partitions can account for civil engineering structures by setting inlets/outlets that are close to them or by comparing the computed and observed discharges when a gauging device is available. A partition also allows us to represent the heterogeneity of the terrain and the meteorological variability over the catchment.

From a computational point of view, finite element models that are applied to sub-units with a similar area involve linear systems with a similar number of unknowns. Thereby, their solution has a comparable computational cost, and the domain decomposition method is more efficient (Kumar & Duffy 2015).

The partition method we propose in the section Area-balanced domain decomposition method aims to:

• 1.

  account for engineered structures as user-defined outlets,

• 2.

  consider small, but homogeneous hydrological sub-units,

• 3.

  ensure a balance of the sub-unit areas,

• 4.

  obtain a prescribed number of sub-units, and

• 5.

  perform simulations with a reasonable number of unknowns per sub-unit.

The first two properties ensure a better quality of modeling. Properties 3 and 4 improve the efficiency of the simulation. The last enables the opportunity for rain-on-grid simulations on very large watersheds. The last three properties are not satisfactorily addressed by previous methods as reported in the section Results.

GIS are precious tools to support hydrological modeling. Software like ArcGIS or QGIS provides graphical user interfaces and programming interfaces with Python. In contrast to this, TopoToolbox (Schwanghart & Kuhn 2010) is a set of Matlab functions dedicated to DEM analysis, with a special focus on hydrology. An overview of the respective performance of these tools is proposed in Schwanghart & Scherler (2014). Among them, we choose TopoToolbox because it is integrated in a scientific computation framework with good efficiency.

TopoToolbox is an object-oriented toolbox provided with a user guide that describes the methods that are necessary to our partition tool. As figured out in listing 1, the DEM is read as a grid object. The flow direction and flow accumulation rasters are then computed to deduce the stream network as an ordered list of cells that drained a minimum number of cells and drainage basins, using or not input information (prescribed points, for instance). The minimum drained cells threshold is denoted by Delta.

State-of-the-art partition methods are based on the split method operated, first, with respect to stream confluences. The partition with respect to a prescribed reach length has been implemented using the distance property of the stream object to identify complementary points. These partition methods are exemplified and analyzed in the section Results.

Targeting rain-on-grid simulations for large-scale watersheds has motivated the design of the new partitioning method. To accommodate parallel computing in the near future, a special emphasis is put on area balancing by considering the flow accumulation raster (2D information) rather than the reach length (1D information) or the confluence information (point-wise information).

Put briefly, for each cell of the flow accumulation raster, we maintain the number of upstream cells it drains. This parameter is, thus, related to the area of its upstream watershed.

The workflow is illustrated in Figure 2. A preprocessing is set up to read the DEM and, possibly, user-defined points that represent the outlets and engineered structures of interest.

The area-balancing watershed decomposition method is organized as follows. The sub-unit with the largest area is selected (imax) and its flow accumulation information is computed. A splitting point is determined by searching for the cell that drains half of that flow (coor(iSplitPt); see Figure 3). This is added to the set of outlet points defining the drainage basins. This process is iterated until the desired stopping criterion is met. This can be based on the number of sub-units and/or on the sub-unit areas. In listing 2, the process stops when all the sub-basins have an area less than a prescribed maxArea.

The sub-watersheds are saved in the shapefile format using the shapewrite function of the Matlab Mapping toolbox such that it can then be used in the HEC-RAS simulations, for instance.

Depending on the actual topography or its digitization, the watershed delineation may produce sub-units that present isolated cells, notably near to the outlet; see Figure 4. Defects with respect to the 4-neighbor connectivity may occur for any partition method or partition software. These are not appropriate to a finite element or a finite volume modeling.

At the end of the partition process, a connectivity check may be performed to correct the decomposition by moving some of the outlet points upward. This can be achieved automatically by inspecting the properties of each sub-unit by means of the regionprops Matlab function, for instance. It is worth noting that the user can also take them into account in a channel, as proposed in the 2D–1D modeling described hereafter.

The RAS Mapper interface is used to model the terrain starting from a DEM as follows.

• 1.

  Sub-unit shapefiles are imported to implement an area-balanced domain decomposition method.

• 2.

  Each sub-unit is meshed to discretize the 2D PDEs, either (1)–(2) or (1) and (3).

• 3.

  A 1D river model is added at the sub-unit outlet to extract the output flow hydrograph (Equation (6)) as proposed in Goodell (2018).

• 4.

  This flow hydrograph serves as an input to the downstream sub-unit following Equation (5).

Note that this 2D–1D modeling of the sub-unit allows for the automation by means of the HEC-RAS Controller (Goodell 2014) that has been designed for 1D problems.

In a nutshell, the HEC-RAS Controller provides a collection of routines that allows us to open a project, to select and to run the so-called plans, and to extract output data.

The DDM sequence and the communications are then automated by means of a Virtual Basic for Applications (VBA) driver written using the HEC-RAS Controller functions.

The VBA script, see Figure 5, organizes the computation sequence by accounting for the natural water flow. It runs sub-unit models from the most upstream to accumulate the stream flow at the watershed outlet. Each sub-unit calculation is carried out in one go, providing all the necessary input flow data to its downstream sub-unit. The number of HEC-RAS plans corresponds to the number of sub-units.

The performances of the partition and domain decomposition methods are compared on two case studies: the Saar watershed in France (1,747 km²) and the Moderbach sub-watershed (89 km²). These are evaluated and analyzed on the generation of nested partitions (Case study 1) and the solution of 2D unsteady rain-on-grid simulations (Case study 2). HEC-RAS version 5.0.7 is used to model the hydrological behavior where the catchments are discretized using a mesh size of 25 m.

Run-times were recorded in minutes on an Intel Xeon 64 bits @3.80 GHz 6-core processor with 16 Gb of RAM memory. The sparse parallel solver of HEC-RAS makes use of the six cores of that processor.

The slightly undulating catchment of the Moderbach stream (89 km², Région Grand-Est, France) is represented in Figure 6. Elevations range from 210 to 305 meters above sea level (m.a.s.l). The DEM raster is issued from BD Alti^®25 m – V2 and is available for research purposes.

Due to the clayey soil, the watershed was converted into a part of the Maginot defense waterline in the 1930s, comprising 11 dams delineating six reservoirs on the Moderbach's tributaries (blue polygons) and five forebays along the main stream delineated by five dams (red lines) to be flooded. These structures that may strongly impact the flow are used to compare guided and unguided area-balanced partition methods.

This HEC-RAS version ignores infiltration. Such a strong hypothesis is plausible for this watershed since soils are quasi-impervious under winter conditions.

Figure 7 plots partitions computed with the two state-of-the-art methods (confluences, then a reach length of 2 km between two outlets) and the unguided area-balanced partitions.

The characteristics of the partitions are reported in Table 1. On the one hand, the confluence method and the reach length method produce very small sub-watersheds. In both cases, the ratio between the maximum and the minimum area is larger than 20. The balance of areas is not satisfactory. On the other hand, the ratio of 2 obtained for the novel partition method (U. Area) demonstrates the balance of the sub-watershed areas when similar mesh sizes are used for the different sub-units.

The computation time for the confluence method corresponds to the time needed to proceed the DEM following the listing 1. The time measured for the unguided area-balanced method comprises those of listing 1 and 2. Note that this difference of a few seconds is very small as compared to the computation time needed for rain-on-grid simulations.

Confluence method. The stream network is a set of connected segments, the number of which is parameterized by means of a threshold Δ on minimum drained cells, see Listing 1. The characteristics of the nested partitions produced with the thresholds are reported in Table 2.

As the number of confluences increases, there is an increasing probability of a very small sub-unit in the partition. Consequently, the ratio between the maximum and the minimum areas increases.

Guided area-balanced method. The guided method is carried out by starting from the seven user-defined outlets that are located at the pond outlets and the watershed outlet. The generated sequence of nested partitions is illustrated in Figure 8, and partition properties are reported in Table 3. One observes that the ratio decreases from 19.3 corresponding to the user-defined partition to 9.1 because the partitions into 7 and 8 sub-units involve very large sub-units.

Unguided area-balanced method. The unguided method only requires the watershed outlet as an input. The nested partitions into a number of sub-units (that can be a power of 2) are plotted in Figure 8. Corresponding properties are reported in Table 4. This shows an increase of the ratio from 1.1 to approximately 2.9 since the first two sub-units have almost the same area.

To summarize, the two area-balanced methods provide similar performances when the number of sub-units is large. It should be mentioned that in these two methods, the minimum drained cell threshold Δ is set to 1,000 allowing us to reach a large number of sub-units. The unguided method outperforms the confluence method in any case.

The catchment is partitioned, then meshed using about 14,5703 elements. The domain decomposition performances are evaluated by simulating a precipitation event of 28 h with an amount of 80 mm/day to saturate the watershed in a reasonable time. An additional 44 h are carried out with no precipitation in order to drain the watershed until reaching a null discharge at the outlet. The Manning's coefficient is set to 0.03. The boundary condition at the 1D river uses a normal depth estimated using the RAS Mapper tool: the slope parameter is set to 0.001. The maximum flood extent is plotted in Figure 9.

Presented in Figure 10(a) are the discharge results computed at the outlet with HEC-RAS and those obtained with the DDMs (workflow of Figure 5) involving nested partitions of 2, 4 and 8 sub-units, respectively. Not shown are very similar results with the confluence partitioning method.

The relative errors displayed in Figure 10(b) exhibit very small discrepancies. These may be due to the value of the slope parameters, topographical effects, and the partition method, since the null flux boundary condition (7) is applied onto the boundaries ∂Ω_side for the catchment and the sub-units, respectively.

The computational efficiency of the DDM is evaluated for a number of partitions carried out with respect to either the confluences or the (unguided) area balancing.

To evaluate the performances of the sparse parallel solver of HEC-RAS and the general behavior of the DDM, an overview of the computation times with respect to the sub-unit areas and the partition methods is provided in Figure 11. In Figure 11(a), points of the same color and of the same shape correspond to sub-units of the same partition. Triangles and filled circles are used to the unguided area-balanced partition method and the confluence partition method, respectively.

As expected, points are almost aligned in this log-log representation. The computation times and the areas are linked by a power law, the slope of which reflects the performance of the HEC-RAS sparse parallel solver. The differences in the solution times for the sub-units of a very similar area indicate that the latter is a key parameter, but not the only one in this nonlinear solver.

Moreover, Figure 11(a) exhibits three colored clusters (filled circles of the area-balanced method) with a very similar area and computing time. In contrast to this, the dispersion of the triangular points is much higher for the confluence method and very small sub-units are created.

The total time and the speed-up are plotted, Figure 11(b) and 11(c), against the number of sub-units. We see that for both partitioning methods, the decrease in computing time is multiplicative. As we have seen, the computation of the linear system generally grows super-linear in the size of the parts; thus, with Jensen's inequality (Jensen 1906), the sum of the times of the linear systems of all parts is much smaller than for the linear system of the original watershed.

Additionally, Figure 11(b) shows regressions of the form ax^b, where b is −0.27 and −0.32, respectively. This clearly shows the advantage of the partition method by a DDM in general, but also that the unguided balanced partition method improves the confluence method by an important factor.

In a first attempt, we designed a HEC-RAS plan modeling the Saar catchment in France with an element size of 25 m, i.e. involving about 2,800,000 elements, to verify the occurrence of a memory overflow. Such an issue has motivated the development of both the area-balanced partition and decomposition method.

Originating in the Vosges Mountains (Région Grand-Est, France), the Saar River (246 km long) drains an area of 1,747 km² before it enters the German territory at the confluence with the Blies River located in Sarreguemines. Elevations range from 800 m.a.s.l. down to 243 m.a.s.l. Twelve rasters issued from BD Alti^®25 m – V2 are concatenated into a unique raster comprising 12 × 10⁶ cells. Precipitation data are from the Kappelkinger station (Météo-France).

Eleven gauging stations are distributed along the main stream, the data of which may be found on http://www.hydro.eaufrance. Among them, we select the stations located in the two main cities, namely Sarralbe (Station A9200100) and Sarreguemines Centre (Station A9311050), to compare computed and measured discharges for the extreme flood event that occurred from 22nd December to 24th December 2010.

A flood map and the related shapefiles were produced by the SERTIT (Service régional de traitement d'images et de télédétection, Strasbourg) from SAR Cosmo-Skymed images (5 m resolution). It is worth noting that these images cause an underestimation in urban and forested areas. Moreover, the permanent waterbodies (river and ponds) and very small water objects (less than 100 m²) were not reproduced. The shapefile of the flood extent on 22nd December is used for comparison with our simulation results.

Computational performances. The catchment is split into 20 sub-basins using the guided area-balanced method (Figure 12) to account for the 11 gauging stations located on the main river (red dots). Other sub-basin outlets are represented using green dots. The sub-basins range from 11 to 164 km² in area (see Table 5), 17,253 to 26,2183 in mesh elements and 15 to 405 min in computational time. The simulation organized following the workflow of Figure 2 was completed in 86 h. As expected, the partition time is several orders of magnitude smaller than the simulation time.

Table 5 also reports the characteristics of an area-balanced partition method involving small sub-basins (an average area of 13.3 km²) and a partition at the confluences, both performed from the same stream network. Here again, the ratio is in favor of the area-balanced method since the confluence method generates very small sub-units.

Hydrological results. The flood event of 22–24 December was preceded by heavy precipitation a few days before (50 mm, 6–9 December) that saturated the soil. This makes case study 2 consistent with the 2D rain-on-grid simulation in HEC-RAS that does not consider the infiltration effect.

In Figure 13, the precipitation histogram (blue bars) and the discharge data (blue and black dash-dotted lines) provided at Sarralbe and Sarreguemines gauging stations are plotted against the discharge results (blue and black solid lines) computed at the same locations. Table 6 provides complementary information on the observed and computed water volumes, Nash–Sutcliffe efficiency (NSE) criterion (Nash 1970), the root mean square errors (RMSE), and the standard deviation (SD) computed for the timeseries.

First, relative error lower than 12% for the water volume comes in support of the quasi-imperviousness assumption during that event. Second, both NSE values are close to 1, which represents the perfect match (Nash 1970). Third, according to Singh et al. (2005), a ratio of RMSE/SD lower than 0.5 as in Table 6 indicates good modeling. These indicators allow to conclude to a very good agreement between the data and the calculated discharge and validate our modeling for that event.

The observed flood extent (FO) and simulated water surface (FS) are compared in Figure 14 by representing the common flood extent (A = FO ∩ FS) in gray. The complementary sets (B = FS  /FO) and (C = FO /FS), plotted in blue and magenta, respectively, highlight the differences. These are notably observed in coniferous/forested areas that were not analyzed with teledetection. The partition points, including the gauging stations, are plotted as dots.

This section discusses advantages, limits, and applicability of the proposed methods from the hydrological and computation points of view. Some outlooks are raised.

Some user-defined outlets may be set to facilitate the account for civil engineering structural or topographical effects in the simulations. Complementary outlets may be positioned in an automatic manner either at confluences or at points determined with respect to particular criteria constraints without, however, fully addressing the area-balancing issue as demonstrated in Table 2.

One of the novelties of the paper is to organize this additional phase in order to delineate sub-units with a similar area using the flow accumulation information. As demonstrated in Tables 3 and 4, both versions (user-guided or unguided) of our automated area-balanced partitioning method perform much better. Moreover, these sub-units can be used with a distributed hydrological model in an efficient manner from both the computational time, see Figure 11, and the memory requirements.

The DDM method allows us to build simulation plans with homogeneous data over the sub-units, but with heterogeneity over a large watershed. In particular, this may facilitate the account for more actual precipitation and soil properties (Manning's coefficient). It is worth mentioning that the Saar catchment is submitted to Cfb climate (Köppen classification) that is with temperate oceanic climate. In particular, winter precipitation events occur at a regional scale. Thus, heterogeneity in the precipitation amount is not implemented to carry out fair hydrological and methodological comparisons.

This HEC-RAS version ignores infiltration. As demonstrated in Figure 13, such a strong hypothesis is plausible for this watershed since soils are quasi-impervious under winter conditions. Such an assumption may also agree with surface sealing of soil in urban areas (although additional factors could be considered) or dry soils. To overcome the imperviousness assumption, a SWAT/HEC-RAS integrated modeling (Zeiger & Hubbart 2021) could be considered, for instance.

A 2D PDE solver is time- and memory-consuming, especially when the domain under study becomes large. This applies to the 2D unsteady rain-on-grid solver recently proposed within the HEC-RAS software. It processes linear systems by means of a sparse parallel solver to take advantage of a unique multi-core processor. This software is, thus, a perfect candidate to assess the benefits that could be brought by an additional partition of the watershed into hydrological sub-units.

The impact of the improved balance in the sub-unit areas on the 2D PDE computation is assessed using HEC-RAS in Figure 11. For our first case study, the DDM workflow described in Figure 5 resulted in a speed-up of 2 when eight sub-units of similar size are considered. Performances are a little lower with the confluence method.

For our second case study, the DDM workflow allows us to run simulation involving a larger number of unknowns (2,800,000) as compared to previous rain-in-grid simulations (David & Schmalz 2021) carried out with HEC-RAS. This is a major improvement.

HEC-RAS plans should be built using the graphical user interface. The multiplication of the sub-domains requires some initial manpower to set up 2D–1D models. Nevertheless, the cost–benefit ratio becomes important as soon as the rain-in-grid model is used in an operational manner and run many times as in a data assimilation process.

Local mesh refinements that could impact the load balance in a parallel implementation will be considered in the near future.

For now, the VBA code developed using the HEC-RAS Controller functions allows us to automate the DDM/HEC-RAS coupling strategy. The implementation on a multi-processor cluster could bring an additional reduction of the time to results. For the moment, this strongly depends on the maximum number of the sub-units in the upstream-to-downstream sequence.

Assuming the availability of the HEC-RAS sources, the DDM/HEC-RAS approach could benefit from two levels of parallelism. First, the sparse parallel solver of HEC-RAS performs the sub-unit computation on a multi-core processor. Second, the distribution of the sub-units on a multi-processor cluster may then further improve computing times.

At a watershed scale, rain-on-grid simulations for the modeling of 2D unsteady flows in response to precipitation input are implemented in the HEC-RAS software by means of a finite volume method, a discretization of the topography, and a sparse parallel linear solver running on a unique processor. Obviously, this may fail due to memory shortage when the watershed of interest and the chosen discretization yield too large linear systems. Such large simulations are good candidates for the design of dedicated partitioning and domain decomposition methods for both the modeling of very large catchments and the reduction of the overall computational time.

Thinking in a load-balanced domain decomposition method, we first propose an area-balanced partition method carried out by considering the flow accumulation information rather than the confluence points and/or the reach length constraint. Second, the domain decomposition method is exemplified for the HEC-RAS software. Special attention has been brought to the set-up of 2D–1D models that combine 2D sub-unit meshes to 1D channels for the extraction of outflow data from upstream sub-units and their propagation to their downstream unit. The workflow is automated by means of the HEC-RAS Controller.

For both case studies, the new partition and decomposition methods have demonstrated a better equilibrium for the sub-basin areas. As a consequence the computational load that is required for the solution of the related sub-problems, respectively, is balanced between the parts. As expected, the total computational time is substantially reduced when the domain decomposition method is used.

Rain-on-grid simulations over the Saar watershed (1,747 km²) with a mesh size of 25 m cannot be carried out using HEC-RAS in a straightforward manner because these would yield 2,800,000 mesh elements, overflowing the 16 Gb of RAM memory. The proposed domain decomposition method allows for the solution of such large rain-on-grid simulations, while also reducing the overall computational time. From a hydraulic point of view, resulting discharge and flood maps show very good agreements with public discharge data and a flood extent map produced from SAR Cosmo-Skymed images, respectively.

The authors are grateful to H. Yésou (SERTIT) for the communication of the flood extent shapefiles.

All relevant data are included in the paper or its Supplementary Information.