Numerical investigation of the performance of engineered barriers in reducing flood risk

In this paper, 2-dimensional, hydro-mechanically coupled finite element analyses are performed to assess the performance of an engineered barrier aimed at reducing flood risk in urban environments. The barrier consists of an unsaturated compacted soil layer and a drainage layer of a coarse granular material, constructed on top of the natural soil, in this case London clay. The barrier is vegetated so that its water storage capacity is renewed after each rainfall event. Sophisticated boundary conditions are used to simulate the effect of precipitation and evapotranspiration. The rainfall water infiltration and the initiation of water run-off during intense precipitation events are simulated. The effect of the hydraulic properties of the unsaturated soil layer on the performance of the system is investigated by means of parametric analyses. The effect of precipitation rate and geometry of the barrier is also discussed. Design recommendations for the properties of the compacted layer and the dimensions of the system are given at the end of the paper.


Introduction
Extreme precipitation events that lead to flooding have become more common in recent years. Brisbane Australia Climate projections for the rest of the century show continued intensification of daily precipitation extremes [1], while continuous urbanization and infrastructure development create more impervious surfaces in cities. The combination of intense rainfall events and limited access of the water to the ground increase the flooding risk significantly.
A potential solution to control the risk of surface run-off is the use of engineered barriers, with enhanced water holding capacity. These systems may also return the stored water back to the atmosphere, by evaporation and plant transpiration, during dry periods. To investigate the effectiveness of this solution and to optimise the design of barriers, sophisticated rainfall infiltration analyses must be performed.
A review of rainfall infiltration analysis was presented in [2]. The use of empirical and analytical models provides tools for better understanding of water infiltration in soil. However, due to the non-linear mechanical and hydraulic behaviour of soils during the process, advanced numerical methods are needed for accurate infiltration simulation.
Hydro-mechanically (HM) coupled finite element (FE) analysis, that considers the effect of rainfall infiltration, is used in slope stability geotechnical analysis (see for example [3]). However, very limited work has been published on the usage of HM-FEM for the simulation of surface run-off initiation. Although mechanical instability is not the predominant risk in surface flooding events, especially in urban environments, the accurate modelling of soil stress evolution during rainfall infiltration leads to successful simulation of the evolving hydraulic properties (e.g. variable permeability due to desiccation cracking).
This paper presents a HM-FE model that simulates the construction of an engineered barrier for flood protection, the infiltration of rainfall water into the barrier, potential evapotranspiration and water run-off initiation when intense rainfall events are applied. Based on these capabilities, different scenarios can be tested, to evaluate the performance of the barrier during rainfall events with different intensities, as well as testing barriers with different material properties and geometries.

Finite element analysis overview
The rainfall infiltration into soil is modelled by means of hydro-mechanically coupled (HM) finite element (FE) analysis. The Imperial College Finite Element Program (ICFEP; [4]) is utilised, where Biot's consolidation theory extended to unsaturated soils has been implemented. The governing equations and their implementation were presented in [5], with subsequent modification detailed in [6] and [7]. Atmospheric boundary conditions (BC) are used to simulate rainfall, evaporation, and transpiration in plants. An outline of ICFEP capabilities in terms of boundary conditions, mechanical, permeability and soil water retention models necessary for the analysis of this type of problems can also be found in [8].
A 2-D plane strain mesh is created, as a column of 92, 8-noded, quadrilateral elements (Fig.1c). Each element has 16 displacement and 8 pore water pressure (PWP) degrees of freedom. The mesh represents a 45metre deep soil column. The column in Fig. 1a represents the first 45 meters of the in-situ soil conditions, until the depth where the bedrock is assumed. It consists of 3 meters of Weathered London Clay (WLC) and 42 meters of London Clay (LC). This column is subjected to a long duration (4-7 years) of soil-atmosphere interaction analysis, to estimate the stable seasonal PWP profile (1 st stage). The atmospheric BC and material properties are explained subsequently.
At Stage 2, the top 1 m of the column is excavated. The 1 m thick engineered barrier is then constructed, consisting of a 0.2 m thick Drainage Layer (DL) at the base, and a 0.8 m thick Compacted Soil layer that reaches the original ground surface (CS; Fig.1b).
At Stage 3, one additional year of soil-atmosphere interaction is simulated with the soil column of Fig. 1b. Three rainfall events of different intensity are then applied: each is 10 hours long, starting from a dry summer day. The water run-off initiation and the rate of surface suction decrease is investigated, as well as the effect of rainfall intensity, hydraulic properties, and geometry of the geo-barrier.
In the FE model of Fig. 1c, the horizontal displacements along the two vertical boundaries are set to zero (ux=0), whereas the material can deform freely in the vertical direction. The two vertical boundaries are assumed impermeable (qf=0). The bottom boundary of the column is fully fixed in terms of displacement degrees of freedom (ux=uy=0), with no change in pore water pressure (Δpw=0). Atmospheric boundary conditions are applied at the top of the column, as explained below.

