A New Approach to Quantify Shallow Water Hydrologic Exchanges in a Large Regulated River Reach

Hydrologic exchange is a crucial component of the water cycle. The strength of the exchange directly affects the biogeochemical and ecological processes that occur in the hyporheic zone and aquifer from micro to reach scales. Hydrologic exchange fluxes (HEFs) can be quantified using many field measurement approaches, however, in a relatively large river (scale > 103 m), these approaches are limited by site accessibility, the difficulty of performing representative sampling, and the complexity of geomorphologic features and subsurface properties. In rivers regulated by hydroelectric dams, quantifying HEF rates becomes more challenging because of frequent hydropeaking events, featuring hourly to daily variations in flow and river stages created by dam operations. In this study, we developed and validated a new approach based on field measurements to estimate shallow water HEF rates across the river bed along the shoreline of the Columbia River, USA. Vertical thermal profiles measured by self-recording thermistors were combined with time series of hydraulic gradients derived from river stages and inland water levels to estimate the HEF rates. The results suggest that the HEF rates had high spatial and temporal heterogeneities over the riverbed, with predicted flux rates varied from +1 × 10−6 m s−1 to −1.5 × 10−6 m s−1 under different flow conditions.


Hydrologic Exchange and River Regulations
Hydrologic exchange is a concept introduced by Harvey and Gooseff [1] which combines surface water-groundwater interaction processes along river corridors at multiple spatiotemporal scales, including hyporheic exchange, bank storage, and regional groundwater discharge and recharge.River water interacts with subsurface water through the hydrologic exchange fluxes (HEFs), which facilitate the nutrient and carbon cycling, organic biodegradation, fish spawning, metal transport, and other key biogeochemical and hydroecological processes in the subsurface region of the river corridor [2][3][4][5][6][7].
Hydrologic exchange dynamics, including the direction, path, magnitude, and residence time of the HEFs, define where and when these aforementioned processes occur in the river bed and bank [8][9][10].The HEF dynamics are controlled by channel geometry, catchment geology, and hydrology in both space and time [11][12][13], and the HEF rates are essentially governed by the permeability of the sediments and local hydraulic gradients.While the overall permeability of the river bed sediment is a physical property that is relatively stable (if not considering the sediment clogging/colmation at the river bed surface), the hydraulic gradient across the river bed is a time-varying variable that is highly dependent on the river stage variations.In natural rivers, the river stage reflects the hydrometeorological and drainage characteristics in the upstream watershed.However, in a dam-regulated river, variations in stage are dramatically altered by upstream dam operations on a wide spectrum of timescales from hourly to monthly [14], and these variations extend to a few hundred kilometers downstream of the dam [15].Globally, over 30,000 large dams were built in the past 50 years, which considerably affected terrestrial water resources [16][17][18] and thus changed the HEF directions and rates as well as the biogeochemical processes in downstream reaches [19].A number of studies have investigated HEFs in regulated river systems based on measurements or modeling at either the point [15] or transect scale [20][21][22].However, monitoring HEFs in a large regulated river reach (>10 3 m spatial scale) with high spatiotemporal resolutions remain less explored.

Approaches to Measuring HEF Rates
The challenges of monitoring the HEFs in a regulated large river reach arise from the size of the domain and dam-induced rapid hydraulic gradient variations that complicate the system dynamics.Current field measurement methods used in a large river system usually include mass balance approaches that integrate the groundwater discharge or river loss over the entire reach.These approaches include longitudinal flow gauging, longitudinal chemistry sampling, or hydraulic head measurements [23].Alternatively, interpolation over a number of intensively measured points (<1 m spatial scale) is a widely used approach to establish spatial distributions across the river bed for a given river reach.This type of measurement can be categorized into three groups: direct measurements, indirect measurements based on Darcy's Law, and indirect measurements based on heat transport equations [24].Direct measurements monitor the water volume flow across the surface water-groundwater interface using seepage meters, which are chambers with the bottom open to the sediment [25][26][27][28].This type of measurement can only estimate the flux at the interface, and the seepage meters themselves disturb the surrounding flow field and therefore affect the flux measurements.
The method based on Darcy's Law estimates HEF rates by solving Darcy's equation: where q is the flux, K is the hydraulic conductivity of the sediment, and dh/dl is the hydraulic gradient.K can be estimated from various experiments such as pumping tests, slug tests, and permeameter tests, while dh/dl can be derived from water level measurements in piezometers at different depths in the riverbed (for vertical flux) or at different locations at the same level beneath the riverbed (for horizontal flux) [29].The concept of the Darcy's Law-based method is simple and equipment installation is straightforward, which makes this method popular in some small-scale applications.However, uncertainties could be introduced in K estimations due to the heterogeneity and variations of the hydraulic conductivity; and the approach may not be suitable for large river systems because of its intensive labor requirements [24].Methods that use heat as a naturally occurring tracer to estimate HEFs have been developed, widely adopted, and improved since the 1960s [4].The underlying concept of such methods is that the temperature distribution in the subsurface zone is determined by the surface water temperature with strong diurnal fluctuations and relatively stable temperature in deep groundwater.Heat transports in the subsurface region through thermal convection and conduction processes, which Water 2017, 9, 703 3 of 18 modifies the distribution and dynamics of temperature following the 3D heat transport equation.The one-dimensional thermal transport equation along the vertical direction can be written as: where T is temperature; t is time; q z is the vertical flux; z is the depth in sediment; ρ w and ρ are the density of water and water-sediment matrix; c w and c are the specific heat of water and river-sediment matrix, respectively; and D is the thermal diffusion coefficient.A number of analytical solutions to the heat transport equation have been provided and applied successfully in specific cases by assuming the system is in steady state with vertical fluxes only, in a homogeneous semi-infinite sediment domain [30][31][32].Transient state HEF rate can be estimated by applying numerical modeling tools (e.g., VS2DH [33]) based on the heat transport equations, with thermal boundary conditions and initial states defined at the river bed and in deep sediment [34][35][36].Another method for transient state estimation is analyzing the amplitudes and phase shift of the damping diurnal temperature signals at different depths in the sediment [37][38][39][40]).Computational codes, such as Ex-Stream and VFLUX, were also developed to facilitate the application of these models [41,42].One of the great advantages of these approaches is that temperature is a reliable tracer and can be relatively easily measured at different depths using one single measuring rod installed in the river bed compared to a piezometer nest for the water head measurements.However, it might be a challenge to apply these approaches in a highly regulated river system, because the temperature signals might be significantly dampened or distorted by the fluxes with rapidly changing directions and magnitudes.

