Lagrangian Complexity Persists with Multimodal Flow Forcing in Compressible Porous Systems

Abstract

We extend previous analyses of the origins of complex transport dynamics in compressible porous media to the case where the input transient signal at a boundary is generated by a multimodal spectrum. By adding harmonic and anharmonic modal frequencies as perturbations to a fundamental mode, we examine how such multimodal signals affect the Lagrangian complexity of flow in compressible porous media. While the results apply to all poroelastic media (industrial, biological and geophysical), for concreteness we couch the discussion in terms of unpumped coastal groundwater systems having a discharge boundary forced by tides. Particular local regions of the conductivity field generate saddles that hold up and braid (mix) trajectories, resulting in unexpected behaviours of groundwater residence time distributions and topological mixing manifolds near the tidal boundary. While increasing spectral complexity can reduce the occurrence of periodic points, especially for anharmonic spectra with long characteristic periods, other signatures of Lagrangian complexity persist. The action of natural multimodal tidal signals on confined groundwater flow in heterogeneous aquifers can induce exotic flow topologies and mixing effects that are profoundly different to conventional concepts of groundwater discharge processes. Taken together, our results imply that increasing spectral complexity results in more complex Lagrangian structure in flows through compressible porous media.

Appendices

Appendix 1: Characteristic Periods of Commensurable M-Spectra

We wish to determine the characteristic period P of a periodic boundary forcing condition F(t) with discrete Fourier spectrum ${\mathcal {F}}$ given by

$$\begin{aligned} {\mathcal {F}}(\omega ) = \textit{g}_{0} \; \mathbf {\delta }(\omega ) + \sum _{m=1}^{M} \textit{g}_{m} e^{i\theta _{m}} \mathbf {\delta }(\omega - \omega _{m}) \end{aligned}$$

(17)

with amplitudes $\textit{g}_{0} \ge 0$, $\textit{g}_{m} > 0$, the phases satisfying $0 \le \theta _{m} \le 2\pi$, the commensurable modes $\omega _{m} > 0$ and M finite. $\mathbf {\delta }$ is the Dirac delta function. Several cases are important.

Phase-split spectra: Consider a 2-spectrum, i.e. a spectrum with two identical frequency modes $\omega _{1} = \omega _{2}$, with unequal phases $\theta _{1} \ne \theta _{2}$. Via phasor algebra, it is easy to show that the sum of these two modes is equivalent to a single mode with well-defined amplitude $\big ( \textit{g}_{1}^{2} + \textit{g}_{2}^{2} + 2 \textit{g}_{1} \textit{g}_{2}\, \text {cos}[\theta _{1}-\theta _{2}]\big )^{1/2}$ and phase $\text {arctan}[(\textit{g}_{1}\, \text {sin}\, \theta _{1} + \textit{g}_{2}\, \text {sin}\, \theta _{2})/(\textit{g}_{1}\, \text {cos}\, \theta _{1} + \textit{g}_{2}\, \text {cos}\, \theta _{2})]$. Thus, by extension, any M-spectrum of modes with identical frequencies can be expressed as a single mode with deterministic amplitude and phase. Hence, phase-split spectra are equivalent to the single-mode forcing scenario modelled earlier (Trefry et al. 2019; Wu et al. 2019) and, as a consequence, in heterogeneous flows subject to boundary conditions represented by phase-split spectra no new topological structures can arise than those previously discussed.

Another consequence of the phasor analysis is that any N-spectrum (phase-split or otherwise) can be reduced (deterministically) to an ordered sum of M distinct modes $(\textit{g}_{m},\theta _{m},\omega _{m})$ for $m = 1,2,\ldots ,M$ where the strict inequalities $\omega _{1}< \omega _{2}< \cdots < \omega _{M}$ hold and $M \le N$. Thus, we can restrict attention to M-spectra where the frequency modes are distinct without loss of generality.

$\underline{\hbox {Indeterminate}\,P:}$ Consider the trivial case where modes are of equal frequency and amplitude, and exactly out of phase, for example a 2-spectrum with $\omega _{1} = \omega _{2}$, $\textit{g}_{1} = \textit{g}_{2}$ and $|\theta _{1} - \theta _{2} | = \pi$. Then, the sinusoidal terms cancel exactly and the spectrum ${\mathcal {F}} = \textit{g}_{0} \; \mathbf {\delta }(0)$. Hence, F(t) is a constant and any real number is formally a period of F. This trivial case is easily detected and avoided in practice.

