Propagation of hydropeaking waves in heterogeneous aquifers: effects on flow topology and uncertainty quantification

Topological flow properties are proxies for mixing processes in aquifers and allow us to better understand the mechanisms controlling transport of solutes in the subsurface. However, topological descriptors, such as the Okubo–Weiss metric, are affected by the uncertainty in the solution of the flow problem. While the uncertainty related to the heterogeneous properties of the aquifer has been widely investigated in the past, less attention has been given to the one related to highly transient boundary conditions. We study the effect of different transient boundary conditions associated with hydropeaking events (i.e., artificial river stage fluctuations due to hydropower production) on groundwater flow and the Okubo–Weiss metric. We define deterministic and stochastic modeling scenarios applying four typical settings to describe river stage fluctuations during hydropeaking events: a triangular wave, a sine wave, a complex wave that results of the superposition of two sine waves, and a trapezoidal wave. We use polynomial chaos expansions to quantify the spatiotemporal uncertainty that propagates into the hydraulic head in the aquifer and the Okubo–Weiss. The wave-shaped highly transient boundary conditions influence not only the magnitude of the deformation and rotational forces of the flow field but also the temporal dynamics of dominance between local strain and rotation properties. Larger uncertainties are found in the scenario where the trapezoidal wave was imposed due to sharp fluctuation in the stage. The statistical moments that describe the propagation of the uncertainty highly vary depending on the applied boundary condition. Deterministic and stochastic scenarios to describe the groundwater flow field under river stage fluctuations during hydropeaking. Propagation of uncertainty of highly transient boundary conditions in the Okubo–Weiss metric. Highly transient boundary conditions can significantly affect mixing potential. Deterministic and stochastic scenarios to describe the groundwater flow field under river stage fluctuations during hydropeaking. Propagation of uncertainty of highly transient boundary conditions in the Okubo–Weiss metric. Highly transient boundary conditions can significantly affect mixing potential.


