Abstract

It is well known that land surface topography governs surface–groundwater interactions under some circumstances and can be separated in a Fourier-series spectrum that provides an exact analytical solution of both the surface and the underlying three-dimensional groundwater flows. We evaluate the performance of the current Fourier fitting process by testing on different scenarios of synthetic surfaces. We identify a technical gap and propose a new version of the approach which incorporates the spectral analysis method to help identify the statistically significant frequencies of the surface to guide the refinement and mesh. Our results show that spectral analysis is the method that can help improve the accuracy of representing the surface, thus further improving the accuracy of predicting the bedform-driven hyporheic exchange flows.

HIGHLIGHTS

  • Identify a research gap in the traditional Fourier fitting approach in the hyporheic exchange model.

  • Propose a new mathematical approach that incorporates spectral analysis to improve the accuracy of the Fourier fitting process.

  • Verify the validity of spectral analysis in the Fourier fitting process with synthetic scenarios.

INTRODUCTION

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.

METHOD

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.

Test the performance of Wörman et al.’s (2006) approach: 3D analytical solution of the surface–groundwater interaction flow for a block of sediment

In Wörman et al.’s (2006) work, they solved the following Fourier series which satisfy the Laplace equation , where ∇ is the Nabla operator and h is the energy potential of the groundwater flow (hydraulic head) (m) as follows:
formula
(1)
where , is the amplitude coefficients (m), is the aerial arithmetic mean value, N is the number of wavelengths in the x- and y-direction, is the wave number, (m) is the wavelength and ε is the depth-to-impermeable surface (m) (Wörman et al. 2006). Here, a generally accepted relationship is adopted: the groundwater surface and its energy potential generally follow land surface topography, i.e. , where H is the phreatic surface and Z is the ground surface elevation (Wörman et al. 2006). Based on Darcy's law, the subsurface velocity field is easily obtained from the head distribution as , where q is the Darcy's velocity vector and K is the hydraulic conductivity. Streamlines and residence time distributions can be calculated from the flow field using the numerical particle tracking technique.

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).

To test the accuracy of Wörman et al.’s (2006) exact 3D analytical solution of the surface–groundwater interaction, we first create a very simple harmonic synthetic surface and apply Wörman et al.’s (2006) analytical solution on this synthetic surface to obtain an analytical solution of both the approximated surface and the boundary fluxes. The synthetic surface used here is a simple harmonic synthetic surface with one frequency in both x- and y-directions shown in Figure 1 as
formula
(2)
where the amplitude of the surface is set as 30 (m), the wavelength of the surface is set the same for both the x- and y-direction as (m), thus the wave number for both directions is . In addition, we also used COMSOL Multiphysics to numerically solve the boundary fluxes of this synthetic surface by taking it as the top boundary condition. It should be noted that we kept every variable the same in both analytical and numerical calculation scenarios with a subsurface flow and transport model configuration under homogeneous and isotropic circumstances similar to Ferencz et al. (2019) (density , porosity = 0.3 and hydraulic conductivity ). The hypothesis is that if Wörman et al.’s (2006) analytical solution is exact, the analytical solution should be very close to the numerical solution of the subsurface flow field obtained from COMSOL Multiphysics for this simple harmonic synthetic scenario.
Figure 1

The method/strategy adopted in the first step of this study: we first create a synthetic surface to test Wörman et al.’s (2006) analytical solution by comparing it with the numerically solved results in COMSOL Multiphysics.

Figure 1

The method/strategy adopted in the first step of this study: we first create a synthetic surface to test Wörman et al.’s (2006) analytical solution by comparing it with the numerically solved results in COMSOL Multiphysics.

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.

Figure 2

(a) A simple harmonic synthetic surface; (b) influence of the number of fitting frequencies used to approximate the real synthetic surface, the vertical coordinate is the mean error of the approximated surfaces; (c) we used eight fitting frequencies to fit the real synthetic surface, the red points are the real frequencies of the synthetic surface; and (d) we used 16 fitting frequencies to fit the real synthetic surface, the red points are the real frequencies of the synthetic surface.

