{'fraction_non_alphanumeric': 0.027063697816974777, 'fraction_numerical': 0.008867782812007851, 'mean_word_length': 5.093002657218777, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
When Do Complex Transport Dynamics Arise in Natural Groundwater Systems?

In a recent paper (Trefry et al., 2019) we showed that the interplay of aquifer heterogeneity and poroelasticity can produce complex transport in tidally forced aquifers, with significant implications for solute transport, mixing and reaction. However, what was unknown was how broadly these transport dynamics can arise in natural groundwater systems, and how these dynamics depend upon the aquifer properties, tidal and regional flow characteristics. In this study we answer these questions through parametric studies of these governing properties. We uncover the mechanisms that govern complex transport dynamics and the bifurcations between transport structures with changes in the governing parameters, and we determine the propensity for complex dynamics to occur in natural aquifer systems. These results clearly demonstrate that complex transport structures and dynamics may arise in natural tidally forced aquifers around the world, producing solute transport and mixing behaviour that is very different to that of the conventional Darcy flow picture. Key Points: * Transient Darcy flows can generate complex transport dynamics in heterogeneous compressible aquifers. * This complex transport can trap dispersing solutes for many years. * Global tidal maps indicate widespread potential for complex transport dynamics in coastal zones.


Introduction
Many natural groundwater systems involve discharge boundaries at rivers, lakes and coasts where subsurface aquifer fluids interact with fresh or saline fluids at the discharge boundary, leading to strong chemical and nutrient gradients and highly localized hotspots of chemical and biological activity. Invariably, these discharge boundaries have transient forcing via e.g. tidal action, seasonal variations, barometric effects or long-range climatic cycles, which interact with the regional aquifer flow to generate complicated mixing phenomena near the discharge boundary.
One important class of such systems is coastal aquifers, where tidal forcing is critical to both saltwater intrusion into the aquifer and a wide range of physical and chemical processes. Whilst much attention has been focused on the role of density-driven mixing in the Ghyben-Herzberg zone, interplay of tidal forcing with regional groundwater flow can also generate complex mixing phenomena, leading to predictions of surprising intrusion dynamics for which no field studies have yet been attempted. Similarly, transient forcing of lacustrine and riverine aquifers can generate complex mixing phenomena between water bodies of different bio-geo-chemical origin, leading to highly localized concentration gradients and complex transport phenomena. The transient nature of these complex flow dynamics means they can be difficult to visualize and understand (even in two-dimensional idealizations), and specialized mathematical tools and techniques are required to elucidate the impacts on solute transport and associated phenomena such as chemical reaction or biological activity.
Recently, in Trefry et al. (2019) (called TLMW henceforward) we uncovered mechanisms that can lead to complex mixing in periodically forced Darcy groundwater systems, where the interplay of tidal forcing, aquifer compressibility and heterogeneity caused important changes in the flow structure and transport dynamics. Using a highly simplified, vertically averaged twodimensional (2D) aquifer model, we showed that the interplay of these properties can result in localized flow reversal (where the groundwater flow vector makes a sweep through all orientations) over the periodic forcing cycle. Such changes in flow orientation reorganize the effective Darcy flow topology away from the conventional flow structures of steady Darcy flow.
We term these changes topological bifurcations of the aquifer transport structure and, as shown below (Figure 2), these bifurcations fundamentally change the organization of fluid trajectories in the aquifer and lead to transport dynamics that are radically different to those of steady Darcy flow.
Unlike steady 2D flows, where the transport dynamics can be directly visualized from the streamlines of the flow, transient 2D flows require different methods to understand the transport structure. The simplest way to visualize this transport structure is via a Poincaré section (Aref, 1984;Ottino, 1989), which records fluid particle positions at integer multiples of the forcing periods, effectively filtering out the "fast" intra-period particles trajectories, leaving only the "slow" time-averaged transport. TLMW showed that flow reversal can lead to two major changes in the slow, time-averaged transport structure. First, flow reversal can lead to a change in the effective transport structure of natural aquifers from open 1D particle trajectories without fixed points (i.e. similar to streamlines for steady Darcy flow), to closed 1D particle trajectories with fixed points. Second, even stronger flow reversal can then lead to a bifurcation of the closed flow structure to chaotic particle trajectories, where the closed 1D particle trajectories (in the periodic sense) break up into a 2D scattering of particle locations that lead to rapid mixing and strongly anomalous transport. We collectively term these profound changes in the transport structure of these flows to closed and chaotic particle orbits as complex transport dynamics in the aquifer.
Complex transport dynamics determine the mixing, dispersion and residence time characteristics of these flows, with significant implications for solute transport, chemical reactions and biological activity (Heiss et al., 2017;Liu et al., 2018;Slomp & Van Cappellen, 2004;Tél et al., 2005).
Such fundamental changes in the transport structure of tidal aquifers are governed by a small number of dimensionless parameters which are related to the aquifer properties, tidal forcing, and regional flow characteristics. In TLMW, we showed that complex transport dynamics can arise in coastal aquifers with physical properties similar to those of several naturally occurring systems.
However, as in that study we considered a single set of governing parameters, the transport and mixing behaviour of these systems across the full dimensionless parameter space is presently unknown. This question is relevant to understanding the nature and ubiquity of complex transport dynamics in tidally forced aquifers and implications for physical processes such as mixing, solute dispersion, seawater intrusion, pollutant transport, chemical reactions and biological activity. Therefore, we scan the parameter space and examine the propensity for complex transport dynamics to occur in field settings.
To address this question, we perform a range of simulations using the periodically forced 2D aquifer model introduced in TLMW to resolve the aquifer transport dynamics over the relevant parameter space. Precisely defined below, these parameters include the strength of tidal forcing relative to that of the regional flow, relative compressibility of the aquifer, and the attenuation of the tidal fluctuation into the aquifer, as well as the statistical parameters of the random hydraulic conductivity field. From these results we determine the key transport characteristics over the flow parameter space and develop a phase diagram for the different transport structure topologies (i.e. open, closed or chaotic). We then use these results to develop mechanistic arguments for how and why transport varies over the parameter space, and how these governing parameters for natural systems relate to these regions of augmented transport. Whilst solute transport is governed by the interplay of fluid advection and molecular diffusion or hydrodynamic dispersion, in this study we focus on the purely advective dynamics imparted by transient forcing of the aquifer. This approach, which uncovers the Lagrangian kinematics of these flows, facilitates a clearer visualization of the transport structures that govern solute transport in periodically forced aquifers. Subsequent studies will consider the interplay of advective transport, solute diffusion and dispersion. Finally, we estimate the impacts of diffusion upon solute transport and the propensity for complex transport to occur in natural coastal aquifers.
The remainder of this paper is organized as follows. In Section 2 we review the model for flow in a heterogeneous tidally forced aquifer and define the key dimensionless parameters. We then examine the governing tidal mixing mechanisms in Section 3 and topological bifurcations of the transport structure (Lagrangian kinematics) over the aquifer control parameter space in Section 4. In Section 5 we estimate the impact of complex transport upon solute diffusion and dispersion, and in Section 6 we consider the propensity for complex transport to occur in natural systems before conclusions are made in Section 7.

Governing equations
To begin we briefly review the model for flow in a heterogeneous tidally forced aquifer introduced in TLMW as this will serve as a basis for the calculations herein. We consider in Figure   1, a vertically averaged 2D linear groundwater model that is driven by a constant regional flow gradient J from right to left in an aquifer with a spatially heterogeneous conductivity field K. We note that the inland (right, x = L) and no flow (bottom, y = 0; top, y = L) boundaries do not correspond to physical features of the aquifer, rather these boundaries are the physical extent of the computational model. Conversely, the tidal boundary (left, x = 0) is a discharge boundary to an open flow system such as oceanic, riverine or lacustrine flow, where tidal, seasonal, barometric or climatic variations give rise to a temporally fluctuating flow potential. The interaction between this oscillating potential at the tidal boundary and the regional flow gradient can give rise to complex transport dynamics within the aquifer. In this linear groundwater model, the Darcy flux q and aquifer porosity are governed respectively by the Darcy flow and linearized porosity-storativity equations where h is the pressure head (with reference value ℎ -./ ), the aquifer porosity (with reference value -./ ), S is the constant aquifer storativity and K(x) is the heterogeneous isotropic saturated hydraulic conductivity which varies with the physical space coordinate x = (x, y). Equations (1) and (2) may be combined via the continuity equation into the linear groundwater flow equation in terms of the pressure head h subject to the no flow boundary conditions at y = 0 and y = L, fixed head condition at the inland boundary x = L and periodically fluctuating head at the tidal boundary x = 0, where g( ) = g : <=> is the periodic tidal forcing function. For simplicity we consider tidal forcing g( ) via a single Fourier mode with constant forcing frequency (and period = 2 / ) and amplitude g : .
With such unimodal forcing, equation ( To visualize the transport structures in the periodically forced aquifer, we heavily utilize particle tracking, where the position x of the fluid particle evolves with time t as per the advection equation where ( , ) is fluid velocity defined as the Darcy flux normalised by the local porosity. TLMW detailed the computational method to accurately compute ( , ) to ensure exact satisfaction of the continuity equation (3), in the Supplementary Material (S1) we describe a method to rapidly solve This conventional linear groundwater flow model is the basis for study of the interplay of tidal forcing, non-zero matrix compressibility and aquifer heterogeneity, leading to complex transport dynamics in periodically forced aquifers. As the advection equation (9) represents a transform from the Eulerian to the Lagrangian frame, we shall make extensive use of particle tracking via (9) to study the Lagrangian kinematics and transport dynamics of periodically forced aquifers.

Key dimensionless parameters
TLMW identified a set of dimensionless parameters that characterize the solutions to the 2D tidal forcing problem introduced above. In addition to these dynamical parameters, there also exists a number of statistical parameters which characterize the heterogeneous hydraulic conductivity field K(x). Whilst generation of complex transport dynamics in periodically forced aquifers requires a heterogeneous hydraulic conductivity field to generate flow reversal, we make the conjecture that this phenomenon persists regardless of the particular statistical autocorrelation model used for the conductivity field. In this study we consider a random log-Gaussian conductivity field K that is statistically determined by its mean conductivity .// , log-variance fgh i j and correlation length , and so the independent set of physical parameters that characterize the statistical conductivity field is then ≡ ( .// , fgh i j , ). In this paper we will study advective flow in the linear groundwater model over the tidal flow parameter space ≡ × × , identify regions of anomalous and chaotic transport, and comment on how and together control these topological transitions and the underlying physical mechanisms. Table 1 summarises the various dimensionless dynamical and heterogeneity parameters and provides a brief physical motivation for each.

Transport structures of periodically forced aquifers
The complex transport structures which arise in periodically forced aquifers require specialized tools for their visualization and understanding, which may be unfamiliar to researchers in groundwater hydrology. As shall be shown, we observe three main types of transport structures (open, closed, chaotic) in periodically forced aquifers over the flow parameter space = × × . In this Section we briefly review these visualization tools, classify and describe these transport structures and show how they change with , , .

Visualisation of transport structures and Lagrangian kinematics
Steady Darcy flow is typified by open streamlines and an absence of stagnation points (Bear, 1972), leading to slow mixing and limited transport dynamics (Dentz et al., 2016). Conversely, unsteady Darcy flows can break these topological constraints due to transient switching of streamlines (Lester et al., 2009;Metcalfe et al., 2010;Trefry et al., 2012), leading to a much richer set of possible transport structures. Particle trajectories in transient flows can become very complicated, so direct plotting does not provide much insight. To circumvent this problem, we use the Poincaré section (Aref, 1984;Ottino, 1989) to resolve the transport structure (or Lagrangian kinematics) of the flow by advecting (via equation (9)) a large number of fluid particles and recording the particle locations at each forcing period P. This filters out the "fast" oscillatory intraperiod particle trajectories, leaving only the "slow" net particle motion. The Poincaré section can contain regular, non-chaotic regions which are comprised of coherent 1D particle trajectories that are non-mixing (in the ergodic sense), and distinct chaotic regions which are comprised of a seemingly random particle distributions which are formed as particles undertake chaotic, space-filling orbits. Chaotic regions are associated with localized strong (exponential) fluid stretching and rapid mixing, whereas the coherent trajectories in regular regions form barriers to transport.
We shall also make extensive use of Residence Time Distributions (RTDs) to characterize transport in the periodically forced aquifers. Although particle trapping is not possible in steady

Open transport structures
As shown in Figure

Closed transport structures
Above a critical value (termed 1) of the relative compressibility (Figure 2b), the aquifer Poincaré section appears to undergo a topological bifurcation from open trajectories to particle trajectories that are now closed. This represents a fundamental change in the transport structure as fluid elements can neither enter nor leave this closed region, and so trapped fluid elements have essentially infinite residence times. In addition to these closed regions, periodic points of the flow also arise (where fluid particles return to the same position after one (k = 1) or more (k > 1) flow periods) which then manifest as fixed (k = 1) or periodic (k > 1) points in the Poincaré section. As shown in Figure 2, these periodic points may be classified as elliptic (E)  Fluid elements near hyperbolic points undergo exponential fluid deformation over each flow period, with stretching (blue) and contracting (red) particle trajectories shown in Figure 2c that are respectively termed the unstable ( ‚ƒ " ) and stable manifolds ( ‚ƒ x ) of the flow. These manifolds play a critical role in organizing transport. If these manifolds connect smoothly (Figure 2b, c), then the associated exponential stretching cancels out and regular particle trajectories result, but more complex transport dynamics arise if they intersect transversely.

Chaotic transport structures
As shown in Figure 2d, further increases in beyond a second critical value 2 (where 2 > 1) leads to breakup of regular, coherent 1D particle trajectories into chaotic, space-filling particle trajectories that manifest as chaotic regions in the Poincaré section. From classical understanding of Hamiltonian chaos (Aref, 1984;Ottino, 1989), these chaotic dynamics arise as the stable and unstable manifolds no longer connect smoothly, resulting in a chaotic saddle within the open flow.
This region is comprised of a fractal tangle of stable and unstable manifolds that imparts exponential stretching to fluid elements and a fractal distribution of residence times, see Tél et al. (2005), Toroczkai et al. (1998) for details. In contrast, elliptic points (E) represent regular, non-mixing regions in the flow. The Kolmogorov-Arnol'd-Moser (KAM) theorem states that the outermost orbits of these regular "islands" (termed KAM islands, see Figure 2d) around elliptic points are the least stable and so will break up into chaotic trajectories with further increases in the compressibility parameter beyond 2. Eventually all of the orbits will become chaotic but the elliptic point itself will persist.

Tidal Emptying Region and Summary
There also exists an inflow/outflow region near the tidal boundary termed the tidal emptying

Parameter Sweep
To determine the range of transport dynamics and Lagrangian kinematics in periodically forced aquifers we compute the linear groundwater model over the aquifer parameter space = × × . We analyse the transport structures outlined in the previous section in order to explore the distribution of Lagrangian kinematics over the parameter space and understand of how the underlying physical mechanisms represented by the dimensionless parameters control these transport structures. Unless stated otherwise, all simulations in this study use a 164 × 164 finite difference grid to resolve flow over the 2D hydraulic conductivity field with mean value .// = 2´10 -4 , log-variance ( fgh … j = 2) and correlation length ( = 0.049). We use the finite-difference and particle mapping methods described in TLMW respectively to provide numerical solutions to the flow field ( , ) and then generate the associated Poincaré sections over the flow parameter space = × × = (1, 10, 100) × (1, 10, 100) × (0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5) in efficient manner. In order to ensure robustness of the particle mapping method to solve (9), we tested two different integration schemes (explicit and implicit Runge-Kutta) and mapping grid resolutions We also test for statistical invariance in all cases by performing simulations for three different realisations of the random hydraulic conductivity field. Of primary interest are the transport dynamics that arise within the periodically forced aquifer, the topological bifurcations between transport structures and the physical implications of the transport dynamics over the aquifer parameter space defined by and . This parametric sweep study elucidates potential ramifications for coastal groundwater systems, e.g. transport, mixing and chemical reactions and uncovers the physical mechanisms that control these processes.

Impact of tidal strength and attenuation on aquifer transport
To begin, we consider the impact of both tidal strength and attenuation on the transport dynamics in the tidally forced aquifer by varying the Townley number and tidal strength at a fixed value ( = 0.1) of the tidal compression ratio. The Poincaré sections generated by the distributed injection protocol are shown in Figure 3, where the purple curves indicate the tidal emptying boundary, and the region between this boundary and tidal boundary (x = 0) is the tidal emptying region which indicates the part of the aquifer for which intense mixing between freshwater and saltwater occurs over a single flow period. We also note that the rightmost boundary (brown) in Figure 3 is the inland boundary, given by the location of particles that have been seeded  and . If the value of ℋ > is large, particles will drift many tidal periods, but their elliptical velocity orbits become too small to generate significant flow reversal. The spatial character ℋ` provides information on how many vortices are likely to occur within the aquifer, so complex transport arises when ℋ` is greater than one. As shown in Figure 3, both ℋ` and the width of the tidally active zone decrease with increasing , and these values increase with increasing tidal strength .
Hence, complex transport structures arise for moderate values of and large values of . These trends persist for different values of the tidal compression ratio beyond that ( = 0.1) used in Figure 3.
The closed regions shown in Figure 3 occur in both the aquifer bulk and in the tidal emptying region, with significant implications for transport of solutes. If trapped regions in the tidal emptying region are initially filled with freshwater, then these trapped regions represent "islands" of freshwater within a "sea" of brackish water in the tidal emptying region, and salt may only enter via diffusion or dispersion. In Section 5.2 we show that these regions persist for long times and strong concentration gradients are expected at their boundaries. Conversely, trapped regions within the aquifer bulk may represent islands of different water composition to that of the fluid entering via the inland boundary (x = L), and the only transport mechanism into or out of these regions is diffusion or dispersion. Hence these augmented dynamics and concentration gradients will have strong impacts upon transport, mixing and chemical reactions within the aquifer.
As the parameter intervals for and in Figure 3

Impact of aquifer compressibility on transport
To understand the role of aquifer compressibility on transport, we now vary the compressibility parameter . Figure 5 shows a series of aquifer Poincaré sections using the distributed injection protocol at various values of in the range = 0 -0.5 for = 10 and = 100. With increasing , the particle drift velocity increases and the tidal emptying region widens.
In TLMW, we predicted that fluid tracer particle trajectories are integrable for incompressible ( = 0) aquifers, which means they propagate in an oscillatory manner along smooth 1D streamlines and so cannot exhibit trapping or chaotic dynamics. Figure 5 confirms that particle trajectories are indeed confined to 1D orbits for = 0, and this behaviour appears to extend to very weakly compressible aquifers ( = 0 -0.002). Closed particle orbits are observed for weakly compressible aquifers ( = 0.005 -0.1), whereas strongly compressible aquifers ( = 0.2 -0.5) undergo a further topological bifurcation to exhibit chaotic transport structures, as indicated by the chaotic saddles highlighted in Figure 5.
As the region of parameter space ( , , ) for which chaotic transport is observed is a subset of that for which closed time-averaged orbits are observed, we hypothesise that as either , or are increased, the basic transport structure may bifurcate from open to closed orbits and admit periodic points that may be classed as elliptic (E) or hyperbolic (H). With further increases in , or , the stable and unstable manifolds associated with the hyperbolic periodic points intersect transversely, corresponding to a bifurcation to chaotic transport dynamics, with important transport implications beyond that of simple trapping in closed orbits (see TLMW for details). Thus, there exists a hierarchy of transport complexity in periodically forced aquifers based on the route to chaos in aquifer flow structures as described above.

Sensitivity to hydraulic conductivity field
All numerical treatments of Darcy flow in porous media depend on the model chosen to represent conductivity heterogeneity. Of primary interest, however, is the essential physics of the studied flow, i.e. those features which are robust to both the class and details of the conductivity model. In TLMW, we used a simple Gaussian autocorrelation function for log-conductivity that is statistically determined by the mean conductivity ( .// ), conductivity log-variance ( fgh i j ) and correlation length ( ) to analyse the transport structure of the associated flow field. The concept was that if complex transport dynamics can occur in these simple conductivity fields, then it is possible that similar dynamics will occur in more sophisticated geo-statistical models.

Solute trapping and transport
In this Section we consider the impact of closed transport regions upon solute transport by considering both the size distribution of closed regions across the parameter space, and the impact of the size of trapped regions upon solute diffusion and dispersion.

Size distribution of closed regions
When closed regions do occur within the aquifer, both the location (i.e. whether in the tidal emptying region or the aquifer bulk) and size distribution of these regions have important implications for solute transport, as solutes will more rapidly disperse out of smaller closed regions. We use tracer particle RTDs under the distributed injection protocol to determine the total of the area (a) of trapped regions in the aquifer, the resultant RTD histogram (probability density function) for = 100, = 100 and = 0.1 is shown in Figure 7a. Particle tracking is performed up to τ/P = 10 5 , most tracer particles leave the domain within 4×10 4 periods, but a significant proportion (4.28%) are confined indefinitely which are deemed to be trapped regions. = 10, = 100, the total trapped area is smaller for =100, = 100 due to the smaller tidal region shown in Figure 3. Figure 7b indicates that trapped regions are typically smaller for small values of and , but with parameter increases the trapped regions increase in size and coalesce into fewer but larger regions. As we show explicitly in the next section, larger trapped regions trap diffusive solutes for much longer time and so represent significant regions of retarded transport.

Solute transport across trapped regions
Whilst the analysis in this study focuses predominantly upon the advection processes that govern these transport structures, solutes are typically transported by advection and diffusion/ dispersion. In the presence of dispersive transport mechanisms, closed regions of the flow become regions of retarded transport. To provide and order-of-magnitude estimate of the transport rate of solutes into or out of these regions, we approximate the trapped region as a disc of radius G (such that G is the disc radius relative to the aquifer dimension L) that is surrounded by fluid with zero solute concentration in the open region of the flow. This estimate considers the dispersive "leakage" of the solute concentration field ( , t) out of such a trapped region subject to homogeneous Dirichlet conditions (c = 0) at the boundary. Under these assumptions, the solute concentration field is governed by the dispersion equation where Ž is the hydrodynamic dispersion tensor.

Diffusive transport
We first consider the case of molecular diffusion, in which case Ž = • , where • is the molecular diffusivity. In this case the total amount of solute ( ) ∈ [0,1] in the trapped region relative to the initial amount (such that (0) = 1) is where ˜ is the n-th zero of the zero-th order Bessel function of the first kind and Pe • = j /( • ) is a Péclet number defined in terms of the forcing period P and dimension L of the aquifer. The molecular diffusivity of water • ~ 10 -9 m 2 /s then yields an estimate of this Péclet number of order Pe • ~ 10 11 for a typical coastal aquifer of width L ~ 2000m subject to a P = 12 hour tidal period. As per equation (12), for all but very short times, the total relative solute amount decays exponentially at a rate very close to − ‚ j /( G j Pe • ) (with ‚ =2.40483), thus 90% of solute has leached out of the trapped regions after time §G% , which is defined as §G% ≈ Pe • G j ‚ j ln(10).
Thus, the residence time of a solute within a given trapped region scales as G j , which is proportional to the trapped region area (a). Hence the histograms of the trapped area size (a) shown in Figure 7b can be used to directly inform how solute transport varies with the aquifer parameters , , . In the supplementary material (S3 and S4) we also show that the relative radius G of trapped regions scales linearly with the correlation length scale at a rate that varies with the Townley number , hence the duration of solute trapping scales quadratically with the correlation length scale. Figure 7b shows that a typical closed region has relative area of a ~ 0.001, and so the corresponding disc has relative radius G ~ 0.0178 and from (13) this takes around §G% ~ 16,000 years to empty. Hence these trapped regions have a profound impact upon transport of diffusive solutes.

Dispersive transport
We quantify solute transport under hydrodynamic dispersion via the Scheidegger parameterization for the longitudinal H component of the hydrodynamic dispersivity tensor D (under the assumption of relatively small transverse dispersion) which is related to the local velocity as where the longitudinal dispersivity for a coastal aquifer H is estimated to be of order 10 m or larger. In the supplementary material (S5) we show that the regional flow « = .// / .// provides a reasonable estimate of the local velocity in (14)

Physical relevance
The parametric studies in the previous section clearly show how complex transport is a natural consequence of matrix heterogeneity and compressibility in formations subject to periodic forcing. Nevertheless, it may not yet be clear when and where to expect complex transport in natural aquifers. TLMW showed that physically plausible values for and yielded Poincaré sections exhibiting elliptic and hyperbolic points at high tidal compressibility = 0.5, and results from the present work have extended the Lagrangian complexity regime to much lower tidal compressibility values. In this section we seek to (i) clarify how realistic coastal aquifer systems relate to the topological bifurcation phenomena, (ii) provide the reader with a simple rule of thumb for assessing the potential for Lagrangian complexity in any particular field setting. In this analysis, we maintain the single-mode spectral approximation for the tidal forcing, although it is recognized that tidal spectra commonly display multimodal forms, and the effects of multimodal forcing is the subject of future studies.

Tidal transition diagram
As discussed in TLMW, published coastal aquifer studies can be used to assess likely parameter ranges for common coastal aquifer types including sediments, limestones, fractured basalts, etc. However, the Lagrangian flow complexity is also coupled with the effective coastal tidal amplitude which is a complicated function of position around the globe. We can decouple this complexity by plotting a small number of published field sites in × × parameter space, superimposing the results of the topological analysis, and seeking zones of overlap for increasing tidal amplitude. We term such a plot the tidal transition diagram as it maps how tidal flows in specific aquifer formations move toward Lagrangian complexity as the effective tidal amplitude increases. The tidal transition diagram in Figure 8 is constructed as follows: i. A rectangular prism is constructed in × j × in log-scale space.
ii. Three-dimensional coloured ellipsoids, with ellipse axes approximating error bars on the parameter estimates, are plotted according to the parameters reported in the published field studies discussed in TLMW.
iii. Shaded zones are superimposed to indicate the parametric extents of the elliptic bifurcation zone (pink) and the zone of chaotic flow (red) found in the present simulations. The phase lines can be computed as exact linear functions of g : from the definitions for and in Section 2. Figure 8 shows that the tidal trajectories are increasing functions of g : and display slightly curvilinear profiles since the axis is not drawn to be exactly logarithmic. Notably, all aquifer types impinge on the pink elliptic bifurcation zone as g : increases, at amplitudes as low as 2m, see notation marks in Figure

Conclusions
TLMW first demonstrated that the interplay of tidal forcing, aquifer heterogeneity and poroelasticity can produce complex Lagrangian transport, which has profound ramifications for solute transport and reaction in transiently forced groundwater systems. What was unknown, however is to what extent these transport dynamics arise in natural groundwater systems and how these dynamics depend upon the aquifer properties and characteristics of the regional and tidal flows. In this study we answered these questions through parametric studies of the governing parameters that characterize the aquifer and tidal properties. This study elucidates the mechanisms that govern complex transport dynamics, bifurcations between structures with changes in the governing parameters, and indicates the propensity of complex dynamics to occur in natural aquifer systems.
To achieve this, we defined two distinct parameter sets; the aquifer parameter set = ( , , ) which controls the dynamics of a periodically forced aquifer, and the hydraulic conductivity parameter set = ( .// , fgh i j , ) which characterises the log-Gaussian conductivity field. We used numerical simulation to consider the transport dynamics of a model tidally forced aquifer over these two parameter sets. We observed a change from open fluid particle orbits (similar to open streamlines characteristic of steady Darcy flow) to closed orbits that can trap fluid particles for indefinitely long periods. With further parametric variation, these closed orbits can then become chaotic, leading to accelerated mixing dynamics and slow "leaking" of particles into and out of these previously closed regions. We find that these three different transport structure types (open, closed, chaotic) arise at particular combinations of the governing parameters, leading to a "phase portrait" of transport dynamics over the aquifer parameter space that we explain via the physical mechanisms that govern flow in periodically forced aquifers.
Stability and robustness of the Gaussian statistical model have also been tested via changing the correlation length ( ), log-variance ( fgh i j ) and log-conductivity realization, and we find that size of closed regions increases linearly with the correlation length ( ), and the propensity for complex transport increases with log-variance ( fgh i j ). Although we chose a simple logconductivity field in this study, these results give confidence that the dynamical parameters will likely generate similar flow characteristics regardless of the aquifer heterogeneity models employed.
We have quantified the impact of closed and chaotic regions upon transport of diffusive and dispersive solutes and find that such trapping impacts solute transport over long timescales at a rate proportional to the size of the closed flow region. We also examine parametric variation of the size distribution of trapped regions with the aquifer, allowing the impacts on solute transport to be directly estimated. Simple calculations show that common coastal aquifer types can display Lagrangian bifurcations to trapping flows for sufficiently high tidal amplitudes; such tidal amplitudes may be prevalent around the world, e.g. in northern Europe, east Africa and the Alaskan-Canadian Pacific coast. These phenomena may profoundly impact on biogeochemical processes in coastal aquifers.
The results in this paper predict that complex transport structures and dynamics arise in natural tidally forced aquifers around the world. When such dynamics do arise, it will lead to mixing and reactive transport patterns that are very different to those predicted by conventional Darcy flow models.

Introduction
Text S1 describes the mapping method for rapid solution of the advection equation.
Text S2 and Figure S1 illustrates the impact of arbitrary realizations for the log-conductivity field.
Text S3 and Figure S2, 3 illustrates the effect of different correlation lengths upon the aquifer transport structures.
Text S4 and Figure S4 shows the distribution of 0 2 and 0 2 /v for various values of the Townley number and tidal strength.
Text S5 provides an estimate of the local velocity of fluid elements in periodically forced aquifers.

Text S1. Mapping Method
We employ a "mapping method" to rapidly advect large numbers of fluid particles based upon solution of the advection equation (13)  In TLMW we show that the fluid velocity is net divergence-free over a single forcing period of the aquifer, hence in principle, this method should yield an area-preserving mapping such that the Jacobian = det(∇ ) , = ( , ) is everywhere unity. However, interpolation of particle trajectories via equation (S1) above introduces small (typically 10 -3 ) deviations in the Jacobian from unity. Whilst these deviations can produce spurious behaviour in e.g. chaotic regions of the flow over very large particle residence times, these effects are relatively minor and do not qualitatively impact the flow topologies in these regions. As expected, large deviations in the Jacobian also occur near the inflow/outflow domain boundaries due to the trapping of particles here, however these errors are highly localised and do not impact the internal transport structures. Once the mapping functions ( , ) are generated, these may be used to rapidly advect large numbers of particles of many flow periods, facilitating rapid generation of the Poincare sections used in this study.
Text S2. Impact of conductivity realization Figure S1 shows three arbitrary realizations of the log-conductivity field (for fixed eff ,     In TLMW, we constructed numerical methods to generate self-consistent numerical solutions to physically constrained flows and accurate advection methods to construct Poincaré sections. It was shown that the accuracy of the method for generating physically constrained flows declines as log-conductivity variance log 2 increases and as the spatial resolution /∆ is reduced (where ∆ is the spatial resolution of the finite difference grid). Thus, while the sections in Figure S2 show a trend of increasing structure and complexity as /∆ decreases, the accuracy of the underlying flow fields increases as /∆ increases. Even so, it is clear from the figure that the calculated time-averaged particle trajectories are smoothly varying for increasing /∆.

Text S4. Distributions of and /
In Section 5.2 of the manuscript, we showed that the residence time of diffusive and dispersive solutes within trapped regions of the flow respectively scales roughly as the squared relative radius ( 0 2 ) of the trapped region, and this area normalized by the local fluid velocity ( 0 2 /v). Figure S4

Text S5. Mean velocity variation
In Section 5.2 we relate the residence time of dispersive solutes to the size of trapped regions and the local fluid velocity. To provide an estimate of this residence time, we need to relate the mean regional flow velocity to the tidally influenced velocity near the tidal boundary. The mean regional groundwater velocity is and the maximum tidal groundwater velocity ̃ can be shown to be ̃= ̅ (1 + (1 − ) taz ) .
We consider some field-scale parameters pertinent to a Safety Bay Sand/Tamala Limestone coastal aquifer near Fremantle (Smith et al., 2012). Representative parameter values are: aquifer width L ≈ 2000 m, dimensionless head gradient J ≈ 0.005, hydraulic conductivity eff ≈ 10 m/d, dimensionless porosity eff ≈ 0.3, dimensionless tidal strength in Fremantle ≈ 0.5 / (JL), tidally affected zone (dimensionless fraction of L) taz ≈ 0.2, dimensionless tidal compressibility ≈ 0.03. Numerical evaluations show that the regional flow velocity ̅ is well within an order of magnitude estimator of the maximum tidally induced velocity ̃ for practical purposes. But for extreme values of and , the tidally induced velocity becomes significantly greater than the regional velocity.