Modeling bioinfiltration surface dynamics through a hybrid geomorphic-infiltration model

drag coefficient

D

median particle size

accumulated depth of infiltrated water

P

pressure

Rep

particle Reynolds number

R

submerged specific gravity

velocity vector

infiltration rate

g

gravitational acceleration

mean flow depth

saturated hydraulic conductivity

Manning's roughness coefficient

time

dimensionless bed transport rate

wetting front suction head

moisture content

density of water

ν

kinematic viscosity

dimensionless shear stress

dimensionless critical shear stress (shield parameter)

density of sediment

Urbanization and the creation of impervious surfaces typically result in increased stormwater runoff, elevated risk of flooding, and instability in urban streams due to substantial bank erosion (Wolman & Schick 1967; Leopold 1968; Hammond et al. 2015; Anim & Banahene 2021). These drastic changes within watersheds and the associated water quality problems in urban streams often have far-reaching consequences such as disruption of the stream's ecosystem and disconnection from the floodplain (McGrane 2016; Waite et al. 2019). Over the past two decades, innovations in urban stormwater engineering have led to the implementation of stormwater management practices (SMPs) such as bioinfiltration systems. Such systems are effective in reducing runoff volume resulting in less impactful stormwater runoff discharge into receiving water bodies (Davis et al. 2012; Lee et al. 2016). Bioinfiltration systems have also been shown to be beneficial for significant water quality improvements for both organic and inorganic pollutants (Davis et al. 2003, 2009; Jeon et al. 2021). The bioinfiltration system's ability to treat total suspended solids (TSS) is particularly effective as a significant quantity of sediment is retained in these systems (Smith et al. 2023). The latter aid municipalities in meeting environmental regulations, resulting in bioinfiltration systems becoming an increasingly important and pervasive feature in the urban landscape. Computational modeling creates an exciting stage for innovation and a deeper understanding of these nature-based systems within the built environment.

The performance of the system over time can be heavily impacted by erosion or deposition of significant amounts of sediment (DelGrosso et al. 2019). In recent years, there has been a growing interest in understanding the impact of sedimentation on the performance of bioinfiltration systems to ensure the sustainable function of the infrastructure (Coustumer et al. 2012; DelGrosso et al. 2019; Homet et al. 2022). Sediment removal can lead to one of the most expensive maintenance costs associated with bioinfiltration systems (Hunt et al. 2020). Clay and silt-sized particles within suspended sediment can lead to clogging the infiltration media and reduce the saturated hydraulic conductivity (k_sat) of the soil (Langergraber et al. 2003; Asleson et al. 2009; Virahsawmy et al. 2014). Many jurisdictions have thus implemented SMP manuals or design standards aimed at mitigating erosion, channelization, and sedimentation within these systems (e.g., Philadelphia Water Department 2021).

Prior studies regarding bioinfiltration soil media have primarily focused on the reduction of k_sat as a result of the introduction of sediment (Gonzalez-Merchan et al. 2010; Coustumer et al. 2012) and not on the surficial movement of sediment. However, the ability to link bioinfiltration to surface morphology is crucial for a more sustainable design and can open the door to targeted maintenance strategies for such systems. This is particularly true for specific geometries of bioinfiltration systems, such as those designed for conveyance. Specifically, understanding the shear stress distribution and sediment flux within a system allows for a more effective prioritization of maintenance of systems with a high risk of failure or underperformance due to sedimentation or erosion.

Existing models for sediment transport, such as the ROWERO (Wu & Meyer 1989) and the surface irrigation model by Mailapalli et al. (2013), are better suited for irrigation furrows, not a vegetated bioinfiltration basin. Currently, geomorphic models do not exist for bioinfiltration systems; however, advancements in fluvial geomorphic modeling (Nelson 2018) provide a starting point for modeling sediment flux in bioinfiltration. There are challenges in applying fluvial geomorphic models to bioinfiltration because they do not account for the high infiltration rates typical in such systems. Infiltration within bioinfiltration systems constitutes a significant volume of water removal (Davis et al. 2009) and, therefore, cannot be ignored when estimating hydrographs or predicting geomorphic changes.