Figure 2

(a) A simple harmonic synthetic surface; (b) influence of the number of fitting frequencies used to approximate the real synthetic surface, the vertical coordinate is the mean error of the approximated surfaces; (c) we used eight fitting frequencies to fit the real synthetic surface, the red points are the real frequencies of the synthetic surface; and (d) we used 16 fitting frequencies to fit the real synthetic surface, the red points are the real frequencies of the synthetic surface.

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.

Explore and test a modified new version of Wörman et al.’s (2006) approach

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?

Two-dimensional discrete Fourier transform

The two-dimensional (2D) discrete Fourier transform (DFT) of a data set consisting of measurements spaced at even intervals and can be written as (Priestley 1981):
formula
(3)
where and are the wave numbers in the x- and y-direction, and m and n are indices in the z array () (Perron et al. 2008). The element at ( in the DFT array corresponds to the two orthogonal frequency components:
formula
(4)
and the ranges of the wavenumbers are and (Perron et al. 2008). To collapse the 2D spectra into one-dimensional (1D) plots, we adopt the radial frequency concept used in Perron et al.’s (2008) work. The radial frequency is defined as
formula
(5)
The power spectrum is estimated by the DFT periodogram which represents how the variance of z varies with frequency:
formula
(6)

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

Preprocessing steps of topographic data: detrend and window function

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).

The Fourier transform also assumes that the input signal is periodic at the edges of the sampled interval (Perron et al. 2008). However, few signals have this property. As a result, it requires many noise (unreal) frequencies to describe the edge discontinuities and thus contaminate the power spectrum (Priestley 1981). This phenomenon is known as spectral leakage and can be mitigated by a window function, W, that has a maximum at its center and tapers smoothly to zero at its edges (Perron et al. 2008). There are many window functions that are suitable for different practical applications. We used a Hann (raised cosine) window in this case:
formula
(7)
formula
(8)
formula
(9)
formula
(10)
formula
(11)
Thus, for any 2D DFT applied by window function , the equation becomes
formula
(12)

Significance levels and significant signals identification

To use the power spectrum to evaluate the significance of these different signal components, we set up a confidence level where we can reject the null hypothesis that an observed periodic signal has occurred by chance in a random topographic surface, so that we can identify the significant frequency signals (Perron et al. 2008). The spectral power corresponding to a confidence level is (Gilman et al. 1963; Torrence & Compo 1998)
formula
(13)
where is the value at which the cumulative distribution function with two degrees of freedom equals . The background spectrum is obtained by creating sets of random topographic surfaces. These random topographic surfaces have the same overall statistical properties as the real topography (the same roughness of the surface – Hurst parameter H, the same total surface variances – the sum of the array), but without a concentration of variance into any particular frequency bands (Perron et al. 2008). We use these topographic data sets and calculate the corresponding power spectra using the same processing approach applied to the real topographic surface and averaged those series of random background surfaces to get the background spectrum (Perron et al. 2008). We set a significance level of 0.05 () to test whether the power at a given frequency exceeds the significance level calculated in the equation, thus identifying the significant surface signals (characteristic spatial scales) (Perron et al. 2008). The 1D power spectra plot helps us to discriminate the differences between the surface signals from the generated background signals, thus answering the question of whether the surface has the characteristic spatial signals. The 2D power spectra plot helps us to identify those characteristic signals of the surface. After obtaining the significant signals of the surface, we use an inverse Fourier transform process to reconstruct the characteristic spatial parts of the surface as well as the background noise part, thus decomposing the surface topography. Figure 3 shows the generic steps of conducting spectral analysis on a surface topography.
Figure 3

The basic steps of applying spectral analysis on surface topography to identify the significant signals of the surface and decompose the surface into characteristic spatial scales and random (background) portion before conducting the Fourier fitting process to approximate the surface.

Figure 3

The basic steps of applying spectral analysis on surface topography to identify the significant signals of the surface and decompose the surface into characteristic spatial scales and random (background) portion before conducting the Fourier fitting process to approximate the surface.

