Systems and Methods for Determining Ground Water-Surface Water Interactions

TECHNICAL FIELD

The present disclosure relates to environmental monitoring and more particularly to surface water and ground water monitoring.

BACKGROUND

Surface water flux across the soil-water boundary in ground water/surface water systems governs a number of ecological functions that are critical, for example, to riverine/lacustrine health and contaminant transport. However, traditional methods for monitoring provide limited, or partial information. For example, flux chambers provide the total mass flux but they do not provide any information on pore velocity or residence time. Vertical temperature arrays provide information on pore velocity and residence time, but do not provide information on total mass flux. Accordingly, improved probes and methods of sensing are needed.

SUMMARY OF THE DISCLOSURE

Systems for determining GW/SW interaction are provided. The systems can include: a sensing assembly comprising sensors for pressure, fluid conductivity, temperature, and transfer resistance; processing circuitry operatively coupled to the sensing assembly and configured to receive data from the sensing assembly and process the data to provide a GW/SW interaction, wherein the data includes pressure, fluid conductivity, temperature, transfer resistance data.

Methods for determining GW/SW interaction are provided. The methods can include: receiving real time data including pressure, fluid conductivity, temperature, and transfer resistance; from at least some of the data received simulating the SW/GW interaction; and fitting the real time data with the simulated data to provide actual SW/GW interaction.

BRIEF DESCRIPTION OF THE FIGURES

FIG. 1 depicts a system according to an embodiment of the disclosure.

FIG. 2A is a depiction of a sensing assembly according to an embodiment of the disclosure.

FIG. 2B is a depiction of another sensing assembly according to an embodiment of the disclosure.

FIG. 3 is a depiction of components of the sensing assembly of FIG. 2B according to an embodiment of the disclosure.

FIG. 4 is a system for determining SW/GW interactions according to an embodiment of the disclosure.

FIG. 5. is another configuration of a system for determining SW/GW interactions according to an embodiment of the disclosure.

FIG. 6 is a portion of a system according to an embodiment of the disclosure.

FIG. 7 is a portion of methods for determining SW/GW interactions according to an embodiment of the disclosure.

FIG. 8 is another portion of methods for determining SW/GW interactions according to an embodiment of the disclosure.

FIG. 9 is another portion of methods for determining SW/GW interactions according to an embodiment of the disclosure.

FIG. 10 is a portion of methods for determining SW/GW interactions according to an embodiment of the disclosure.

FIG. 11 is modeled data for use in testing the systems of the disclosure.

FIG. 12 is modeled data for use in testing the systems of the disclosure.

FIG. 13 is fit data using the modeled data according to embodiments of the disclosure.

FIG. 14 is an experimental apparatus and sensing assembly configuration according to an embodiment of the disclosure.

FIG. 15 is another depiction of the experimental apparatus of FIG. 14.

FIG. 16 is data acquired using systems according to embodiments of the disclosure.

FIG. 17 is data acquired using systems according to embodiments of the disclosure.

FIG. 18 is data acquired using systems according to embodiments of the disclosure.

FIG. 19 is data acquired using systems according to embodiments of the disclosure.

FIG. 20 is another sensing assembly configuration used in combination with an experimental apparatus.

FIG. 21 is data acquired using systems according to embodiments of the disclosure.

FIG. 22. is data acquired using systems according to embodiments of the disclosure.

FIG. 23 is data acquired using systems according to embodiments of the disclosure.

FIG. 24 is data acquired using systems according to embodiments of the disclosure.

FIG. 25 is data acquired using systems according to embodiments of the disclosure.

FIG. 26 is data acquired using systems according to embodiments of the disclosure.

FIG. 27 is data acquired using systems according to embodiments of the disclosure.

FIG. 28 is data acquired using systems according to embodiments of the disclosure.

DESCRIPTION

Disclosed are probes and methods for sensing parameters at an interface between groundwater and surface water. Groundwater/surface water flux monitoring approaches that estimate either transient pore fluid velocity or mass flux rate at the surface water/groundwater interface are limited in usefulness. Described herein are multi-sensor probes and methods for sensing that continuously monitors the spatial (e.g., vertical) distribution of pore fluid conductivity, temperature, pressure, and bulk electrical conductivity. Combined with fluid conductivity, bulk electrical conductivity can be used to estimate the vertical distribution of porosity, which links pore fluid velocity to mass flux. We developed the capability to simulate all time-series data from the probe using PFLOTRAN-E4D, and a joint Occam's inversion for estimating the simplest vertical distribution of porosity, permeability, and dispersivity that honor the data. We have also developed a companion stochastic analysis to investigate uncertainty in parameter estimates and corresponding flux rates. Once parameters are estimated, transient pore fluid velocity and mass flux can be monitored using only the pressure sensors located at the top and bottom of the sensor probe to determine dynamic pore fluid velocity and mass flux.

The systems and methods of the present disclosure were developed to address capability limitations of current sensing techniques to measure the dynamic flux of water across the GW/SW boundary. Surface water flux across the soil-water boundary in ground water/surface water system governs a number of ecological functions that are critical to riverine/lacustrine health and contaminant transport. There are currently two predominant methods to monitor surface water flux; flux chambers and vertical temperature arrays. Flux chambers provide the total mass flux but they do not provide any information on pore velocity or residence time. Vertical temperature arrays provide information on pore velocity and residence time, but do not provide information on total mass flux. The systems and/or methods of the present disclosure can be used to monitor dynamic (i.e. time varying) pore fluid velocity, residence time, and total mass flux, with the added benefits of estimating the in situ porosity and hydraulic conductivity of the sediment adjacent to the probe.

Embodiments of the present disclosure can combine electrical geophysical measurements (collected using ring electrodes on the outer shell of the sensing assembly) with traditional pressure, temperature, and fluid conductivity measurements (collected using sensors installed on the inner spine, inside of the shell of the sensing assembly). The geophysical measurements are sensitive to porosity, which provides the parameter necessary to convert pore fluid velocity to mass flux.

The systems and methods of the present disclosure can utilize a thermal contrast between surface water and groundwater as well as contrasts in fluid electrical conductivity between surface water and groundwater. Temperature based flux measurements only provide pore fluid velocity, and not mass flux, and assume steady state pressure boundary conditions. The systems and methods of the present disclosure can provide both pore fluid velocity and mass flux under dynamic pressure boundary conditions. Temperature based flux measurements do not provide information about vertical sediment structure. The systems and methods of the present disclosure can provide the vertical distribution of sediment porosity and permeability.

The systems and methods of the present disclosure can use geophysical measurements to estimate in situ sediment porosity and permeability which allows for the conversion boundary pressure measurements to Darcy's flux and pore velocity. These natural tracer flux measurements can then used to calibrate differential pressure to use the Darcian approach to measure dynamic flux even when natural tracers cannot be used. In most cases where there is a contrast in fluid conductivity between GW and SW these measurements can be used to determine if the flux of water is truly GW, SW or a mix of GW and SW which can all occur under dynamic stage conditions.

The sensing assembly can include a vertical array of sensors (electrodes, thermistors, fluid specific conductance probes, and pressure sensors) that are installed into the sediment at the lower boundary of a water body (e.g. stream, river, lake or tidal zone). The sensors can be deployed on a linear, non-metallic conduit (approximately 2-5 cm in diameter and 50-300 cm long). The conduit can include a series of ring electrodes that can be used to monitor the bulk electrical conductivity time-series of the sediment along the longitudinal axis. A conduit can also be equipped with miniature sensors that monitor pressure, temperature, and fluid conductivity along the vertical axis of the tool.

Referring first to FIG. 1, a sensing assembly 16 is shown within a GW/SW environment 10. As depicted ground water exists in earth 14 and below surface water 12 and GW/SW barrier 15. Sensing assembly 16 can include sensors configured to sense changes in multiple physical parameters 20. These physical parameters can include pressure, fluid conductivity, temperature, and transfer resistance. Sets of these sensors can be provided for each of the physical parameters. These sets of sensors can be spaced apart predetermined distances and these distances can be associated with a depth below the GW/SW barrier. Accordingly, the systems and methods can be used to determine vertical GW/SW interactions which in turn can be used to determine porosity, permeability and/or dynamic flow processes.

In accordance with at least one embodiment of the disclosure and with reference to FIGS. 2A, 2B, and 3, sensing assemblies are depicted that include both inner 22 and outer rods 24 of the flux tool. Inner rod 22 can include upper 25a and lower 25b pressure sensors and fluid conductivity sensors 26. Pore fluid within the groundwater may accesses the fluid conductivity sensors via ports installed in the outer rod. Outer rod 24 can include embedded thermistors 27 (temperature sensors) and bulk conductivity electrodes 28. Referring to FIG. 3, sets of these sensors can be provided along the length of the sensing assembly a predetermined distance from one another; for example, pressure sensors 25a-b, fluid conductivity sensors 26a-g, temperature sensors 27a-g, and bulk conductivity electrodes 28a-u.

In accordance with example implementations, the rod can be installed vertically into a river-bottom or other groundwater/surface water interface. Data from each sensor can then be collected.

