Improving Groundwater Flow Representation in a Hydrological Model for Châteauguay Watershed

1. Improving groundwater flow representation in a semi-distributed hydrological model using hillslope Boussinesq and regional groundwater flow equations, with applications to the Châteauguay Watershed Ph.D. candidate: Stefan Broda, Centre de recherche pour l’Étude et la Simulation du Climat à l’Échelle Régionale, Département des Sciences de la Terre et de l’Atmosphère Université du Québec à Montréal Director: Marie Larocque, Centre de recherche pour l’Étude et la Simulation du Climat à l’Échelle Régionale, Département des Sciences de la Terre et de l’Atmosphère Université du Québec à Montréal Co-Director: Claudio Paniconi, Institut national de la recherche scientifique Centre Eau, Terre & Environnement, Québec

2. Improving groundwater flow representation in a semi-distributed hydrological model using hillslope Boussinesq and regional groundwater flow equations, with application to the Châteauguay watershed 1 Context This thesis is part of a larger project entitled “Improved models for surface water – groundwater interactions at the watershed and subcatchment scales: Châteauguay river basin, southwestern Québec”. The general goal of the project is to improve the understanding of hydrological processes at the watershed and subwatershed scales to address water resource management issues within the broader framework of surface water – groundwater interactions in a climate change context. The thesis concerns simulation of these interactions at the river basin (regional watershed) scale based on the HYDROTEL model (Fortin et al., 2001). In testing the coupled model to be developed, simulations at smaller scales (hillslope and subcatchment) will also be involved. 2 Problem Many recent studies indicate that coupling surface and subsurface flow is necessary to improve quantitative and qualitative watershed modeling. Groundwater dynamics have a large influence on accurate evaluation of streamflow kinematics and vice versa. In some conditions, baseflow can represent up to 50% of total runoff. HYDROTEL, a semi-distributed hydrological model, has been applied to simulate river flow rates in the Châteauguay River watershed. In HYDROTEL, the groundwater flow behavior is approximated using a linear reservoir. This approach does not allow the simulation of groundwater levels and provides an incorrect estimate of baseflows. This limitation can become problematic under hydrologic stress conditions where the interaction between surface and groundwater can change, such as expected with modified land uses or global warming. For a better understanding and assessment of present and future conditions, it is essential to improve the groundwater flow representation in HYDROTEL, in order to have a comprehensive insight into surface water – groundwater interactions. 3 Objectives The main objective of this thesis is to develop a complete representation of surface water flow, interflow, and groundwater flow incorporated in the HYDROTEL model. The modified model should be able to simulate surface water and groundwater interaction at the basin scale, providing for the first time on the Châteauguay River basin, a complete and accurate simulation of the aquifer and watershed systems. Specific objectives are as follows: 1) to develop a representation of hillslope groundwater flow linked both to the unsaturated zone and to the confined aquifer with two-way interaction potential; 2) to include deep groundwater flow and bidirectional links with subsurface flow and river flow; 3) to develop a computational efficient and applicable model for the watershed scale; 4) to use the improved HYDROTEL model to simulate steady state and transient state past and current conditions on the Châteauguay River watershed