The International River Interface Cooperative Software (iRIC) (Nelson 2018) is a public domain software platform for numerical simulation of river flow and morphodynamics. The iRIC platform hosts a number of two- and three-dimensional solvers for flow, sediment transport, and bed evolution modeling, and has the advantage of both pre- and post-processing of input and results, respectively (Nelson et al. 2010). The flow and sediment transport with morphologic evolution of channels (FaSTMECH) (Nelson et al. 2003), one of the river modeling packages hosted on the iRIC platform, has been successfully implemented in a variety of studies including predicting morphological alterations in backwater zones in rivers (Hosseiny & Smith 2019), quantifying and visualizing uncertainty associated with flood inundation maps (Zarzar et al. 2018), explaining spatial patterns of mussel beds (May & Pryor 2015), investigating chute formation along meandering streams (Harrison et al. 2015), and modeling the effects of flow and sediment transport on fish spawning habitat (McDonald et al. 2010; McKean & Tonina 2013). This model was chosen for modeling the bioretention due to the available solvers and its ability to represent flow and sediment transport in two dimensions.

Here, we present a novel framework for analyzing sediment dynamics within a bioinfiltration system. This new framework combines the Green-Ampt infiltration equation (Green & Ampt 1911) with the FaSTMECH sediment transport model for an estimation of both infiltration and sediment transport. The advantage of using a Green-Ampt infiltration model is its wide acceptability as a suitable infiltration model for bioinfiltration systems, with easily obtainable input parameters compared with other infiltration models (Bedient et al. 2008).

This study introduces a predictive framework for identifying potential high shear stress (erosion) and high sediment deposition (sedimentation) areas within a bioinfiltration system, which will allow for prioritizing areas in design and maintenance. By integrating different mechanistic models iteratively, a representation of bioinfiltration surface dynamics can be produced (Chappell et al. 2006; Lastra et al. 2008). The increasing availability of open-access sediment models and several well-established infiltration models makes this approach computationally less expensive and more feasible (Garcia 2007). The objective of this study is to create a framework to examine morphology within a bioswale bioinfiltration system to identify areas at risk of sedimentation or erosion through combining an existing fluvial geomorphic modeling tool with a widely accepted infiltration model.

Flow into SMP A was monitored by a Blue Siren area/velocity flow sensor installed at all three inlets. Internal flow depth was monitored by a CS451 pressure transducer installed upstream of each weir. Flow measurements were recorded at 5-min intervals from the beginning of each storm event. The combined flows from the N10 inlet and the downstream weir were used in the simulation (Figure 1). Rainfall data were collected from a tipping bucket rain gage located at the southern end of the site. The study used data collected between May 2017 and April 2019 for simulation and model verification. For the purposes of hydraulic simulation, each hydrograph was treated as a separate event. An event was considered to have ended when the flow of the hydrograph goes to zero. The next event is considered as the next non-zero period. This discretization was necessary for computational efficiency considering the hydrographs have such a high temporal resolution. A total of 619 events between May 2017 and April 2019 were analyzed. The number of events was trimmed down to 100 to exclude events with peak discharge less than 0.0002 m³/s because they did not produce a shear stress that was geomorphically significant (these events are summarized in Table S.1 of the Supplemental Information). This threshold was established by estimating the critical shear stress for the median soil (particle d₅₀) and estimating the corresponding critical discharge. Assumptions made in establishing this threshold are later explained.

This high resolution of DEM allows an accurate estimation of the hydraulic parameters required for a reliable evaluation of morphology within the bioinfiltration system. The DEM was modified to include the location of the two outlet structures within the model domain. Computation in FastMECH was performed on a 0.3 m-by-0.4 m grid, which was determined to be a computationally efficient grid size for this study (smaller grids were not feasible due to instabilities in the model).

The initial volumetric water content of the soil (θ_i) and the saturated volumetric water content (θ_s) are required to determine ⁠. Values of θ_s are often assumed equal to the soil porosity, f. In this study, values of f and were initially estimated based on recommendations by Rawls et al. (1983) for loamy sand. Porosity values were confirmed through comparison with measurements from soil cores collected from the field. When data on the initial water content of the soil prior to each storm event were unavailable, θ_i was assumed equal to the field capacity (the of the soil after excess water has drained). Average assumed values for f, field capacity, and were 0.4, 0.14, and 6.13 cm, respectively.