Accordingly, sensors can record the pressure, fluid conductivity, and temperature time-series along the length rod, including the upper and lower boundary conditions. The ring electrodes can be used to make transfer resistance measurements in a 4-electrode configuration, for example. The thermistors (temperature sensors) can be calibrated to operate between 5-30 degrees C. The fluid conductivity cells can use four-electrodes built into a flexible circuit board that can be installed around the circumference of the inner rod (see FIGS. 2-3). When the inner rod is inserted into the outer rod, the O-rings on the inner rod seal to the inside surface of the outer rod, thereby bounding the electrical field within the measurement chamber. Vertical slits on the outer rod allow water to pass through the sensor. These slits can be located at the top and bottom of the sensor cell and were 1 cm long and 0.5 mm wide and can be located 10 cm apart from top to bottom of the cell and also located on opposite sides of the sensing assembly. The cell constant for the fluid conductivity sensors can be calibrated using 4 sodium chloride solutions ranging from 0.0044 to 0.298 S/m. Nylon tubes can be used between the pressure ports and the differential pressure transducers located in the sensing assembly head.

Referring to FIG. 4, this data can be processed to determine GW/SW interaction. In accordance with example implementations, some fundamentals of SW/GW relationships are relied upon. For example, Porosity can be estimated using Archie's law (Archie, 1942). For a saturated system, Archie's law relates the bulk conductivity (σ_bulk)t to the fluid conductivity (ø_fluid) using the formation factor (F) seen in Equation 1.

σ_Bulk=σ_fluidF (1)

σ_Bulkcannot be measured directly, but transfer resistance measurements can be inverted to obtain or σ_Bulk.

Furthermore, F is a function of porosity (n) and the cementation exponent (m) by the relationship seen in Equation 2.

F=n
^−m (2)

For granular media, particle shape is a major factor affecting m. For spherical particles m is about 1.35, but can be as high as 1.65 for prismatic and angular shaped particles.

Pore water velocity can be estimated using time series data from different locations along the flow path of a conservative tracer. In most GW-SW interaction zones, electrical conductivity (EC) can be used as a natural conservative tracer due to the differences between the electrical conductivity of GW and SW. Freshwater lakes and streams have specific conductance values ranging from 2-100 μS/cm, GW, wetlands and bogs generally range from 50-50,000 μS/cm, and ocean water is generally around 50,000 μS/cm. According to EPA, river water in the United States ranges in conductivity from 50-1500 μS/cm.

For a 3-D system, bulk electrical conductivity of the subsurface can be imaged by inverting time lapse electrical resistivity tomography (ERT) measurements. This provides spatial and temporal conductivity estimates at locations where the system is sensitive to ERT measurements. This time lapse data can then be used to identify breakthroughs of a conductive tracer passing through the system. Transfer resistance values were directly used to measure breakthroughs of either GW or SW by using their difference in electrical conductivity as a natural tracer.

For both gaining and losing steams, the temperature profile can be used as a natural tracer to estimate pore water velocity. This can be utilized wherein there is a temperature contrast between GW and SW, which is generally true for streams and rivers that have diurnal and/or seasonal temperature swings. To date, no analytical solution was found to use temperature as a natural tracer in dynamic systems with unpredictable pressure boundary conditions, thus PFLOTRAN, a parallel multi-physics code used to simulate reactive flow and transport through porous materials, was used to predict temperature at 4 of the 5 thermistor locations for comparison with observed data.

Water passing through porous media experiences a resistance to flow. This resistance to flow is quantified using Darcy's law that relates the specific discharge also known as Darcy's flux (q) to the hydraulic gradient

$(\frac{dh}{dl})$

using the hydraulic conductivity (K), which is considered a property of the porous media (see Equation 3).

$\begin{matrix} q = - K (\frac{dh}{dl}) & (3) \end{matrix}$

For a gaining or loosing stream, quasi steady state Darcy's flux can be estimated using the temperature profile below the river bed. The temperature amplitude and phase shift at different depths driven by diurnal temperature fluctuations of surface water to calculate a vertical flux value can be used for steady state Darcy's flux. This approach may be applied easily in the field. The analytical solution to a one-dimensional steady-state heat advection-diffusion equation can be used to calculate the flux values as applied in Schmidt et al. (2007) (see Equation 4).

$\begin{matrix} q_{z} = - \frac{K_{fs}}{(ρ_{f} C_{f} Z)} \ln (T_{Z} - T_{L}) / (T_{0} - T_{L}) & (4) \end{matrix}$

where, T_zis the streambed temperature at depth z (° C.), T_Lis the lower temperature boundary (° C.), T₀is the temperature at Z=0 (° C.), ρ_fC_fis the volumetric heat capacity of the fluid (J m⁻³K⁻¹), and K_fsis the thermal conductivity of the solid-fluid system (J s⁻¹m⁻¹K⁻¹). K_fsis calculated using the geometric mean of the thermal conductivity of sand K_sand the thermal conductivity of fluid K_fas can be seen in Equation 5, where n is the porosity.

K
_fs
=K
_s
^(1-n)
K
_f
ⁿ (5)

K_sis assumed to be between 7.7 to 8.4 J s⁻¹m⁻¹K⁻¹which is the thermal conductivity of solid quartz. However, the thermal conductivity of saturated sand was best fit to 2.86 J s⁻¹m⁻¹K⁻¹where the relationship of the average grain size d_m(micron) to the thermal conductivity of solid-fluid system was calculated using Equation 6.

K
_s=0.215 ln d_m+1.93 (6)

Referring next to FIG. 4, a system 40 is depicted that includes sensing assembly 16 operationally coupled to processing circuitry 42 which is configured to provide GW/SW interaction data 44. In accordance with example implementations, processing circuitry 42 can be operationally coupled to sensing assembly 16 via cables, wirelessly, and/or integrated as part of sensing assembly 16.

Referring next to FIG. 5, system 40 is shown receiving pressure (25), fluid conductivity (26), temperature (27), and transfer resistance (28) data from the sensing assembly. Each of these parameters can be associated with a time and/or a depth below the SW/GW barrier.

Referring next to FIG. 6, and example embodiment of processing circuitry 42 is shown with which the sensing assembly can operationally interfaced. In one example, a computing environment such as shown in FIG. 6 can be used to log data from the sensing assembly, to analyze data from the sensing assembly and/or to perform simulations based on sensing assembly data, as well as fit real time and/or recorded/stored data with simulated data. In accordance with example implementations, the processing circuitry is configured to provide simulated SW/GW interaction and real time SW/GW interaction, and process those interactions to fit the simulated SW/GW interaction with the real time SW/GW interaction to provide a modified SW/GW interaction. The modified SW/GW interaction can be provided continually in an iterative process.

For example, processing circuitry 42 can be used to simulate all time-series data from the sensing assembly using PFLOTRAN-E4D, and a joint Occam's inversion for estimating the vertical distribution of porosity, permeability, and dispersivity from the data received. Furthermore, a companion stochastic analysis can be performed with the processing circuitry to determine uncertainty in parameter estimates and corresponding flux rates. Once parameters are estimated, transient pore fluid velocity and mass flux can be monitored using only the pressure sensors located at the top and bottom of the sensing assembly.

The time series data can be jointly inverted as described below to estimate the vertical distribution of porosity and permeability, pore fluid velocity and vertical mass flux time-series.

The combined monitoring data are processed using analysis techniques developed specifically for the systems and methods of the present disclosure. The techniques involving the use of U.S. Pat. No. 9,772,423 may be utilized, the entirety of which is incorporated herein by reference.

With reference to FIG. 6, an example system for implementing some embodiments includes a general-purpose computing device in the form of a computer 210. Components of computer 210 may include, but are not limited to, a processing unit 220 (which is not limited to CPUs, but can comprise GPUs), a system memory 230, and a system bus 221 that couples various system components including the system memory to the processing unit 220. The system bus 221 may be any of several types of bus structures including a memory bus or memory controller, a peripheral bus, and a local bus using any of a variety of bus architectures. By way of example, and not limitation, such architectures include Industry Standard Architecture (ISA) bus, Micro ChannelArchitecture (MCA) bus, Enhanced ISA (EISA) bus, Video Electronics Standards Association (VESA) local bus, and Peripheral Component Interconnect (PCI) bus also known as Mezzanine bus. Memory and programs described herein be deployed in corresponding portions of FIG. 6.

Computer 210 typically includes a variety of computer readable media. Computer readable media can be any available media that can be accessed by computer 210 and includes both volatile and nonvolatile media, removable and non-removable media. By way of example, and not limitation, computer readable media may comprise computer storage media and communication media. Computer storage media is different from, and does not include, a modulated data signal or carrier wave. It includes hardware storage media including both volatile and nonvolatile, removable and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or other data. Computer storage media includes, but is not limited to, RAM, ROM, EEPROM, sash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical disk storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information and which can be accessed by computer 210. Communication media typically embodies computer readable instructions, data structures, program modules or other data in a transport mechanism and includes any information delivery media. The term “modulated data signal” means a signal that has one or more of its characteristics set or changed in such a manner as to encode information in the signal. By way of example, and not limitation, communication media includes wired media such as a wired network or direct-wired connection, and wireless media such as acoustic, RF, infrared and other wireless media. Combinations of any of the above should also be included within the scope of computer readable media.