Introduction
Mixing plays a critical role in describing solute transport in aquifers (de Anna et al. 2014b; Rolle and Le Borgne 2019). Understanding mixing-limited reactions in the subsurface is particularly relevant to recognize biogeochemical transformations (Boisson et al. 2013;Kang et al. 2019;Pinay et al. 2015), operate and design engineering remediation techniques (Cho et al. 2019;Dentz and Carrera 2005;Mays and Neupauer 2012;McCarty and Criddle 2012;Neupauer et al. 2014), and understand hyporheic processes (Boano et al. 2014). A variety of methods have been developed to understand the transport dynamics in the subsurface and the effect that heterogeneous hydraulic properties have on spreading, dilution, and reactive mixing (Dentz and Carrera 2005;Valocchi et al. 2019). However, transport simulations are often computationally expensive and the quantification of uncertainties may result in a very time-consuming exercise (Lykkegaard et al. 2021;Smith 2013). The relation between topological flow properties and mixing processes in aquifers de Barros et al. 2012) and porous media de Anna et al. 2014a;Engdahl et al. 2014;Wright et al. 2017) provides an interesting alternative to the solution of the transport problem. One advantage of investigating such relations is that the calculation of topological features of the flow field requires only the solution of the flow problem, which is much cheaper from the computational point of view than the solution of the flow and transport equations.
A topological quantity known as the Okubo-Weiss metric (Okubo 1970;Weiss 1991) was shown to be a good proxy for mixing potential de Barros et al. 2012;Wright et al. 2017). This metric is commonly used in geophysics to identify filament from vortex structures (Casella et al. 2011;Roullet and Klein 2010) and characterize them in terms of dominant forces of the flow field, such as vorticity, shear strain, and normal strain (de Barros et al. 2012;Wallace et al. 2021). Still, the quantification of such topological descriptors of the flow field is affected by the uncertainty that is caused by the heterogeneous nature of the aquifer (Geng et al. 2020;Valocchi et al. 2019). Significant efforts have been made to quantify the uncertainty affecting the predictions of solute concentration values caused by the generally unknown hydraulic conductivity field (Moslehi and de Barros 2017;Nowak et al. 2010). However, in this work, we assume the hydraulic conductivity field as properly characterized and well known in order to focus on a different source of uncertainty, which did not receive comparable attention in the literature, i.e., highly transient boundary conditions.
In fact, transient boundary conditions can also be uncertain and affect the estimation of the topological properties of the flow field and, consequently, the understanding of mixing and transport processes in aquifers (Hester et al. 2021;Ziliotto et al. 2021). This is particularly the case of the aquifer area close to surface water bodies (Dudley-Southern and Binley 2015;Merchán-Rivera et al. 2021;Santizo et al. 2020;Singh et al. 2020). In this work, we focus on a river reach that is affected by hydropeaking, i.e., sudden changes in the hydraulic head of the river caused by the operation of hydropower plants. Such fluctuations display some typical periodicity (Pérez Ciria et al. 2020) and modify the natural hydrological behavior and hydraulic conditions of the streams (Hauer et al. 2017;Meile et al. 2011), which can impact the hyporheic zone (Sawyer et al. 2009;Singh et al. 2019) and propagate to the groundwater (Francis et al. 2010;Song et al. 2020). Moreover, since hydropeaking may depend on hydrological conditions (Li and Pasternack 2021) and the dynamic behavior of the energy market (Chiogna et al. 2018;Pérez Ciria et al. 2019;Wagner et al. 2015), the stream head fluctuations entail uncertainty related to the peak amplitude and the temporal occurrence of the event.
The question that we aim at answering in this work is to what extent the shape and the uncertainty of hydropeaking waves affect the topology of the groundwater flow field quantified through the Okubo-Weiss parameter. To achieve our aim, we define one single realization of a two-dimensional heterogeneous aquifer and build modeling scenarios based on four typical settings for the stream fluctuations of the boundary conditions: a triangular wave, a sine wave, a complex wave (realized as the superposition of two sine waves), and a trapezoidal wave (Ferencz et al. 2019;Li and Pasternack 2021;Sawyer et al. 2009). Moreover, we apply polynomial chaos expansion (Xiu and Karniadakis 2002) to quantify the uncertainty due to the oscillatory boundaries and quantify the mean and standard deviation of the temporal and spatial values of the Okubo-Weiss metric.
The paper is structured as follows. Section 2 presents the synthetic case study, the deterministic and stochastic modeling scenarios, the polynomial chaos expansion method, and the topological metric that we use to describe the flow field. In Sect. 3, we present and discuss the results and findings related to the groundwater flow and the Okubo-Weiss metric and the quantification of the spatiotemporal uncertainty. We conclude this work in Sect. 4 by restating major findings and discussing the environmental implications of our results.

Groundwater flow equation
The governing equation of the two-dimensional transient groundwater flow in a heterogeneous, isotropic, and unconfined aquifer can be written as where h is the hydraulic head [L], K x 1 and K x 2 are values of the hydraulic conductivity along the x 1 and x 2 coordinate axis [LT −1 ], s y is the specific yield of the porous medium [-], and Q describes the volumetric flux from source and sink terms [LT −1 ] (Anderson et al. 2015;Bear 1979). Dirichlet and Neumann boundary conditions can be denoted as respectively, where g(x, t) is a continuous function, {x, t} are spatiotemporal dependencies, and n is the outward normal of the boundary. While Dirichlet boundaries D are used to define hydraulic heads h from the elevation of the water level [L], the normal derivative of the head in the Neumann boundary conditions N is used to define a known specific discharge [LT −1 ], including no-flow boundaries (Bear and Cheng 2010;Cheng and Cheng 2005;Liu 2018).