Study Objectives
In this study, we aimed to synthesize the HEF measuring methods reviewed above and develop an approach that is suitable for a large river with highly regulated discharge.Such an approach has to fulfill the following criteria: (1) using less, easily accessible data to infer HEF rates in a relatively large domain; and (2) providing continuous HEF rates at high temporal resolutions (e.g., daily or sub-daily).Note that the river bed in our study reach has small slopes at both longitudinal (1/4000) and lateral (1/40) directions, so that the major HEF direction that penetrates the river bed is vertical, which is used to represent the total HEF rate in our analysis.By revisiting monitoring data from previous studies in our study reach [15,43], we developed and validated our approach by empirically relating the point vertical HEF rates with river hydrologic conditions and inland water table levels.Such an approach could provide point HEF rates for large scale groundwater modeling evaluations [44,45], and provide guidance for the field studies to identify strong HEF hotspots and active biogeochemical reaction areas at the river bed [46].Although hydrologic exchanges occur almost everywhere at the river bed, numerical simulations and field measurements indicated that the exchange fluxes in the shallow water near the river banks were stronger than those in the center of the channel because of the greater pressure gradient between the inland water table and the river stage [21,47].Therefore, in this study, we only focus on the shallow water area near the bank along the river reach.
To summarize, the objectives of this study are two-fold: (1) to demonstrate a new approach that combines field measurement data and regression analysis to infer long term river bed HEF rates in a highly dynamic river reach; and (2) to examine the spatial and temporal distributions of HEFs in the shallow water area along the river reach.

Study Reach
Our study reach is a 5 km long, 800 m wide, nearly straight river segment of the Columbia River near the 300A Area of the U.S. Department of Energy Hanford Site in southeastern Washington State (Figure 1).The Columbia River flows from north to south through the study reach over an unconfined Water 2017, 9, 703 4 of 18 aquifer on top of the impermeable Columbia River Basalt Group.The total thickness of the aquifer ranges from 40 m to 60 m in this area, with a highly permeable Pleistocene flood gravel layer of the Hanford Formation on top and a more consolidated and less permeable Ringold Formation at the bottom.At the top of the river bed, there is a thin layer of alluvium with varying thicknesses ranging from one to three meters [48].Measurements used in this study were all collected from sensors installed in the alluvial layer.Several sandbar islands were deposited in the center of the river, forming a deep primary channel with a meandering thalweg on one side of the islands and a relatively shallower secondary channel on the opposite side of the islands (Figure 1).Priest Rapids Dam, a hydroelectric dam located ~80 km upstream of the study reach, controls the stage for hydropower generation.Based on historical records over the past 40 years, the range of stage spans 105 m to 109 m based on the North American Vertical Datum of 1988 (NAD88), with a maximum daily fluctuation of 2 m.In this study, the river stage observations were recorded by a pressure transducer logger (SWS-1) installed inside our study domain since the year 2001.Nearly 100 monitoring wells were constructed over the past 40 years around the 300A Area to record the inland water levels and solute concentrations for environmental monitoring purposes [49].In this study, one monitoring well (Well 2-3), about 150 m from the river, was selected as the reference point for the inland water level.Here we define "inland" as a contrast of the "river", where the water table does not show sub-daily variations in response to the dam regulations.According to the Koppen climate classification, our study area is located in a semi-arid-desert catchment in the Columbia River basin.The annual precipitation in this region is less than 200 mm, and the evapotranspiration to precipitation ratio is about 1-1.09 [50].Under such conditions, recharge to the groundwater is close to ~2 mm year −1 [51] and precipitation over the river surface can be ignored compared to the volume of water in the river channel.Therefore, the fluctuations observed in the river stage time series are mainly caused by the dam operations.
Water 2017, 9, 703 4 of 17 bottom.At the top of the river bed, there is a thin layer of alluvium with varying thicknesses ranging from one to three meters [48].Measurements used in this study were all collected from sensors installed in the alluvial layer.Several sandbar islands were deposited in the center of the river, forming a deep primary channel with a meandering thalweg on one side of the islands and a relatively shallower secondary channel on the opposite side of the islands (Figure 1).Priest Rapids Dam, a hydroelectric dam located ~80 km upstream of the study reach, controls the stage for hydropower generation.Based on historical records over the past 40 years, the range of stage spans 105 m to 109 m based on the North American Vertical Datum of 1988 (NAD88), with a maximum daily fluctuation of 2 m.In this study, the river stage observations were recorded by a pressure transducer logger (SWS-1) installed inside our study domain since the year 2001.Nearly 100 monitoring wells were constructed over the past 40 years around the 300A Area to record the inland water levels and solute concentrations for environmental monitoring purposes [49].In this study, one monitoring well (Well 2-3), about 150 m from the river, was selected as the reference point for the inland water level.Here we define "inland" as a contrast of the "river", where the water table does not show sub-daily variations in response to the dam regulations.According to the Koppen climate classification, our study area is located in a semi-ariddesert catchment in the Columbia River basin.The annual precipitation in this region is less than 200 mm, and the evapotranspiration to precipitation ratio is about 1-1.09 [50].Under such conditions, recharge to the groundwater is close to ~2 mm year −1 [51] and precipitation over the river surface can be ignored compared to the volume of water in the river channel.Therefore, the fluctuations observed in the river stage time series are mainly caused by the dam operations.