The system memory 230 includes computer storage media in the form of volatile and/or nonvolatile memory such as read only memory (ROM) 231 and random-access memory (RAM) 232. A basic input/output system 233 (BIOS), containing the basic routines that help to transfer information between elements within computer 210, such as during startup, is typically stored in ROM 231. RAM 232 typically contains data and/or program modules that are immediately accessible to and/or presently being operated on by processing unit 220. By way of example, and not limitation, FIG. 6 illustrates operating system 234, application programs 235, other program modules 236, and program data 237.

The computer 210 may also include other removable/nonremovable volatile/nonvolatile computer storage media. By way of example only, FIG. 6 illustrates a hard disk drive 241 that reads from or writes to non-removable, nonvolatile magnetic media, and an optical disk drive 255 that reads from or writes to a removable, nonvolatile optical disk 256 such as a CD ROM or other optical media. Other removable/non-removable, volatile/nonvolatile computer storage media that can be used in the exemplary operating environment include, but are not limited to, magnetic tape cassettes, sash memory cards, digital versatile disks, digital video tape, solid state RAM, solid state ROM, and the like. The hard disk drive 241 is typically connected to the system bus 221 through a nonremovable memory interface such as interface 240, and optical disk drive 255 are typically connected to the system bus 221 by a removable memory interface, such as interface 250.

Alternatively, or in addition, the functionality described herein can be performed, at least in part, by one or more hardware logic components. For example, and without limitation, illustrative types of hardware logic components that can be used include Field-programmable Gate Arrays (FPGAs), Program-specisc Integrated Circuits (ASICs), Program-specisc Standard Products (ASSPs), System-on-a-chip systems (SOCs), Complex Programmable Logic Devices (CPLDs), etc.

The drives and their associated computer storage media discussed above and illustrated in FIG. 6, provide storage of computer readable instructions, data structures, program modules and other data for the computer 210. In FIG. 6, for example, hard disk drive 241 is illustrated as storing operating system 244, application programs 245, other program modules 246, and program data 247. Note that these components can either be the same as or different from operating system 234, application programs 235, other program modules 236, and program data 237. Operating system 244, application programs 245, other program modules 246, and program data 247 are given different numbers here to illustrate that, at a minimum, they are different copies.

A user may enter commands and information into the computer 210 through input devices such as a keyboard 262, a microphone 263, and a pointing device 261, such as a mouse, trackball or touch pad. Other input devices (not shown) may include a joystick, game pad, satellite dish, scanner, or the like. These and other input devices are often connected to the processing unit 220 through a user input interface 260 that is coupled to the system bus, but may be connected by other interface and bus structures, such as a parallel port, game port or a universal serial bus (USB). A visual display 291 or other type of display device is also connected to the system bus 221 via an interface, such as a video interface 290. Video interface 290 can comprise a graphics card having a GPU. The GPU be used for computations. In addition to the monitor, computers may also include other peripheral output devices such as speakers 297 and printer 296, which may be connected through an output peripheral interface 295.

The computer 210 is operated in a networked environment using logical connections to one or more remote computers, such as a remote computer 280. The remote computer 280 may be a personal computer, a hand-held device, a server, a router, a network PC, a peer device or other common network node, and typically includes many or all of the elements described above relative to the computer 210. The logical connections depicted in FIG. 6 include a local area network (LAN) 271 and a wide area network (WAN) 273, but may also include other networks. Such networking environments are commonplace in offices, enterprise-wide computer networks, intranets and the Internet.

When used in a LAN networking environment, the computer 210 is connected to the LAN 271 through a network interface or adapter 270. When used in a WAN networking environment, the computer 210 typically includes a modem 272 or other means for establishing communications over the WAN 273, such as the Internet. The modem 272, which may be internal or external, may be connected to the system bus 221 via the user input interface 260, or other appropriate mechanism. In a networked environment, program modules depicted relative to the computer 210, or portions thereof, may be stored in the remote memory storage device. By way of example, and not limitation, FIG. 6 illustrates remote application programs 285 as residing on remote computer 280. It will be appreciated that the network connections shown are exemplary and other means of establishing a communications link between the computers may be used.

Referring next to FIG. 7, an example process 71 used by processing circuitry 42 for processing transfer resistance data is shown. Accordingly transfer resistance data is received by processing circuitry 42 at 70. Using processing circuitry 42, at 72 an ERT inversion of the data of 70 is executed using E4D Software. (E4D software is available from the Pacific Northwest National Laboratory).

E4D is a parallel ERT forward modeling and inversion code. Data are collected along the sensing assembly using a four-electrode Wenner-alpha configuration, whereby current is injected between two outer source and a sink electrodes, and the resulting electrical potential is measured between a positive and negative potential electrode, each lying between the current electrodes. Measurements proceed down the sensing assembly using 4 adjacent electrodes in sequence to produce thirteen measurements at a given time step.

The direct-current electrical potential arising from a given current source is given by the Poisson equation (7)

∇·σ_b∇φ=−∇·J, (7)

with bulk conductivity σ_b, electrical potential ϕ, and source current density J. Given a_b, J, and appropriate boundary conditions, the objective of the forward modeling is to solve Equation (7) for ϕ, which is measured during an ERT experiment. Simulation of the voltage across potential electrodes normalized by the injected current, namely the transfer resistance, is summarized as

R=A
_e(σ_b), (8)

with transfer resistance R, and ERT forward operator A_e. In time-lapse ERT surveys, both the transfer resistance and bulk conductivity are functions of time. E4D solves Equation (7) using the finite element method on an unstructured tetrahedral mesh. The ring electrodes and rod of the flux tool are modeled in true dimension using the non-point source electrode modeling method described by Johnson and Wellman [2015]. At 74, using processing circuitry 42, a vertical bulk electrical conductivity profile is provided.

Referring next to FIG. 8, an example process used by processing circuitry 42 for processing pressure, temperature, and fluid conductivity boundary data over time is shown. Accordingly, the data is received by processing circuitry 42 and PFLOTRAN (Open Source) simulation is performed at 82 to prepare a simulated time series of pressure temperature and fluid conductivity from data 20.

PFLOTRAN is a parallel multi-physics code that simulates reactive flow and transport in porous materials. PFLOTRAN simulates liquid phase flow and solute transport based on the following conservation equations of mass, energy, and solute transport

$\begin{matrix} \frac{\partial}{\partial t} (φ s ρ) + \nabla \cdot (ρ q) = 0, & (9) \end{matrix}$

$\begin{matrix} \frac{\partial}{\partial t} (φ s ρ U + (1 - φ) ρ_{r} c_{r} T) + \nabla \cdot (ρ qH - κ \nabla T) = 0, & (10) \end{matrix}$

$\begin{matrix} \frac{\partial}{\partial t} (φ sC) + \nabla \cdot (q - φ sD \nabla) C = 0, & (11) \end{matrix}$

with porosity φ, saturation s, water density ρ, Darcy velocity q, fluid internal energy U, rock density ρ_r, rock heat capacity C_r, water temperature T, fluid enthalpy H, bulk thermal conductivity κ, solution concentration C, and dispersion coefficient D. The Darcy velocity (or Darcy flux) q is given by

$\begin{matrix} q = - \frac{k}{μ} \nabla (P - ρℊ z), & (12) \end{matrix}$

with permeability k, viscosity μ, water pressure P, acceleration of gravity g, vertical position z. The Darcy velocity q gives discharge per unit area, which is not the actual fluid speed traveling through the porous medium. Instead, the actual pore fluid velocity (or seepage velocity) v is computed as

$\begin{matrix} v = - \frac{q}{φ} . & (13) \end{matrix}$

For SW/GW interactions in this paper, we assume incompressible flows (i.e. water density is a known constant value) and the following relations

$\begin{matrix} U = H - \frac{P}{ρ} & (14) \end{matrix}$

$\begin{matrix} dH = c_{p} dT + V (1 - α T) dP & (15) \end{matrix}$

with water heat capacity c_p, system volume V, and thermal expansion coefficient α. Equation (14) is obtained from an equation of state for pure water; Equation (15) describes the dependency of enthalpy as a function of temperature and pressure.

There can be 6 scalar unknowns (water pressure P, 3 components of Darcy velocity q, water temperature T, and solute concentration C), which can be obtained by solving 6 scalar equations (conservation of mass, 3 components of Darcy’ Law, conservation of energy, and conservation of solute transport). PFLOTRAN employs backward Euler time discretization and the Newton-Raphson method to solve for flow, heat, and transport state variables over time. In the 1D mass flux monitoring and inversion at the GW-SW interface, the unknown model parameters are the vertically distributed permeability and porosity. Given the boundary and initial conditions of pressure, temperature, and fluid conductivity at the top and bottom of the system, the hydrogeologic modeling of dynamic temperature and fluid conductivity can be summarized as

(T,σ_w)=A_h(k,φ), (16)

with fluid conductivity σ_w, and hydrogeologic forward operator A_h. The hydrogeologic forward operator in Equation (16) represents the process of solving the partial differential equations (9) to (11) with the relation Equation (12) and appropriate boundary and initial conditions.

At 86 a combination of PFLOTRAN and E4D mapping is performed using processing circuitry 42. PFLOTRAN-E4D is a parallel ERT module for PFLOTRAN to simulate the evolution of subsurface resistivity and corresponding transfer resistance measurements arising from a particular process [Johnson et al., 2017]. The code receives the evolution of subsurface states from PFLOTRAN, converts those states to bulk electrical conductivity, and simulates a specified set of transfer resistance data.