3. Improving groundwater flow representation in a semi-distributed hydrological model using hillslope Boussinesq and regional groundwater flow equations, with application to the Châteauguay watershed 4.2 Future representation of groundwater flow In the course of this project, the BV3C module will be modified by eliminating the third layer and keeping the two top layers which are representing the vadose zone and interface between atmosphere and unconfined aquifer. Underneath a third layer will be added, in which the hillslope-storage Boussinesq equation (hsB) will simulate unconfined rapid interflow (Hilberts et al., 2005). Below, an additional layer will be implemented to represent unconfined and confined, saturated, regional groundwater flow. This will be accomplished by application of the Analytic Element based GFLOW model (Hatjema, 1995). The application of the hsB and GFLOW model guarantees low dimensionality and adaptability to complex profile and plan form morphologies. On the Châteauguay basin, this will allow the simulation of groundwater flow in parts of the basin, where bedrock is very close to the surface and provide large contributions to river baseflow. This configuration will be useful in other aquifers where shallow and deep groundwater flow is present. Figure 1 illustrates the proposed representation, including the planned interaction potentials. 4 Research methodology 4.1 Current presentation of HYDROTEL In the current version of HYDROTEL, the subsurface processes are simulated once on each subwatershed, also named Relatively Homogeneous Hydrological Unit (RHHU). This soil column is subdivided into three homogenous layers, called BV3C. The surface layer is thin, to control infiltration and to represent the evaporation out of the soil column. The second layer and the third layer represent interflow and baseflow respectively. Exchange fluxes between layers 1-2 and 2-3 are accomplished by solving the one-dimensional Richards equation. Discharge from the third layer is calculated using water content and an empirical baseflow recession coefficient. At the end of each time step, runoff, interflow and baseflow are summed and become available for the surface routing scheme. The current subsurface representation in HYDROTEL is unable to simulate groundwater flow itself due to the absence of groundwater flow equations. It cannot either simulate its interaction with surface water flow. Fig. 1: future representation of groundwater flow in HYDROTEL within one RHHU

4. 1 day 5 days 10 days Improving groundwater flow representation in a semi-distributed hydrological model using hillslope Boussinesq and regional groundwater flow equations, with application to the Châteauguay watershed 5 Preliminary results Leakage was calculated according to Darcy's law, and hydraulic heads were scanned in the centre of the y-direction at an upslope (refers to A in Figure 2), midslope (B), and downslope (C) position. Results indicate that leakage is highly dependent on soil parameterization, but also exhibits strong dependencies on hillslope inclination (Figure 3, straight plan form geometry) and geometry (Figure 4). Generally, observations have shown leakage percolating in both directions. In particular, geometry is driving the partitioning of the hillslope in up- and downwelling leakage areas (Figure 4), with convergent hillslopes containing the largest portions of upward directed leakage. Besides, hillslope inclination tends to overshadow retention effects of medium to low conductive aquifer materials, leading to quick upslope dry out, and contribution of the confined aquifer to outflow. 5.1 Evaluation of leakage The leakage process was evaluated in a hillslope hydrological context, e.g. dependencies on soil and aquitard parameterization, hillslope parameters (geometry, slope), and boundary conditions (Dirichlet, seepage face, mixed boundaries at the downslope face) were evaluated using the 3D Richards equation solver CATHY (CATchment HYdrology; Bixio et al., 2000), a model based on the finite-element scheme in space and therefore applicable to catchments of arbitrary geometrical shape. Tests were conducted on laboratory hillslopes of straight, convergent, and divergent plan form shape. Figure 3: Darcy leakage (lines; primary axis) calculated at downslope transect and outflow rate (symbols; secondary axis) as function of slope angle Figure 2: Spatial discretization of the uniform hillslope, location of nodal scans (a), convergent (b), and divergent hillslope (c), initial vertical discretization (d); Figure 4: dimensionless leakage for straight, convergent, divergent hillslope plan form shape at t = 10 days 5.2 Coupling of the hsB and GFLOW-model An iterative coupling scheme was developed, calculating leakage which satisfies both, the hsB aquifer as well as the GFLOW aquifer. A medium conductive soil (K=1e-4 m/s) and a three orders of magnitude lower aquitard conductivity were applied in order to create a realistic semi-pervious unit. Figure 5 depicts calculated water tables for the coupled hsB-GFLOW model (in black) and is compared to the hsB model stand alone (grey lines), showing a gently lower head for the coupled model indicating somewhat low leakage towards the deep aquifer. The red line depicts the calculated water table in the GFLOW aquifer which indicates a pretty realistic low gradient. Calculated outflow rates indicate consistency with the results for the water tables between the coupled and hsB stand alone model. Figure 5: calculated water tables for the uncoupled hsB-model, coupled hsB-GFLOW model (both for heads in unconfined unit), and heads in the confined unit (GFLOW)