Data Revisiting
Given that the pressure gradient is the major driving factor of the HEF, we hypothesize that the HEF rates could be empirically related to the head difference between the river stage and the inland

Data Revisiting
Given that the pressure gradient is the major driving factor of the HEF, we hypothesize that the HEF rates could be empirically related to the head difference between the river stage and the inland water table.We tested this hypothesis by revisiting a published data set [43].The data was Water 2017, 9, 703 5 of 18 collected at half-hourly time steps at Spring 9 (Figures 1 and 2), a site inside our study domain, from August 2004 to June 2005.It included temperature time series for surface water and a location 19 cm deep in the riverbed (Figure 3A) as well as a vertical flux time series at the riverbed (Figure 3B) derived using Darcy's Law from measured Vertical Hydraulic Gradient (VHG, by piezometers) and hydraulic conductivity estimated by slug tests at this location.We also extracted the inland water level from the monitoring well mentioned above and the corresponding river stage data for the same time frame (Figure 3C).Strong linear correlations were detected between head difference and the Darcy's Law-based vertical fluxes at both half-hourly (r = 0.971) and daily (r = 0.979) time steps (Figure 4).The scatter plots also showed that the fitted linear model based on the daily time series can be used to describe a similar relationship at the half-hourly time scale with high goodness-of-fit values.This indicates that the linear relation based on daily fluxes and the corresponding head difference can be used to estimate sub-daily fluxes if the continuous sub-daily head difference data are available.Further statistical analyses indicated that at the daily time scale, the fitted sum square of HEF rates based on head difference is 2.61 × 10 −10 , with a residual sum square of 1.28 × 10 −11 , and R 2 of 0.958.The corresponding F-test yields a nearly zero p-value of 2.58 × 10 −202 , rejecting the null hypothesis that the relationship is insignificant, which confirms the significant linear relationship.This finding inspired the method we developed for HEF rate estimations (described below).
Water 2017, 9, 703 5 of 17 water table.We tested this hypothesis by revisiting a published data set [43].The data was collected at half-hourly time steps at Spring 9 (Figures 1 and 2), a site inside our study domain, from August 2004 to June 2005.It included temperature time series for surface water and a location 19 cm deep in the riverbed (Figure 3A) as well as a vertical flux time series at the riverbed (Figure 3B) derived using Darcy's Law from measured Vertical Hydraulic Gradient (VHG, by piezometers) and hydraulic conductivity estimated by slug tests at this location.We also extracted the inland water level from the monitoring well mentioned above and the corresponding river stage data for the same time frame (Figure 3C).Strong linear correlations were detected between head difference and the Darcy's Lawbased vertical fluxes at both half-hourly (r = 0.971) and daily (r = 0.979) time steps (Figure 4).The scatter plots also showed that the fitted linear model based on the daily time series can be used to describe a similar relationship at the half-hourly time scale with high goodness-of-fit values.This indicates that the linear relation based on daily fluxes and the corresponding head difference can be used to estimate sub-daily fluxes if the continuous sub-daily head difference data are available.Further statistical analyses indicated that at the daily time scale, the fitted sum square of HEF rates based on head difference is 2.61 × 10 −10 , with a residual sum square of 1.28 × 10 −11 , and R 2 of 0.958.The corresponding F-test yields a nearly zero p-value of 2.58 × 10 −202 , rejecting the null hypothesis that the relationship is insignificant, which confirms the significant linear relationship.This finding inspired the method we developed for HEF rate estimations (described below).The insets of the plots showed the locations of the measurements in zoomed in 1:1 aspect ratio plots.
Water 2017, 9, 703 5 of 17 water table.We tested this hypothesis by revisiting a published data set [43].The data was collected at half-hourly time steps at Spring 9 (Figures 1 and 2), a site inside our study domain, from August 2004 to June 2005.It included temperature time series for surface water and a location 19 cm deep in the riverbed (Figure 3A) as well as a vertical flux time series at the riverbed (Figure 3B) derived using Darcy's Law from measured Vertical Hydraulic Gradient (VHG, by piezometers) and hydraulic conductivity estimated by slug tests at this location.We also extracted the inland water level from the monitoring well mentioned above and the corresponding river stage data for the same time frame (Figure 3C).Strong linear correlations were detected between head difference and the Darcy's Lawbased vertical fluxes at both half-hourly (r = 0.971) and daily (r = 0.979) time steps (Figure 4).The scatter plots also showed that the fitted linear model based on the daily time series can be used to describe a similar relationship at the half-hourly time scale with high goodness-of-fit values.This indicates that the linear relation based on daily fluxes and the corresponding head difference can be used to estimate sub-daily fluxes if the continuous sub-daily head difference data are available.Further statistical analyses indicated that at the daily time scale, the fitted sum square of HEF rates based on head difference is 2.61 × 10 −10 , with a residual sum square of 1.28 × 10 −11 , and R 2 of 0.958.The corresponding F-test yields a nearly zero p-value of 2.58 × 10 −202 , rejecting the null hypothesis that the relationship is insignificant, which confirms the significant linear relationship.This finding inspired the method we developed for HEF rate estimations (described below).

