A mathematical approach to improving the representation of surface water–groundwater exchange in the hyporheic zone Open Access

Rivers are among the most fascinating freshwater bodies on Earth, integrating many forms of ecological processes that influence the water quality (Brunke & Gonser 1997; Boulton et al. 1998; Boano et al. 2014). Recently, much progress has been made in understanding the river corridor system; however, there is one critical limitation: river modules cannot be treated as a closed conduit and connectivity between rivers and their surroundings needs to be investigated because of their ability to trigger complicated forms of biogeochemical functions that can influence the water quality. For example, vertical exchange flow through the riverbed sediments and the lateral exchange flow within the surrounding riparian and floodplain increases the residence time in reactive environments and their associated dissolved, suspended materials, and thus increases the opportunities for biological and geochemical processing (Jones & Mulholland 1999; Battin et al. 2008). Hence, the river corridor must be studied as a holistic system instead of an individual ecological element (National Research Council, Division on Earth and Life Studies, Board on Environmental Studies and Toxicology, Water Science and Technology Board, Committee on Riparian Zone Functioning and Strategies for Management 2002; Boano et al. 2014; Harvey & Gooseff 2015).

Speaking of those various forms of connectivity between the river corridor and its surroundings, the hyporheic exchange flow (HEF) is one of the key elements. It is critical to control the water quality because it is the river's continuous exchange of water, solutes and energy with the riverbed sediments that leads to complex biogeochemical functions (Dudley Williams & Hynes 1974; Winter 1998; Jones & Mulholland 1999; Sophocleous 2002). Characteristic biogeochemical processes include denitrification, and carbon and nitrogen cycling. For example, the stream denitrification rate that controls the possibility of eutrophication is closely related to the residence time of water and uptake rates of the oxygen, controlled by the hyporheic transportation rate (Harvey et al. 2013). Thus, HEFs are critical to investigating the role of river corridor in the natural ecological system.

Many studies have been conducted to understand the physical mechanisms that drive the HEFs by unit geomorphic bedform features (Buffington & Tonina 2009; Tonina & Buffington 2009; Deb et al. 2019). Elliott & Brooks (1997a, 1997b) combined their theoretical and experimental studies to highlight the two mechanisms that drive water and its constituents back and forth from the water column to the bed sediment with the presence of the bedform: ‘pumping’ and ‘turnover’. The pumping model is the most fundamental component in the current state-of-the-art multiscale hyporheic exchange model system (Wörman et al. 2006; Stonedahl et al. 2010; Gomez-Velez & Harvey 2014) and serves as the basis to develop upscaled hyporheic exchange models. In addition to the physical mechanisms of the HEFs associated with the unit geomorphic bedform features, the modulating effects of the HEFs, such as the regional groundwater flow (Cardenas & Wilson 2006, 2007b; Boano et al. 2009; Mojarrad et al. 2019) and the dynamic stream flow conditions (Wondzell & Swanson 1996; Wroblicky et al. 1998; Boano et al. 2007, 2013; Sawyer et al. 2013; Singh et al. 2019), are also investigated.

Mathematical models that contain key hydrological components are ideal solutions for practical applications because computational fluid dynamics simulation is computationally expensive and impractical to apply on a large scale. Since the groundbreaking research on groundwater flow patterns from Toth (1963) to Vaux (1968), who first developed a physical-based model for the hyporheic flow zone, fruitful achievements have been made using mathematical models to understand both the physical mechanisms and the biogeochemical implications associated with the hyporheic exchange system (McClain et al. 2003; Mulholland et al. 2008; Alexander et al. 2009; Kennedy et al. 2009; Krause et al. 2009). Some of these efforts have been made using physical-based analytical models to study the individual stream channel features like ripples and dunes (Thibodeaux & Boyle 1987; Elliott & Brooks 1997b; Packman & Brooks 2001; Marzadri et al. 2005; Cardenas & Wilson 2007a; Wörman et al. 2007; Cardenas 2008, 2009b), and some of them have focused on a larger scale such as bars and stream meanders (Harvey & Bencala 1993; Boano et al. 2006; Cardenas 2009a, 2009b; Gomez et al. 2012).