Fig. 3 Precipitation and potential evapotranspiration for London area.
Rainfall is simulated with the use of the precipitation BC in ICFEP, which is a dual boundary condition: a flow rate qn, equal to the rainfall intensity, and a pore water pressure condition, pfb, are prescribed in every increment for which this BC is active. At the beginning of such increments, the algorithm compares the pore water pressure at each of the boundary nodes to pfb; if found to be more tensile, then the flow rate of qn is applied (rainfall intensity); otherwise the constant pressure pfb is applied. For more details about the development and implementation of the precipitation boundary condition in ICFEP the reader is referred to [5] and [9].
Evaporation and transpiration are simulated in a simplified manner, as a combined evapotranspiration outflow source in the continuity equation. This is possible by applying the vegetation BC presented in [10]. The potential evapotranspiration rate, ܶ is an input for each increment of the analysis. Based on a root water uptake model (Fig. 2), the maximum rate of water outflow is first calculated as: where rmax is the maximum root depth (input); bellow this depth Smax=0. In the analysis presented here rmax=0.5 m. The actual outflow rate is calculated by multiplying Smax by a suction-dependent parameter α, according to [11]. The variation of α with suction is illustrated in Fig. 2. Suctions S1 (anaerobiosis point), S2, S3 and S4 (wilting point) are input parameters and were taken as 0, 5, 50 and 1500 kPa, respectively. The fact that Sact depends on suction via α makes the root water uptake model nonlinear, as was highlighted in [10]. The input precipitation and potential evapotranspiration rates (mm/day) for a typical year in London are presented in Fig. 3. These are 29-year average (1972-2000) monthly rates, measured at Greenwich weather station, London. The evapotranspiration rates are calculated based on the FAO Penman-Monteith method [12]. These values are assigned for each increment of Stage 1, to compute the stabilised seasonal PWP profile. Synthetic data for Stage 2 are constructed for daily averages of London meteorological data as explained in Section 3.4. Table 1 Constitutive parameters for Mohr-Coulomb model.

Parameter
Value Parameter Value Table 2 Constitutive parameters for non-linear elastic small strain overlay model [13].   [13], for all three materials. The material parameters for the latter model are presented in Table 2.
The calibration of the MC and small-strain elastic models for weathered and unweathered London clay was presented in [3] and is followed here. The same parameters are assumed for the DL for simplicity in this study.
The mechanical behaviour of the unsaturated compacted soil (CS) in the engineered barrier is simulated by a constitutive model detailed in [14]. This is an extended and modified version of the Barcelona Basic Model [15], first implemented in ICFEP by [16]. The mechanical material properties in the current paper are the same as those employed in [17] and are summarised in Table 3, corresponding to a silty soil tested by [18]. The over-consolidation ratio of the compacted soil is assumed at 1.5 in this analysis.
The hydraulic conductivity properties for WLC, LC and DL are summarized in Table 4. The hydraulic parameters of London Clay are well described in the literature and have been employed in various studies [19]. A variable permeability model for LC is employed here, which is a function of mean effective stress [4]: where α is a fitting parameter. As far as the WLC is concerned, the variation of permeability due to desiccation crack opening (dry months) and closing (wet months) is simulated as: where σΤ is the current tensile total principal stress, and σT1 and σT2 are two tensile total stress limits, in between which the logk varies linearly. The permeability of the drainage layer is assumed constant and significantly higher than that of London Clay ( Table 4). The logarithm of the permeability in the compacted soil varies linearly with suction; an initial saturated value, ksat, corresponding to suction SA is assumed and a final value kmin, is computed when suction equals SB (Table 5). A non-hysteretic soil-water retention curve (SWRC) is employed, described as: where seq=s-sair; sair is the suction value at air-entry value; s0 the suction value at zero degree of saturation; ad,w is a fitting parameter for the primary drying (d) and wetting (w) curves. The parameters are summarized in Table 5.
3 Numerical Analysis