Using Temperature Records to Infer Vertical HEF Rates
The method we used to quantify vertical HEF rates is based on temperature time series measured in both surface water and in the riverbed.Rapid changes in river stage lead to quick variations in magnitude and frequent directional changes in the vertical fluxes across the river bed.Because no single method was capable of estimating vertical fluxes from such a complex system, we applied two different methods for time periods with different observed subsurface thermal characteristics.
The first method is a simple analytical solution of a 1D heat transfer equation described by [32], which assumes that the system is in a steady-state condition and that the direction of HEFs is constantly upward: where Ks and Kf are thermal conductivities of fluid and solid, respectively; n is porosity; ρfCf is fluid volumetric heat capacity; z is the depth of sensor; and T0, TZ, and TL represent the surface water temperature, the temperature at the sensor location, and groundwater temperature (assumed to be constant).This method is only applicable when there is constant upward flux.The second method is the Local Polynomial method with a Maximum Likelihood estimator (LPML) developed by Vandersteen et al. [39] and has been extensively tested with time-series data [52].The LPML model transforms temperature signals from the temporal domain to the frequency domain and finds the best vertical flux rate to fit equations that describe the subsurface signal frequency as a response function of the surface signal frequency.The 1D heat transport Equation ( 2) can be written as:

Using Temperature Records to Infer Vertical HEF Rates
The method we used to quantify vertical HEF rates is based on temperature time series measured in both surface water and in the riverbed.Rapid changes in river stage lead to quick variations in magnitude and frequent directional changes in the vertical fluxes across the river bed.Because no single method was capable of estimating vertical fluxes from such a complex system, we applied two different methods for time periods with different observed subsurface thermal characteristics.
The first method is a simple analytical solution of a 1D heat transfer equation described by [32], which assumes that the system is in a steady-state condition and that the direction of HEFs is constantly upward: where K s and K f are thermal conductivities of fluid and solid, respectively; n is porosity; ρ f C f is fluid volumetric heat capacity; z is the depth of sensor; and T 0 , T Z , and T L represent the surface water temperature, the temperature at the sensor location, and groundwater temperature (assumed to be constant).This method is only applicable when there is constant upward flux.
The second method is the Local Polynomial method with a Maximum Likelihood estimator (LPML) developed by Vandersteen et al. [39] and has been extensively tested with time-series data [52].The LPML model transforms temperature signals from the temporal domain to the frequency domain and finds the best vertical flux rate to fit equations that describe the subsurface signal frequency as a response function of the surface signal frequency.The 1D heat transport Equation (2) can be written as: Water 2017, 9, 703 7 of 18 where ρ and ρ w are densities of the fluid-sediment matrix and the fluid; c and c w are the specific heat of the fluid-sediment matrix and the fluid; D is the effective thermal diffusion coefficient that can be estimated from the bulk thermal conductivity (K) through D = K/ρc; and α, β, and γ are parameters in the frequency response function (FRF) that converts the surface temperature signal to the signal in the sediment at depth z in the frequency domain.These parameters can be fitted using a maximum likelihood estimator based on observations and therefore can be used to determine the vertical flux: The LPML model can incorporate the full or any portion of the temperature time series at daily to seasonal time steps.It can also deal with transient temperature signals by separating them into periodic, non-periodic, and additive noise portions using the local polynomial (LP) method.However, similar to other frequency-based methods such as those described by Hatch et al. [37] and Keery et al. [38], LPML works well when both surface and subsurface temperature time series share the same frequency, especially for the fundamental frequency range (e.g., diurnal cycles).This method, therefore, is more applicable when the subsurface temperature variations follow a pattern comparable, and have a fundamental frequency similar to that of the surface water temperature.
To identify periods suitable for the two methods, we employed the Dynamic Harmonic Regression (DHR) model [53] to decompose the observed surface water and subsurface time series into signals that represent the general trend and fundamental frequency.The DHR model has been successfully applied to detect the amplitude and phase of thermal signals in other studies (e.g., [38,42]).It treats observed temperature time series (y t ) as the sum of a general trend (T t ), a fundamental signal and its associated harmonics (C t ), and a white noise component (e t ): where a i,t and b i,t are time-varying parameters, ω 1 is the fundamental frequency, and ω i is the harmonics of ω 1 with ω i = i × ω 1 .The decomposition of the time series using the DHR model was conducted by applying the CAPTAIN toolbox using an auto-regression technique [54].
The decomposed general trend and amplitude of the fundamental frequency (i.e., diurnal) from surface and subsurface locations were used to first identify periods with clear upward groundwater fluxes when the analytical solution method is applicable.The surface water trend line follows a clear seasonal cycle that gradually changes between 2 • C and 24 • C within a year, while the deep groundwater temperature at this location is fixed at about 16 • C in our study reach.In summer and winter seasons, the temperature difference between surface and groundwater could reach up to 10 • C. As a result, each period with significant disparities between the two trend lines from the water and river bed temperatures indicates a groundwater discharge event and can be considered a clear upward flux period.Here we define "significant departure" as periods when the difference between surface and subsurface trend lines is more than three times the standard deviation of the fundamental (diurnal) signal.During these periods, the shallow subsurface temperature is highly affected by the upwelling of deep groundwater and the diurnal signal from the solar radiation is dampened or even disappeared (when the upward flux is very strong), which makes the frequency-based method impractical.Therefore, in these periods of strong upward fluxes, we assume that the system remains at a quasi-steady state at the daily scale and is suitable for the analytical solution (Equation ( 3)) for the daily flux estimations.
Given that the temperature at the bottom of the subsurface remains nearly constant and the diurnal temperature signal is forced by solar radiation heating the river water, the amplitude of the diurnal signal in the subsurface should be always less than that in the river water.However, the diurnal Water 2017, 9, 703 8 of 18 temperature amplitude of the subsurface decomposed by the DHR model could be overestimated at the edge of strong upwelling events because of the rapid temperature variations.Here, we simply compared the diurnal amplitudes from surface water and subsurface and discarded the data with greater subsurface amplitudes to reduce the uncertainties caused by the edging effects.The data that passed the filtering were then identified as inputs to the LPML model.Although data were screened at a daily time step, the LPML model was applied at the original 10-min interval over a three-day moving window for each day of interest.The workflow of data screening and model application are shown in Figure 5.
Water 2017, 9, 703 8 of 17 Given that the temperature at the bottom of the subsurface remains nearly constant and the diurnal temperature signal is forced by solar radiation heating the river water, the amplitude of the diurnal signal in the subsurface should be always less than that in the river water.However, the diurnal temperature amplitude of the subsurface decomposed by the DHR model could be overestimated at the edge of strong upwelling events because of the rapid temperature variations.Here, we simply compared the diurnal amplitudes from surface water and subsurface and discarded the data with greater subsurface amplitudes to reduce the uncertainties caused by the edging effects.The data that passed the filtering were then identified as inputs to the LPML model.Although data were screened at a daily time step, the LPML model was applied at the original 10-min interval over a three-day moving window for each day of interest.The workflow of data screening and model application are shown in Figure 5. Data screening clearly introduces discontinuities to the inferred riverbed HEFs at the daily time scale.However, the data points would be sufficient to establish the linear relations between the daily fluxes and the corresponding head difference between the river stage and the inland water table.As we demonstrated in Section 2.2, the same relationship can also be further used to estimate sub-daily fluxes at the riverbed.
We tested the approaches on the Spring 9 temperature data set (Figure 2) and compared the estimated fluxes with the published data based on Darcy's Law.In data filtering, out of the total of 293 days, 122 days were selected to apply the analytical method, 103 days were selected for to apply the LPML model, and 68 days were discarded.Figure 6 demonstrates how the periods were selected based on the criteria described above.To incorporate the uncertainties associated with sediment physical properties, we randomly generated 100 sets of parameter values within the ranges based on sediment properties provided by Ma et al. [55] (Table 1).The estimated daily vertical fluxes from the 100 realizations were plotted against the daily head difference to fit a linear model, and the range defined by the maximum and minimum slopes generally covered the fitted line derived from the observations (Figure 7).We did a linear regression test for these data points and the p-value of the Ftest was 1.22 × 10 −12 , indicating the linear relationship was significant at α = 0.05.Then we applied the 100 ensemble models to the head differences obtained from the river stage and inland water level Data screening clearly introduces discontinuities to the inferred riverbed HEFs at the daily time scale.However, the data points would be sufficient to establish the linear relations between the daily fluxes and the corresponding head difference between the river stage and the inland water table.As we demonstrated in Section 2.2, the same relationship can also be further used to estimate sub-daily fluxes at the riverbed.
We tested the approaches on the Spring 9 temperature data set (Figure 2) and compared the estimated fluxes with the published data based on Darcy's Law.In data filtering, out of the total of 293 days, 122 days were selected to apply the analytical method, 103 days were selected for to apply the LPML model, and 68 days were discarded.Figure 6 demonstrates how the periods were selected based on the criteria described above.To incorporate the uncertainties associated with sediment physical properties, we randomly generated 100 sets of parameter values within the ranges based on sediment properties provided by Ma et al. [55] (Table 1).The estimated daily vertical fluxes from the 100 realizations were plotted against the daily head difference to fit a linear model, and the range defined by the maximum and minimum slopes generally covered the fitted line derived from the observations (Figure 7).We did a linear regression test for these data points and the p-value of the F-test was 1.22 × 10 −12 , indicating the linear relationship was significant at α = 0.05.Then we applied the 100 ensemble models to the head differences obtained from the river stage and inland water level observations and compared the envelope defined by the 5th and 95th percentiles of the ensemble with the Darcy's Law-based flux (Figure 8).The comparisons revealed that the temperature-inferred vertical flux range adequately captures the direction and magnitude of the HEFs at both daily and half-hourly observations and compared the envelope defined by the 5th and 95th percentiles of the ensemble with the Darcy's Law-based flux (Figure 8).The comparisons revealed that the temperature-inferred vertical flux range adequately captures the direction and magnitude of the HEFs at both daily and half-hourly time scales.The same modeling framework was then used to estimate the vertical fluxes in the shallow water along the river, as introduced in Section 2.4.observations and compared the envelope defined by the 5th and 95th percentiles of the ensemble with the Darcy's Law-based flux (Figure 8).The comparisons revealed that the temperature-inferred vertical flux range adequately captures the direction and magnitude of the HEFs at both daily and half-hourly time scales.The same modeling framework was then used to estimate the vertical fluxes in the shallow water along the river, as introduced in Section 2.4.