Model description
In our system, we consider a two-dimensional unconfined aquifer with lognormal heterogeneous isotropic hydraulic conductivity field (x) = ln[K (x)] defined by the geometric mean μ = 1 × 10 −3 m/s, the geometric standard deviation σ = 1.5 m/s and the correlation length λ = 10 m, equivalent to porous medium formed by sands and gravels (Coduto 1999). The area of the squared domain is L 1 × L 2 = 10λ × 10λ with cell size 0.1λ × 0.1λ (see Fig. 1b). The distance between the ground surface and the bottom of the aquifer is z = d(−z/2, z/2), the reference datum is the middle point z 0 = 0. The initial groundwater level conditions h 0 are set to a uniform water level in the domain, which matches the reference datum, so that h 0 = z 0 (see Fig. 1c) and any simulation output h represents the relative movement of the groundwater head with respect to z 0 . The dimensionless time is defined as τ/T , where T represents the simulation time and τ is the size of a single time step that corresponds to 1/(80 f ), where f is the wave frequency.
The discontinuous release of water from hydropeaking events tends to follow a periodicity. Hence, the stage fluctuation in the stream can be described by a periodic function y(t + F) = y(t), where F = 1/ f is a nonzero value defined as the period [T]. We introduce the periodic waves as specified-head boundary conditions, i.e., Dirichlet boundary conditions (Bear and Cheng 2010), by fixing head values at the left border of the domain. The wave-shaped time series of heads represent the stream stage and the values at the boundary are time-dependent and updated as the simulation progresses. This setup assumes that the stage changes for the entire river reach considered simultaneously (i.e., no hydraulic model is used to describe the wave propagation along the reach) and strong connectivity between surface water and groundwater. These assumptions are consistent, for example, with the study of Sawyer et al. (2009) focusing on the Colorado River in Texas. The top, bottom and right boundaries are defined as no-flow boundaries. We define the following periodic functions (see Fig. 1a) to represent the effect of the hydropeaking on the left boundary conditions: • Triangular wave: • Sine wave • Complex wave: • Trapezoidal wave Then, y V (t), y S (t), y C (t) and y Z (t) represent the time-variant specified-head boundaries for the triangular, sine, complex and trapezoidal wave scenario, respectively, t is the evaluation time step, A is the amplitude of the wave, p is the phase shift and f is the frequency. We shaped the trapezoidal wave using the stepwise function y Z (t) to remove the values that exceed the minimum and maximum limits defined by A, which depends on the outcomes of a triangular function y * V (t). A coefficient b Z is introduced to this triangular function to extend the amplitude A, which will determine the interval between minimum and maximum stage to be equal to 4/b Z given the slope of the trapezoid legs with base angles θ = arctan(4b Z A f ). The complex waveform is the result of combining two different sine waves (y α S (t) and y β S (t)) shaped by two different amplitudes (A α and A β ), two phase shifts ( p α and p β ) and two frequencies ( f α and f β ).

Deterministic problem
As starting point, we want to observe the responses of the aquifer and the flow field topology assuming that all model inputs are known. Hence, there are four deterministic scenarios of one single deterministic transient flow problem, each applying one of the wave-shaped boundary conditions. The complex wave is defined by A β = 2 A α /5, p β = 3 p α and f β = 3 f α to satisfy symmetry relations between the periodic waves to match maximum (i.e., wave crest), minimum (i.e., wave through), temporal axis interceptions, and lag between two events. The trapezoidal wave is shaped considering b Z = 6, so that the interval between minimum and maximum is equivalent to 1/12 f .
The spatial and temporal distribution of the groundwater heads and the Okubo-Weiss are analyzed with two-dimensional arrays. We are also interested in the variability of the groundwater head and the Okubo-Weiss metric at a certain distance x 1 /L 1 from the transient boundary conditions. Hence, we compute the expected value and variance in time at each discrete cell, and then the arithmetic mean of the expected value and the arithmetic mean of the variance relative to the distance x 1 /L 1 . For simplicity we call c i, j t the output of interest (i.e., groundwater heads or Okubo-Weiss metric) at a discrete cell with row i, and column j, at the time t, and follow where μ i, j c is the mean of the output of interest at a discrete cell {i, j} along all the evaluation time T . Then, the arithmetic mean of the variances at the distance x 1 /L 1 can be computed along column j ≈ x 1 /L 1 , such that, being n i the number of the discrete rows. Hence, from now on, the bar notation is used to represent the transverse spatial average of the different quantities of interest at a specific distance from the transient boundary conditions.