Previous studies have shown that k_sat values play the most significant role in modeling infiltration rates of stormwater control measures using the Green-Ampt method, with less sensitivity to other input parameters such as ψ or Δθ (Weiss & Gulliver 2015). In addition, the presence of vegetation and activities of other organisms within the soil are known to impact in situ infiltration rates and may result in a higher performance of the system. Therefore, the field k_sat and influence of vegetation were investigated through field infiltration testing. Spot-infiltration tests were performed with an automated SATURO dual head infiltrometer (METER Group, Inc. USA). Infiltration tests were performed at multiple locations throughout the SMP along the basin centerline to obtain a representative geometric mean k_sat for the basin. Values of k_sat from field spot-infiltration tests can be highly spatially variable and typically are log-normally distributed, supporting the use of a geometric mean of the measured k_sat values as the overall representative k_sat for the basin (Ebrahimian et al. 2019). Press (2019) recommended that, for testing of bioinfiltration sites with infiltrometers such as the SATURO, spot-infiltration tests should be performed at a minimum of five locations to obtain a geometric mean k_sat representative of the basin. Thus, the geometric mean k_sat values in this study were based on the results of infiltration tests performed at at least five locations throughout the basin. Each of the k_sat values measured in the field was temperature corrected. The geometric mean k_sat for the field spot-infiltration measurements was 5 cm/h.

In addition, ponding level was continuously monitored in the field with a level sensor located near the B2 outlet structure (Campbell Scientific CS451). The sensor had a total error band of ±2 mm. Monitoring occurred from May 2017 to April 2019. For all storm events during that period that produced a response from the level sensor, a curve fitting function was applied to the 5-min ponding data to calculate an average recession rate during the recession limb when rainfall was not occurring. The average recession rate was 6 cm/h. However, the measured water levels near the gage at the B2 outlet structure did not necessarily represent the average water level throughout the SMP. This information was used to calibrate the hydraulics of the model.

To allow for comparison with the results of the field-measured k_sat from the infiltration testing, hydraulic conductivity tests were performed in the laboratory on 3.8-cm-long soil cores collected from the field. The laboratory k_sat was measured under constant-head conditions in general accordance with ASTM D2434. The average hydraulic gradient during testing was 5. For all of the laboratory tests, permeation continued until at least four consecutive values of the ratio of outflow to inflow were between 0.75 and 1.25, and at least four consecutive k_sat values were within ±25% of the mean k_sat value of the consecutive measurements. The mean k_sat based on the laboratory testing was 3 cm/h, which was lower than the field k_sat of 5 cm/h. However, the lower k_sat value observed for the laboratory testing was expected, as laboratory measurements on smaller samples do not account for the presence of roots and other macropore effects present in the field. Thus, as the field infiltration measurements were considered more representative of actual infiltration conditions, a k_sat value of 5 cm/h was used in the analysis.

For each storm, the result of the Green-Ampt infiltration rate was converted to a volume rate by multiplying with the surface area of the study area (240 m²). Using the estimated volumetric infiltration rate, the storm events were filtered further to only include the significant storm events (events greater than the estimated minimum volumetric infiltration rate, 0.003 m³/s). 66 such significant storm events were identified in this study. The excess runoff was computed as the difference between infiltration and inflow hydrograph, which was then used for the geomorphic simulation (shown in Figure 3).

A boundary condition was defined for the upstream and downstream ends of the modeling domain. The downstream boundary condition was set using ponding depth measurements from the B1 gage. The upstream boundary condition was specified as a hydrograph. As FaSTMECH only allows for a single upstream boundary, the two upstream flows were combined into one total inflow hydrograph.

To quantify the accuracy of the prediction framework, two statistics were computed for each comparison of observed and predicted ponding depth: the Nash–Sutcliffe efficiency (NSE) coefficient and the correlation coefficient (CC) (Supplementary Material, Table S1). The NSE (Nash & Sutcliffe 1970) is commonly used for assessing the predictive power of hydrologic models. The NSE has a range of −∞ to 1 with a value of 1 representing a perfect match. A negative NSE value indicates that the observed mean is a better predictor than the model. The CC measures the linear correlation between the observed and predicted results.