Monitoring near Riverbed Temperature along the River
We deployed iButton ® temperature sensors (Maxim Integrated Products, Model number DS1922L) at five locations along the west bank of the river (Figure 1) to record temperature time series in the river water and the riverbed.The iButton ® sensor is a cylindrical, wireless device 17 mm in diameter and 6 mm thick.Its data storage capacity is 4096 values with 0.0625 °C measuring resolution and 0.5 °C accuracy.The elevations for all the monitoring locations were low enough to make sure the sensors would not be exposed to the air even in low-flow conditions.These locations were also selected based on the local hydrodynamic conditions and river bed properties and were labeled alphabetically from south to north.Sites A and B were located in the secondary channel that featured relatively low-flow velocities and small sediment sizes on the river bed, while sites D and E were located in the primary channel that featured relatively high flow velocities and large grain sizes in the alluvial sediment layer.Site C was in the transition zone between Sites B and D (Figure 2).
Temperature time series from the riverbed top and certain depth in the sediment were required for estimating the flux across the riverbed.Therefore, at each location, two sensors were secured in separate open drill holes along a 1 in.diameter, 1.5 m long solid plastic rod half planted in the river bed to record surface and subsurface temperature time series.All sensors were programmed to record every 10 min from 2 March to 30 March 2016.The coordinates and elevations of each sensor  We deployed iButton ® temperature sensors (Maxim Integrated Products, Model number DS1922L) at five locations along the west bank of the river (Figure 1) to record temperature time series in the river water and the riverbed.The iButton ® sensor is a cylindrical, wireless device 17 mm in diameter and 6 mm thick.Its data storage capacity is 4096 values with 0.0625 • C measuring resolution and 0.5 • C accuracy.The elevations for all the monitoring locations were low enough to make sure the sensors would not be exposed to the air even in low-flow conditions.These locations were also selected based on the local hydrodynamic conditions and river bed properties and were labeled alphabetically from south to north.Sites A and B were located in the secondary channel that featured relatively low-flow velocities and small sediment sizes on the river bed, while sites D and E were located in the primary channel that featured relatively high flow velocities and large grain sizes in the alluvial sediment layer.Site C was in the transition zone between Sites B and D (Figure 2).
Temperature time series from the riverbed top and certain depth in the sediment were required for estimating the flux across the riverbed.Therefore, at each location, two sensors were secured in separate open drill holes along a 1 in.diameter, 1.5 m long solid plastic rod half planted in the river Water 2017, 9, 703 11 of 18 bed to record surface and subsurface temperature time series.All sensors were programmed to record every 10 min from 2 March to 30 March 2016.The coordinates and elevations of each sensor are listed in Table 2.The method described in Section 2.3 was then applied to the monitored temperature records at each location to obtain the local HEF rates at the river bed.Note that here we used the surface water sensors to approximate the temperature at the top of the river bed and the vertical HEF fluxes were estimated based on the distances between the subsurface sensors to the river bed.