We use Archie's Law [Archie, 1942] to describe the petrophysical relationship between the fluid conductivity and bulk conductivity as

$\begin{matrix} σ_{b} = \frac{1}{a} φ^{m} s^{n} σ_{w}, & (17) \end{matrix}$

with tortuosity factor a, cementation exponent m, and saturation exponent n. The formation factor F is defined as the ratio of the fluid conductivity to the bulk conductivity of saturated sediments (s=1)

$\begin{matrix} F = \frac{a}{φ^{m}} = \frac{σ_{w}}{{(σ_{b})}_{s = 1}} . & (18) \end{matrix}$

With predetermined parameters of a, m, and n, the petrophysical transform in Equation (17) can be denoted by

σ_b=A_p(σ_w,φ). (19)

For the saturated applications described herein, s=1 so s and n do not contribute to solution uncertainty. Conversely, a and m contribute to uncertainty in the inversion solution as described in the forthcoming discussion. The hydrogeophysical modeling of the data collected at the SW/GW interface can be represented as

(T,σ_w,R)=A(k,φ). (20)

where A denotes the hydrogeophysical forward operator that integrates operators A_h, A_p, and A_e.

A(k,φ)=(A_h(k,φ),A_p(σ_w,φ),A_e(σ_b)) (21)

There can be three steps in the hydrogeophysical forward modeling, (1) PFLOTRAN simulation of temperature and fluid conductivity, given the boundary and initial conditions of pressure, temperature, and fluid conductivity measured at the top and bottom of the flux tool (Eq. 13); (2) computation of bulk conductivity using Archie's Law after mapping porosity and fluid conductivity from PFLOTRAN to the E4D mesh (Eq. 8); and (3) E4D simulation of time-lapse transfer resistance data (Eq. 9).

We further simulate the hydrogeophysical data (i.e. temperature, fluid conductivity, and transfer resistance) at N_r^sumlocations (total number of thermistors, fluid conductivity sensors, and Wenner-alpha arrays, sensors used for boundary conditions excluded) and N_ttime steps. Using discrete data and model parameters, the 1D hydrogeophysical modeling can be written in matrix form as

d=A(m), (22)

where d=[T^T,σ_w^T,R^T]^Tdenotes the N_dlength data vector, N_d=N_r^sumN_t, A denotes the discretized hydrogeophysical operator A in Eq. (21), and m=[k^T,φ^T]^Tdenotes the N_mlength model vector of layered permeability and porosity, N_m=2N_l. Superscript T indicates vector transpose.

At 88, time series of bulk conductivity is received and processed using E4D to prepare a time series of transfer resistance data at 92. At 94, simulated time series data for pressure, temperature, fluid conductivity, and transfer resistance is provided.

Referring next to FIG. 9, and iterative process 100 for fitting observed data with simulated data is provided. At 102 observed data is received from the sensing assembly. At 104, the hydro parameter (porosity and permeability) layer boundaries provided in step 74 of process 71 are provided. At 106, permeability and porosity layering is determined from the processed data. At 108, the permeability and porosity layering of 106 is updated. At 110, the simulated time series data provided in step 94 is received.

At 112, the observed time series data and the simulated time series data are jointly inverted at 114. Joint inversion can include application of a Jacobian Matrix.

The joint inversion estimates layered hydrogeologic model parameters by minimizing

Φ=∥W_d(A(m)−d_obs)∥²+λ∥W_m(m−m_apr)∥² (23)

with data weighting matrix W_d, observed data vector d_obs, regularization parameter A, model weighting matrix W_m, and apriori model vector m_apr. Two different approaches can be implementated, the conjugate gradient (CG) method, which belongs to the group of deterministic inversion algorithms, and the ensemble Kalman filter (EnKF), which belongs to the group of stochastic inversion algorithms, to solve the optimization problem Eq. (23).

The nonlinear conjugate gradient algorithm as below in Table 1.

TABLE 1

Nonlinear conjugate gradient algorithm.

Given initial guess m₀, predicted data vector A(m₀), Jacobian matrix F_m,

and tolerance tol_r and tol_b, set at first r₀= A(m₀) − d_obs, δ_r=

∥W_dr₀∥², r_δ = ∥W_dd_obs∥², and k = 0.Then:

1
while δ_r> tol_r²r_δ

2
b = F_m^TW_d^TW_dd_obs+ λW_m^TW_mm_apr

3
s_k= F_m^TW_d^TW_dr_k+ λW_m^TW_m(m_k− m_apr)

4
δ_k= ∥s_k∥², b_δ = ∥b∥², p_k= s_k

5
while δ_k> tol_b²b_δ

6
q_k= F_mp_k

7

α_{k} = \frac{δ_{k}}{‖ W_{d} q_{k} ‖^{2} + λ ‖ W_{m} p_{k} ‖^{2}}

8
m_k+1 = mk − α_kp_k

9
r_k+1 = r_k− α_kq_k

10
s_k+1⁼ F_m^TW_d^TW_dr_k+1 + λW_m^TW_m(m_k+1 − m_apr)

11

δ_{k + 1} = ‖ s_{k + 1} ‖^{2}, p_{k + 1} = s_{k + 1} + \frac{δ_{k + 1}}{δ_{k}} p_{k}

12
k = k + 1

13
end

14
Compute A(m_k) and F_m

15
r_k= A(m_k) − d_obs, δ_r= ∥W_dr_k∥²

16
end

with Jacobian matrix F_m, data residual r_k, steepest ascent direction s_k, conjugate gradient p_k, and iteration number k. For 1D inversions F_mcan be computed using a finite difference approximation. There are two loops in the CG algorithm. The main outer loop updates predicted data and Jacobian matrix; the inner loop solves a localized linear system iteratively.

The nonlinear ensemble Kalman filter algorithm is described below with introduction of the following matrix notations