However, the current stage of the limited knowledge of the connection between those individual river streambed features and the hyporheic zone cannot satisfy our need to optimize effective management of the multiscale river corridor system. Even though we have built many exciting and insightful hyporheic exchange models, we still have little predictive, transferable, confident understanding of the system so far (Ward & Packman 2019). This is because the controlling processes and mechanisms (Fehlman 1985; Elliott 1990; Elliott & Brooks 1997a, 1997b) used in all of these current models were based on simplified and empirical theories under certain assumptions and simplifications of natural bedform characteristics (Boano et al. 2014; Ward & Packman 2019). Even the most up-to-date model developer Susa acknowledged that the simplified representation of the streambed topography does not adequately characterize patterns and rates of hyporheic exchange (Stonedahl et al. 2010). Thus, we want to step further to investigate how bedform characteristics control the hyporheic exchange system in this proposed study.

In this study, the accuracy of Wörman et al.'s (2006) analytical solution was tested on a very simple harmonic synthetic surface, the performance of the approach on various synthetic surface scenarios was evaluated with different complex shapes, and a new version of the approach that can help improve the accuracy of both approximating the surface and the boundary fluxes was proposed. Through these analyses, we can answer three research questions: (1) Is Wörman et al.’s (2006) approach accurate to evaluate the surface–groundwater interaction flows? (2) How does Wörman et al.’s (2006) approach vary on different synthetic bedforms with different complex shapes? (3) Can spectral analysis help improve the accuracy of Wörman et al.’s (2006) approach in both approximating the surface and the prediction of the boundary fluxes?

By answering these research questions, here a novel version of the Fourier fitting process is developed. The novelty of this approach is incorporating the spectral analysis as a tool to guide us to fit the surface topography more accurately and reliably, and thus can further improve the prediction accuracy of hyporheic exchange metrics such as boundary fluxes in this study. We also demonstrate its ability to improve the accuracy of representing the surface and predicting the boundary fluxes, especially for those surfaces with characteristic scales and significant frequency signals.

In this study, the research activities were carried out in three stages. We began by testing the performance of Wörman et al.’s (2006) analytical approach to evaluate the surface–groundwater interactions on a very simple synthetic surface by comparing their analytical solution with the numerical solution solved in COMSOL Multiphysics. The numerical solution was used here as a ‘benchmark’ to test the validity of the analytical solution. The reason for using COMSOL Multiphysics is that it can be used to approximate the solutions for Partial Differential Equations (PDEs) based on the very fine discretization. These discretization methods approximate the PDEs with numerical model equations, which can be solved using numerical methods. The solution to the numerical model equations is, in turn, an approximation of the real solution to the PDEs with very high accuracy that can serve as a ‘benchmark’. Based on these considerations, we decide to use COMSOL Multiphysics to do the calculations. The purpose of doing this experiment is to verify the applicability and accuracy of their analytical solution in three-dimensional (3D) hydrological applications. Then, more scenarios of synthetic surfaces with different complex shapes were proposed to test how the Wörman et al. (2006) approach's performance varies with different complex shapes of synthetic surfaces. Last, we explored and proposed a new modified version of Wörman et al.’s (2006) approach by using the spectral analysis to identify the statistically significant frequencies of the surface to guide the refinement and mesh of the Fourier fitting process and test how this modified version of the Fourier fitting process improves the accuracy of both approximating the top surface and predicting the boundary fluxes.

In order to use the Fourier fitting to fit the surface topography, the amplitude coefficients are evaluated for a pre-assigned spectrum of and ⁠, where is wavelength (Wörman et al. 2006). Thus, the analytical solution can represent several topography observations in the form of ⁠, where C is a coefficient matrix containing the pre-assigned harmonics of (1), H contains the amplitudes and Z is an vector containing the observed -values in M locations (Wörman et al. 2006). A least-square estimation of the amplitude vector is obtained as (Wörman et al. 2006).

In order to examine the performance of Wörman et al.’s (2006) approach on different synthetic surfaces with different complexities, we created a series of synthetic surfaces with different numbers of known frequencies in the x- and y-direction, and applied Wörman et al.’s (2006) analytical solution on each of those synthetic surfaces. The series of synthetic surfaces created here are those harmonic synthetic surfaces with 1, 2, 4, 8, 16 and 32 frequencies in the x- and y-direction.