Results
The temperature data from the five iButton monitoring sites (Figure 9) showed that the river water temperature measurements were nearly identical across the sites, featuring a clear diurnal cycle over the period and a moderate raising trend with a daily mean temperature of about 5.5 • C on 2 March and 7 • C on 30 March, while the subsurface temperature differed among the sites.At Sites A and B, which are located in the secondary channel, the subsurface temperature had a similar general trend but a damped magnitude in diurnal variations compared to that in the surface water.At the sites located in the primary channel (D and E), strong variations in subsurface temperature and great departures (over 10 • C) from the surface temperature were evident, indicating strong upwelling events.Site C is located at the transition zone between the primary and the secondary channels, which has small subsurface temperature variations, with the magnitude between the primary channel and the secondary channel.Head differences at these sites were extrapolated from the reference river gauge observations at SWS-1 and inland well readings at Well 2-3 (Figure 10) based on the river stage and inland water table gradients established based on historical data.The river stage during the monitoring period ranged from 104.9 m to 105.9 m, which was quite representative because it covered about 80% of the full spectrum of the historical stage in the past 30 years that has a 5th percentile at 104.8 m and a 90th percentile at 106 m.For each site, we performed a linear regression analysis to test whether a significant linear relationship existed between the estimated HEF rate and the head difference.The p-values of the F-tests were 4.95 × 10 −14 , 6.65 × 10 −14 , 1.86 × 10 −13 , 4.29 × 10 −12 , and 5.89 × 10 −13 , for the sites A-E, respectively, indicating that the linear relationship were all significant.Then, the linear model for esimating flux rates as a function of head differences between the river stage and inland water level was established at each iButton location with 100 realizations.The corresponding time series of estimated HEF rate envelopes from the 5th and 95th percentile of the ensembles are shown in Figure 11.Generally speaking, both upwelling and downwelling events were observed at all sites and extreme values are shown in response to maximum and minimum flow stages.However, the estimated flux rates show strong spatial and temporal variability and differ significantly between the sites in the primary and secondary channels.At Sites A and B the fluxes were at magnitudes of about ±5 × 10 −7 m s −1 and ±4 × 10 −7 m s −1 , while at sites D and E, the values were nearly one magnitude greater (i.e., up to ±1.8 × 10 −6 m s −1 ).The flux ranges at D and E were found to be comparable to fluxes at Spring 9 (about 2 × 10 −6 to −4 × 10 −6 m s −1 ) derived using Darcy's Law.
Three flow conditions from the iButton monitoring period were selected for further comparison.The three conditions represented the high flow on 10 March, median flow on 28 March, and low flow on 8 March (Figure 12).We connected the median value of the predicted flux rates across the five monitoring locations under different flow conditions and found that at high and median flow conditions most of the fluxes were predicted to be negative, indicating aquifer recharging conditions; while at low flow, the aquifer was discharging at all locations.At Site B, the predicted range of flux rates were overlaping with each other across the zero flux line, showing that the uncertainty of the prediction may lead to opposite flux directions at certain locations with relatively small flux magnitudes.At Sites D and E, the predicted flux rates varied greatly (from +1 × 10 −6 m s −1 to −1.5 × 10 −6 m s −1 ) with respect to different flow conditons.The results also suggested that sites located in the primary channel (i.e., D and E) had remarkably higher flux magnitudes (up to 6-9 times higer) compared to those in the secondary channel (i.e., A and B), indicating strong spatial heterogenerity induced by geomorphological features.