Initialisation of stress field
At the beginning of the analysis, soil stresses are initialised based on the unit weight of London Clay (LC), i.e. 19.1 kN/m 3 , above and below the ground water table (GWT). The coefficient of earth pressure at rest, K0, is 2.1 at the ground surface, reducing to 0.6 at 15m below the ground surface. The GWT is 1 m deep and the pore water pressure profile is hydrostatic, with suction developing above the GWT. In Stage 1, a long duration analysis was performed to compute the stabilised seasonal pore water pressures (PWP) in the soil profile, with the assigned mechanical and hydraulic properties discussed in Section 2.3. Fig. 4 presents the simulated PWP, for the first 20 meters of the 45-meter deep soil column. It is concluded that 4 years of analysis is sufficient to achieve the stabilised PWP profile for this set of soil properties, as the filled circular points (4 years) are shown to coincide with the solid lines (6 years). This result is consistent for both September (dry month) and March (wet month). It is of interest to note a different trend in simulated PW pressures in the WLC compared to the LC. This is attributed to the differences in the reference permeability, but also to the different evolution of permeability in the two layers. The reference permeability of the WLC is 1.55x10 -4 m/hour (4.3x10-8 m/s), one order of magnitude higher than the reference permeability of the LC (1.33x10 -5 m/hour). Fig. 5 shows the evolution of WLC and LC permeabilities, during the first year of the long-term analysis (Stage 1), at 1.5 and 3.5 meters, respectively.

Seasonal pore water pressure
The mean effective stress dependent permeability of LC (Eq. 2) reaches a peak of 1x10 -5 m/hour during the first year of analysis, whereas the desaturation dependent permeability of WLC reaches a peak of 5.5x10 -3 m/hour (Fig. 5). The increasing suction leads to an increase of the tensile total stress, causing cracks to open when the limit value of the tensile stress is mobilised and in turn increasing the overall permeability of this layer. The

Construction of the engineered soil barrier
At the 2 nd Stage, the engineered soil barrier is constructed. The process is simulated by applying the excavation and construction boundary condition in the FE analysis [4]. The elements of the top 1 m of the soil column in Fig. 1a are removed and the same elements are constructed again, but with DL and CS properties and become again a part of the active mesh [4]. The process occurs in 4 time steps, of 't=1 day each. The initial suction in the CS upon construction is an input in the simulation and 50 kPa was assumed. A degree of saturation consistent with the SWRC was also input.

Daily rainfall and evapotranspiration rates
For a more accurate estimation of the initial PWP profile, before applying the intense rainfall events, one additional year of soil-atmosphere interaction is simulated. The rates during this year are applied in daily averages. Each increment of the analysis represents 1 day, and thus 365 increments are applied to the FE model.
Daily average meteorological data, consistent with the monthly average meteorological data for London (Fig.  3), were not available to the authors at present. Since the main goal is to develop a conceptual numerical procedure, representative synthetic daily data was derived. This data is generated by the monthly average dataset, by assuming a non-uniform distribution of the rainfall within a month. The distribution is defined by analysing daily data from another area of the UK, Newcastle, based on the dataset presented in [20]. It should be noted that this process may not reflect with absolute accuracy the distribution of meteorological phenomena in London, within a month; however, the monthly average precipitation and potential evapotranspiration rates that are applied match exactly those of Fig. 3. Daily meteorological data for the London area will be included in this analysis in future work.   6 shows the evolution of PWP at the surface of the engineered barrier over the period of one year, starting on April 1 st . After the first 5 days, when the construction of the soil barrier is taking place, the initial suction is set to 50 kPa, as explained earlier. Increasing surface suction is computed during the first six dry months (approximately 180 days to end of September), and then it decreases during the precipitation dominant wet months. The intense rainfall events are applied during a dry month, as described below.