As an example shown here in Figure 2, we created a simple harmonic synthetic surface (panel (a) in Figure 2) with one wavelength in the x- and y-direction, and tested the current Fourier fitting method. Panels (c) and (d) of Figure 2 show that we used 8 and 16 fitting frequencies in the x- and y-direction to fit the real synthetic surface, respectively. Blue points in panels (c) and (d) of Figure 2 are pre-assigned fitting frequencies to fit the synthetic surface and the red points are the real frequencies. From panel (d) of Figure 2, we can clearly see that increasing the number of fitting frequencies helps approach the real frequencies (more blue points can approach the red points more closely), and the trade-off is incorporating many unreal frequencies (bringing in many blue points to approach the red points). Panel (b) of Figure 2 shows that although increasing the number of fitting frequencies improves the accuracy of approximating the real surface (red line in panel (b) of Figure 2), the accuracy of predicting the boundary fluxes decreases (blue line in panel (b) of Figure 2) since the fluxes are the derivative of the hydraulic head (Darcy's law) and can be damaged by those unreal fitting frequencies (many blue points in panels (c) and (d) of Figure 2). Thus, we hypothesize that it is the absence of the bedform characteristic signals in the current Fourier fitting process that leads to the absence of the characteristic hyporheic exchange patterns in the current multiscale model system.

Perron et al. (2008) illustrated that for natural surfaces, they might have characteristic spatial scales. They showed the method of using spectral analysis to identify the characteristic spatial scales from the background signals (parallel to the example shown above is to identify the red points before assigning too many blue points to approximate the surface). Hence, we believe that the spectral analysis can be incorporated into the traditional Fourier fitting process to capture the characteristic bedform signals, further improving the approximation accuracy of bedform surface, hydraulic head and hyporheic exchange.

In this section, a new version of the Fourier fitting process was proposed by incorporating the spectral analysis method shown in Perron et al.’s (2008) work to evaluate its ability to improve the accuracy of current Wörman et al.’s (2006) approach in both fitting the surface and predicting boundary fluxes. In Perron et al.’s (2008) work, they used the spectral analysis to identify the characteristic spatial scales of landscape and showed that landscapes are not always scale-invariant or random surfaces. Here, we adopt this idea and combine this method with Wörman et al.’s (2006) approach to better capture the real surface's information in a more efficient manner. For any arbitrary surface topography, before performing the Fourier fitting process, the spectral analysis method is used to evaluate the surface to obtain two critical pieces of information: (1) Does this arbitrary surface have characteristic spatial scales? (2) If the surface is not a scale-invariant or random surface, what are the significant frequencies of this arbitrary surface?

We evaluate the significance of the signals of the topography data by assessing the of the DFT array.

The Fourier transform treats the input signal (surface topography data) as if it is stationary, with the same mean, variance and frequency content throughout the sampled interval (Perron et al. 2008). However, few natural signals are strictly stationary (Weedon 2005); the power spectrum can only provide a useful description of a signal if its mean and variance are roughly constant (Priestley 1981). To remove the spatial trends in the mean of a topographic data set, a linear function in x and y (i.e. a plane) is fit to the input signal and then subtracted from z. It is called the detrend function which reduces the effects that can contaminate the power spectrum (Perron et al. 2008).

For any arbitrary surfaces, results obtained from the spectral analysis step are used to guide us in the following Fourier fitting process. Figure 4 shows the generic steps of how spectral analysis is incorporated into the original Fourier fitting process: if the surface is identified as a random surface that has no characteristic spatial scales, we directly apply the Fourier fitting process and assign the random frequencies to fit the surface; if the surface is identified as having characteristic spatial scales, we first retrieve and assign those characteristic signals of the surface to describe the surface, only apply the Fourier fitting process to fit the background (noise) part of the surface and combine two parts to generate final results.

The traditional Fourier fitting approach is applied to a series of synthetic surfaces. The influence of the number of Fourier fitting frequencies on synthetic surfaces with different complexities is investigated. If the traditional Fourier fitting process is optimized, we should always expect more accuracy in approximating the synthetic surface when increasing the number of Fourier fitting frequencies.

The spectral method (Wörman et al. 2006) is applied to the one-frequency harmonic synthetic surface. In the meantime, we use the COMSOL Multiphysics to numerically solve the boundary fluxes of the exact same one-frequency harmonic synthetic surface. We compare the analytical solution with the numerical solution using various quantitative metrics to evaluate the accuracy of the spectral method. As mentioned above, if Wörman et al.’s (2006) approach is exact, we should expect that the analytical solution should be very close to the numerical solution.

We used six synthetic surface scenarios with different complexities to investigate the influence of the number of Fourier fitting frequencies on approximating those synthetic surfaces. The complexity of those synthetic surfaces is represented by the number of frequencies in the x- and y-direction. From Figure 5, two interesting phenomena are shown. First, the accuracy of approximating the synthetic surfaces using the traditional Fourier fitting process does not always increase with the increasing number of incorporated Fourier fitting frequencies. From Figure 5, we clearly see the accuracy of approximating synthetic surface decreases with the increasing number of incorporated Fourier fitting frequencies in approximating less complex synthetic surface scenarios such as synthetic surfaces with one, two and four frequencies in the x- and y-direction. However, the accuracy of approximating synthetic surface increases with the increasing number of incorporated Fourier fitting frequencies when we approximate more complicated synthetic surfaces such as surfaces with 8, 16 and 32 frequencies in the x- and y-direction. This indicates that the influence of the number of Fourier fitting frequencies on approximating surface topography varies and depends on shapes and complexities of the surfaces. In other words, we hypothesize that there is a certain complexity threshold of the surface that controls the relationship between the number of incorporated Fourier fitting frequencies and Fourier fitting process accuracy: the increasing number of Fourier fitting frequencies can bring more accuracies only when approximating surfaces with a larger complex index than that specific complexity threshold. Second, the absolute mean error of approximating synthetic surfaces decreases with the increasing complexity of the synthetic surfaces. This indicates that the traditional Fourier fitting approach is more effective in approximating more complicated surface signals.

From these analyses, the conclusion is that Wörman et al.’s (2006) analytical approach can provide an exact solution under certain assumptions, indicating that the current hydrological model has the theoretical foundation to be further improved.

A simple synthetic example was used to compare the surface and boundary fluxes approximation accuracy between the traditional Fourier fitting approach and our new proposed approach. We first created a synthetic surface, which is composed of a pre-defined sinusoidal function (characteristic part) and some red noises (noise part) to imitate a random natural bedform surface. We first preprocessed our synthetic surface data and applied 2D DFT to transform the surface data from the Cartesian domain to the frequency domain. Then, a set of random topographic surfaces were created with the overall statistical properties, the same as the real synthetic surface, but without a concentration of variance into any particular frequency bands (Perron et al. 2008). We used a confidence level (⁠⁠) at which we can reject the null hypothesis that an observed periodic signal has occurred by chance in a random topographic surface to identify those significant frequencies of the real synthetic surface. Given the background power spectrum of those random fractal surfaces, this confidence level corresponds to the probability that the spectrum of a random surface will exceed the background by the observed amount at that frequency. This confidence level was used to identify the significant frequencies of the surface, retrieved those significant frequencies and used inverse Fourier transform to reconstruct the characteristic part of the original synthetic surface. Meanwhile, COMSOL Multiphysics was used to numerically simulate the boundary fluxes corresponding to the characteristic part of the surface and the noise part of the surface. We verified that the total boundary fluxes of the whole surface can be perfectly summed up by those two components of the boundary fluxes corresponding to the two components of the surface. For the characteristic part of the surface, the Fourier fitting process was skipped because we already have the exact mathematical expression and thus can calculate a very accurate analytical solution of the boundary fluxes, capturing the characteristic hyporheic exchange patterns; for the noise part of the surface, the traditional Fourier fitting method was used to approximate the noise part of the surface as well as the boundary fluxes. We summed up these two boundary flux components to form the integrated result. Figure 7 indicates that this new proposed approach successfully improved the approximation accuracy of both surface and boundary fluxes (red line is lower than the blue line in both graphs). It should be noted that although the improvement is not tremendous from the absolute value standpoint, it does decrease the mean error around 10% from a perspective of percentage for both fitting surface and boundary flux predictions. Thus, we believe that this new proposed approach has the theoretical foundation to be effective when dealing with natural bedform surfaces.

From the results and analysis above, we preliminarily verified our hypothesis that spectral analysis can help improve the accuracy of Wörman et al.’s (2006) traditional Fourier fitting approach in both fitting surfaces and predicting boundary fluxes. We also want to illustrate that the limitation of this study is that we only examine the validity of our new Fourier fitting strategy with very simple synthetic surface scenarios. While the validity of our new mathematical approach in real natural surface scenarios needs to be further investigated, we believe that synthetic surface scenarios can verify the validity of this new approach in this case since the proposed new fitting process might significantly improve the fitting accuracy in one natural surface but might have very trivial improvement in another natural surface; however, that does not necessarily indicate the invalidity of this new mathematical model. Besides, since the natural surface is extremely varied and very hard to control, we choose the synthetic surfaces which are very easy to control to illustrate and highlight our point that the traditional Fourier fitting approach's performance varies and depends on shapes and complexities of the targeting surfaces with the idea that our new proposed approach in this work can address this problem. As a result, although this study only examined very simple synthetic scenarios, we believe that this new mathematical approach is of significance to the hydrological modeling community in terms of all kinds of scales and purposes. The future work incorporates this new mathematical approach in hydrological models of all kinds of scales and further verifies its effectiveness in improving the validity and accuracy of the model's output and its implications for various purposes. Besides, the connection between the effectiveness of this mathematical approach and the geomorphology of real natural surfaces is also a research gap that needs to be further investigated and addressed.

We believe that the spectral analysis-incorporated new mathematical approach proposed in this study has the foundation to improve the accuracy and validity of hydrological exchange models of all scales: from river corridor to watershed. To the best of our knowledge, this new proposed method is the first systematical spectral-analysis integrated approach that incorporates the analysis of information on significant bedform features and has the ability to provide more accurate first-order approximation of hyporheic exchange information under various complex streambed topographies either with characteristic scales or not. Our finding aims to break through basic science and modeling barriers in river hydrological science and serves as the basis to understand larger scale drivers of change with features of energy-water systems that modulate the influence of these drivers of hyporheic exchange flux and their biogeochemical implications. The proposed new mathematical approach is expected to produce more accurate and reliable results for global-scale researchers who are interested in hyporheic exchange modeling and its implications.

In this study, a new Fourier fitting strategy that incorporates the spectral analysis as a tool to identify the characteristic scale of the surface and then applies the traditional Fourier fitting approach only on the random part of the surface was proposed. We also used a simple synthetic surface to verify that this new fitting strategy will not only improve the surface approximation accuracy but also improve the prediction accuracy of boundary fluxes. Our overarching goal in this work is to highlight the validity of spectral analysis in identifying the significant signals of surface and improve the surface fitting accuracy, as well as calculating boundary fluxes in the hydrological model. We acknowledge that this new proposed fitting method might have very few effects or no improvements on some of the natural surfaces since not all surfaces in the natural world have characteristic scale, but this should not conceal the importance of this new proposed mathematical approach that can improve the fitting accuracy of those surfaces with characteristic scales and significant frequency signals. We conclude that our proposed mathematical approach in this study has important implications for improving the overall accuracy of the current hydrological exchange model. The specific findings in this study are listed here:

• (1)

  Wörman et al.’s (2006) analytical approach can provide an exact solution to boundary fluxes under the circumstance of known hydraulic head distribution which can be approximated by bedform surface topography.

• (2)

  The accuracy of approximating the synthetic surfaces using the traditional Fourier fitting process does not always increase with the increasing number of incorporated Fourier fitting frequencies. Instead, the influence of the number of incorporated Fourier fitting frequencies on approximating surface varies and depends on the shape and complexity of the surfaces to be approximated.

• (3)

  It is the absence of the bedform characteristic signals in the current Fourier fitting process that leads to the absence of the characteristic hyporheic exchange patterns in the current multiscale model system. The absence of the bedform characteristic signals can be retrieved and identified by spectral analysis.

• (4)

  Our new proposed approach is effective in improving the accuracy of representing the surface and predicting the hyporheic exchange flux. The new proposed method will be more efficient when dealing with regulated surface topography than random surface topography.

Data cannot be made publicly available; readers should contact the corresponding author for details.