Discussion and Conclusions
In this study, we introduced an approach for inferring HEF rates based on temperature

Discussion and Conclusions
In this study, we introduced an approach for inferring HEF rates based on temperature measurements using the stage-flux relationship.This approach allows us to obtain long-term, high spatiotemporal resolution hydrologic exchange dynamics in a large river reach based on the river

Discussion and Conclusions
In this study, we introduced an approach for inferring HEF rates based on temperature measurements using the stage-flux relationship.This approach allows us to obtain long-term, high spatiotemporal resolution hydrologic exchange dynamics in a large river reach based on the river flow and inland groundwater conditions with relatively easy field installations.Our predicted HEF rates vary between +1 × 10 −6 m s −1 to −1.5 × 10 −6 m s −1 under various flow conditions in the shallow water region along the river bank.This results are comparable with previous studies performed in other regulated rivers (e.g., [20,21]), in terms of the order of magnitude of the flux rates and the mechanism between the HEFs and the dam operations.The strength of this approach is that it remains robust under high-frequency river stage variations imposed by flow regulations (i.e., dam operations).Such variations rarely occur in a natural river except for during an extreme flooding event, which means it may not be generally applicable to fluvial systems with constant groundwater discharge from a confined aquifer.In small highland rivers, local geomorphological features such as riffles and pools might be a greater dominant factor than the head difference in controlling HEF rates [56,57].
Uncertainties in this approach could come from field measurements, parameterization of physical properties, model assumptions, modeling resolution, and model structure [58].These include (1) the accuracy of the iButton sensors and the lags on temperature signal recording due to the thermal skin effects of the measuring rod might affect the inferred HEF rates [59]; (2) river bed geomorphic features (e.g., dunes) may create local perturbations on the HEF patterns and thermal processes; (3) the riverbed temperature was approximated by measurements ~20 cm above the river-bed surface due to site and equipment limitations.The water layer between the sensor and the river bed surface may create temperature signal attenuations and lead to uncertainties in flux estimates.Our future studies will evaluate and consider the attenuation effects of the water in reducing the uncertainties; (4) the ranges of the physical property values used in the LPML model parameterizations (Table 1) were from the literature, which may not fully cover the uncertainties of these parameters; (5) the approach we used to identify the time period for different models were based on a frequency analysis (DHR), which is arbitrary and might introduce errors to the analysis.In fact, the choice of applying DHR for data-screening is site-specific, which is only required in locations with strong upwelling fluxes periodically diminishing or even removing the diurnal signals from the temperature data; (6) fluxes at the horizontal and lateral directions, which are smaller in magnitude compared to the vertical fluxes, may still affect the vertical flux estimations [60]; (7) the size of the moving window and the temporal resolution used in the LPML model may introduce uncertainties to the modeling results [39].Another limitation of this approach is that it is only valid at locations where the HEFs are dominated by the stage gradient between river and inland water table.In other words, the hydrostatic driver controls the HEF dynamics [44].However, in some areas in the river, the HEFs are controlled by the hydrodynamic driver, which is determined by the local hydraulic conditions.For example, at the upstream tip of the island, where the hydrodynamic pressure is always high due to the speed of the river flows, the HEF rate may always be negative (i.e., downwelling).
The results of this study suggest that the magnitude of vertical fluxes in the primary channel could be about 6-9 times greater than in the secondary channel.Given that the head differences between the river stage and inland water table are similar across these locations, the major controlling factor that leads to this difference is the permeability of the riverbed sediment.We measured the flow velocity at the iButton ® locations and found that the mean velocity at the primary channel sites (0.4 m/s) was about two-fold higher than at the secondary channel sites (0.2 m/s).According to the Hjulstrøm Curve [61], which demonstrates the sediment erosion, transport, and deposition with regard to the flow speed and grain size, the averaged grain size of the riverbed sediment between the primary and secondary channels might differ by about three times, which then leads to approximately nine times the difference in permeability according to a number of empirical models such as Krumbein and Monk's equation [62], or the Kozeny-Carman equation [63].This is an overly simplified explanation of the flux differences in the two channels and requires further development and validation, but linking the HEF rate and the surface flow speed might be an informative way of estimating vertical flux in deep waters where the field measurements are not available and are hard to obtain.For example, the center of the primary channel in the study reach is as deep as 15 m under low-flow conditions, which poses a safety issue for installing iButtons ® or any other instruments.
To conclude, in this study we successfully inferred the sub-daily vertical HEF rates at five shallow water sites along the shoreline based on temperature profile measurements and relationships between the riverbed hydrologic exchange rate and river flow conditions.The results reveal that the HEF rate had high spatial and temporal heterogeneities over the riverbed, and a magnitude of fluxes 6-9 times higher in the primary channel than in the secondary channel.This approach can be easily employed in other river reaches to facilitate large scale river corridor studies or to inform biogeochemical and ecological studies in highly dynamic large river reaches.

Figure 1 .
Figure 1.The study river reach.The iButton sites A-E are marked as red triangles.The Spring 9 location is marked as a blue triangle.

Figure 1 .
Figure 1.The study river reach.The iButton sites A-E are marked as red triangles.The Spring 9 location is marked as a blue triangle.

