Flood analysis comparison with probability density functions and a stochastic weather generator

Case study

Morelia, the largest and most populous urban area in Michoacán, Mexico, is characterized by its significant population and complex fluvial network. The city has a fluvial network of fifteen perennial tributaries, as shown in Fig. 1, and is the object of analysis in this work. The Rio Grande is the city’s main river channel, collects the runoff generated in the area, and is regulated upstream by the Cointzio dam.

Morelia is prone to frequent floods from overflowing rivers, which directly impact houses and disrupt residents’ daily lives. For instance, in 2018 a severe storm caused injuries, inundations in homes, landslides, and wrecked vehicles as shown in Fig. 2. Today, floods are perceived as having only negative consequences and as disasters that inflict harm, sometimes irreparable. These impacts are aggravated by the expansion of cities that are naturally prone to flooding, such as near rivers and streams.

Methodology

The proposed methodology utilizes maximum runoff results for various return periods within the hydraulic model of Morelia City’s drainage network. It also incorporates scenarios created with the daily MASVC (Hernández-Bedolla et al., 2023). Moreover, the HEC-RAS model was adjusted to obtain the flood areas based on topographic studies and changes in the Manning roughness coefficients (Demir & Keskin, 2020; Yalcin, 2020). The methodology is illustrated in Fig. 3.

Surface runoff