Stochastic problem
We consider the amplitude A and phase p as uncertain parameters and we assume them as mutually independent random variables. On one hand, A represents the maximum water level in the stream, which fluctuates due to the variable discharges from the power plant, which in turn depends on market and seasonal conditions. The uncertainty is defined from a minimum and maximum stage fluctuation that follows a uniform where a A = 0.9 and b A = 1.1. On the other hand, p represents the shift in the stage signal. This random variable then introduces the temporal uncertainty due to changes in the gate management, turbine control and discharge duration. The random variable p is uniformly distributed, such that p ∼ U(a p , b p ). The parameters a p and b p describe the phase difference and relates the offset with f using a factor o p , so that a p = −1/(o p f ) and b p = 1/(o p f ). This arrangement simplifies the application at multiple hydropeaking scale events (e.g., sub-daily, daily, and weekly) because any shift in the phase is a ratio of the periodicity. In our problem setup, we set o p = 8, which is equivalent to a phase difference of 1/8 f . In the case of the complex wave, the problem increases to 4 stochastic dimensions given that it is formed by the superposition of two sine waves y α S (t) and y β S (t). Like in the deterministic problem, we keep the proportional relations between the two waves y α S (t) and y β S (t) that form the complex wave and the random variables that are assumed mutually independent. Hence, two random variables In this work, we represent the dynamic system of the groundwater heads and flow topology as stochastic processes with uncertain boundary conditions using generalized polynomial chaos expansions. The statistical moments of the output of interest are estimated from the polynomial chaos coefficients using the pseudospectral collocation approach. Analogous to Eq. (10), the arithmetic mean is also used to analyze these statistical moments. A summary of the specific parameters used in the deterministic and stochastic scenarios can be found in Table 1.

Stochastic formulation
Let ( , F, P) be a probability space, where is a sample space, F is a σ -algebra on , and P is a probability measure on . Consider a function u(t, x, ) on the probability space ( , F, P), where = [ 1 , . . . , d ] : → R is a random vector with a finite set of d mutually independent random variables with marginal probability density functions ρ i (ϕ i ), i = 1, . . . , d , and {t, x} represent the deterministic temporal and spatial dependencies with a finite temporal horizon t ∈ [0, T ] within the spatial domain D ⊂ R 2 formed by fixed grid points x = (x 1 , x 2 ). Since each parameter i (ω) : → R is associated to a density ρ i (ϕ i ) and ω ∈ is a realization in the underlying probability space, we reformulate the problem in the image probability is the sample space for the range of i , B( ) is the Borel σ -algebra on , and ρ (ϕ) is the joint density associated with , described by The output function is then a random process u(t, x, ) : [0, T ] × D × → R with a finite variance. Following the generalized Cameron-Martin theorem (Cameron and Martin 1947), we can represent it as an infinite series expansion of polynomials, which can be truncated to order K , such that where u κ (t, x) are deterministic expansion coefficients, κ ( ) represent the multivariate orthogonal polynomial basis function, and κ = {κ 1 , . . . , κ d } ∈ N d 0 is a multi-index of non-negative integers of size d to identify the degree of the polynomials for the input variable i . To achieve n order of polynomials, K can be optimally defined by [(n + d)!/n!d!]−1 (Smith 2013;Xiu 2010), or by using experimental designs, such as the empirical rule K = (d − 1) × (n + 1) (Sudret 2008), to reduce the computational demand of the experiment. In this research, we define K = 9 and n = 3 for the triangular, sine and trapezoidal stochastic scenarios, leading to a total of 100 realizations for each scenario. The expansions in the complex scenario are calculated considering K = 4 and n = 2, which leads to 625 realizations.
The orthogonal basis κ must be accordingly specified to ρ i (ϕ i ) (Xiu and Karniadakis 2002). In this work, the uncertain input parameters associated to the boundary conditions are considered random variables uniformly distributed in the interval [a, b], denoted by i ∼ U(a, b). An appropriate basis is formed by the family of Legendre polynomials, which are an orthogonal basis with respect to the weight function ρ i (ϕ i ) = 1/2 for all normalized ϕ i ∈ [−1, 1]. Given the assumption that the random variables are mutually independent, the multivariate Legendre polynomial basis function κ ( ) can be defined as the tensor product of the associated univariate orthogonal polynomials ψ κ i , such that GEM -International Journal on Geomathematics (2022) 13:11 which satisfies orthonormality conditions given that where ·, · denotes the inner product of the sequence of two polynomials ψ κ i , ψ τ i of degree κ i and τ i in the i th variable, δ κ i τ i represents the Kronecker delta. The deterministic coefficients u κ (t, x) can be approximated by exploiting the orthonormality of the basis function and projecting u(t, x, ) onto each basis function k ( ) to obtain the representation