Figure 2 .
Figure 2. Transect profiles at T-A, T-C, T-spring9, and T-E at 1:10 aspect ratio (vertical : horizontal).The insets of the plots showed the locations of the measurements in zoomed in 1:1 aspect ratio plots.

Figure 3 .
Figure 3. Temperature data collected from Spring 9 (A); estimated cross-bed flux rates based on piezometer data and Darcy's Law (B); and river stages and associated inland well water levels during the same period (C).

Figure 2 .
Figure 2. Transect profiles at T-A, T-C, T-spring9, and T-E at 1:10 aspect ratio (vertical : horizontal).The insets of the plots showed the locations of the measurements in zoomed in 1:1 aspect ratio plots.

Figure 2 .
Figure 2. Transect profiles at T-A, T-C, T-spring9, and T-E at 1:10 aspect ratio (vertical : horizontal).The insets of the plots showed the locations of the measurements in zoomed in 1:1 aspect ratio plots.

Figure 3 .
Figure 3. Temperature data collected from Spring 9 (A); estimated cross-bed flux rates based on piezometer data and Darcy's Law (B); and river stages and associated inland well water levels during the same period (C).

Figure 3 .
Figure 3. Temperature data collected from Spring 9 (A); estimated cross-bed flux rates based on piezometer data and Darcy's Law (B); and river stages and associated inland well water levels during the same period (C).

Figure 4 .
Figure 4. Vertical fluxes versus the head difference between the river stage and inland water table level at Spring 9 at daily and half-hourly scales.

Figure 4 .
Figure 4. Vertical fluxes versus the head difference between the river stage and inland water table level at Spring 9 at daily and half-hourly scales.

Figure 5 .
Figure 5.The framework employed to estimate hyporheic fluxes.

Figure 5 .
Figure 5.The framework employed to estimate hyporheic fluxes.
scales.The same modeling framework was then used to estimate the vertical fluxes in the shallow water along the river, as introduced in Section 2.4.Water 2017, 9, 703 9 of 17

Figure 6 .
Figure 6.Demonstration of the data filtering procedure based on Spring 9 temperature data.

Figure 7 .
Figure 7. Relationship between model-estimated vertical flux rates and the head gradients between the river stage and inland water level based on 100 realizations of randomly selected combinations of sediment properties at Spring 9.

Figure 6 .
Figure 6.Demonstration of the data filtering procedure based on Spring 9 temperature data.

Figure 6 .
Figure 6.Demonstration of the data filtering procedure based on Spring 9 temperature data.

Figure 7 .
Figure 7. Relationship between model-estimated vertical flux rates and the head gradients between the river stage and inland water level based on 100 realizations of randomly selected combinations of sediment properties at Spring 9.

Figure 7 .Figure 8 .
Figure 7. Relationship between model-estimated vertical flux rates and the head gradients between the river stage and inland water level based on 100 realizations of randomly selected combinations of sediment properties at Spring 9.

Figure 8 .
Figure 8.The range of fluxes (i.e., the 5th and 95th percentiles) estimated by the proposed model compared with the Darcy's Law-based fluxes at Spring 9 at daily (A) and half-hourly (B) time scales.

Table 1 . 6 2. 4 .
Ranges of physical parameters used to generate parameter sets in the model.Parameters Range (from Ma et al. 2012) Porosity 0.15~0.18Solid density (kg/m 3 ) 2650~2760 Fluid-specific heat capacity (J/(kg K)) 4186 Soild-specific heat capacity (J/(kg K)) 715~920 Solid thermal conductivity (W/(m K)) 1.2~2.2Fluid thermal conductivity (W/(m K)) 0.58 Thermal diffusivity (m 2 /s) 5.19 × 10 −7 ~2.24 × 10 −Monitoring near Riverbed Temperature along the River with each other across the zero flux line, showing that the uncertainty of the prediction may lead to opposite flux directions at certain locations with relatively small flux magnitudes.At Sites D and E, the predicted flux rates varied greatly (from +1 × 10 −6 m s −1 to −1.5 × 10 −6 m s −1 ) with respect to different flow conditons.The results also suggested that sites located in the primary channel (i.e., D and E) had remarkably higher flux magnitudes (up to 6-9 times higer) compared to those in the secondary channel (i.e., A and B), indicating strong spatial heterogenerity induced by geomorphological features.

Figure 9 .
Figure 9.Time series of observed temperature at the five iButton locations.

Figure 10 .
Figure 10.Time series of river stage and inland water level during the iButton monitoring period.

Figure 9 .
Figure 9.Time series of observed temperature at the five iButton locations.

Figure 9 .
Figure 9.Time series of observed temperature at the five iButton locations.

Figure 10 .
Figure 10.Time series of river stage and inland water level during the iButton monitoring period.Figure 10.Time series of river stage and inland water level during the iButton monitoring period.

Figure 10 .
Figure 10.Time series of river stage and inland water level during the iButton monitoring period.Figure 10.Time series of river stage and inland water level during the iButton monitoring period.

Figure 11 .
Figure 11.time series of estimated fluxes during the monitoring period.

Figure 12 .
Figure 12.Estimated vertical fluxes at the five iButton locations under three representative flow conditions.

Figure 11 . 17 Figure 11 .
Figure 11.Time series of estimated fluxes during the monitoring period.

Figure 12 .
Figure 12.Estimated vertical fluxes at the five iButton locations under three representative flow conditions.

Figure 12 .
Figure 12.Estimated vertical fluxes at the five iButton locations under three representative flow conditions.

Table 2 .
Coordinates (based on NAD 1983 StatePlane Washington South system) and depths of the iButton sensors installed at the monitoring sites.