Mapping of uplift hazard due to rising groundwater level during floods

SPECIAL ISSUE Mapping of uplift hazard due to rising groundwater level during floods TomášJulínek | David Duchan | Jaromír Říha Faculty of Civil Engineering, Institute of Water Structures, Brno University of Technology, Brno, Czech Republic Correspondence TomášJulínek, Faculty of Civil Engineering, Institute of Water Structures, Brno University of Technology, Veveri 331/95, 602 00 Brno, Czech Republic. Email: [email protected] Funding information FAST-S19-5714 Probabilistic Assessment of Soil Instability Due to Seepage in Earth Structures and Their Foundations; TH04030087 Tools for Optimisation of the Management of Levee Systems Abstract European Directive 2007/60/EC only briefly mentions the problem of hazard arising due to groundwater flooding, and techniques for the mapping of hazard occuring due to rising groundwater have not yet been scientifically developed. The groundwater-related threats that occur during floods may include concentrated leakage of groundwater behind levees, heave, or potential uplift of the topsoil layer at the protected area. The hazard corresponding to rising groundwater level depends on a number of factors related to the flood course, groundwater regime, geology, and topology of the protected area. The limit state approach is applied to the assessment and mapping of hazard induced by rising groundwater level in the area behind flood protection barriers, and the contributing factors are discussed, quantified, and incorporated into the limit state condition for topsoil layer uplift (UPL). An over-design factor is expressed as a function of the spatial co-ordinates (x,y). Data from geological and hydrogeological surveys and groundwater flow modelling are used to evaluate individual terms in the limit state condition. Uncertainties in the input data are expressed via partial factors. Data collection and their geographic information systems analysis completed with hydraulic modelling are crucial techniques in the hazard mapping of potential UPL during floods. The article includes a case study featuring a flood protection scheme for a shopping centre in the city of Brno, Czech Republic. KEYWORDS European Directive 2007/60/EC, flooding, groundwater, hazard mapping, limit state, uplift 1|INTRODUCTION Floods are characterised by a growth in river flow rates, increased water stages in streams, overbanking, and the inundation of floodplains. Increased water stages during the flood event may significantly affect the groundwater flow regime in the aquifer adjacent to the river. In the case of the construction of flood control works, their subsurface parts can affect the natural groundwater regime during periods when flooding is absent, at which times groundwater usually flows from higher-lying areas towards streams which drain adjacent aquifers (Fetter, 2001). Therefore, to maintain free communication between the river and aquifer (Directive, 2006), for example, seepage barriers such as cut-off walls or injection walls should be designed to be partially penetrable. Received: 12 June 2019 Revised: 4 January 2020 Accepted: 27 January 2020 DOI: 10.1111/jfr3.12601 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. Journal of Flood Risk Management published by Chartered Institution of Water and Environmental Management and John Wiley & Sons Ltd. J Flood Risk Management. 2020;e12601. wileyonlinelibrary.com/journal/jfr3 1of13 https://doi.org/10.1111/jfr3.12601 As a consequence of the rise in the water level in a stream during a flood, surface water temporarily infiltrates the banks of the river and propagates into the aquifer in the opposite direction to that of normal flow. The groundwater table (or piezometric head) rises and may cause one of the reported types of groundwater flood events (CIRIA, 2013; Conant, Robinson, Hinton, & Russell, 2019; Fleckenstein, Krause, Hannah, & Boano, 2010; Robins & Finch, 2012; Sophocleous, 2002): •a true “groundwater flood”in which the water table elevation rises above the ground elevation, •a“groundwater induced flood,”which occurs when intense groundwater discharge via bourne springs and highly permeable shallow horizons discharges to surface waters and causes overbank flooding, •inundation of subsurface infrastructure. If impermeable topsoil covers the aquifer, excess uplift may act on the topsoil base and cause its rupture, followed by uncontrolled concentrated leakage at the protected area behind flood protection measures (FPM). All of the mentioned cases may be related to the local flooding of alluvial lowlands behind flood levees or floodwalls, which is caused when the groundwater table or piezometric head rises above the terrain. This phenomenon is essentially controlled by the local river stage (MacDonald, Bloomfield, et al., 2008; Robins & Finch, 2012), but may also be affected by groundwater inflows from higher geological formations like alluvial terraces. Concentrated or diffuse leakage manifests itself via the ponding of groundwater on the surface (JACOBS, 2006). Various aspects of groundwater flooding in the context of morphology, the geological composition of the subbase, and so on have been a subject of great concern to authors such as Adams et al. (2010), MacDonald, Hughes, et al. (2008), and others. A groundwater flood occurring alongside streams during a flood affects (in relationship to the duration of the flood) the area behind the flood protection line by increasing the stages of the groundwater table. Moreover, uncontrolled seepage may result in a hazard both to flood protection elements and to the subsurface parts of structures like cellars, subsurface garages, and so on in the territory protected against flooding (MacDonald, Dixon, Newell, & Hallaways, 2012). The threat is posed by the pressure of subsurface water on the inundation and subsurface infrastructure (Abboud, Ryan, & Osborn, 2017; Kreibich & Thieken, 2008), namely underground parts of civil structures, and by the high hydraulic gradients along the foundations of flood protection elements and at places of groundwater leakage onto the ground surface. Seepage into the protected area is also unfavourable during a flood as it increases the amount of “internal waters”that must be pumped back into the stream. However, up to now, losses caused by rising groundwater levels have been neglected in flood risk studies (Kreibich & Thieken, 2008). According to Directive (2007), the scope of flood hazard maps related to groundwater flooding shall be limited to extreme event scenarios. However, such maps are not routinely prepared as a part of flood risk management plans. Mathematical modelling has been employed in the analysis of the seepage flow regime for decades. Current hydraulic modelling methods (Bear & Verruijt, 1992) and the existing software enable the effective evaluation of groundwater regime changes due to changes in hydrological conditions and boundary conditions. The methods provide effective assessment of the impact of groundwater on foundation soils, the foundations of flood protection structures, buildings in protected areas and the seepage of groundwater onto terrain (Mansour et al., 2013; Šoltész & Baroková, 2014). Hydraulic models linked to computeraided design (CAD) or geographic information systems (GIS) systems provide spatial data for the evaluation and quantification of groundwater-related hazard and its representation using computer mapping techniques (Büchele et al., 2006; Hughes et al., 2001; MacDonald, Bloomfield, et al., 2008; Morris, Cobby, Zaidman, & Fisher, 2015; Naughton, Johnston, McCormack, & Gill, 2015; Sommer, 2007). Such techniques are routinely applied in the development of flood risk management plans which take into account surface flooding (Arrault et al., 2016; CDS, 2013; Directive, 2007). However, there is a lack of studies focused on the mapping of the hazard induced by groundwater, namely due to the uplift on the topsoil layers during flood events. The evaluation of the uplift (UPL) limit state is described in (Bond & Harris, 2008; Eurocode 7, 2004; Frank et al., 2004). In the case of UPL at an area protected against flooding, all variables entering the limit state condition have a spatial character. The over-design factor (Frank et al., 2004) was introduced as a function of spatial coordinates for the assessment of hazard. In this study, all of the concerned variables (including the overdesign factor [OF]) are expressed using thematic maps processed by GIS techniques. 2|MATERIALS AND METHODS The method of evaluating groundwater hazard due to uplift during floods consists of the following steps: •identification and conceptual description of hazard arising from UPL scenarios, 2of13 JULÍNEK ET AL. •limit state definition, •input data collection, verification, and analysis, •hydraulic modelling, •quantification of hazard due to UPL via flood mapping. 2.1 |Identification of hazard During a flood, the water level in the stream increases as the piezometric head of the groundwater in the aquifer rises. In Figure 1A,B, a scheme of groundwater flow is shown for two scenarios. Figure 1A shows the state just before a flood, when the water level in a river begins to rise and starts to infiltrate into the adjacent aquifer. In Figure 1B, a scheme of the flood peak is shown. The piezometric head rises significantly above the terrain (Jandora & Říha, 2008; CIRIA, 2013). Figure 1A,B illustrate a situation where a protected area is geologically formed by an upper topsoil layer composed of low permeable materials (alluvial loams) overlaying a permeable (e.g., gravel-sand) aquifer. The area behind the flood protection barrier is thus endangered by the increase in the piezometric head above the terrain and by the uplift of the less permeable topsoil layer. 2.2 |The application of the limit state procedure In general, the internal erosion process can be divided into four stages. These are initiation, continuation, progression, and failure (Fell & Fry, 2007). In standard studies on the safety of levees, the initiation phase should be avoided and is subject to the (recently recommended) performance of an assessment using the limit state approach according to Eurocode 7, 2004. In our case, the initiation phase is represented by the rupture of the topsoil layer at the “dry”side of FPM due to upward water pressure in the aquifer. This type of internal erosion failure is classified as uplift (Bond & Harris, 2008; Eurocode 7, 2004). The general form of the limit state condition holds: Fdst,d≤Fstb,d,ð1Þ where F dst,d is the design value of the destabilising force and F stb,d is the design value of the stabilising force acting on the topsoil. In the safety assessment, it is necessary to deal with uncertainties in the input data which in the limit state method are incorporated via partial factors (sometimes referred as reliability coefficients; Eurocode 7, 2004). The design values of the forces are obtained by multiplying their characteristic values by the corresponding partial factors. For UPL assessment, the design values in Equation (1) acting on the topsoil may be specified as follows: Fdst,d=Fupl,d=γuplFupl,k,ð2Þ Fstb,d=Fg,d=γgFg,k,ð3Þ where F upl,d is the design value of the destabilising force due to the vertical water pressure (uplift), γ upl is the partial factor of variable load from the pore water pressure in the aquifer, FQ upl,kis the corresponding characteristic value of the uplift force, F g,d is the design value of the force from the total weight of the topsoil layer, γ g is the partial factor expressing uncertainties in the topsoil layer weight, and FG g,kis the characteristic value of the force from the weight of the topsoil layer. In Equations (2–6) all forces are expressed per unit area. When introducing the importance factor γ 1 (expressing the importance of the structure and the potential hazard to which it is subject) and substituting the design values from Equations (2) and (3) to Equation (1), one obtains: (A) (B) FIGURE 1 The topsoil uplift process, (A) just before a floodand (B) during a flood JULÍNEK ET AL.3of13 γ1γuplFupl,k≤γgFg,k:ð4Þ Using the notation in Figure 2, the characteristic values of the forces in Equation (4) may be expressed as follows: Fupl,k=γwhmax −LBS ðÞ,ð5Þ Fg,k=γsLT−LBS ðÞ=γsb,ð6Þ where γ w is the specific weight of the water and h max is the maximum piezometric head reached during the flood period (Section 2.4), L BS is the level of the topsoil base, γ s is the total specific weight of soil in natural conditions, L T is the level of the terrain, and bis the topsoil layer thickness. Variables L T ,L BS ,L BA ,b, and h max are functions of spatial coordinates (x, y). When implementing Equations (5) and (6) into the limit state relation (4) for the point (x, y) one obtains: γ1γuplγwhmax x,yðÞ−LBS x,yðÞðÞ≤γgγsbx,yðÞ ð7Þ Condition (7) is used for the limit state assessment performed for the area behind the FPM (see Section 2.5). 2.3 |Data The input data as well as hazard quantifiers vary in space (x, y), while the piezometric head in the aquifer during the flood varies with time (t). In order to evaluate groundwater flood hazard arising due to uplift, the set of input data has to be aquired. The necessary geological and hydrogeological data are related to the topsoil layer, hydraulic conductivity, and the thickness of the aquifer (Figure 2). A typical geological composition alongside streams consists of a relatively impermeable aquifer base (e.g., Neogene sublayers) covered by Quaternary aquifer layers. Poorly permeable topsoil layers often overlay aquifer soils. The basic information about the geological composition of the area is taken from geological and hydrogeological maps, as well as from the central database of geological works (CGS, 2019a), which is a cheap source providing all registered geological surveys in the Czech Republic. A more detailed onsite geological survey is usually necessary for the design of FPM as well as for modelling purposes (Section 2.4). The data are usually obtained in the form of a point dataset with irregular placement. The UPL assessment employs the following data: •layer thickness is identified in the boreholes and is specified by the level of the base of the topsoil L BS and the level of the base of the aquifer L BA , •geotechnical properties of soils such as specific weight, porosity, grain size distribution, and so on are obtained from geotechnical analysis carried out in the laboratory on the samples taken from boreholes, •the hydrogeological properties of the soil, such as the hydraulic conductivity kand storage Sof aquifer soils, are obtained from pumping or infiltration tests, •the groundwater levels for a given time are taken from monitoring boreholes; these data should be completed by the corresponding water stage in adjacent watercourses. Geographic data include: •basic raster maps (including digital terrain model [DTM]) and aerial photos, •geodetic survey data in the form of point datasets or vector maps specifying the spatial characteristics of the studied area and related objects. Hydrological and hydraulic data for surface streams are used for the specification of boundary and initial FIGURE 2 Scheme for conditions (4) and (7) 4of13 JULÍNEK ET AL. conditions for seepage modelling. Data representing the relation between surface water in streams and the groundwater level in hydrological monitoring boreholes are particularly important. These data, which are provided by the Czech Hydrometeorological Institute, are usually available from the state groundwater monitoring network (CHMI, 2019). The description of existing and newly designed FPM, which may significantly influence the groundwater regime, is taken from the design of individual structures (levees, floodwalls) and their components (slurry walls, cut-off trenches, drainage, etc.). All spatial data must be georeferenced. The application of CAD systems and/or GIS is a common procedure during the analysis, preparation, and interpretation of input data and results. 2.4 |Groundwater flow modelling Mathematical modelling is a common tool for predictions regarding a given groundwater flow regime and the changes that may affect it due to hydrological and hydrogeological conditions, flood course, human actions, and other factors. For spatial UPL assessment, the piezometric head h(x, y, t) in the area behind FPM has to be determined using the appropriate groundwater flow model. Some simplifications of reality due to the complexity of hydrogeological, hydrological, and topological conditions have to be accepted (Bear & Verruijt, 1992). Dimensional simplifications are applied to the model with respect to the hydrogeological conditions and groundwater flow regime (direction of seepage, flow regime, etc.). For the determination of piezometric head within large aquifers, a general three-dimensional seepage flow problem may be reduced to a two-dimensional (2D) unsteady flow model in a horizontal plane. In this case, the results are obtained in terms of the piezometric head h(x,y,t), which changes with location and with time. For the assessment, the most unfavourable situation is described by its maximum value h max (x,y) reached during the flood. For 2D, unsteady horizontal flow in confined aquifer the governing equation holds (Bear & Verruijt, 1992): ∂ ∂xT∂h ∂x  +∂ ∂yT∂h ∂y  −S∂h ∂t=0, ð8Þ where Tis aquifer transmissivity (T=k.b A ), kis hydraulic conductivity, b A is aquifer thickness, Sis the storage coefficient, and his the piezometric head. For the boundary it holds that: htðÞ= htðÞ,ð9Þ where  htðÞ is the known water level in the adjacent stream during the flood. At a boundary with prescribed flux it holds that: k∂h ∂xnx+k∂h ∂yny=qð10Þ where n x ,n y are directional cosines related to the outer normal vector to the boundary with prescribed flux q(per unit width of the boundary). Equation (10) may be applied at the Neumann “no flow”boundary, where q=0. The initial condition expresses the known piezometric head h 0 over the flow domain at the beginning of the flood (t= 0). This may be taken from the calibrated steady state solution: hx,y,0ðÞ=h0x,yðÞ,ð11Þ The maximum piezometric head h max at the given point of the domain is taken from the resulting piezometric head time course at the given location: hmax x,yðÞ= max hx,y,tðÞ fg ,ð12Þ The model domain usually follows local streams with known water level or piezometric head determined by groundwater level monitoring (hydrogeologic studies). To determine the boundary condition in Equation (9), the time course of the discharge and corresponding water stage in rivers has to be modelled. The flood hydrograph may be theoretical or taken from a real flood event. The numerical solution was done using the finite element method implemented in the code HPV2D (2015) developed at the Institute of Water Structures, Faculty of Civil Engineering, Brno University of Technology. 2.5 |Quantification of hazard The uplift stability assessment for the protected area is performed using Equation (7) with h max determined by Equation (12). For further analysis, the OF is introduced (Frank et al., 2004) as the ratio of the right and left sides of the Equation (7): OF x,yðÞ=γgγsbx,yðÞ γ1γuplγwhmax x,yðÞ−LBS x,yðÞðÞ ð13Þ JULÍNEK ET AL.5of13 Safety against uplift is achieved if OF ≥1. In the area where OF < 1 measures have to be adopted, for example, raising the terrain level, or using relief wells. The values of partial factors may be derived using various methods. General guidance is provided by Eurocode 7 (2004), where partial factors are assigned according to the characteristic design situation (see Eurocode 7, 2004, annex A. 3.1): •Permanent favourable action γ g = 0.9 •Variable unfavourable action γ upl = 1.5 In practical applications, the values of partial factors are frequently determined using expert opinion based on previous experience. In the case of a more extensive survey, where a sufficient amount of data is available, probabilistic methods may be used. 3|HAZARD MAPPING Input data and the results of modelling and assessment can be efficiently analysed, manipulated, and displayed using GIS mapping techniques. As a result, the following maps are gradually generated for UPL assessment in the flood protection context: 1. maps related to the input data: •map of the terrain level (L T ) compiled using a DTM and geodetic survey data, FIGURE 3 Flowchart for spatial uplift assessment via mapping 6of13 JULÍNEK ET AL. •map of the topsoil layer base (L BS ), •aquifer base (L BA ), •aquifer thickness b A =L BS -L BA , •monitored groundwater levels (h) for given measured water stages in surface streams, 2. maps related to destabilising action containing hydraulic modelling results (Section 2.4) showing: •the maximum piezometric head h max during flooding over all modelled scenarios (Equation (12)), •map of h max −L BS differences denoting the uplift pressure head at the base of the topsoil, 3. map related to stabilising action displaying the topsoil layer thickness b = LT −LBS, 4. hazard map displaying the over-design factor OF. The flowchart in Figure 3 summarises the analysis procedures. It includes the processing of input data from surveys, monitoring, and hydraulic modelling. Based on discrete spatially oriented data for individual parameters, the raster layers are processed using standard GIS tools for interpolation, calculation, and logical operations. In Figure 3, the procedure is divided into two parts for both the stabilising and destabilising action calculations. The destabilising action is interpreted by the thematic map of maximum piezometric head h max and the map of pressure head h max −L BS . The stabilising action mapping involves the analysis of land surface and geological survey data and provides a map of parameter b(x,y) determined from the map of the terrain and topsoil base level (L T −L BS ). 4|CASE STUDY 4.1 |Description of the study area The procedures described above shall now be demonstrated in connection with the flood protection scheme for the Olympia shopping centre in the city of Brno. The area of interest is located in the left-bank floodplain between the Svratka River and the north-west highway between Brno and Bratislava. Along its western side, the area is protected by flood levees, which are substituted in places by floodwalls. From the eastern floodplain related to the Svitava River and the Ivanovicky stream, the area is divided by the highway (Figure 4). The geological composition generally corresponds to the schemes in Figures 1, 2, and 5. The impervious base FIGURE 4 Area of interest, groundwater flow domain, and geological profile JULÍNEK ET AL.7of13 is about 6.4 –8.4 m below the terrain (aquifer base L BA ) and is composed of Neogene clays. The base is overlaid by permeable Quaternary fluvial gravels with an aquifer thickness of 3.6 –4.7 m, with hydraulic conductivity ranging from k= 3.010 −4 –4.210 −4 m/s and with a storage coefficient for the confined aquifer of 2.10 −4 [−]. The aquifer is covered by the topsoil layer, which is composed of relatively impermeable fluvial silty clays of irregular thickness in the range b=3–5 m. The composition can be seen in Figure 5. 4.2 |Data analysis, groundwater flow modelling, and mapping The first phase of the analysis includes data collection and analysis. General data on geology and hydrogeology were obtained from the map server operated by the Czech Geology Service (CGS, 2019b). Geological data were compiled from about 40 historical boreholes included in the CGS (2019a) database and from 10 additional boreholes drilled down to the aquifer base in the protected area. The corresponding levels and soil FIGURE 5 Geological profile FIGURE 6 Map showing digital terrain model and land survey points (L T ) FIGURE 7 Map of topsoil layer base L BS (CGS, 2019a) 8of13 JULÍNEK ET AL. characteristics were identified during the survey and subsequent geotechnical laboratory testing. Aquifer parameters such as hydraulic conductivity and storage were determined by constant-rate pumping tests carried out at two boreholes during the geological survey. The hydrograph related to a 50-year flood was provided by the Czech Hydrometeorological Institute (CHMI, 2019). It was transformed into a time-series of water stages in the Svratka River in order to define the boundary condition for the model (Equation (9)). The mapping of the ground level L T (x,y), topsoil base L BS (x,y), and aquifer base L BA (x,y) involved the spatial analysis of the input data from the survey. A DTM of the studied area was compiled using a combination of laser scanning data available for the whole territory of the Czech Republic and a specially conducted geodetic land survey, which covered the area in more detail (Figure 6). The map of the topsoil layer base (Figure 7) and aquifer base was interpreted from the historical and new boreholes, which were reasonably well-scattered over the area. For further processing, a map of the topsoil layer thickness b(x,y) was generated using two previously mentioned maps b(x,y)=L T (x,y)−L BS (x,y; see Figure 8). The accuracy of the spatial distribution of individual characteristics (L T ,L BS ,L BA ) depends on the accuracy, amount, and location of input data. The ground level obtained from laser scanning has an expected accuracy of ± 0.2 m, and when supplemented by a surface geodetic land survey the error drops to the magnitude of single centimetres. In the case of L BS and L BA , the error of the reading during the borehole drilling may be less than 0.1 m. More significant error may arise during the interpolation process in the case of scarce trial boreholes. In our case, due to the reasonably dense network of boreholes the expected error is in the order of single decimetres (about 0.2 m). The generation of maps depicting L BS and L BA involved data filtering and interpolation techniques. As mentioned above, the accuracy of the thematic maps strongly influences the interpolation method. Since the borehole data cover the area of interest reasonably well, the “natural neighbour”interpolation method (Bobach & Umlauf, 2007) was applied (Figure 7). Its advantage (Dumitru, Plopeanu, & Badea, 2013; Ledoux & Gold, 2005) is that it uses only a subset of samples that surround a query point. The results of interpolation do not produce anomalies such as peaks, ridges, and so on unless they are represented by the original input data. The degree of uncertainty involved in the interpretation of all layers should be taken into account and FIGURE 8 Map of topsoil layer thickness b(CGS, 2019a) FIGURE 9 Example of calibration results for a no-flood period JULÍNEK ET AL.9of13