Spatial Variability Pattern of Hyporheic Exchange in a braided River

The streambed flux is variable in space; the spatial variability results in part from bedforms, but few works on streambed fluxes in channels with strongly abrupt varying bedforms are carried out. Heat as a tracer to delineate the streambed flux pattern has been widely adopted in numerous fields. In this paper, a braided channel with complicated topography was selected as study site, where the temperature was monitored. One-dimensional (1-D) analytical method based on the amplitude attenuation (Ar) and 1-D numerical method were used to interpret the temperature. As a result, streambed fluxes of a total of 50 sites in the braided channel are obtained. From the results we can know the magnitude and direction of streamed flow velocity are spatially variable, even within a 1-m distance. Then, this study summarizes five bedform-driven flux patterns: 1 downward flow driven by the head difference between groundwater and stream, 2 downward flow related to a meter-scale pool, 3 a transition from upward to downward flow associated with a centimeter-scale riffle, 4 horizontal flow in braided bars and 5 upward flow driven by vegetation roots. Overall, multiple physical mechanisms together contributed to the complex streambed flow system, which reflected great challenges for the scaling up of point-in-space seepage flux.


Introduction
Hyporheic zone is a hot spot as a biogeochemical response to hyporheic exchange. Studies for this zone play an important role in regulating water resources and maintaining ecosystem function. Understanding and quantifying flux exchanges between streams and groundwater is very crucial in the prevention of contamination, the restoration and maintenance of ecosystem and the management of coupled water resources.
Hyporheic exchange fluxes are varied in space and are controlled by both natural and anthropogenic factors, e.g. asynchronous oscillation of stream water and groundwater stage (Wroblicky et al. 1998;Wu et al. 2015), river geomorphology and stream curvature (Boano et al. 2006), the depth and permeability of bed sediments (Anderson et al. 2005;Buffington and Tonina 2009), etc. The spatial variability in streambed fluxes is ubiquitous. Conant (2004) and Schmidt et al. (2006) took snapshots of heterogeneous groundwater exfiltration in a river. Even, the rates of exchange between surface water and groundwater differ dramatically in magnitude and direction over a meter scale (Lautz and Ribaudo 2012). In a view of overall watershed, the potential for surface water/groundwater exchange flow decreases along the continuum from headwater to mid-order mountain streams (Anderson et al. 2005). The heterogeneity in K is ubiquitous (Kalbus et al. 2009), which results in the spatial difference of water fluxes. Conant (2004) observed high discharge areas that were located in areas with high hydraulic conductivity deposits in the streambed. Besides, a great number of literatures (e.g. Hester and Doyle 2008; Wu et al. 2016) reported tremendous impacts of various in-stream geomorphic structures on hyporheic exchange. Winter (1998) pointed out that the direction of seepage flux across the streambeds is commonly related to abrupt changes in the slope of the streambed or to meanders in the stream channel. The stream usually loses water at the upstream end of a riffle or a channel bend, and gains water at their downstream ends (Winter 1998). Lautz et al. (2010) found the beaver dam or log dam altered previous vertical hydraulic gradient, influencing the stream/groundwater interchange upstream and downstream of these obstacles. Lu et al. (2012) modeled that groundwater discharge was disproportionally distributed along the cross section of a stream, and mostly occurred on the part of streambed in the deeper-water areas. To date various conceptual models of hyporheic exchange under different topography conditions have been developed and to some extent revealed topography-dependent streambed flow characteristics (Conant 2004). In this paper, we used these flow models to generalize the flow characteristics of a braided stream.
A range of methods are available to determine surface water/ groundwater exchange rates, such as seepage meters, differential discharge measurements (Essaid et al. 2008), hydrograph separation, etc. In terms of spatial scale, resolution and accuracy, these traditional methods face great challenges. Their applicability and limitations have been summarized in detailed (Hatch et al. 2006;Kalbus et al. 2006).
Temperature as a naturally occurring tracer can be exploited to infer the movement of water between aquifers and surface-water bodies. Though using heat as a proxy to estimate bulk groundwater recharge/discharge may not be practical because of naturally exiting problems with heterogeneities, horizontal flow components, etc., it may still provide enriched information on spatial and temporal variability in flux rates on a point scale (Ferguson and Bense 2011). Even scaling up of point-in-space heat tracing has been initially developed to determine seepage flux where temperature measurements were made (Conant 2004;Lautz and Ribaudo 2012). Temperature measurements over time of both stream water and its underlying sediments in conjunction with the hydraulic head measurements at depth has become a routine procedure (Essaid et al. 2008) to determine spatial patterns and the magnitude of surface water/groundwater interactions, which benefits from significant improvements in computational techniques and the accessibility of high-quality temperature acquisition. The benefits of heat transport method have been recognized for decades. The earliest step of utilizing one-dimensional flow and heat transport equation to computer infiltration rates was made by Stallman (1965). A decade ago, Hatch et al. (2006) and Keery et al. (2007) extended the method by deriving an analytical solution to the equation and enabled the automation to remote and continuous monitoring of seepage rates. That is to say, one benefit of this seepage estimation method is the ability to model how seepage changes over time. Furthermore, information on vertical variability in water flux can be detected using heat transport method. Besides, numerical approaches based on temperature data have become a useful tool to calculate seepage rates (Skinner 2006) with the development of computer technology. Though more data are needed (i.e., initial condition, boundary condition, etc.), numerical modeling probably offers more accurate results compared to analytical methods because it is closer to real world scenarios (i.e., dimensionality and heterogeneity inherent to natural near-surface systems) (Rau et al. 2014). In recent years, a device employing heat as an artificial tracer was developed to measure hyporheic flow direction and velocity in three dimensions at a scale of a few centimeters (Angermann et al. 2012). There is no doubt that the burgeoning interest of utilizing diel temperature fluctuation to calculate water flow velocities will advance near-surface hydrological process understanding, albeit with several unresolved issues.
In terms of software development and application of experimental apparatuses, analytical models derived from heat transport equation have been greatly improved in particular during the past decade. But the validity of using the 1-D heat transport model under different field circumstances is still ongoing; the influence of bed heterogeneity, extreme flow events or violated boundary conditions on water flow though porous media has not yet been fully unraveled and require further research. Ferguson and Bense (2010) modeled the effect of K heterogeneity on 1-D heat flow analysis to estimate streambed flux. Lautz (2012) conducted sand column experiment to explore the effect of storm events on the 1-D analytical solution. However, the influence analysis of non-ideal conditions based on field data has rarely been made. In this study, temperature array sites with different heterogeneities and flow conditions (horizontal to oblique flow) were chosen. The detailed insight into the spatial variability of vertical hyporheic flux in such a complex scenario will be revealed in next sections.