Incorporating the initial moisture content, the framework also accounts for the antecedent soil moisture content within the soil leading to a more accurate estimate of the soil infiltration capacity prior to a storm event. The framework is also adaptable as it allows for soil properties to be modified depending on the soil type. This flexibility allows for this modeling framework to be extended to systems with soil types other than loamy sand.

Quantifying evapotranspiration (ET) is important in assessing the volume reduction potential of a bioinfiltration system (Hess et al. 2015; Wadzuk et al. 2015). ET, however, occurs on a much larger timescale compared with infiltration and fluvial geomorphology of a bioinfiltration system. For example, in their study of two bioretention systems over a 2-year span of time, Wadzuk et al. (2015) measured an average ET of 6.1 mm/day. In comparison, the median in situ saturated hydraulic conductivity rate in this study over a 3-year time span was 8 cm/h (1,900 mm/day). For this reason, ET was ignored in the development of the current framework but should be explored in future work for different types of bioinfiltration installations.

As seen in Table 1 (or Supplemental Material Table S1), although there were a few instances where the model performed poorly (represented by a negative minimum NSE and low minimum absolute CC), the overall performance of the model on the 31 events was very good (median NSE = 0.93 and median CC = 0.97).

This high level of prediction accuracy can be attributed to assigning the correct boundary conditions in this study. Specifically, an accurate prediction of the system's ponding depth largely depends on setting the right upstream and downstream boundary conditions. In this study, using water surface elevation data from the B1 gage as the downstream boundary condition was critical for the accuracy of the model. This observation highlights the importance of having a good knowledge of the configuration of the outlet structure when using this framework. It was also notable that the prediction accuracy was better for events where the maximum ponding depth exceeds 0.2 m, such as events 282, 302, 311, and 327 in Figure 7. Saturated soil conditions associated with ponded water allow the system to function more like a river channel with minimal infiltration and with hydraulic conditions controlled by the downstream boundary condition. This was tested using depth data from the pressure transducer at gage B2, located between the inlet and B1 (Figure 2).

The results show that for a typical significant event (0.03 m³/s) in SMP A, there is not substantial erosion within the system. More likely, there will be deposition in a vast region within the system. The significance of these shear stress distribution results is the ability to identify regions with a high potential for erosion either during design or planning for maintenance of bioinfiltration systems. Sediment flux provides an insight into the rate of mobilization of sediment within the bioinfiltration system. Such knowledge is crucial in predicting sediment erosion and deposition, which can educate long-term maintenance and rehabilitation needs of a bioinfiltration system (DelGrosso et al. 2019; Smith et al. 2021; Wadzuk et al. 2021).

SMP A currently shows minor spot erosion and micro-channelization around vegetation but no major erosion within the system except for scouring downstream of the weirs due to relatively higher velocities. It is important to highlight that one limitation of the current framework is the inability to model localized erosion around individual clusters of vegetation and hydraulic structures such as the weirs. It is expected that the vortices around individual vegetation or clusters of vegetation will lead to increased scouring. However, modeling this scenario is outside the scope of this study and may be better suited for more advanced unsteady 3-D hydrodynamic or morphodynamic modeling.

Despite the high level of accuracy achieved comparing the observed hydraulic depth in the pond and predicted pond depth of the bioinfiltration system, this framework has some limitations. As highlighted earlier, the FaSTMECH sediment transport model utilized in this study allows for a single inflow boundary. This presents a challenge when evaluating systems with multiple inlets, like ours. Future studies on this framework must incorporate other transport models with more flexibility to accommodate additional inflow locations.

Another important factor to note is that the current framework assumes that the bioinfiltration media is homogenous. This may be true for engineered system such as a bioinfiltration system. However, caution should be exercised when the framework is being applied to non-homogenous systems or systems where the soil bed has been significantly altered by the introduction of different soil types. Use of the Green-Ampt method allows for input parameters that are easily obtainable from the literature. The limitations of this method, however, are that it may not fully address the spatial variability typical of a bioinfiltration system and only provides an approximate solution to Richard's Equation. The assumptions of an initially uniform moisture content profile and a sharp wetting front may introduce additional errors associated with the use of the Green-Ampt method. Another limitation is the fact that the results only represent the median particle sizes (d₅₀) and thus may not fully capture the dynamics of the entire range of particle sizes within the soil media.