Pseudospectral collocation approach
In this section, we explain the use of the pseudospectral approach as a solution technique to estimate the multidimensional integral that describes u κ (t, x) and the derivation of the statistical moments of interest to describe the propagation of uncertainty. The pseudo-spectral approach is employed due to the possibility of decoupling the stochastic and deterministic dependencies of the problem, the relatively low computational effort and the flexibility of the method to apply different basis functions to represent different random variables, which can be particularly useful for the application of the method in real case scenarios. It is a discrete collocation method that relies on quadrature techniques to calculate u κ (t, x) at selected quadrature nodes g q = g q 1 1 , . . . , g q d d ∈ R defined on with the associated weights w q = w q 1 , . . . , w q d ∈ R. Therefore, the deterministic solvers that describe the groundwater flow and the topological responses of the aquifer are not modified (i.e., non-intrusive spectral projection) because it is only required to evaluate u(t, x, g q ) at the given g q . We employ Gaussian quadrature rules (Golub and Welsch 1968) over a full tensor product grid to distribute g q according to the probability density functions ρ i (ϕ i ). Using Q to represent the quadrature integration, the extension of the univariate Gaussian quadrature yields to the summation over all possible combinations over m i nodes that can be recast for the sake of simplicity to the multi-index approximation where the total number of grid points is M = m d given that we opt for m 1 = · · · = m d = m.
The expected value of u(t, x, ) can be estimated from the polynomial chaos coefficients u 0 (t, x), given that Similarly, the variance σ 2 = V[u(t, x, )] and the standard deviation σ = 2 √ V[u(t, x, )] are quantified following Finally, Sobol' indices can be formally obtained by exploiting the polynomial expansion representation to perform a global sensitivity analysis (Le Maitre and Knio 2010). The use of the expansion coefficients reduces significantly the computational burden that is associated to typical methods to quantify the indices (e.g., Monte Carlo schemes) because the variance decomposition can be obtained directly from the available expansion coefficients. We refer the reader to Formaggia et al. (2013) and Sudret (2008) for further detail related to the computation of polynomial chaos-based Sobol' indices.

Okubo-Weiss
We consider a flow deformation metric based on a two-dimensional velocity gradient tensor, where v(x 1 , x 2 , t) is the velocity at the space coordinates x 1 and x 2 at time t and it is estimated with the calculated head gradient by solving Eq. (1). The Okubo-Weiss function (Okubo 1970;Weiss 1991) is defined by which in the horizontal plane with coordinates x 1 and x 2 is written as We take the definition used by Okubo (1970) for stretching deformation α, vorticity ω, and shear deformation σ , to be By substituting Eq. (23) into Eq. (22), and following de Barros et al. (2012), in which for a two-dimensional transport scenario ∂v x 1 ∂ x 1 = − ∂v x 2 ∂ x 2 , and therefore stretching deformation α = 2 ∂v x 1 ∂ x 1 , the deformation tensor can be then rewritten as and the Okubo-Weiss function ξ [1/T 2 ] is calculated by Positive Okubo-Weiss values, ξ > 0, correspond to regions where shear and stretching forces dominate, and are associated to mixing hotspots (Engdahl et al. 2014;Wright et al. 2017). On the other hand, negative values, ξ < 0, correspond to regions dominated by vorticity and local mixing potential is low.