Fourier fitting process based on the results of the spectral analysis

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.

Figure 4

Here, we show an example of the new version of the Fourier fitting method by incorporating the spectral analysis method to first determine whether the surface has characteristic spatial scales. If the surface is a random one which has no significant signals, we directly apply the Fourier fitting process to approximate the surface. If the surface has characteristic spatial scales, we first use the spectral analysis to identify its characteristic spatial scales and use inverse Fourier transform to reconstruct and decompose the surface into two parts: characteristic spatial scales part and random (noise) part. We hypothesize that by evaluating the two parts separately, we get more accurate results in both approximating the real surface and predicting the boundary fluxes.

Figure 4

Here, we show an example of the new version of the Fourier fitting method by incorporating the spectral analysis method to first determine whether the surface has characteristic spatial scales. If the surface is a random one which has no significant signals, we directly apply the Fourier fitting process to approximate the surface. If the surface has characteristic spatial scales, we first use the spectral analysis to identify its characteristic spatial scales and use inverse Fourier transform to reconstruct and decompose the surface into two parts: characteristic spatial scales part and random (noise) part. We hypothesize that by evaluating the two parts separately, we get more accurate results in both approximating the real surface and predicting the boundary fluxes.

RESULTS AND DISCUSSION

The performance of the Fourier fitting process and the accuracy of Wörman et al.’s (2006) 3D analytical solution

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.

Surface comparison

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.

Figure 5

A series of synthetic examples shows the influence of the traditional manner of directly increasing the number of fitting frequencies to approximate surface topography varies and depends on shapes and complexities of the surfaces. (N is the number of frequencies of synthetic surfaces in both x- and y-directions, representing the complexities of synthetic surfaces.)

Figure 5

A series of synthetic examples shows the influence of the traditional manner of directly increasing the number of fitting frequencies to approximate surface topography varies and depends on shapes and complexities of the surfaces. (N is the number of frequencies of synthetic surfaces in both x- and y-directions, representing the complexities of synthetic surfaces.)

Boundary fluxes comparison (analytical vs. numerical)

Here, the quantitative comparison between the boundary fluxes calculated by both the analytical solution and numerical solution in COMSOL Multiphysics is shown. From Figure 6(a), we see that the analytical solution captures the numerical solution very well. We also show the spatial distribution of the mean error of the boundary fluxes (panel (b)). It indicates that for the most part of the location, the mean error is less than 5% which is small. The mean error is calculated as:
formula
(14)
where A is the analytical solution value and N is the numerical solution value of the boundary fluxes, respectively. The boxplot in Figure 6(c) shows that the analytical solution of the boundary fluxes has a very similar distribution to the numerical solution of the boundary fluxes. Figure 6(d) shows that the mean error of the analytical solution along the x-direction is generally less than 5%, corresponding to panel (b) again.
Figure 6

(a) The boundary fluxes comparison between the analytical solution and the numerical solution; (b) the spatial distribution of the mean error of the analytical solution of the boundary fluxes using the numerical solution as the benchmark; (c) the boxplot of the boundary fluxes of two solutions (all position); and (d) the boxplot of the mean error of the analytical solution according to the x-direction.

Figure 6

(a) The boundary fluxes comparison between the analytical solution and the numerical solution; (b) the spatial distribution of the mean error of the analytical solution of the boundary fluxes using the numerical solution as the benchmark; (c) the boxplot of the boundary fluxes of two solutions (all position); and (d) the boxplot of the mean error of the analytical solution according to the x-direction.

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.

Surface and boundary fluxes comparison between traditional Fourier fitting process and the new proposed approach

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.

Figure 7

A simple synthetic example shows that spectral analysis can identify the characteristic signals of the bedform surface and improve the prediction accuracy of hydraulic head and boundary fluxes. (a) The influence of the number of fitting frequencies on approximating surface by the original Fourier fitting process and the new proposed approach; (b) the influence of the number of fitting frequencies on predicting boundary fluxes by the original Fourier fitting process and the new proposed approach.

Figure 7