Study site
Our study site is located in Zhongtian River, southern part of Jiangsu Province, China. Zhongtian River flowing from south to north is a typical river in hilly areas, with a river gradient of 92.7 ‰ and a basin area of 47.43 km 2 . It is one of the three most important rivers entering Lake Tianmuhu, which is the primary drinking water source for 700,000 people in the city of Liyang. The Zhongtian River basin with a subtropical monsoon climate is the main source of water conservation for Lake Tianmuhu. It has hot summers and cold winters with an annual average temperature of 17.5 ℃. The area receives approximately 1150 mm of precipitation annually, almost as rain, but its monthly average precipitation differs greatly throughout the year. The amount of precipitation in summer almost reaches 50% of annual precipitation, but in contrast winter only has a proportion of 12.7% (according to average data of 1971-2010). The basin is underlain by shallow carbonate rocks, and some hard rocks emerge from the land surface. Nowadays, in the catchment area, the development of irrigated agriculture, especially as tea plantation expanded, led to the rise of nitrogen concentration in the groundwater. At the same time, surface water body will be threatened by contaminated groundwater exfiltration.
Stream environments range from steep, narrow, bedrock chutes, to wide alluvial reaches. The field experiment is conducted in the middle reaches of Zhongtian River. This reach is unconfined and developed into a braided channel, with bundles of bushes gathered around it. Our study area is around 30 m length and 10 m width, with a great local heterogeneity; the middle of the stream or narrow channels between braided bars comprise mixed sand and gravel, with a few larger boulders of poorly consolidated sand stone, whereas, in contrast, the braided bar itself distributed within the stream primarily consists of fine-grained sand, with occasionally silt or clay. Braided rivers occur when a threshold level of slope is reached and a threshold slope was experimentally determined to be 0.016 (m/m) for a stream with poorly sorted coarse sand (Schumm and Khan 1972), which is consistent with the channel in this paper. A total of 4 braid bars were formed in the channel as shown in Figure 1. Braid bars are usually highly mobile, with the river layout often changing significantly during flood events (Hickin and Sichingabula 1988). Relatively these bars in this field site are hard to erode and can remain stable due to the consolidation imposed by bushes growing on them.

