The ability to model groundwater and surface water flows as two interacting components of a single resource is highly important for robust catchment management. Existing methods for spatially-distributed numerical modelling of flow in connected river-aquifer systems treat rivers and aquifers as separate sub-domains, with different governing equations for the flow in each. Mass-fluxes exchanged between the sub-domains are modelled using one of several coupling methods, which do not accurately represent the physics of the flow across the interface between the surface and subsurface flows. This can be problematic for model stability and mass conservation. This thesis investigates the feasibility of modelling interacting surface water and groundwater flows in a single domain, using a single system of equations. It is shown that the governing equations in existing numerical models for river and aquifer flow can be derived from the Navier-Stokes Equations. A time- and space-averaged form of Navier-Stokes Equations, the Double-Averaged Navier-Stokes (DANS) Equations, can be used to model both groundwater and surface water flows. The volume- averaging process allows the porous medium to be represented as a continuum. A novel two-dimensional numerical model is developed from the DANS Equations to simulate flows in connected groundwater and surface water systems. The DANS equations are solved using the finite-volume method. The model simulates two-dimensional flow in a vertical slice.

This allows the horizontal and vertical velocity components and pressure to be modelled over the depth of a stream and the underlying aquifer or hyporheic zone. The model does not require the location of the interface between surface and subsurface flows to be specified explicitly: this is determined by the spatial distribution of hydraulic properties (permeability and porosity). The numerical model handles the transition between laminar and turbulent flows using an adaptive damping approach to modify the terms in a single-equation turbulence model, based on a locally-defined porous Reynolds number, Rep. This approach removes the need to specify a priori whether flows in any part of the domain are laminar or turbu lent. Turbulent porous media flows can be simulated. The model is verified for porous-media and clear-fluid flows separately, before being used to simulate coupled groundwater - surface water flow scenarios. For porous-media flows with low Rep the numerical model results agree exactly with Darcy’s Law. The value of Rep at which the model results begin to deviate from Darcy’s Law is consistent with published values. For turbulent clear-fluid flows the time-averaged velocity and turbu lent kinetic energy (TKE) results from the numerical model are ver ified against a RANS model and published data. A good match is achieved for both velocity and TKE. Energy grade-line slopes for free-surface flows simulated in the numerical model are a reasonably good match to equivalent results to the one-dimensional hydraulic model HEC-RAS. Idealised river-aquifer interaction experiments are conducted in a lab- oratory flume to provide verification data for the numerical model. An innovative combination of optical flow measurement and refractive- index-matched transparent soil is used to measure two-dimensional velocities and turbulent statistics in laboratory flow scenarios that

simulate flow in both losing and gaining streams, and the underlying connected porous layer. The “gaining stream” laboratory scenario is replicated using the numerical model. The model simulates the key features of the mean flow well. Turbulent statistics deviate substantially from the laboratory results where vertical velocities across the surface-subsurface interface are high, but are a better match elsewhere. The “losing stream” laboratory results are unable to be reproduced with the numerical model. Results for a similar scenario with lower outflow velocities are presented. These results are qualitatively consistent with the laboratory results. The numerical model is expected to perform better in simulations of field-like conditions that involve less extreme gradients than the laboratory scenarios.