Harmonic spectra: Consider the commensurable case where the modes are ordered so that $\omega _{1}< \omega _{2}< \cdots < \omega _{M}$ and the fundamental mode $\omega _{1}$ is an exact divisor of all modes. In this case, all spectral modes are harmonically related to $\omega _{1}$ and the characteristic (minimum) period of F(t) is $P = 2\pi /\omega _{1}$.

Anharmonic spectra: Consider the commensurable case where the modes are ordered so that $\omega _{1}< \omega _{2}< \cdots < \omega _{M}$ and $\omega _{1}$ is not an exact divisor of all modes. In this case, the spectral modes are anharmonically related, but a characteristic period can still be determined as follows (see Raykh (2017)). Form the list of modal periods $P_{modal} = \{2\pi /\omega _{1},2\pi /\omega _{2},\ldots ,2\pi /\omega _{M}\}$ and express each list member in rational form p/q for integer values p, q. Multiply the $q_{m}$ to form the product $\varDelta$ and determine the greatest common denominator GCD of the entries in the list $M = \varDelta \times P_{modal}$. Calculate the least common multiple LCM of the list $M/\text {GCD}$. The characteristic period P of F(t) is then given by $(\text {GCD}/\varDelta ) \times \text {LCM}$. As an example, a 3-spectrum with modes $\omega _{1}=2\pi , \omega _{2}=2.1\pi , \omega _{3}=4\pi$ yields $P = 20$.

Appendix 2: Comments on Spectral Scaling

We consider here the problem of scaling two different tidal spectra ${\mathbb {S}}_{M}$ and ${\mathbb {S}}_{N}$ so that the effective problem parameters ${\mathcal {C}}_{\text {eff}},{\mathcal {G}}_{\text {eff}}$ (as defined in Table 1) are comparable between the two spectra. Noting that both ${\mathcal {C}}_{\text {eff}}$ and ${\mathcal {G}}_{\text {eff}}$ are both defined in terms of modal amplitudes $\textit{g}_{m}$, it is natural to pursue a scaling of modal amplitudes. For example, if M-spectrum ${\mathbb {S}}_{M}$ contains modes $(\alpha _{1},\theta _{1},\omega _{1}),(\alpha _{2},\theta _{2},\omega _{2}), \ldots ,(\alpha _{M},\theta _{M},\omega _{M})$, then the sum of amplitudes for ${\mathbb {S}}_{M}$ is $a_{{\mathbb {S}}_{M}} = \alpha _{1} + \alpha _{2} + \cdots + \alpha _{M}$. Similarly, if ${\mathbb {S}}_{N}$ contains modes $(\beta _{1},\phi _{1},\upsilon _{1}),(\beta _{2},\phi _{2},\upsilon _{2}),\ldots , (\beta _{N},\phi _{N},\upsilon _{N})$, then the corresponding sum of amplitudes is $a_{{\mathbb {S}}_{N}} = \beta _{1} + \beta _{2} + \cdots + \beta _{N}$. In this picture, the rescaled spectrum ${\mathbb {S}}_{N}^{*}$ with modes $(\beta _{1} a_{{\mathbb {S}}_{M}}/a_{{\mathbb {S}}_{N}},\phi _{1},\upsilon _{1}),(\beta _{2} a_{{\mathbb {S}}_{M}}/a_{{\mathbb {S}}_{N}},\phi _{2},\upsilon _{2}), \ldots , (\beta _{N} a_{{\mathbb {S}}_{M}}/a_{{\mathbb {S}}_{N}},\phi _{N},\upsilon _{N})$ would yield peak amplitudes (and hence ${\mathcal {C}}_{\text {eff}},{\mathcal {G}}_{\text {eff}}$ values) equivalent to those generated by ${\mathbb {S}}_{M}$.

While the above amplitude scaling method is both appealing and simple, the kinematics generated from ${\mathbb {S}}_{M}$ and ${\mathbb {S}}_{N}^{*}$ contain a systematic bias, as illustrated in the first column of Fig. 3. In the figure, it is seen that as the amplitude of the first harmonic increases (down the A simulations column) the size of the velocity locus (solid curve) declines steadily, even though amplitude scaling is employed to keep ${\mathcal {C}}_{\text {eff}}$ and ${\mathcal {G}}_{\text {eff}}$ constant. The reason for this is that the waveform is increasingly being distorted by the first harmonic, reducing the number of peak amplitudes per period (Fig. 2). It is possible to use a spectral power scaling technique instead, where the metric of interest is the power spectral density $p_{{\mathbb {S}}_{M}}$ defined by