Figure 1The
location of the study reach in the Zhongtian River catchment that is nested in the Lake Tianmuhu Watershed is shown together with the numbered temperature sensor array used throughout. The relative bed surface elevation is surveyed on 22 November 2010 at the beginning of the study period and remained stable thereafter.

Temperature data monitoring
Temperature sensors were used to monitor the temperature of stream water and groundwater, respectively. One temperature sensor was installed in a nearby well that is about 100 m away from the river bank, and the other on the streambed surface. The monitoring campaign was from 22 November to 2 December, lasting for a 10-days period, during which surface water was relatively cold compared to the shallow groundwater temperature with approximately 18 ℃. The deep groundwater even had a lower temperature of 16.5 ℃ (Provided by Jiangsu Provincial Water Resources Bureau). The stream water temperature signal had large diurnal amplitude swings, and even approach linear without minor fluctuation for 24th, 25th and 30th, November, which was subject to precipitation events ( Figure 2). Still, the temperature signal can be considered as a quasi-sinusoidal oscillation, with a 1-day period, maxima generally occurring mid-day. Surface-water temperature followed a decreased trend over time, so did shallow groundwater ( Figure 2). According to the preliminary survey, the streambed flux varied spatially, occasionally under upward conditions. In order to capture the amplitude attenuation or phase lag of temperature signal propagation at depth, PVC (Polyvinyl chloride) pipes (2.0 cm outer diameter) were slotted normal to pipe length at close 0.02 or 0.04-m spacing ( Figure 7). The thin, flat, circular iButtons were strapped tightly on the slots using nylon cords, and covered with layers of water-proof PTFE (polytetrafluoroethylene film) tape to prevent water from damaging sensors. This PTFE material conducts heat actively to help maintain good thermal contact between the sensors and sediments, and reduce thermal skin effects (Cardenas 2010).Resolution and accuracy of sensors were quantified prior to deployment using a reference sensor with resolution and accuracy of 0.005℃.The pipes, each along with five sensors, were driven vertically into the sediments using a rubber sledge so that the iButtons could be positioned at 0.02, 0.04, 0.06, 0.08 and 0.12 cm below the sediment-water interface, respectively. The temperature sensors were programmed to record temperature every 10 min, which was sufficient to accurately characterize the heat wave signal in the streambed (Hatch et al. 2006). Four cross sections were arranged along the stream channel, with an interval of 5 or 10 m (shown in Figure  1). Two to four pipes were installed at each transect to derive temperature profile across the stream.
As an attempt to obtain the head gradient in the streambed, PVC apparent pipes were vertically driven into the desired depth; thereby, the head gradient can be recorded roughly through observations. However, studies from previous literatures showed that it can be challenging to obtain high-quality head measurements in streambeds where head differences are probably on the order of a few centimeters, even less in coarse sediment  2004) found it difficult to obtain significant head gradients over the 0.5-1 m vertical distances at which heads were measured. Considering large uncertainties associated with streambed head measurements and spatial variability of head gradients, these head gradient data were not used to determine the groundwater recharge/discharge relation. Differential discharge measurements were taken to estimate the exchange between streams and groundwater. They were conducted along each transect between 14:00 and 15:00, 22 November, 2014, during which we assumed that variations in discharge were minimal. Detailed procedures for measuring discharge can be seen in Bartram and Ballance (1996). The discharge results at each cross section were presented in Table 1. As braided bars had much lower K compared to channels, the influence of braided bars on overall discharge contribution was ignored. Thus from Table 1 we can see in our study area the stream water lost water into the groundwater.