Algorithm implementation
We use MODFLOW-2005 as groundwater flow equation solver, and the time-variant specified-head boundary package (i.e., CHD package) was used to set the Dirichlet boundaries (Harbaugh 2005). The model was built using FloPy (Bakker et al. 2016) and the polynomial chaos expansion approach was implemented using the Chaospy library (Feinberg 2019). The Okubo-Weiss calculations, random field generator, wave functions, postprocessing scripts, and the code for coupling the polynomial expansions and the MODFLOW-2005 model were written in Python 3. The scripts and results are available in the online repository of the research (Merchán-Rivera et al. 2022). To validate the application of the gPCE, we include a comparison with results obtained from 1000 realizations from quasi Monte Carlo (MC) samples using Halton sequences (Halton 1964;Smith 2013) in the Supplementary Material of this document. Figure 2 shows the spatial effect of the waveform on the groundwater heads at specific time steps τ/T ∈ {60, 80, 86, 100}. Steeper hydraulic gradients are observed in the trapezoidal wave, which significantly impact the flow field magnitude. These gradients occur in the trapezoidal wave due to two reasons. First, this wave exposes longer intervals of constant head at the wave crest and wave through. Second, the trapezoidal wave presents a sharp fluctuation from minimum to maximum stage. Furthermore, the behavior of the sine scenario is very similar to the triangular one. As expected, the porous medium acts as a damper that gradually moderates the propagation of groundwater head signals, converting all of them to sinusoidal patterns after travelling a certain distance and later vanishing them (see Fig. 3). The dampening primary depends on the value of the diffusivity of the aquifer which depends on the hydraulic conductivity, specific yield, and saturated aquifer thickness, in accordance with the analytical solution for the head response in a semi-infinite aquifer presented by Singh (2004) and Sawyer et al. (2009) for the case of a homogeneous porous medium.

Deterministic scenarios
The outcomes from the deterministic scenarios show some remarkable differences in the mean variance of the groundwater heads and the Okubo-Weiss. In Fig. 4a, we observe a larger σ 2 h from the trapezoidal wave scenario, followed by the complex wave scenario. These large values of the mean variance can be explained by the spread from the mean that occurs when the groundwater heads rapidly fluctuate between minimum and maximum values. For instance, the trapezoidal boundary conditions introduce a quasi-bimodal forcing characterized by short periods with heads in between or close to the mean. A similar behavior with abrupt changes is also observed in the complex wave. In the sine and triangular wave, the spread is lower because the head values are closer to a uniform distribution in time. We also observe that the behavior of the four scenarios is very similar after x 1 /L 1 = 40 and that σ 2 h → 0 after x 1 /L 1 = 60. Figure 4b shows the mean variance in the results of the Okubo-Weiss metric. We see that ξ may vary by several orders of magnitude depending on the wave used as boundary condition. Similar to σ 2 h , we see larger variations in σ 2 ξ in the trapezoidal and complex scenarios.
We show the spatial and temporal pattern of the Okubo-Weiss metric in Fig. 5. We observe regions where the Okubo-Weiss metric is very high ξ > 0.1e −6 , when the GEM -International Journal on Geomathematics (2022) 13:11 Spread of deterministic results at different distances x 1 /L 1 from the boundary conditions over the whole simulation period: a mean variance of the groundwater heads σ 2 h , and b logarithm of the mean variance of the Okubo-Weiss σ 2 ξ flow is dominated by stretching and strain and very low ξ < −0.1e −6 , when vorticity dominates. These regions correspond to areas of high hydraulic conductivity and high conductivity contrasts (see Fig. 1b), where also flow focusing may occur. Hence, the location of these spots is fully controlled by the configuration of heterogeneous field in all the scenarios. In temporal terms, the most remarkable discrepancies in ξ among the four scenarios occur during the sharp ramp upwards and the sharp drop of the trapezoidal wave. In a lower magnitude, this is also visible in the complex wave. Highest positive and lower negative values of ξ are found in these two scenarios. Furthermore, it is possible to observe cells changing from positive ξ to negative ξ , and vice versa, in all the scenarios. This can be observed before and after the apexes of the trapezoidal wave (Fig. 5b, c), the peak of the triangular wave (Fig. 5c, d), the GEM -International Journal on Geomathematics (2022) 13:11 local maximum and local minimum of the complex wave (Fig. 5b, c). This swap of dominance is transitory and can be repeatedly observed in all the scenarios at critical points of the wave-shaped boundaries, such as stationarity points (i.e., constant value), inflection points, local maxima, and local minima. This could occur due to flow reversal caused by the deacceleration of the transient boundary signal into the aquifer. Overall, this behavior gives evidence of the waveform's role in the temporal dynamics of the topology of the flow field.