Simulation of intense rainfall events
Three precipitation events of different intensity are applied in Stage 4, starting from mid-July, thus from a dry month with initial surface suction (positive PWP) of 120 kPa. The duration of each precipitation event is 10 hours. The rate of surface suction decrease is a function of the hydraulic and mechanical parameters of the simulated layers. The point where surface suction becomes zero is the run-off initiation point in time, thus, a potential initiation of a flooding event. Although in practice a flooding event depends on the geometry of the ground surface and the available volume of surface run-off water, the initiation point as a conservative method to estimate potential flooding.
In Fig. 7, the simulations of suction decrease during the qn=5, 25 and 60 mm/hour rainfall events are presented. These three rainfall magnitudes are thought to correspond to moderate, heavy, and violent events respectively [21]. It is concluded that the numerical model can simulate the effect of rainfall intensity for the intended run-off prediction. A moderate event of 5 mm/hour does not lead to run-off initiation, thus, there is no risk of run-off flooding during a 10-hour period. However, the heavy and violent intensity events lead to run-off.  In Fig. 8, the heavy intensity event is applied to the engineered barrier, in which three different reference permeability values are assigned to the Compacted Soil (CS). It is concluded that a permeability of magnitude 3E-07 m/s significantly impacts the infiltration process. The ability of water to infiltrate is significantly decreased, when compared with the k=3E-06 m/s CS. However, the difference in water infiltration is not as large, when CS with k=3E-06 m/s and k=3E-05 m/s is compared. This leads to the conclusion that there is an optimum permeability range, above which increasing the permeability does not effectively improve the performance of the barrier.

Parametric analyses
Finally, two barriers, with compacted soil (CS) layers of different depth, are tested and compared. The first one is the 0.8-metre-deep CS layer of Fig. 1b. The new one is a 0.5-metre-deep CS layer, with the results presented in Fig. 9. It can be observed that the initial surface suction is higher in the 0.5-meter-deep layer case, even though the Stage 4 of analysis starts from the same day in July. This is due to the difference in the stabilised seasonal pore water pressure (PWP) profiles for the two cases, that were computed in Stage 3 of the analysis. However, the rate of surface suction decrease is higher for the shallower CS layer, since it is primarily the volume of the CS material that accommodates the volume of the infiltrating rainfall water.

Conclusions
In this paper, a 2-dimensional, hydro-mechanically coupled finite element (FE) model was built to perform rainfall water infiltration analyses, in in-situ soil and in a constructed engineered barrier. Τhe purpose of the constructed barrier is to minimise the surface run-off risk. The excavation and construction of the barrier, as well as its performance during intense precipitation events were simulated. The FE model accounts for the non-linear mechanical and hydraulic behaviour of the simulated soils.
The numerical model demonstrated numerical stability in all cases, even when applying high precipitation rates compared to the CS layer's permeability. The run-off point in time is computed, as well as the rate of surface suction decrease, that is a function of rainfall intensity and soil properties. It is demonstrated that the rainfall intensity has a significant effect at the point of surface run-off initiation in time. A violent rainfall event (60 mm/hour) led to run-off initiation in less than four hours, whereas a moderate rainfall event (5 mm/hour) did not lead to surface run-off initiation within a 10 hour event. Those events were applied during a dry month (July) with initial surface suction measuring 120 kPa.
It is also demonstrated that increasing the permeability of the upper soil layer of the barrier has significant effect in minimising the risk of run-off, until an optimum permeability is reached. Increasing the permeability above optimum has only a marginal effect on the barrier performance.
Finally, it is concluded that a reduced thickness of the engineered barrier will increase the initial surface suction compared to a thicker layer during dry summer months. However, the rate of surface suction decrease will MATEC Web of Conferences 337, 04003 (2021) PanAm-UNSAT 2021 https://doi.org/10.1051/matecconf/202133704003 be significantly higher, and thus run-off is expected earlier than in the thicker layer.