Sampling
After the temperature collection campaign was completed, a total of 22 streambed sediment cores were collected at pipe installation points at the shallow depth of 0.05-0.15 m, manually using a shovel. Field observations told us that the sediment compositions differ greatly between samples; the midstream cores consist of coarse sand and gravels, and the cores near banks silt or clay. These samples were taken to the state key laboratory of hydrology-water resources and hydraulic engineering, China, for gain size analysis. These fine sediment samples were analyzed with a laser diffraction particle size analyzer (LS13320, Beckman Coulter), and for the other samples, the standard sieving analysis method was used to separate grain-size grades. This sieving analysis was made through seven sieves ranging in size from 0.075 to 2 mm in aperture diameter; thereby grain size distribution were determined, and the class and sorting of sediments were presented in Table 2.

Heat transport modeling
Previous hydrogeologists primarily studied two zones in the subsurface: the surficial zone and the geothermal zone; the former is influenced by diurnal/seasonal heating and cooling of the land surface, but the latter is isolated from them (Parsons 1970). In this study, stream hyporheic zones are a part of the surficial zone, which are subject to diurnal fluctuation of stream water temperature. In general, stream water temperature may vary with a quasi-sinusoidal waveform, while groundwater temperature remains constant. Surface water temperature signal propagates downwards into underlying sediment matrix, by conduction and advection when under the downwelling condition (losing river), while the downward propagation is transported by conduction only and is restrained by the opposite advection when the direction of flow is reversed (gaining river). Natural temperature variation in streams propagating downward into the streambed is the driver to track the streambed water flux. Early in 1960s, the quantitative study of utilizing heat to estimate groundwater flow velocity has been conducted. The governing equation for one-dimensional conductive and advective heat transport, analogous to the equation of solute and water movement, was presented by Suzuki (1960) and Stallman (1965): It has been assumed that the sediment properties do not change between the two points. In this study, the one-dimensional heat transport model presented by Hatch et al. (2006) was used to quantify the flux rate of water through the sediment matrix. Briggs et al. (2012) pointed out that relative temperature accuracy is more important than absolute accuracy. Given that the data of sensors along temperature profiles were not logged synchronously, interpolation between temperature data is done.

Analytical modeling
Raw temperature time series of stream water and its underlying streambed are typically not pure sinusoids, and may be perturbed, more or less by numerous processes including precipitation, snowmelt, upstream watershed characteristics, sensor accuracy. Thus to isolate the daily temperature signal of interest from raw signals with multiple frequencies, signal processing techniques must be used. Keery  rates between sensor pairs. In this study, vertical flux rates of every 1.5-h were derived during the study period using the one-dimensional heat transport model described above. Table 3 presents parameter values used in the model.

Numerical modeling
Although analytical models are straightforward to implement, especially owing to an emergence of the graphical user interface (Swanson and Cardenas. 2011), there exist to be uncertainties or errors when solving problems in the presence of surface-water fluctuations (e.g. induced by precipitation, water release by dams) or other complicated forcing mechanisms (Voytek et al. 2014). In this study, 1DTemProwas utilized to analyze vertical 1Dtemperature profiles (Voytek et al. 2014). 1DTemPro is a graphical user interface to the numerical model VS2DH (Healy and Ronan 1996) that numerically solves the flow and heat transport equations in porous media. 1D vertical temperature profiles are analyzed by fitting modeled temperature data to observed data. The calibration commonly focuses on adjustment of flux or, if head data are available, K, which is facilitated by immediate display of model results against measured data with each model run. For the field data, the time-varying temperature data recorded by temperature sensors at 0.02 m and 0.12 m depth were applied to the top and bottom boundary conditions. The 0.06 m depth temperature data were fit through manual adjustment of flux and K. The level of fit can be judged by mean absolute error (MAE). The same thermal parameters used in VFLUX were employed for the numerical model domain. Homogenous material condition was set throughout the vertical profile. Due to the influence of precipitation processes and subsequent stage fluctuation, the measured time series were broken up into 1-d segments, each for calculating streambed flux rates, to accommodate potential daily variation in vertical flux.