M
₀=[m₀+u⁽¹⁾, . . . ,m₀+u^(N^s⁾]=[m₀⁽ⁱ⁾],

D
_obs=[d_obs+v⁽¹⁾, . . . ,d_obs+v^(N^s⁾]=[d_obs⁽ⁱ⁾],

M
_apr=[m_apr+w⁽¹⁾, . . . ,m_apr+w^(N^s⁾]=[m_apr⁽ⁱ⁾],

where M₀is the initial ensemble matrix of the size N_c×N_swith mean m₀and covariance Q, D_obsis the observed data matrix of the size N_d×N_swith mean d_obsand covariance R=(W_d^TW_d)⁻¹, M_apris the apriori model matrix of the size N_c×N_swith mean m_aprand covariance S=(λW_m^TW_m)⁻¹, and u⁽ⁱ⁾, v⁽ⁱ⁾, and w⁽ⁱ⁾are Gaussian distributed random vectors in the ith column, i=1, . . . , N_s. Vectors m₀⁽ⁱ⁾, d_obs⁽ⁱ⁾, and m_apr⁽ⁱ⁾, are therefore m₀, d_obs, and m_aprplus random vectors of u⁽ⁱ⁾, v⁽ⁱ⁾, and w⁽ⁱ⁾, respectively.

TABLE 2

Nonlinear ensemble Kalman filter algorithm

Given initial guess M₀= [m₀⁽ⁱ⁾], predicted data vector A(m₀), matrix

D₀= [d₀⁽ⁱ⁾] = [A(m₀⁽ⁱ⁾], and tolerance tol_r, set at first δ_r=

∥W_d(A(m₀) − d_obs)∥², r_δ = ∥W_dd_obs∥, k = 0. Then:

1
while δ_r> tol_r²r_δ

2

Δ M_{k} = [m_{k}^{(i)} - \frac{1}{N_{s}} (M_{k} e_{N_{s} \times 1}) e_{1 \times N_{s}}]

Δ D_{k} = [d_{κ}^{(i)} - \frac{1}{N_{s}} (D_{k} e_{N_{s} \times 1}) e_{1 \times N_{s}}]

P_{k} = \frac{Δ D_{k} Δ D_{k}^{T}}{N_{s} - 1} + R

M_{k + \frac{1}{2}} = M_{k} - \frac{Δ M_{k} Δ D_{k}^{T}}{N_{s} - 1} P_{k}^{- 1} (D_{k} - D_{o b s})

Δ M_{k + \frac{1}{2}} = [m_{k + \frac{1}{2}}^{(i)} - \frac{1}{N_{s}} \sum_{i = 1}^{N_{s}} m_{k + \frac{1}{2}}^{(i)}]

P_{k + \frac{1}{2}} = \frac{Δ M_{k + \frac{1}{2}} Δ M_{k + \frac{1}{2}}^{T}}{N_{s} - 1} + S

M_{k + 1} = M_{k + \frac{1}{2}} - \frac{Δ M_{k + \frac{1}{2}} {ΔM}_{k + \frac{1}{2}}^{T}}{N_{s} - 1} P_{k + \frac{1}{2}}^{- 1} (M_{k + \frac{1}{2}} - M_{a p r})

m_{k + 1} = \frac{1}{N_{s}} (M_{k + 1} e_{N_{s} \times 1})

10
k = k + 1

11
Compute A(m_k) and D_k= [d_k⁽ⁱ⁾] = [A(m_k⁽ⁱ⁾]

12
r_k= A(m_k) − d_obs, δ_r= ∥W_dr_k∥²

13
end

where e denotes the matrix of all ones of indicated size. The EnKF algorithm has only one loop (i.e. the main outer loop) that updates predicted data vector and matrix. The inner loop collapses because direct method of matrix inversion is evolved, for example,

$\begin{matrix} M_{k + \frac{1}{2}} = M_{k} - \frac{Δ M_{k} Δ D_{k}^{T}}{N_{s} - 1} {(\frac{Δ D_{k} Δ D_{k}^{T}}{N_{s} - 1} + R)}^{- 1} (D_{k} - D_{obs}) . & (24) \end{matrix}$

The regularization can be treated as a second data set to be fit in the EnKF algorithm (lines 6 to 8 in Table 2). This so-called two-stage EnKF is based on the principle that assimilating two independent observations jointly gives the same result as assimilating them sequentially. Alternatively, one can arrange the measurements and regularization into an augmented data, and then perform assimilation on this new data in a similar way. In such a case, an equivalent formulation of (24) is worth considering for problems with increased data size

$\begin{matrix} M_{k + \frac{1}{2}} = M_{k} - Δ {M_{k} (Δ D_{k}^{T} R^{- 1} Δ D_{k} + (N_{s} - 1) I_{N_{s}})}^{- 1} Δ D_{k}^{T} R^{- 1} (D_{k} - D_{obs}), & (25) \end{matrix}$

where I denotes the identity matrix of indicated size. To achieve an efficient inversion, Equation (24) is suggested when N_d<<N_s, and Eq. (25) is recommended when N_s<<N_d. Cholesky decomposition can be further used for a more economical matrix inversion. For completeness, the CG and EnKF algorithms with regularization are provided above. In the case of overdetermined systems (e.g. 1D inversion of the flux tool data), regularization can be skipped by setting α=0. The EnKF provides uncertainty of parameter estimates, which is drawn from the ensemble variance after data assimilation. However, the uncertainty becomes trivial in overdetermined problems, especially for inverting the flux tool data since redundant measurements are heavily introduced by time-lapse monitoring.

Let C_m=ΔM_kΔM_k^T/(N_s−1) denote the model covariance prior to data assimilation at iteration number k. By representing ΔD_k≈F_mΔM_k, we can rewrite Eq. (24) as the following data space Gauss-Newton formula

$\begin{matrix} M_{k + \frac{1}{2}} = M_{k} - C_{m} {F_{m}^{T} (F_{m} C_{m} F_{m}^{T} + R)}^{- 1} (D_{k} - D_{obs}) . & (26) \end{matrix}$

Similarly Equation (25) can be rearranged into the model space Gauss-Newton formula

$\begin{matrix} M_{k + \frac{1}{2}} = M_{k} - {(F_{m}^{T} R^{- 1} F_{m} + C_{m}^{- 1})}^{- 1} F_{m}^{T} R^{- 1} (D_{k} - D_{obs}), & (27) \end{matrix}$

where F_mis the Jacobian of operator A with respect to m_k. Equations (26) and (27) indicate that the EnKF algorithm is a stochastic implementation of the Gauss-Newton (GN) method. As a result, the EnKF inherits both advantages and disadvantages of the GN method with respect to convergence, e.g., superlinear convergence rate near the solution, local minima due to improper initial guess, step length adjustment for nonlinear inversions, etc. There are two major differences in the implementation of the GN and the EnKF algorithms. First, the prior model covariance C_mdoes not change in the GN iteration but is continuously updated in the EnKF. Second, the EnKF does not require an explicit form of the Jacobian matrix. For example, it approximates F_mC_mF_m^T, by ΔD_kΔD_k^T/(N_s−1) without evaluating F_m.

Referring again to FIG. 9, at 116 convergence of the simulated and observed data is determined using the processing circuitry. Accordingly, converging the data determines if the observed and simulated data match to an accuracy consistent with the noise observed in the time series data. If yes, the process continues to 118 where layered hydro parameters are provided and then the SW/GW interaction provided at 120. If no, the process returns to 106.

In accordance with example implementations, a joint inversion framework can be orchestrated within processing circuitry 42. Referring to FIG. 10, a joint inversion framework is shown. The joint inversion framework for the system is summarized in FIG. 10, which depicts of three main parts, layering analysis, model update per outer iteration, and hydrogeophysical modeling. In the layering analysis, transfer resistance data are inverted with E4D to derive vertical bulk conductivity distribution. The bulk conductivity profile is then analyzed to choose the number and thickness of layers in the inversion domain. In the model update per outer iteration, either the CG or EnKF algorithm (lines 2 to 13 in Table 1 or lines 2 to 10 in Table 2) is used to update the hydrogeologic model of permeability and porosity. In the hydrogeophysical modeling, PFLOTRAN is first used for forward flow simulations with the transient boundary and initial conditions of pressure, temperature, and fluid conductivity measured by the flux tool. After mapping the simulated fluid conductivity and estimated porosity from PFLOTRAN to E4D mesh, transient bulk conductivity is computed using Archie's law, and E4D is called to simulate the time-lapse transfer resistance data. A time-interpolation approach can be used to accelerate time-lapse transfer resistance simulation. The model update per outer iteration and hydrogeophysical modeling are repeated until the difference between the predicted and observed flux tool data is less than the noise level.

Computational time of the CG and EnKF algorithms greatly depends on the efficiency of the hydrogeophysical modeling, especially the time-lapse ERT simulation. The CG needs N_mruns of hydrogeophysical modeling to compute the Jacobian matrix, while the EnKF requires N_sruns of hydrogeophysical modeling to calculate the predicted data of ensemble, per outer iteration in an inversion. Each hydrogeophysical modeling calls PFLOTRAN once and E4D N_ttimes. Even with 1D models, the time-lapse ERT simulation is 3D due to the use of ring electrodes and the cavity caused by the flux tool measuring within the sediment. This results in time-consuming modeling using E4D. The hydrogeologic modeling is 1D on the other hand.

A time-interpolation technique to speed up the 3D time-lapse ERT simulation for 1D models is provided. The time-interpolation of fluid resistivity (reciprocal of fluid conductivity) time series at depth z

ρ_w(z,t_j)=Σ_k=1^N^bc_jkρ_w(z,t_k),j=1, . . . ,N_t, (28)

with fluid resistivity ρ_w, number of interpolation basis N_b, and interpolation coefficient c_jkof the kth basis for the jth datum. The set of basis time steps t_k, k=1, . . . , N_b, is a subset of the entire time steps t_j, j=1, . . . , N_t. Equation (28) holds because the solution (i.e. PFLOTRAN simulated fluid conductivity) of the 1D advection-diffusion equation is the superimposition of a set of error functions under transient boundary conditions. Using the linear relationship between the fluid resistivity and bulk resistivity (reciprocal of bulk conductivity) according to Archie's Law, we can represent bulk resistivity time series at depth z using a similar time-interpolation

ρ_b(z,t_j)=Σ_k=1^N^bc_jkρ_w(z,t_k),j=1, . . . ,N_t, (29)

with bulk resistivity p_b, and N_b, C_Jk, and t_kthe same as those in Equation (28).

As the fluid resistivity contrast between the GW and SW is usually small (i.e. less than one order of magnitude), the related bulk resistivity varies in a small range over time. For this small range of temporal variations in bulk resistivity, the transfer resistance time series at depth z is approximately a linear function with respect to the bulk resistivity time series at the same depth. We further reach the following time-interpolation of transfer resistance time series

R(z,t_j)=Σ_k=1^N^bc_jkρ_w(z,t_k),j=1, . . . ,N_t, (30)

where R(z,t_k)=A_e(σ_b(z,t_k)) denote E4D simulated transfer resistance at depth z and basis time steps t_k. Once N_b, t_k, and c_jkare known, we can first run E4D to compute R(z,t_k), and then interpolate R(z,t_k) to R(z,t_j) using Eq. (30).

Given the PFLOTRAN-E4D simulated bulk conductivity on E4D mesh over the entire time steps, we use Equation (29) to find N_b, t_k, and c_jk. The overall procedure is summarized as follows. (1) Generate an auxiliary 1D mesh z. This is done by discretizing the spatial extent of Wenner-alpha arrays into a uniform grid of N_zcells, with nodes at the electrode positions. The spacing of the auxiliary 1D mesh therefore represents the vertical resolution of the ERT survey. (2) Build bulk resistivity matrix ρ_b=[ . . . , ρ_b(z,t_j), . . . ]. We interpolate the bulk resistivity from the E4D mesh to the auxiliary 1D mesh at every time step, and arrange the interpolated bulk resistivity into a matrix of the size N_z×N_t, with N_z<<N_t. (3) Find basis time steps t_k. QR decomposition with pivoting can be used to find the basis columns of the bulk resistivity matrix. We take N_b=N_zand the basis time steps t_kfrom the indices of basis columns. (4) Calculate interpolation coefficients c_jk. This is done by plugging N_band t_kinto Eq. (23) and solving for c_jk.

Using the above time-interpolation process, the time-lapse ERT simulation needs N_btimes of E4D runs rather than N_t, where N_b=N_zis close to the number of spatial transfer resistance measurements. The time-interpolation has enough accuracy when the fluid resistivity contrast between the GW and SW is small (i.e. within one order of magnitude).

In accordance with example implementations, the system can be configured to assess vertical SW/GW transport under conditions of transient flow and heterogeneous sediments (e.g. high-frequency stage fluctuations within the hyporheic zone). The sensing assembly, installed vertically at the SW/GW interface, collects bulk electrical resistance (i.e. transfer resistance) through time at multiple depths. It also records corresponding time series of vertical pressure gradient, streambed temperature, and fluid conductivity. Under vertical flow conditions, analysis of sensing assembly data as described can provide: (1) vertical component of pore velocity, (2) vertical distribution of porosity, (3) vertical Darcy velocity, and (4) vertical distribution of hydraulic conductivity (or the relevant permeability).

Accordingly, a joint inversion method to estimate porosity, permeability, pore velocity and Darcy velocity is provided. The joint inversion can exploit all available data sets towards an integrated hydrogeological model and provide dynamic SW/GW interactions.

Referring to FIG. 11, layered hydrogeologic model (a) and vertical ERT inverted bulk conductivity profile (b) is shown. The model depicts of five layers between z=−1 to 0 m, with orange squares the depth of fluid conductivity and temperature sensors, and blue diagonal crosses the center of Wenner-alpha arrays. The vertical bulk conductivity profile was conservatively interpreted as three layers between z=−0.85 and −0.1 m.

Referring next to FIG. 12, EnKF results of hydrogeologic modelling (a and b), Darcy velocity (c), temperature (d), and fluid conductivity (f), with ensemble size N_s=30 and a relative data misfit of 7.7% are shown. CG result is similar except variations in permeability.

Referring next to FIG. 13, transfer resistance data fit at different depths using the EnKF algorithm. The largest misfit occurs at the depth between z=−0.36 and −0.32 m, which coincides with that of the estimated porosity. CG result is similar.

TABLE 3

CG and EnKF recovered models.

CG
EnKF

algorithm
algorithm
True model

Data misfit
7.66%
7.72%
—

Porosity
top
15.30%
15.58%
15%

layer

mid
43.78%
42.88%
b/t 25%-45%

layer

bot
20.43%
20.19%
20%

layer

Permea-
top
2.60 × 10⁻¹¹
1.45 × 10⁻¹⁰
1.18 × 10⁻¹⁰

bility
layer

[m²]
mid
5.84 × 10⁻¹¹
4.95 × 10⁻¹²
b/t 1.18 × 10⁻¹⁰-

layer

5.91 × 10⁻¹⁰

bot
1.13 × 10⁻¹¹
2.56 × 10⁻¹⁰
1.18 × 10⁻¹⁰

layer

Average
1.88 × 10⁻¹¹
1.84 × 10⁻¹¹
1.92 × 10⁻¹¹

permeability

[m²]

SW/GW interactions were simulated in the laboratory using a 1.5 m tall and 0.6 m diameter test cell. The test cell consisted of five sand layers with an accumulated thickness of 1.2 m, as shown in FIG. 14. Exterior (a) and interior (b) of the test cell. The test cell was 1.5 m high and 0.6 m diameter, including five sand layers with a total thickness of 1.2 m. The sensor assembly and system collects time series of vertical pressure gradient, temperature, fluid conductivity and transfer resistance between z=−0.725 and −0.125 m at the positions indicated. (ER: electrical transfer resistance; EC: electrical fluid conductivity; TM: temperature)

Tap water of different temperature and salinity were used for the surface water and groundwater. The SW reservoir had a temperature of ˜10° C. and a conductivity of ˜0.011 S/m, while the GW reservoir had a temperature of ˜20° C. and a conductivity of ˜0.045 S/m. FIG. 15 illustrates the setup of the flow experiment configured for downward flow. During the experiment, surface water is circulating through a chilled reservoir to maintain a temperature of ˜10° C. The test cell and the surface water reservoir were coated with foam insulation to minimize the room temperature drift in the measurements.

The sensor assembly was vertically installed across the streambed down to a depth of 0.85 m. Within this depth range, three layers of two types of sand (LM30 and LM70) were encountered. There are 16 ring electrodes between z=−0.8 and −0.05 m, with a vertical spacing of 0.05 m, for the time-lapse ERT survey using Wenner-alpha configuration. Time series of temperature and fluid conductivity are collected at 5 depths between z=−0.725 and −0.125 m, with a vertical spacing of 0.15 m. Besides, the flux tool records pressure difference between z=−0.725 and −0.125 m over time. The top and bottom temperature and fluid conductivity, together with the vertical pressure gradient data, were used to set the transient boundary conditions in the simulation. The rest of data were jointly inverted to estimate the hydrogeologic parameters of the sand column.

The LM30 and LM70 samples (in 5 cm scale) had hydraulic conductivity of 0.148 and 0.0117 cm/s, respectively, which were obtained through constant head tests. With flux, pressure gradient, and length of sample column measured from the constant head tests, hydraulic conductivity was calculated using Darcy's Law. Upscaling of the sample measurements gave an average hydraulic conductivity of 0.0514 cm/s (equivalent to permeability of 5.25×10⁻¹¹m²) for the 0.6 m height sand column next to the flux tool sensors, based on the geometric mean K_G=K_LM30^35/60K_LM70^25/60.

Laboratory measured porosity φ of different layers were given in FIG. 14. These values were calculated by φ=1−ρ_bulk/ρ_particle, with the saturated bulk density ρ_bulkmeasured during the packing, and the particle density ρ_particleassumed to be 2.65 g/cm³. Porosity of the top three layers (where the sensor assembly was located) were similar and had an average value of 42% for the sediments between z=−0.725 and −0.125 m, based on the arithmetic mean

$φ_{A} = (\frac{22.5}{60}) φ_{1} + (\frac{25}{60}) φ_{2} + (\frac{12.5}{60}) φ_{3},$

with φ₁=40%, φ₂=44%, and φ₃=43%.

Experiment at variable flow rate between −0.06 m/hr and 0.06 m/hr (positive for upward flow), with observed data plotted in FIG. 16. The experiment was (1) started at a downward flow rate of ˜0.06 m/hr for 16 hours, (2) switched to a slower downward flow rate of ˜0.03 m/hr for 1 hours, (3) changed to zero flow rate for 1 hour, (4) restarted at an upward flow rate of ˜0.03 m/hr for 2 hours, and (5) switched to a faster upward flow rate of ˜0.06 m/hr for 20 hours. A total of forty hours of flux tool data were collected at a time-sampling rate of 1 min.

The bulk conductivity profile obtained from the ERT inversion suggested one layer between z=−0.725 and −0.125 m, which agreed with the fact that the measured porosity within this depth range did not show obvious layering. Time series of temperature, fluid conductivity, and transfer resistance were then jointly inverted for a homogeneous hydrogeologic model between z=−0.725 to −0.125 m. The joint inversion recovered an average permeability 6.60×10⁻¹¹m²and porosity 42.08%, with a relative data misfit of 9.4%. The system data fit and predicted pore velocity are shown in FIG. 17. Compared to the measured values, the joint inversion recovered an exact average porosity but higher permeability (25.71% higher), indicating fast flowpaths in the measurements. The fast flowpaths are possibly caused by fluid passing through the fluid conductivity sensor ports. This is a design flaw that will be addressed in the next generation of the tool design.

The sampling volumes of the fluid conductivity and bulk conductivity sensors are different. The pore velocity related to the breakthrough of the fluid conductivity time series represents a local fluid transport, which may be affected by fast flowpaths. In contrast, the breakthrough of the transfer resistance time series characterizes the fluid seepage within a larger volume, which is dominated by flowpaths in the saturated sediments. The forward modeling in the joint inversion does not account for fast flowpaths through the fluid conductivity sensor ports. This may lead to a predicted pore velocity reflecting local fast flowpaths, and thus may result in biased estimation of permeability and Darcy velocity. For example, a better data fit of the fluid conductivity than transfer resistance time series in terms of breakthrough arrival time was found in FIG. 17b and FIG. 17c. This coincided with larger estimated pore velocity and permeability mentioned above.

An additional inversion for a more accurate estimation of permeability was performed. With porosity fixed as 42%, we converted the measured transfer resistance data to the bulk fluid conductivity and used the “converted bulk conductivity” as the boundary conditions in the PFLOTRAN simulation. Temperature and fluid conductivity were not inverted in this test. The inversion recovered an average permeability of 5.60×10⁻¹¹m², with a relative data misfit of 8.6%. The predicted transfer resistance and pore velocity can be found in FIG. 18. This time, the estimated permeability and pore velocity are close to the actual values, and the transfer resistance data fit in terms of the breakthrough arrival time were improved. This inversion can be used as a supplementary step of the joint inversion if the fluid conductivity measurements introduce significant local fast flowpaths (i.e., incompatible data fit of fluid conductivity and transfer resistance breakthrough curves).

It is possible that there is no SW/GW breakthrough in the fluid conductivity and transfer resistance time series, e.g., for constantly upwelling or downwelling streams. For this scenario, the fluid conductivity and transfer resistance data contain only subsurface porosity information. Bulk conductivity obtained from ERT inversion can be correlated with fluid conductivity measurements to calculate porosity through Archie's Law. With porosity estimated and fixed in the hydrogeologic model, the permeability and Darcy velocity can then be derived using the temperature data collected by the embodiments of the systems of the present disclosure. The temperature data can be inverted for permeability, with porosity fixed as 42%. The inversion recovered an average permeability of 6.11×10⁻¹¹m², with a relative data misfit of 4.1%. The temperature data fit and predicted pore velocity are given in FIG. 19.

The systems and method of the present disclosure were laboratory tested in a 152 cm tall, 60 cm in diameter insulated soil column. Calculations showed that a large diameter column was necessary to minimize boundary effects on the electrical field produced by the transfer resistance measurement system. In accordance with another experiment, the layering packed around the sensing assembly can be seen in FIG. 20. These layers were packed in 2.5 cm lifts. The bottom two layers of sand labelled 12/20A and LM70 were designed to distribute the flow of water evenly in the column by having a high hydraulically conductive zone, 12/20A, on the bottom and a low conductivity zone, LM70, above it to provide back pressure and distribute flow. The 3rd and 5th layers from the bottom were packed with sand labelled LM30. The fourth layer from the bottom was packed with LM70 sand. A diagram of the probe and layering of each sand type in this test can be seen in FIG. 20.

Three Environmental Tracer experiments were performed using three flow rates (50, 100, and 150 mL/min). For all three of these experiments GW was pumped from the bottom of the soil column into the GW reservoir, pulling SW down through the soil column. This was done until a full breakthrough of SW was observed. Once a full breakthrough was observed the flowrate was reversed at the same rate by pumping GW from the GW reservoir back into the bottom of the soil column.

The 100 mL/min Environmental Tracer experiment was repeated after the fluid conductivity sensor inlets and outlets were modified so that they were separated by 4 cm instead of 10 cm. The only difference between this experiment and the previous 100 mL/min flow rate experiment was the flow was incrementally changed from 100 mL/min downward, to 50 mL/min downward, to zero flow, to 50 mL/min upward and then to 100 mL/min upward where it remained until a full breakthrough of the GW was achieved. Each of the 4 steps was about an hour except there were 2 hours between the 50 and the 100 mL/min upward instead of 1 hour.

For these experiments the temperature corrected bulk resistance measurements were calibrated from bulk conductivity data for the same sand as packed at a similar bulk density and porosity as the sand packed in the soil column. The calibrated bulk conductivity and fluid conductivity peaks of these experiments were used to estimate the porosity.

The temperature corrected transfer resistance measurements were converted to normalized breakthrough curves for the location between the middle electrodes of each 4 electrode resistance measurement. The bulk resistance measured was assumed to increase during the time that the SW front moves between these two electrodes for each of the 13 measurement locations. Therefore, the center of the interval was chosen to extract SW travel times. The breakthrough time was defined as the time at which 0.5 times the total breakthrough magnitude was observed. There may be some error in using the transfer resistance due to changes in the current flow path associated with changes in fluid conductivity within the pore water. However, this should be minimized in this experiment by the similar geometry and porosity within the soil column at various depths canceling out these errors for times between breakthrough curves.

The thermal breakthrough curves of the repeated 100 mL/min experiment were also compared to PFLOTRAN simulated breakthrough curves using the theoretical pore water velocities within the soil column during Experiment 4. The observed and predicted values were evaluated using the coefficient of determination and the root mean square error between the observed and the predicted data.

A Steady-State Thermal Experiment was performed where the SW and GW had the same electrical conductivity but the temperature of the SW remained at 10° C. and GW at room temperature (approx. 21° C.). For this experiment, an upward flow was initiated until thermal equilibrium was reached than the flowrate was increased until steady-state temperature values were collected at 100 and 150 mL/min.

The measured pore water velocities of Experiments 1-3 were all faster than the theoretical values calculated for plug flow through a homogenous porous media. Also, fluid conductivity breakthrough curves had significantly faster pore water velocities than the bulk resistance measurements, indicating that there was a fast flow path through these sensors that had less of an effect on the bulk measurements. The thermistors, being co-located with these fluid conductivity sensors, also may have been affected by this fast flow path. Once a fast flow path was realized, the fluid conductivity sensor port locations were moved from 10 cm apart to 4 cm apart for Experiment 4, which had the same flow rate as Experiment 2. Moving these ports closer together reduces the amount of soil that water passing through the sensor can potentially bypass. The resulting data for Experiment 4, showed measured pore water velocities much closer to the theoretical for both the bulk resistance and fluid conductivity derived values, but even though there is a large improvement in the fluid conductivity derived pore water velocity, a fast flow path is still observed through this sensor. However, with the fast flow path greatly reduced in Experiment 4, the transfer resistance derived pore water velocities are not significantly different than the theoretical. Thus, more weight should be given to pore water velocity values derived using bulk resistance measurements than those derived using fluid conductivity measurements unless the fluid conductivity sensors are improved. Improvement have already been made in the next version which minimize excess flow through these sensors while still providing enough flow to quickly measure changes in fluid conductivity.

The pressure transducers in all of the 4 experiments appeared to have considerable noise. This noise is thought to be from the AC power supply and would most likely be greatly reduced during field deployment where the system will run off of battery and a solar panel. This noise in the pressure transducers measurements was observed to increase when other instruments were running on the same circuit. For the purpose of this study, the pressure measurements were assumed to be constant at a given flowrate so the pressure at each flowrate was averaged. Despite the AC noise, the hydraulic conductivity values from the differential pressure transducer were reasonable and repeatable. This AC noise should not be present when installed in the field as long as clean DC power is used.

TABLE 4

Pore water velocity from breakthrough times of the transfer

resistance measurements and the fluid conductivity measurements.

Exp. 1
V_p_avg (m/hr)
q (m/hr)
% Error

Theoretical
0.0277
0.0117
—

Measured using Transfer Resistance
0.0328
0.0138
18.2%

Measured using fluid conductivity
0.0438
0.0184
47.3%

Exp. 2

Theoretical
0.0555
0.0233
—

Measured using Transfer Resistance
0.0652
0.0274
17.4%

Measured using fluid conductivity
0.0842
0.0353
42.1%

Exp. 3

Theoretical
0.0832
0.0349
—

Measured using Transfer Resistance
0.1023
0.0430
22.9%

Measured using fluid conductivity
0.1192
0.0500
39.7%

Exp. 4

Theoretical
0.0555
0.0233
—

Measured using Transfer Resistance
0.0566
0.0238
1.9%

Measured using fluid conductivity
0.0701
0.0294
25.3%

Referring to FIG. 21, experimental data collected during Experiment 1 with a flowrate of 50 mL/min is provided where (a.) is the differential pressure between a depth of 12.5 and 72.5 cm below the top of the sediment, (b.) is the transfer resistance measurements of the 13 measurements, (c.) are measurements form temperature sensors TM(1), TM(2), TM(3), TM(4), and TM(5), and (d.) are measurements from fluid conductivity sensors EC(1), EC(2), EC(3), EC(4) and EC(5), which are located at depths of 12.5, 27.5, 42.5, 57.5, and 72.5 cm below the soil surface, respectively.

Referring to FIG. 22, experimental data collected during Experiment 2 with a flowrate of 100 mL/min is provided where (a.) is the differential pressure between a depth of 12.5 and 72.5 cm below the top of the sediment, (b.) is the bulk resistance measurements of the 13 measurements, (c.) are measurements form temperature sensors TM(1), TM(2), TM(3), TM(4), and TM(5), and (d.) are measurements from fluid conductivity sensors EC(1), EC(2), EC(3), EC(4) and EC(5), which are located at depths of 12.5, 27.5, 42.5, 57.5, and 72.5 cm below the soil surface, respectively.

Referring to FIG. 22, experimental data collected during Experiment 3 with a flowrate of 150 mL/min is provided where (a.) is the differential pressure between a depth of 12.5 and 72.5 cm below the top of the sediment, (b.) is the bulk resistance measurements of the 13 measurements, (c.) are measurements form temperature sensors TM(1), TM(2), TM(3), TM(4), and TM(5), and (d.) are measurements from fluid conductivity sensors EC(1), EC(2), EC(3), EC(4) and EC(5), which are located at depths of 12.5, 27.5, 42.5, 57.5, and 72.5 cm below the soil surface, respectively.

Referring to FIG. 24, experimental data collected during Experiment 4 which was a repeat of Experiment 2 after the fluid conductivity sensors were modified and insulation was installed on the outside of the soil column is provided. (a.) is the differential pressure between a depth of 12.5 and 72.5 cm below the top of the sediment, (b.) is the bulk resistance measurements of the 13 measurements, (c.) are measurements form temperature sensors TM(1), TM(2), TM(3), TM(4), and TM(5), and (d.) are measurements from fluid conductivity sensors EC(1), EC(2), EC(3), EC(4) and EC(5), which are located at depths of 12.5, 27.5, 42.5, 57.5, and 72.5 cm below the soil surface, respectively.

In FIGS. 21-24, GW is initially pumped out of the bottom of the test column pulling SW down past the TEFLUX. This is why the pressure starts out negative in FIGS. 21a-23a. Once a full breakthrough of SW was observed, the flow is reversed to pump GW back into the bottom of the column until a full breakthrough of GW is observed. In FIG. 24a the flow is changed in stages to better simulate flow changing directions in a steam bed. In FIGS. 21b-23b the transfer resistance breakthrough curves are closer together than expected because of an observed fast flow path through the fluid conductivity sensors, which was corrected in 24b. This same fast flow path is even more pronounce in FIG. 21d-23d, but was greatly improved in FIG. 24d due to the modifications to the system.

The porosity was estimated using Archie's law parameters where the slope of the fluid electrical conductivity and bulk electrical conductivity were used to calculate the formation factor. The results of these calculations can be seen in Table 5. It is also important to note that the fast flow path in Experiments 1-3 did not negatively affect these porosity measurements because they were only calculated from data where just the SW and GW are observed, and not during breakthroughs of GW or SW. The cementation exponent was assumed to be 1.4 for sand grains close to spherical. The estimated values of porosity based on measured parameters range from −14 to 15 percent of the values calculated from the packing parameters. However, most of the system derived porosity values have absolute errors below 10%, with the exception of the upper most measurement (12.5 cm depth). The upper most measurements may be skewed to a higher porosity because of their proximity to the soil surface. Also between Experiment 3 and 4 the experimental system was removed from the test column for modifications and reinstalled, which may have disturbed some of the packed sediment, thus potentially increasing the error in porosity during experiment 4.

TABLE 5

Porosity estimations using transfer resistance measurements inverted for bulk conductivity

and the measured fluid conductivity with an assumed cementation exponent of 1.4.

Porosity values calculated from packing parameters were calculated for approximately

a 7.5 cm of sand centered on the fluid conductivity locations.

Porosity

Depth
Calculated Porosity
from

(cm)
Exp. 1
% Err
Exp. 2
% Err
Exp. 3
% Err
Exp. 4
% Err
Packing

12.5
0.431
8%
0.453
15%
0.442
11%
0.443
12%
0.394

27.5
0.429
3%
0.399
−5%
0.411
−2%
0.410
−2%
0.411

42.5
0.392
−9%
0.406
−6%
0.419
−2%
0.370
−14%
0.421

57.5
0.407
−12%
0.424
−8%
0.459
2%
0.427
−7%
0.451

72.5
0.409
−4%
0.416
−2%
0.425
1%
0.391
−8%
0.416

Average
0.414
−3%
0.420
−1%
0.431
2%
0.408
−4%
0.419

The saturated hydraulic conductivity values calculated for each experiment using the theoretical and the measured differential pressure between the upper and lower pressure ports can be seen in Table 6. These values seem reasonable considering that LM70 and LM30 sand had laboratory measured hydraulic conductivity values of 0.0117 and 0.148 cm s⁻¹, respectively. Using the geometric mean for the two above values, the calculated effective vertical hydraulic conductivity was calculated to be 0.0514 cm s⁻¹, which appears be close to the measured hydraulic conductivity values seen in Table 6.

TABLE 6

Measured Hydraulic Conductivity for Exp. 1-4

Experiment
Measured Ksat (cm/s)

Exp. 1
0.0444

Exp. 2
0.0440

Exp. 3
0.0442

Exp. 4
0.0480

PFLOTRAN was used to predict the temperature breakthrough curves for Experiment 4. Parameters used to predict temperature data include porosity and the thermal boundary conditions. The observed and predicted temperatures at 4 thermistor locations [TM(1), TM(2), TM(3), and TM(4)] (seen in FIG. 25) were evaluated using the coefficient of determination and the RMSE between the predicted and observed values. The coefficient of determination was 0.947, 0.966, 0.982, and 0.967 and the RMSE was 0.0908, 0.0244, 0.0015, and 0.0037, respectively. Referring to FIG. 25, observed data from Exp. 4 (solid) and predicted data using PFLOTRAN (dashed). TM(1) is the upper most thermistor and TM(4) is the second thermistor from the bottom are provided. A supplemental thermistor at the top of the sand and the lowest thermistor (TM(5)) were used by PFLOTRAN as boundary conditions for the predicted data. As can be seen in FIG. 25, the observed temperature is generally higher than the predicted for all 4 locations which may be from insufficient insulation around the soil column, but also could have been caused by slight differences between the model parameters and the actual. However, the observed and predicted data appeared to have a good correlation for all 4 thermistor locations indicating that the pore water velocity used in the predicted values must have been close to the actual pore water velocity.

The calculated pore water velocity and Darcy's flux can be seen in Table 7. These values were calculated by fitting a K_fsof 3.28 W m-1 K-1 to the data collected and is within the range ranges (2.5-3.4 W m-1 K-1) described by Midttomme and Roaldset (1998) for the saturated sand used these experiments.

TABLE 7

Resulting pore water velocities and Darcy’s flux values calculated

using the TEFLUX. These values were calculated using the Turcotte and

Schubert (1982) model as discuss in Schmidt et al. (2007).

Flowrate @ 100 mL/min
Flowrate @ 150 mL/min

V_p(m/hr)
q (m/hr)
V_p(m/hr)
q (m/hr)

Measured
0.0564
0.0231
0.0834
0.0333

Theoretical
0.0555
0.0228
0.0832
0.0333

As has been shown, porosity can be calculated using measured Archie's law parameters and can be used to calculate vertical flux from pore water velocity. Most of the porosity values calculated had errors less than 10%, except the upper most porosity values. It is thought that the upper most values may be skewed to a higher porosity because of their proximity to the sand surface boundary. Unless significant laying is observed only one porosity value will be assigned by averaging the measured values. When all 5 values are averaged, the porosity error goes down to less than 5%.

FIG. 26 shows example laboratory scale testing data collected the systems and methods of the present disclosure. FIG. 26A shows temperature time-series on each thermistor, FIG. 26B shows the fluid specific conductance time-series, FIG. 26C shows the bulk electrical resistance time-series, and FIG. 26D shows the pressure differential between upper and lower boundary pressure sensor.

To simulate the flux tool survey at the groundwater/surface water interface, we consider a 1D hydrogeological model with vertically varying permeability, k(z), and porosity, φ(z). The 1 m deep synthetic model consists of five layers, with permeability ranging from 1.2×10⁻¹²to 6×10⁻¹⁰m², and porosity from 0.15 to 0.45, as shown in FIG. 27. A survey configuration using Wenner Alpha array is deployed for the system ERT measurement, whereby the two potential measurement ring electrodes are placed in between the two current ring electrodes, with an equal spacing between any two adjacent ring electrodes. There are a total of 13 Wenner Alpha arrays in the synthetic model measuring time series of transfer resistance, R(z,t), where the array centers are evenly distributed between the depths of 0.175 m and 0.775 m with a vertical spacing of 0.05 m. In addition to the ERT measurement, there are 3 pairs of hydrological sensors measuring time series of temperature, T(z,t), and fluid conductivity, σ_w(z,t), at depths of 0.275 m, 0.425 m, and 0.575 m. The entire simulation spans 2.5 days with a time interval of 15 min. Finally, we assume the system records the boundary as well as initial conditions of the pressure, temperature, and fluid conductivity at the top (z=0) and bottom (z=−1 m) of the modeling domain.

The simulated time series of temperature, fluid conductivity, and transfer resistance, along with the boundary and initial conditions, are used as the input data of the joint inversion. The E4D inversion can be used to obtain the 1D layered mesh, and found three bulk layers between the depths of 0.1 m and 0.8 m, except the thin layer between the depths of 0.3 m and 0.35 m and the bottom layer below the depth of 0.9 m, due to the limit of vertical resolution and survey coverage. The next joint inversion of the hydrological and ERT data estimated permeability and porosity of the three bulk layers. Both the CG and EnKF algorithms produced good data fittings between the predicted and synthetic temperature and fluid conductivity, and accurate estimations of the Darcy velocity, as shown in FIG. 28.

Due to the missing thin layer, the recovered porosity model between the depths of 0.3 m and 0.6 m has certain offset away from the true value. The error in the model space then propagates into the data space, resulting a static shift between the predicted and observed data at corresponding depths throughout the entire simulation time, as shown in FIG. 28. As the center of the Wenner Alpha array moving away from the thin layer and its neighbors, e.g., above the depth of 0.3 m or below the depth of 0.6 m, the static shift goes away and we again obtain good data fitting between the predicted and synthetic transfer resistance.

Unlike the porosity, the missing thin layer in the recovered permeability model does not affect the data fitting. Indeed, the survey design of the sensor rod allows an estimate of the average permeability instead of the permeability variations of the 1D column. Comparison of the CG and EnKF inversions in FIG. 27 indicates that the recovered permeability of the layers is non-unique. However, the average permeability of the two inversion results are similar and close to that of the true model. This can explain why both inversions can fit the overall data at the same level. To confirm the average permeability hypothesis, we ran another test of the joint inversion by using three layers for porosity but one layer for permeability in the 1D mesh. This inversion produced the same data fitting as before, and the recovered permeability of the single bulk layer was close to the average permeability of the true model.

	Number	Date	Country
Parent	16438312	Jun 2019	US
Child	17988704		US

Systems and Methods for Determining Ground Water-Surface Water Interactions

Information

Publication Number

Date Filed

Date Published

Inventors

Original Assignees

CPC

International Classifications

Abstract

Description

Claims

CROSS REFERENCE TO RELATED APPLICATION

STATEMENT AS TO RIGHTS TO INVENTIONS MADE UNDER FEDERALLY-SPONSORED RESEARCH AND DEVELOPMENT

Provisional Applications (1)

Continuations (1)