Stochastic scenarios
We summarize the results from the global sensitivity analysis in Fig. 6, which shows the total Sobol indices computed using the relations reported in Table 1b. The variations of the sensitivity along the simulation period indicate that the influence of amplitude and phase follow a periodic dynamic. The phase p has a primary influence in the triangular and sine wave. In the complex wave, p α has a primary influence, p β and A α have secondary influence, and A β has negligible influence. On the trapezoidal wave the influence of the A and p on the groundwater head outputs depends on the time step of evaluation. The latter occurs in the trapezoidal boundaries due to the low influence of p when the heads are constant, and the low influence of A during the periods of transition between the minimum and maximum values. The uncertainty in the amplitude and phase of the waves propagates in the groundwater head following different patterns (see Fig. 7). The results of μ h are similar to the results of h in the deterministic scenarios. Regarding the standard deviation, in Fig. 7a, we see at τ/T = 60 that all scenarios present similar snapshots, with slightly higher σ h close to the left boundary for the triangular wave. In contrast, in Fig. 7b, a significant difference in σ h can be observed at τ/T = 80 in the complex and trapezoidal wave as compared to the triangular and sine waves. This time step corresponds to the change between low and high river stage.
We also computed the probability density functions from the output expansions of the groundwater heads. The results are shown in Fig. 8. We observe small uncertainties  Fig. 8a, c). Larger uncertainties are observed when μ h ≈ 0, when τ/T = 80 (see Fig. 8b). These behaviors occur in all the scenarios. However, two peaks of high probability are observed in Fig. 8b in the trapezoidal scenario due to the rapid fluctuation of the heads in the transient boundary conditions. The trapezoidal scenario also shows the smallest uncertainty at τ/T = 60 and τ/T = 100, because of the low influence of the phase uncertainty in the points where the heads in the boundaries are constant (i.e., minimum and maximum).
The influence of the uncertain transient boundary conditions is also reflected in the Okubo-Weiss values shown in Fig. 9. As expected, the spots with large uncertainty appear in the regions with high hydraulic conductivity contrast and large hydraulic conductivity. Specifically, we allocate large uncertainties in the complex and trapezoidal scenarios at τ/T = 80 (Fig. 9b). This is a consequence of the large uncertainty in the groundwater heads observed previously (Fig. 7b), which occur due to the effect of the phase shift over the sharp movements in the boundary conditions. The steepness of the slopes in these waves creates a wide range of variability when we introduced the offset uncertainty. Moreover, the magnitude of the μ ξ values vary depending on the scenario, finding larger values in the trapezoidal and complex scenarios. Similar to the outcomes of the deterministic scenarios, the results of μ ξ also reveal spots with variable dominances of the deformation and rotational forces of the flow field.
The responses of both h and ξ , as well as the statistics that define their uncertainty, follow periodic patterns. To evaluate the average behavior of the statistics of h and ξ at a certain distance from the boundaries and their differences, we chose a relatively close distance to the stream, at distance λ/2 (i.e., x 1 /L 1 = 5), where the propagation of the signal is clear. We computed the arithmetic means of the uncertainty statistics at a distance x 1 /L 1 = 5, which are shown in Fig. 10. We see that μ h is highly fluctuating in the trapezoidal and the complex scenario. Moreover, according to the interval [μ ξ − σ ξ , μ ξ + σ ξ ], the trapezoidal scenario exhibits the highest uncertainty in ξ , followed by the complex wave. The results also indicate that it is more likely to find rotation properties dominating in the flow field under the trapezoidal scenario conditions than to find them under the conditions of the other scenarios. On the other  hand, we see smaller variability of σ h and σ ξ in the sine wave scenario, showing a similar spread in the outputs along the whole simulation period. While the behavior of the triangular wave is similar to the sine wave, the complex wave is comparable with the trapezoidal wave.
Overall, our results from the deterministic and stochastic scenarios show that waveshaped boundary conditions can influence not only the magnitude of the deformation and rotational forces of the flow field (i.e., shear, stretching, and vorticity) but also the temporal dynamics of dominance between local strain and rotation properties. Although our results show that their location is determined by the areas with high hydraulic conductivity contrast, as can be seen in Fig. 1b, we provide evidence that the mixing potential in these areas is significantly affected by highly transient boundary conditions. This occurs due to the variety of hydraulic gradient responses as a consequence of the highly fluctuating head boundary conditions. To observe in detail the temporal variation and the two-dimensional distribution of the groundwater heads, the Okubo-Weiss values, and the statistical moments that describe the uncertainty, we refer to a series of videos included as part of the Supplementary Material of this research.