Temperature profile
In this study, temperature measurements were conducted at 5 depths for 12 locations (Figure 1). Vertical temperature profiles can reveal many interesting and valuable attributes of streambed flow condition, even before the 1-D flux calculation was completed. In this study site, temperature profiles differed greatly from one location to another, which indirectly revealed significant differences in vertical hyporheic exchange between various locations. Under downwelling conditions, fluid flux transports stream water temperature signals into sediments along with the inherent heat conduction, so that the temperature front can reach a relatively deep depth, whereas, under upwelling conditions, the temperature can only conduct downwards against the opposite flow and it is restricted at a shallow depth. The penetration depth of temperature signals can tell us the qualitative information on hyporheic fluxes. For example, the high temperature of 18.5℃at around 4 pm penetrated into the deepest depth more than 0.12 m for the D1L2, while for the D1L3it was only constricted near the bed/stream interface, which indicates D1L2 and D1L3 were under different flow conditions ( Figure 3). Interestingly, this remarkable contrast can be observed over a lateral distance of no more than 3 m within the streambed. Thus from the temperature profiles we can find the flow conditions differed between various locations ( Figure 3). Also, the diel vertical thermal envelops can serve as a qualitative indicator of vertical fluid flux through near-surface sediment systems. The narrower the lateral width of the thermal envelop is, the smaller the diurnal fluctuation is. When groundwater flows upwards towards the stream, the diurnal thermal envelope in the streambed converges on the groundwater temperature at a shallower depth than when the diurnal signal is partially carried down into the streambed by downwelling water. For example, the contrast in the envelope convergence between D1L2 and D1L3 also verified the distinct flow conditions at these two locations (Figure 4). Spurious low-temperature or high-temperature discrete patches was detected in Figure 3 and non-uniform changes of temperature with depth were found in Figure 4, which together indicates that advection of heat and associated vertical hyporheic flux were not uniform with depth. Most temperature profiles had only low-temperature and did not have high-temperature discrete patches except for D2L2. The monitored data showed there was a 0.01-0.02 ℃ difference between these low temperature patches and their surrounding portions. This small difference was probably attributed to either noise ratios of sensors or subtle flow change with depth due to the heterogeneity of the bed medium. In contrast, the D2L2 profile had both low-temperature and high-temperature discrete patches. Different from other profiles, both the low and high temperature discrete patches for the D2L2 have a relatively great temperature difference; the low temperature patch was 0.2 ℃ different from the area around it, but the difference in temperature for the high temperature patch reached a maximum of 1.4 ℃. As the D2L2 was located in a braided bar with low permeability, complex flow condition may convectively alter the vertical temperature distribution.  No temperature data is available for the 12 cm depth at the locations D2L1 and D3L1 due to data logger failure in the field. As no prior knowledge of the penetration depth of diurnal temperature signal was obtained, the thermal envelop was in part presented and no points of convergence could be seen in each panel above.

Spatial variability of streambed flux
Streambed temperature profiles from all locations were modeled with the diurnal signal-based 1D analytical method and with 1D numerical model. The temperature data at 0.02 m and 0.12 depths were input into the model as the upper and lower boundary, respectively. 1DTempPro was utilized to simulate the 0.06 m depth data for every location in 1-day sub-sections. Specific discharge as a hydrologic constraint was adjusted to match the model-derived temperature data with the observed data. The performance of fit was accessed using the root mean square (RMS). The RMS values were nearby the intrinsic system noise (0.02 ℃) so good calibration had been made ( Figure 5), indicating the accuracy and reliability of 1DTempPro-derived fluxes. Temperature time series at 0.02 m and 0.12 depths also was applied to the analytical method, which yielded the mean flux value between 0.02 m and 0.12 depths.
Overall, there is a moderate agreement between the two methods in that linear fits to the modeled data indicated the a marked trend with increase in exchange flux during the period of investigation, except for the declining trend at locations D4L2 and D4L4. Thus it can be inferred that the 1D analytical method can successfully predict the relatively low or high flow condition and the trend in fluxes over time. Nevertheless, a small discrepancy between these two methods was detected. Compared to the numerical method, the analytical method yielded much smaller flux values, consistently less than 0.6 m/d, and variationsin fluxes over time was subdued ( Figure 5). The streambed flux varied spatially in both magnitude and direction. The fluxes calculated showed, during the investigation period, most locations were primarily under the losing condition, but the fluxes differed from one location to another. It can be inferred that this study site had a great spatial heterogeneity in seepage rates.

