## 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

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

*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 aswhere 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.

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.

### 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

*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: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

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

*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:

#### Significance levels and significant signals identification

*et al.*2008). The spectral power corresponding to a confidence level is (Gilman

*et al.*1963; Torrence & Compo 1998)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.

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

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

#### Boundary fluxes comparison (analytical vs. numerical)

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

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.

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