As has been demonstrated in this study, the framework requires a downstream boundary condition to be defined in order for the simulation to proceed. In fact, the accuracy of the results obtained from this framework is largely a function of correctly defining the downstream boundary condition. Additionally, this geometry may not be representative of typical bioinfiltrations systems and as such may not be applicable in every case. Future studies may need to examine the applicability of this model to other systems. Understanding this limitation is important during data collection for using this framework.

This study introduces a novel framework for spatial prediction of sediment erosion and deposition within a bioswale. By identifying regions of potential high shear stress (erosion) and sediment deposition (sedimentation) maintenance requirements and long-term performance of bioinfiltration systems can be more focused. This has the potential to aid in dynamic maintenance (Wadzuk et al. 2021). This study highlights an important, yet often ignored aspect of low impact development design, which is understanding the evolution of the soil media. Substantial soil loss within soil media impedes the performance of the system by limiting the storage volume capacity of the system and may overwhelm downstream infrastructure with eroded soil material. A good understanding of the dynamics of the soil media morphology helps in predicting and prioritizing future maintenance needs as well as SMP design to minimize impact from excessive erosion or deposition. For large systems, understanding the shear stress distribution and rate of soil loss can have major cost saving implications by knowing which storm events can significantly impact the integrity of the infiltration bed. As nature-based solutions play an increasingly significant role in urban stormwater management, leveraging existing modeling tools for natural systems provides an opportunity for advancing the science and engineering in tangential fields. Enabling the ability to identify areas where deposition and erosion are likely to occur will allow for more efficient and sustainable designs, and aid in identifying the drivers of underperformance or non-performance of systems, since the depositional areas are likely to correspond with areas of reduced infiltration rates.

The current study has introduced a unique method of modeling the hydraulics and morphology of a bioinfiltration system with a high degree of accuracy. This was achieved by combining an existing model, FaSTMECH, with an infiltration model (Green-Ampt). This novel approach allows for an accurate accounting of the high infiltration rates typical of bioinfiltration systems while modeling sediment transport within the system. The proposed framework is unique as it repurposes an existing river model for a bioinfiltration system and provides an inexpensive method for modeling the hydraulics and morphology of systems with high amounts of infiltration using easily accessible open access tools.

In developing the current framework, the most important data requirements include a DEM, k_sat, d₅₀ of the soil media, inflow hydrograph, and the ponding depth. Of the hydrologic data required for utilizing this modeling framework, downstream ponding depth information appears to be the most critical input. This is because ponding is largely controlled by the configuration of the hydraulic structures. The accuracy of a model, therefore, is highly dependent on a good understanding of the ponding elevations.

While some limitations were identified with FaSTMECH, there are other solvers within the iRIC platform that may have the potential to represent systems more accurately with multiple inlets. This needs to be investigated in future studies. As these tools evolve it will enable testing a range of storm events and infrastructure geometries, providing insight into how sediments move through these systems, potentially impacting their function.

The authors have presented the first known attempt to model morphology and bioinfiltration soil loss in order to identify potential problematic areas within the system where significant erosion or channelization is likely to occur. The bed shear stress and sediment flux outputs of the proposed framework will be a crucial tool for low impact development (LID) engineers and planners in targeting and prioritizing specific regions within a system where failure or underperformance is likely to occur. This framework will also be an effective tool for evaluating design alternatives for minimizing soil loss within a bioinfiltration system.

The authors would like to thank the Pennsylvania Department of Transportation (PennDOT) for their support and funding. The opinions presented in this publication are those of the authors and do not necessarily express the opinions of the PennDOT. Reference in this report to any commercial product, process, or service, or the use of any trade, firm, or corporation name is for general informational purposes only and does not constitute an endorsement or certification of any kind by the authors. This project is a research initiative of the Villanova Center for Resilient Water Systems. The authors would also like to thank their colleagues on this project, faculty, graduate research assistants, and undergraduate research assistants, in this project for their invaluable insights during the preparation of the manuscript.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.