Discussion
In this section we will analyze the influence of forcing mechanisms (e.g. topography, heterogeneity, vegetation, etc.) on the streambed fluid flux. This study with such diversity is conducted in a peculiar site of braided rivers that has riffle/pool bedforms, tremendously heterogeneous sediments, sporadically covered vegetation, etc.
Generally, the hyporheic flux varies with depth. Vogt et al. (2010) utilized a temperature profiler with an improved spatial resolution to derive flux rates that were found to be vertically variable at the centimeter scale. Though our temperature profile was not that high-resolution, it can still infer spatial information on water flow velocity in the streambed. In our study, the departure of envelop lines with depth from the linear as well as low/high-temperature discrete patches along temperature profile simplicitly indicates the non-uniformity of streambed water flow with depth; a spectrum of horizontal to oblique flow may exist around the temperature pipes, which was probably caused by the great heterogeneity of sediments.
The vertical flux rates varied spatially; magnitude or even direction in flux differed greatly between close locations ( Figure 5). Large spatial variability in flux has been reported at the scales of 10 s of meters and even down to 1 s of meters by numerous researchers. Although these locations in our study were not located far away from each other, they may be controlled by different physical mechanisms. An interesting finding is that the D1L2 and D1L3 locations had an absolutely opposite direction in flux. The upward flow at Location D1L3 may be, to a great extent, influenced by its upstream bushes. The schematic flow line around D1L3 can be seen in Figure 7.In contrast, the D1L2 location with downward movement of water was located at the tail of a pool, where the water flux was controlled by local bedform change associated with the pool. Generally, stream water flows into sediments at the tail of the pool (Winter 1998;Tonina and Buffington 2009). Likewise, the stream flux at the D1L1 was influenced by the pool. Interestingly, the flux across the streambed surface at Location D2L1 was impacted by a riffle, and this riffle was small, nested within this larger-scale pool (Figure 1  and 7).The riffle can exert a strong influence on the flow when they are partially exposed or only shallowly submerged ), so the D2L1 flux moved upward during the early days of the investigation, which was controlled by the shallowly submerged riffle. As flow increased, so does water depth, the bedforms exert less influence on the streambed flow (Tonina  and  Buffington  2009). Thus, the precipitation-induced stream stage rise in the latter days made the D2L1 streambed flux much less influenced by bedforms; the hydraulic gradient drove stream water into the sediments (Figure 7). Though the Location D3L1 was situated at the head of the pool, this bedform change did not impose great influence on the D3L1 flux. The streambed flux at Location D3L1 may be dominantly influenced by the head difference between stream water and groundwater (Figure 7). Likely, the D4L1 downward flow was attributed to this downward hydraulic gradient (Figure 7). Considering very coarse-grained sediments and the low submergence, the flow condition at the D4L2 and D4L3locationswas probably unstable, and the turbulent flow may increase the difficulty of understanding the exchange flux in the near-surface streambed.

Conclusion
The combination of A r -based analytical solution to 1-D heat transport equation and temperature time series data was used to interpret the flux rates of water across the groundwater/surface water interface during the 10-days period with three precipitation events. 1-D numerical method was applied to demonstrate the applicability of the A r analytical method under non-ideal conditions. Though the analytical method based on the amplitude decay was subject to a great number of limitations, on the whole, the A r method yielded acceptable results and the results can distinguish relatively high and low rates of water flux and detect the temporal trend in streambed flux.
Though situated in a 20×10 m area, local bedform change or vegetation cover may create unique small flow systems nested within the relatively large study area. The Location D2L1 was dominated by a small riffle and had water flux moving upward when the stream flow was slow and water depth was shallow at the early days of the investigation. As stream stage increased due to precipitations, the riffle was less influential. The downward hydraulic gradient drove stream water into the sediments, producing downward flow at Location D2L1. Likewise, the D1L3 flux was influenced by a small-scale physical mechanism. The roots submerged in upstream sediments of D1L3 greatly altered the flow system around the location D1L3, contributing to the upward movement of fluxes at Location D1L3. Besides, the pool in our study area imposed great influence on the D1L1 and D1L2 locations. Other locations (D3L1, D4L1, D4L2 and D4L4) were less influenced by the bedforms and vertical hydraulic gradient drove stream water moving into the streambed. Note that D4L2 and D4L4 were situated in channels with a great number of cobbles and may be potentially dominated by turbulent flow. The bedforms, sediments and even some plants grown in near-bank channels, together, contributed to the spatial and temporal variability of fluxes in this study area. Therefore, a few spot measurements will not adequately capture the spatial flow characteristics, which provide great challenges for the scaling up of point-in-space seepage flux in complex bedform conditions.