Conclusions
We studied the effect of the periodic stage conditions due to hydropeaking events on the groundwater flow topology in terms of the Okubo-Weiss metric. We imposed Dirichlet boundary conditions in the form of wave-shaped specified-heads with four types of waveforms: triangular, sine, complex, and trapezoidal. The formulation of the system considers a deterministic solution of the heterogeneous hydraulic conductivity field for all the scenarios. The first part of our analysis was done assuming no input uncertainty over the four waves that define the transient boundary conditions. The second part of the study approached the problem as a stochastic system with uncertainty in the parameterization of the transient boundary conditions. Here, the wave amplitude and phase were considered uncertain and treated as mutually independent random variables. These variables introduced the uncertainty related to unknown fluctuations in the discharge volume and discharge duration and temporal uncertainty due to energy market demands and powerplant management. The application of polynomial expansions and pseudo-spectral collocation method allowed us to estimate the statistical moments (i.e., mean, and standard deviation) of the outputs of interest (i.e., groundwater heads and Okubo-Weiss metric). The method was convenient to extract the required spatial and temporal detail with low computational effort.
One of the main messages that the Okubo-Weiss metric can provide us is the identification of reaction hotspots. The spatial distribution of the Okubo-Weiss responses is fundamentally controlled by the hydraulic conductivity. In accordance, our results show that their location is determined by the areas with high hydraulic conductivity contrast. However, we also provide evidence that the mixing potential in these areas is significantly affected by the highly transient boundary conditions. The magnitude and temporal behavior of this topological indicator of mixing significantly vary according to the imposed boundary conditions. Different highly transient boundary conditions influenced in different degree the temporal dynamics of dominance between local strain and rotation properties and the magnitude of the deformation and rotational forces of the flow field. Therefore, given the dynamic responses of the flow field to the time-variant head boundary conditions, the detailed temporal characterization of this metric is important to reliably predict, for instance, mixing-driven reactions.
The evaluation of hydropeaking impacts on subsurface flow requires to characterize the management of the surface water system and the intensity of the impact (e.g., shape, amplitude, and periodicity of the wave). Hence, we think it is essential to estimate hydropeaking effects on flow and transport processes in aquifers using a stochastic approach, not only due to the essential uncertainty in the aquifer heterogeneity but also due to the uncertain stream stages. The statistical moments that describe the propagation of the uncertainty show a periodic behavior and a varying degree of uncertainty depending on the applied wave-shaped boundary. Further work should focus on the characterization of real hydropeaking events to explicitly acknowledge the inherit uncertainty of these systems and its effect in the estimation of topological descriptors.