Surface runoff was calculated using two types of methodology. One method was the probability density functions with the Soil Conservation Service Curve Number (PDF-SCS-CN), which is used to estimate runoff from small-to-medium-sized watersheds SCS-CN method (AL-Hussein et al., 2022; Ansori, 2023; Sathya, Thampi & Chithra, 2023). The other method was the stochastic weather generator, which applied the MASVC-SCS-CN approach (Hernández-Bedolla et al., 2023). The MASVC is a multivariate stochastic model that uses lag one autoregressive multivariate parameters. It consists of two modeling phases. The first phase models the occurrence of precipitation (wet-dry), and the second phase estimates the amount of precipitation on a daily scale. Then, the maximum precipitation values are extracted, and the runoffs are computed using the SCS-CN method. The observed daily precipitation was obtained from the National Water Commission (CONAGUA). The weather stations used were 16,022, 16,055, 16,114, and 16,247 (Fig. 1). The historical period of record was from 1980 to 2009, and it is available for the four rainfall stations (https://smn.conagua.gob.mx). In addition, hourly precipitation for October 22, 2018, from the Automated Weather Stations (AWS) database, available at https://smn.conagua.gob.mx/es/observando-el-tiempo/estaciones-meteorologicas-automaticas-ema-s, was used. The hydrographs for the modeling were derived from the maximum flows suggested by Hernández-Bedolla et al. (2023).

Flood inundation

Hydrodynamic model HEC-RAS is used to model rivers (one-dimensional), flood areas (two-dimensional), channels, drains, and dams (Açıl et al., 2023; Bharath et al., 2021; Bush et al., 2022; Goswami, Prasad & Kumar, 2023; Namara, Damisse & Tufa, 2022; Ongdas et al., 2020). HEC-RAS calculates the depth of inundation and velocity based on floodwater discharge hydrographs (Hydrologic Engineering Center, 2009). HEC-RAS was applied for the two-dimensional hydraulic model and was solved by the 2-D shallow water equations. The solution method for both equations was by volume finite differences. The unsteady differential form of the mass conservation equation is Eq. (1). (1)${\displaystyle r=\frac{\partial H}{\partial t}+\frac{\partial \left(hu\right)}{\partial x}+\frac{\partial \left(hv\right)}{\partial y}}$

where r is the source/sink term, H is the surface elevation (m); t is the time, hu and hv are the flow in x and y (m²s⁻¹); HEC-RAS recently incorporates fully 2-D shallow equations (Costabile et al., 2020). The shallow water equations are presented in Eqs. (2) and (3).

(2)${\displaystyle \frac{\partial p}{\partial t}+\frac{\partial }{\partial x}\left(\frac{p^2}{h}\right)+\frac{\partial }{\partial y}\left(\frac{pq}{h}\right)=-\frac{n^{2}pg\sqrt{p^2+q^2}}{h^2}-gh\frac{\partial H}{\partial x}+pf+\frac{\partial }{\rho \partial x}\left(h{\tau }_{xx}\right)+\frac{\partial }{\rho \partial y}\left(h{\tau }_{xy}\right)}$ (3)${\displaystyle \frac{\partial q}{\partial t}+\frac{\partial }{\partial x}\left(\frac{pq}{h}\right)+\frac{\partial }{\partial y}\left(\frac{q^2}{h}\right)=-\frac{n^{2}pg\sqrt{p^2+q^2}}{h^2}-gh\frac{\partial H}{\partial y}+qf+\frac{\partial }{\rho \partial x}\left(h{\tau }_{xy}\right)+\frac{\partial }{\rho \partial y}\left(h{\tau }_{yy}\right)}$

where p and q are hu and hv, (m²s⁻¹); n is the Manning’s roughness coefficient (s m^−1/3), g is the gravity acceleration (ms⁻²), ρ is the water density (kg m⁻³), τ_xx, τ_yy, and τ_xy are the components of the stress tensor and f is the Coriolis parameter (s⁻¹). We propose the diffusive wave algorithm to select boundary conditions and response times. In this case, the inertial terms of the equation are neglected.

The Eulerian Shallow Water equation (EM-SWE) was the calibration solution method. This method utilizes the momentum-conservative discretization assuming local conservation of momentum about control volume centered on all cell face $\left(V\cdot \nabla \right)u_N$ (Hydrologic Engineering Center, 2021). (4a)${\displaystyle {\left[\left(V\cdot \nabla \right){\mathrm{u}}_N\right]}_f\approx {\left(\frac{1}{h}\left[\nabla \cdot \left(hV{u}_N-u_N\nabla \cdot hV\right)\right]\right)}_f}$ (4b)${\displaystyle {\left[\left(V\cdot \nabla \right){\mathrm{u}}_N\right]}_f\approx \frac{{\alpha }_f^L}{{\overline{h}}_f{A}_L}\sum _{k\in K\left(L\right)}{Q}_{i,k}^{-}\left(V_k^u\cdot n_f-u_{N,f}\right)+\frac{{\alpha }_f^R}{{\overline{h}}_f{A}_R}\sum _{k\in K\left(L\right)}{Q}_{R,k}^{-}\left(V_k^u\cdot n_f-u_{N,f}\right)}$

where $h_i={\Omega }_i/A_i^w$; $Q_{i,k}^{\text{\_}\text{\_}}$ is the minimum value of inflow at face k to cell i s_i,k, Q_k, 0;  A_k is the face vertical area; A_L is the left cell horizontal area, A_R is the right cell horizontal area, V_j is the cell average current velocity vector $V_k^u$ is V_L for u_N,k major than cero and is V_R for u_N,k less than zero; n_f is the face-normal unit vector.

DEM and land use

The quality of flood modeling depends on the topographic information. A hybrid mesh was used, combining a gridded DEM dataset with five m of resolution from the National Institute of Statistics and Geography (INEGI, https://www.inegi.org.mx/temas/topografia/) and more detailed 0.5 m data from a specific topographic study to improve the accuracy in rivers. The land use data was obtained from INEGI (https://www.inegi.org.mx/temas/usosuelo/), and field visits to the rivers were conducted to determine the Manning coefficient (n).

Boundary conditions

Two-dimensional modeling was done with flexible and triangular meshes. The meshes gave more detail in the rivers due to a high-accuracy topographic survey (Muñoz et al., 2022; Ontowirjo et al., 2023). For the flood areas, a regular grid with a 5-meter resolution was proposed for roads, as well as for urban and rural areas. Boundary conditions describe how water behaves at the model domain’s boundaries, including different conditions at the model edges (Hydrologic Engineering Center, 2009). Upstream boundary conditions are needed at the upstream end of all reaches that do not connect to other reaches or storage areas. Hydrographs corresponding to various return and historical periods were assigned upstream of each sub-basin studied. For the downstream boundary condition, the normal depth option was chosen for the HEC-RAS model (Hydrologic Engineering Center, 2021).

Simulation and calibration

The process of calibrating the hydraulic model is based on the different steps described below. (1) Simulation of the entire study area to identify the simulation times and potential flooding zones near the streams and rivers. (2) Simulation of the diffusive wave algorithm to generate the first flood from historical rainfall data. (3) Calibration using the EM-SWE by changing the Manning roughness coefficient from 0.01 to 0.75, and field visits to perform mesh refinement. (4) Validation of the results from the actual flooding for a specific date. In this case, validation was conducted by juxtaposing the simulation results with an orthophoto that documented the flooding event which occurred on October 22, 2018 (Roblero-Escobar et al., 2025). (5) Generation of floods from different return periods and results from MASVC-SCS-CN and PDF-SCS-CN. For the MASVC, we generated 1,000 series of the same length as the observed rainfall. For PDF, we consider the annual maximum rainfall. Subsequently, we estimated the surface runoff based on curve number “CN” and calibrated de hydraulic model changing n’s Manning.