A simple synthetic example shows that spectral analysis can identify the characteristic signals of the bedform surface and improve the prediction accuracy of hydraulic head and boundary fluxes. (a) The influence of the number of fitting frequencies on approximating surface by the original Fourier fitting process and the new proposed approach; (b) the influence of the number of fitting frequencies on predicting boundary fluxes by the original Fourier fitting process and the new proposed approach.

Discussion

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.

CONCLUSIONS

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 AVAILABILITY STATEMENT

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

REFERENCES

Alexander
R. B.
,
Böhlke
J. K.
,
Boyer
E. W.
,
David
M. B.
,
Harvey
J. W.
,
Mulholland
P. J.
,
Seitzinger
S. P.
,
Tobias
C. R.
,
Tonitto
C.
&
Wollheim
W. M.
2009
Dynamic modeling of nitrogen losses in river networks unravels the coupled effects of hydrological and biogeochemical processes
.
Biogeochemistry
93
(
1–2
),
91
116
.
Battin
T. J.
,
Kaplan
L. A.
,
Findlay
S.
,
Hopkinson
C. S.
,
Marti
E.
,
Packman
A. I.
,
Newbold
J. D.
&
Sabater
F.
2008
Biophysical controls on organic carbon fluxes in fluvial networks
.
Nature Geoscience
1
(
2
),
95
100
.
Boano
F.
,
Camporeale
C.
,
Revelli
R.
&
Ridolfi
L.
2006
Sinuosity-driven hyporheic exchange in meandering rivers
.
Geophysical Research Letters
33
(
18
),
L18406
.
Boano
F.
,
Revelli
R.
&
Ridolfi
L.
2007
Bedform-induced hyporheic exchange with unsteady flows
.
Advances in Water Resources
30
(
1
),
148
156
.
Boano
F.
,
Revelli
R.
&
Ridolfi
L.
2009
Quantifying the impact of groundwater discharge on the surface–subsurface exchange
.
Hydrological Processes: An International Journal
23
(
15
),
2108
2116
.
Boano
F.
,
Revelli
R.
&
Ridolfi
L.
2013
Modeling hyporheic exchange with unsteady stream discharge and bedform dynamics
.
Water Resources Research
49
(
7
),
4089
4099
.
Boano
F.
,
Harvey
J. W.
,
Marion
A.
,
Packman
A. I.
,
Revelli
R.
,
Ridolfi
L.
&
Wörman
A.
2014
Hyporheic flow and transport processes: mechanisms, models, and biogeochemical implications
.
Reviews of Geophysics
52
(
4
),
603
679
.
Boulton
A. J.
,
Findlay
S.
,
Marmonier
P.
,
Stanley
E. H.
&
Valett
H. M.
1998
The functional significance of the hyporheic zone in streams and rivers
.
Annual Review of Ecology and Systematics
29
(
1
),
59
81
.
Cardenas
M. B.
&
Wilson
J.
2006
The influence of ambient groundwater discharge on exchange zones induced by current–bedform interactions
.
Journal of Hydrology
331
(
1–2
),
103
109
.
Cardenas
M. B.
&
Wilson
J. L.
2007a
Dunes, turbulent eddies, and interfacial exchange with permeable sediments
.
Water Resources Research
43
(
8
),
W08412
.
Cardenas
M. B.
&
Wilson
J. L.
2007b
Exchange across a sediment–water interface with ambient groundwater discharge
.
Journal of Hydrology
346
(
3–4
),
69
80
.
Dudley Williams
D.
&
Hynes
H.
1974
The occurrence of benthos deep in the substratum of a stream
.
Freshwater Biology
4
(
3
),
233
256
.
Elliott
A. H.
1990
Transfer of Solutes Into and out of Streambeds
.
PhD Thesis
,
California Institute of Technology
,
Pasadena
,
USA
.
Elliott
A. H.
&
Brooks
N. H.
1997a
Transfer of nonsorbing solutes to a streambed with bed forms: laboratory experiments
.
Water Resources Research
33
(
1
),
137
151
.
Elliott
A. H.
&
Brooks
N. H.
1997b
Transfer of nonsorbing solutes to a streambed with bed forms: theory
.
Water Resources Research
33
(
1
),
123
136
.
Fehlman
H. M.
1985
Resistance Components and Velocity Distributions of Open Channel Flows Over bed-Forms
.
PhD Thesis
,
Colorado State University
,
Fort Collins
,
USA
.
Gilman
D. L.
,
Fuglister
F. J.
&
Mitchell
J. M.
Jr.
1963
On the power spectrum of ‘red noise’
.
Journal of the Atmospheric Sciences
20
(
2
),
182
184
.
Gomez
J. D.
,
Wilson
J. L.
&
Cardenas
M. B.
2012
Residence time distributions in sinuosity-driven hyporheic zones and their biogeochemical effects
.
Water Resources Research
48
(
9
),
W09533
.
Harvey
J. W.
,
Böhlke
J. K.
,
Voytek
M. A.
,
Scott
D.
&
Tobias
C. R.
2013
Hyporheic zone denitrification: controls on effective reaction depth and contribution to whole-stream mass balance
.
Water Resources Research
49
(
10
),
6298
6316
.
Jones
J. B.
&
Mulholland
P. J.
1999
Streams and Ground Waters
.
Elsevier
,
Netherlands
.
Kennedy
C. D.
,
Genereux
D. P.
,
Corbett
D. R.
&
Mitasova
H.
2009
Spatial and temporal dynamics of coupled groundwater and nitrogen fluxes through a streambed in an agricultural watershed
.
Water Resources Research
45
(
9
),
W09401
.
Krause
S.
,
Heathwaite
L.
,
Binley
A.
&
Keenan
P.
2009
Nitrate concentration changes at the groundwater-surface water interface of a small Cumbrian river
.
Hydrological Processes: An International Journal
23
(
15
),
2195
2211
.
Marzadri
A.
,
Vignoli
G.
,
Bellin
A.
&
Tubino
M.
2005
Effect of bar topography on hyporheic flow in a gravel-bed river
. In
River, Coastal and Estuarine Morphodynamics: Proceedings of the 4th IAHR Symposium on River, Coastal and Estuarine Morphodynamics
.
CRC Press
,
Urbana, IL
,
USA
.
McClain
M. E.
,
Boyer
E. W.
,
Dent
C. L.
,
Gergel
S. E.
,
Grimm
N. B.
,
Groffman
P. M.
,
Hart
S. C.
,
Harvey
J. W.
,
Johnston
C. A.
,
Mayorga
E.
,
McDowell
W. H.
&
Pinay
G.
2003
Biogeochemical hot spots and hot moments at the interface of terrestrial and aquatic ecosystems
.
Ecosystems
6
(
4
),
301
312
.
Mojarrad
B.
,
Riml
J.
,
Wörman
A.
&
Laudon
H.
2019
Fragmentation of the hyporheic zone due to regional groundwater circulation
.
Water Resources Research
55
(
2
),
1242
1262
.
Mulholland
P. J.
,
Helton
A. M.
,
Poole
G. C.
,
Hall
R. O.
,
Hamilton
S. K.
,
Peterson
B. J.
,
Tank
J. L.
,
Ashkenas
L. R.
,
Cooper
L. W.
,
Dahm
C. N.
,
Dodds
W. K.
,
Findlay
S. E. G.
,
Gregory
S. V.
,
Grimm
N. B.
,
Johnson
S. L.
,
McDowell
W. H.
,
Meyer
J. L.
,
Valett
H. M.
,
Webster
J. R.
,
Arango
C. P.
,
Beaulieu
J. J.
,
Bernot
M. J.
,
Burgin
A. J.
,
Crenshaw
C. L.
,
Johnson
L. T.
,
Niederlehner
B. R.
,
O'Brien
J. M.
,
Potter
J. D.
,
Sheibley
R. W.
,
Sobota
D. J.
&
Thomas
S. M.
2008
Stream denitrification across biomes and its response to anthropogenic nitrate loading
.
Nature
452
(
7184
),
202
205
.
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
Riparian Areas: Functions and Strategies for Management
.
National Academies Press
,
Washington, DC
,
USA
.
Packman
A. I.
&
Brooks
N. H.
2001
Hyporheic exchange of solutes and colloids with moving bedforms
.
Water Resources Research
37
(
10
),
2591
2605
.
Perron
J. T.
,
Kirchner
J. W.
&
Dietrich
W. E.
2008
Spectral signatures of characteristic spatial scales and nonfractal structure in landscapes
.
Journal of Geophysical Research: Earth Surface
113
(
F4
),
F04003
.
Priestley
M. B.
1981
Spectral Analysis and Time Series
.
Elsevier
,
Netherlands
.
Sawyer
A. H.
,
Shi
F.
,
Kirby
J. T.
&
Michael
H. A.
2013
Dynamic response of surface water-groundwater exchange to currents, tides, and waves in a shallow estuary
.
Journal of Geophysical Research: Oceans
118
(
4
),
1749
1758
.
Singh
T.
,
Wu
L.
,
Gomez-Velez
J. D.
,
Lewandowski
J.
,
Hannah
D. M.
&
Krause
S.
2019
Dynamic hyporheic zones: exploring the role of peak flow events on bedform-induced hyporheic exchange
.
Water Resources Research
55
(
1
),
218
235
.
Stonedahl
S. H.
,
Harvey
J. W.
,
Wörman
A.
,
Salehin
M.
&
Packman
A. I.
2010
A multiscale model for integrating hyporheic exchange from ripples to meanders
.
Water Resources Research
46
(
12
),
W12539
.
Thibodeaux
L. J.
&
Boyle
J. D.
1987
Bedform-generated convective transport in bottom sediment
.
Nature
325
(
6102
),
341
343
.
Tonina
D.
&
Buffington
J. M.
2009
Hyporheic exchange in Mountain Rivers. I: mechanics and environ-mental effects
.
Geography Compass
3
(
3
),
1063
1086
.
Torrence
C.
&
Compo
G. P.
1998
A practical guide to wavelet analysis
.
Bulletin of the American Meteorological Society
79
(
1
),
61
78
.
Toth
J.
1963
A theoretical analysis of groundwater flow in small drainage basins
.
Journal of Geophysical Research
68
(
16
),
4795
4812
.
Vaux
W.
1968
Intragravel flow and interchange of water in a streambed
.
Fishery Bulletin
66
(
3
),
479
489
.
Ward
A. S.
&
Packman
A. I.
2019
Advancing our predictive understanding of river corridor exchange
.
Wiley Interdisciplinary Reviews: Water
6
(
1
),
e1327
.
Weedon
G. P.
2005
Time-series Analysis and Cyclostratigraphy: Examining Stratigraphic Records of Environmental Cycles
.
Cambridge University Press
,
Cambridge
,
UK
.
Winter
T. C.
1998
Ground Water and Surface Water: A Single Resource
.
DIANE Publishing Inc.
,
PA, USA
.
Wondzell
S. M.
&
Swanson
F. J.
1996
Seasonal and storm dynamics of the hyporheic zone of a 4th-order mountain stream. I: hydrologic processes
.
Journal of the North American Benthological Society
15
(
1
),
3
19
.
Wörman
A.
,
Packman
A. I.
,
Marklund
L.
,
Harvey
J. W.
&
Stone
S. H.
2006
Exact three-dimensional spectral solution to surface-groundwater interactions with arbitrary surface topography
.
Geophysical Research Letters
33
(
7
),
L07402
.
Wörman
A.
,
Packman
A. I.
,
Marklund
L.
,
Harvey
J. W.
&
Stone
S. H.
2007
Fractal topography and subsurface water flows from fluvial bedforms to the continental shield
.
Geophysical Research Letters
34
(
7
),
L07402
.
Wroblicky
G. J.
,
Campana
M. E.
,
Valett
H. M.
&
Dahm
C. N.
1998
Seasonal variation in surface-subsurface water exchange and lateral hyporheic area of two stream-aquifer systems
.
Water Resources Research
34
(
3
),
317
328
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).