$$\begin{aligned} p_{{\mathbb {S}}_{M}} = \frac{1}{P} \int _{0}^{P} |\alpha _{1} e^{i (\omega _{1} t + \theta _{1})} + \alpha _{2} e^{i (\omega _{2} t + \theta _{2})} + \cdots + \alpha _{M} e^{i (\omega _{M} t + \theta _{M})}|^{2} dt \end{aligned}$$

(18)

where P is the characteristic period of ${\mathbb {S}}_{M}$. For harmonic spectra, where all frequency modes are integer multiples of the fundamental $\omega _{1}$, the power spectral density reduces to $p_{{\mathbb {S}}_{M}} = (\alpha _{1}^{2} + \alpha _{2}^{2} + \cdots + \alpha _{M}^{2})/2$. The power-scaled spectrum ${\mathbb {S}}_{N}^{*}$ may be formed from modes $(\beta _{1} p_{{\mathbb {S}}_{M}}/p_{{\mathbb {S}}_{N}},\phi _{1},\upsilon _{1}),(\beta _{2} p_{{\mathbb {S}}_{M}}/p_{{\mathbb {S}}_{N}},\phi _{2},\upsilon _{2}),\ldots , (\beta _{N} p_{{\mathbb {S}}_{M}}/p_{{\mathbb {S}}_{N}},\phi _{N},\upsilon _{N})$; this spectrum displays the same power spectral density as ${\mathbb {S}}_{M}$ but a different waveform. In particular, peak amplitudes will differ between the two spectra ${\mathbb {S}}_{M}$ and the power-scaled ${\mathbb {S}}_{N}^{*}$, leading to variations in ${\mathcal {C}}_{\text {eff}},{\mathcal {G}}_{\text {eff}}$ values. However, as seen from the dashed grey loci in the first column of Fig. 3, the velocity loci maintain approximately constant size as the amplitude of the first harmonic increases.

It is not possible to deduce a scaling procedure that ensures that two different spectra generate identical tidal problem parameters ${\mathcal {T}}_{\text {eff}},{\mathcal {C}}_{\text {eff}},{\mathcal {G}}_{\text {eff}}$ and hence are directly comparable in terms of their Lagrangian structures. This difficulty makes comparison of simulation cases in Table 2 and Fig. 3 an imprecise process. Nevertheless, as we have discussed here, amplitude and power scaling both lead to topologically similar velocity loci at points internal to the problem domain, albeit with different magnitudes. The similarity of the velocity loci between the two scaling approaches provides confidence that our focus on comparing topological characteristics between simulation cases is reasonable. Apart from the limited results for power-scaled harmonic spectra presented in Fig. 3, amplitude scaling is employed by default in all simulations and results throughout this paper.

Appendix 3: N-step Poincaré Mapping Method

In Wu et al. (2019), we described an efficient mapping method for advecting large numbers of fluid particles over one period P of the flow. This method involved advecting a uniform square grid of tracer particles distributed over the entire aquifer domain for the duration of a single forcing period. The resultant particle locations were then used to construct the spline interpolation functions $(f_x,f_y)$ for the x- and y-coordinates at the end of a forcing period $(x_{n+1},y_{n+1})=(x((n+1) P),y((n+1) P))$ as a function of the coordinates $(x_n,y_n)=(x(n P),y(n P))$ at the start of the forcing period:

$$\begin{aligned} x_{n+1}=f_x(x_{n},y_{n}), \quad y_{n+1}=f_{y}(x_{n},y_{n}). \end{aligned}$$

(19)

We use a similar method in this study to advect particles subject to multimodal forcing, with the exception that for cases with long characteristics periods we employ an N-step mapping method such that mapping is performed over a fraction P/N (for some integer $N>1$) of the characteristic period N to improve temporal resolution of RTDs, e.g. for anharmonic simulation case C5. In this case, the (x, y) coordinates at time $t=n P/N$ are denoted as $(x_{n},y_{n})=(x(n P/N),y(n P/N))$, and so the mapping to coordinates $(x_{n+1},y_{n+1})=(x(n P/N),y(n P/N))$ at time $t=(n+1) P/N$ is

$$\begin{aligned} x_{n+1}=f_{x,n}(x_{n},y_{n}), \quad y_{n+1}=f_{y,n}(x_{n},y_{n}), \end{aligned}$$

(20)

where the interpolated functions $(f_{x,n},f_{y,n})$ are constructed by integrating particles from the uniform square grid from time $t=n P/N$ to $t=(n+1) P/N$. Construction of N such maps (where $n=1,2 \ldots N$) then allows rapid particle mapping over arbitrary time periods. The accuracy of the N-step map process was compared to the fully periodic case ($N=1$) and was found to be accurate to within the error of the fully periodic map.