Daniel R Lester
email: daniel.lester@rmit.edu.au

Marco Dentz

Tanguy Le Borgne

Felipe P J De Barros

Fluid Deformation in Random Steady Three Dimensional Flow

Keywords: Fluid deformation, steady flow, porous media, stochastic modelling

The deformation of elementary fluid volumes by velocity gradients is a key process for scalar mixing, chemical reactions and biological processes in flows. Whilst fluid deformation in unsteady, turbulent flow has gained much attention over the past half century, deformation in steady random flows with complex structure -such as flow through heterogeneous porous mediahas received significantly less attention. In contrast to turbulent flow, the steady nature of these flows constrains fluid deformation to be anisotropic with respect to the fluid velocity, with significant implications for e.g. longitudinal and transverse mixing and dispersion. In this study we derive an ab initio coupled continuous time random walk (CTRW) model of fluid deformation in random steady three-dimensional flow that is based upon a streamline coordinate transform which renders the velocity gradient and fluid deformation tensors upper-triangular. We apply this coupled CTRW model to several model flows and find these exhibit a remarkably simple deformation structure in the streamline coordinate frame, facilitating solution of the stochastic deformation tensor components. These results show that the evolution of longitudinal and transverse fluid deformation for chaotic flows is governed by both the Lyapunov exponent and power-law exponent of the velocity PDF at small velocities, whereas algebraic deformation in non-chaotic flows arises from the intermittency of shear events following similar dynamics as that for steady two-dimensional flow.

Introduction

Whilst the majority of complex fluid flows are inherently unsteady due to the ubiquity of fluid turbulence, there also exist an important class of steady flows which possess complex flow structure. Such viscous-dominated flows typically occur in complex domains through natural or engineered materials, including porous and fractured rocks and soils, granular matter, biological tissue and sintered media. The heterogeneous nature of these materials can impart complex flow dynamics over a range of scales [START_REF] Cushman | The Physics of Fluids in Hierarchical Porous Media: Angstroms to Miles[END_REF]). These heterogeneities (whether at e.g. the pore-scale or Darcy-scale) arise as spatial fluctuations in material properties (e.g. permeability), and so the velocity field structure of the attendant steady flows is directly informed by these fluctuations. These flows play host to a wide range of physical, chemical and biological processes, where transport, mixing and dispersion govern such fluid-borne phenomena. An outstanding challenge is the development of models of mixing, dilution and dispersion which are couched in terms of the statistical measures of the material properties such as hydraulic conductivity variance and correlation structure. As fluid deformation governs transport, mixing and dispersion, a critical step in the development of upscaled models of physical phenomena is the determination of fluid deformation as a function of these material properties.

Whilst fluid deformation in unsteady, turbulent flows has been widely studied for over half a century [START_REF] Cocke | Turbulent hydrodynamic line stretching: Consequences of isotropy[END_REF][START_REF] Wmtkerstein | Alignment of vorticity and scalar gradient with strain rate in simulated navier-stokes turbulence[END_REF][START_REF] W09407 | Material-element deformation in isotropic turbulence[END_REF][START_REF] Charles | Lagrangian dynamics and models of the velocity gradient tensor in turbulent flows[END_REF][START_REF] Simonkrstulovic | Turbulent pair dispersion as a continuous-time random walk[END_REF], much less attention has been paid to deformation in random steady flows, which have a distinctly different character. The steady nature of these flows imposes an important constraint upon the evolution of fluid deformation in that one of the principal stretches of the deformation gradient tensor must coincide with the local velocity direction [START_REF] Tabor | Stretching and Alignment in General Flow Fields: Classical Trajectories from Reynolds Number Zero to Infinity[END_REF]. Furthermore, the magnitude of this principal deformation is directly proportional to the local velocity, hence there is no net fluid stretching along the flow direction. This constraint, henceforth referred to as the steady flow constraint, means that fluid deformation of material lines and surfaces can be highly anisotropic with respect to the flow direction, leading to e.g. different mixing dynamics longitudinal and transverse to the mean flow direction. This difference in deformation dynamics is illustrated by the evolution of material elements in the mean translational flow shown in Figure 1. Here a continuously injected material line ℋ 1D 1 (green line) evolves with the mean flow direction 〈v〉 to form the steady 2D material sheet ℋ 2D , the cross-section of which ℋ 1D

(3) (blue line) can evolve into a complex lamellar shape due to fluid deformation transverse to the mean velocity field. In contrast, if we consider the material line ℋ 1D 1 to be instantaneously injected at the inlet plane, this line evolves into the 1D material line ℋ 1D

(2) (red line) which resides within the 2D surface ℋ 2D .

The line ℋ 1D

(2) deforms due to fluid deformation parallel to the local flow direction, leading to the fingering pattern observed in Figure 1, whereas ℋ 1D (3) evolves solely due to stretching transverse to the mean flow direction.

We denote fluid deformation in the longitudinal and transverse directions as Λ l , Λ t respectively, which directly govern longitudinal and transverse mixing and dispersion. In steady 2D flow, deformation transverse to the local flow direction (Λ t ) only fluctuates with time, whilst longitudinal deformation (Λ l ) can grow without bound [START_REF] Attinger | Exact transverse macro dispersion coefficients for transport in heterogeneous porous media[END_REF]).

Due to the steady flow constraint, fluid deformation evolves via different mechanisms depending upon whether it is parallel or transverse to the local velocity direction. These differences mean that in general, fluid deformation in steady 3D flows cannot be described in terms of a single scalar quantity (such as mean stretching rate).

As such, models of fluid mixing longitudinal and transverse to the mean flow must reference appropriate components of the fluid deformation tensor, and these can have very different evolution rates. For example, Le [START_REF] Borgne | Stretching, coalescence, and mixing in porous media[END_REF][START_REF] Borgne | The lamellar description of mixing in porous media[END_REF] show that longitudinal deformation Λ l governs mixing of a pulsed tracer injection in steady 2D Darcy flow. Conversely, [START_REF] Lester | Chaotic mixing in three-dimensional porous media[END_REF] illustrate that transverse fluid deformation Λ t governs mixing of a continuously injected source in a steady 3D flow at the pore scale, and Cirpka et al. (2011) show that the same holds for continuously injected sources in 2D at the Darcy scale.

In contrast, the unsteady nature of turbulent flows means that the deformation dynamics are largely isotropic with respect to the flow direction, whilst quantities such as vorticity tend to align with principal strain directions (Meneveau 2011). As a result, stochastic modelling of material deformation has focused on the measurement and development of stochastic models for orientation-invariant measures of deformation, typically in terms of invariants of the strain rate and velocity gradient tensors [START_REF] Charles | Lagrangian dynamics and models of the velocity gradient tensor in turbulent flows[END_REF][START_REF] Simonkrstulovic | Turbulent pair dispersion as a continuous-time random walk[END_REF]. Hence these models are not capable of resolving the constraints associated with steady flows. Conversely, whilst deformation in steady flow has been the subject of a number of studies [START_REF] Adachi | Calculation of strain histories in Protean coordinate systems[END_REF][START_REF] Adachi | A note on the calculation of strain histories in orthogonal streamline coordinate systems[END_REF][START_REF] Finnigan | A streamline coordinate system for distorted two-dimensional shear flows[END_REF][START_REF] Finnigan | Streamline coordinates, moving frames, chaos and integrability in fluid flowTopological Fluid Mechanics[END_REF][START_REF] Tabor | Stretching and Alignment in General Flow Fields: Classical Trajectories from Reynolds Number Zero to Infinity[END_REF], the development of stochastic models of deformation in these flows is an outstanding problem. To begin we first outline the stochastic modelling approach utilised in this study in §2, and identify a moving streamline coordinate frame, termed a Protean frame, to couch this model.

In §3 we consider evolution of fluid deformation in the Protean frame and derive expressions for evolution of the deformation tensor components in this frame. In §4 we derive and solve the coupled CTRW model for fluid deformation in the Protean frame coordinates and apply the stochastic model to several model steady flows. In §5 the velocity and velocity gradient statistics of several model flows in the Protean frame are examined, and model predictions are compared with direct computations. These results are discussed and then conclusions are made in §6.

Stochastic Modelling of Fluid Deformation in Steady Flows

Steady flows, whether 2D or 3D, involve topological and kinematic constraints [START_REF] Tabor | Stretching and Alignment in General Flow Fields: Classical Trajectories from Reynolds Number Zero to Infinity[END_REF] which impact the deformation of material elements. To develop a general stochastic modelling framework which honours these constraints, it is necessary to resolve the deformation tensor in a frame relative to the local velocity field. By utilising a coordinate frame which aligns with the streamlines of the flow, the topological constraints upon fluid deformation associated with steady flow are automatically imposed, and the stretching and shear contributions to deformation are clearly resolved. Specifically, we seek to develop stochastic models for the evolution of Λ t and Λ l in various flow regimes and across timescales which observe these kinematic constraints.

There exist several choices for streamline coordinate frames in steady flows, ranging from moving Cartesian coordinate frames to curvilinear streamline coordinate systems and convected coordinates. [START_REF] Adachi | Calculation of strain histories in Protean coordinate systems[END_REF][START_REF] Adachi | A note on the calculation of strain histories in orthogonal streamline coordinate systems[END_REF] presents a rotating Cartesian coordinate frame, termed a Protean coordinate frame, to study fluid deformation in 2D planar and axisymmetric flows. Similar to other streamline coordinate frames (such as Frenet-Serret coordinates [START_REF] Finnigan | Streamline coordinates, moving frames, chaos and integrability in fluid flowTopological Fluid Mechanics[END_REF] or curvilinear streamline coordinates), this frame automatically imposes the topological and kinematic constraints outlined in §1 that are inherent to steady flows.

In contrast to these coordinate systems, the moving Protean coordinate frame renders the velocity gradient tensor upper triangular. As shown in Appendix B, orthogonal streamline coordinate systems do not possess this property, and moreover it appears difficult to define a consistent streamline coordinate system for chaotic flows [START_REF] Finnigan | Streamline coordinates, moving frames, chaos and integrability in fluid flowTopological Fluid Mechanics[END_REF]). The upper triangular form of the velocity gradient in the Protean frame greatly simplifies solution of the deformation gradient tensor, which may now be expressed in terms of definite integrals of the velocity gradient components [START_REF] Adachi | Calculation of strain histories in Protean coordinate systems[END_REF][START_REF] Adachi | A note on the calculation of strain histories in orthogonal streamline coordinate systems[END_REF]. This representation then facilitates the development of stochastic models of fluid deformation, based upon solution of these integrals as a random walk process along fluid streamlines. Dentz et al. (2016b) used

the Protean frame as a the basis for the development of a stochastic model of fluid deformation in random steady 2D flows. In this study we are concerned with the extension of this approach to random steady 3D flow.

Unlike steady 2D flow, the moving streamline coordinate frame is not unique for steady 3D flow due to the additional degree of freedom associated with orientation of the frame about a streamline. As shall be shown in §3.2, the velocity gradient tensor is also no longer upper triangular for any arbitrary orientation about a streamline. However, we extend the Protean coordinate frame of [START_REF] Adachi | Calculation of strain histories in Protean coordinate systems[END_REF][START_REF] Adachi | A note on the calculation of strain histories in orthogonal streamline coordinate systems[END_REF] to arbitrary steady 3D flows by choosing the orientation angle that renders the velocity gradient upper triangular, recovering the integral solutions for the deformation tensor components. In some sense this 3D Protean frame for calculating fluid deformation in steady flow is analogous to continuous QR decomposition methods [START_REF] Lucavan | On the error in QR integration[END_REF][START_REF] Marcolester | Coupled continous time random walks for fluid stretching in two-dimensional heterogeneous media[END_REF]) for autonomous dynamical systems, which automatically recover constraints such as regularity and volume preservation of phase space.

In contrast to 2D flow, 3D steady flow also admits the possibility of chaotic fluid trajectories and exponential fluid deformation. In this context, an appropriate framework is needed which can also correctly capture either the Lyapunov exponents (exponential stretching rates) associated with chaotic flows or the power-law indices associated with algebraic fluid stretching in non-chaotic flows (such as isotropic Darcy flow). In §3.2 we show that the ensemble average of the diagonal components of the velocity gradient in the Protean frame correspond directly to the Lyapunov exponents of the flow. With respect to non-chaotic flows, special care is also required to suppress spurious exponential stretching in stochastic models, even for steady 2D flows. For example, [START_REF] Duplat | A nonsequential turbulent mixing process[END_REF] discuss a series of stretching mechanisms to explain observed sub-exponential (algebraic) stretching across a broad range of flow. Rather, algebraic stretching is inherent to steady 2D flow as a direct result of the topological constraints of these flows (associated with the Poincaré-Bendixson theorem). Dentz et al. (2016b) show that the use of streamline coordinates for 2D flows automatically recovers these constraints without resort to any specialised stretching mechanism. In Appendix B we show that certain components of the velocity gradient tensor are rendered zero for non-chaotic 3D flows such as complex lamellar [START_REF] Finnigan | A streamline coordinate system for distorted two-dimensional shear flows[END_REF][START_REF] Kelvin | Reprint of Papers on Electrostatics and Magnetism[END_REF] or zero helicity density [START_REF] Moffatt | The degree of knottedness of tangled vortex lines[END_REF]) flows, automatically recovering the kinematic constraints associated with these steady flows. Hence the Protean frame automatically adheres to the topological and kinematic constraints associated with chaotic and non-chaotic steady flows, so seems well-suited to the development of stochastic models of fluid deformation in these flows.

In steady random flows, Lagrangian velocities along streamlines fluctuate on a scale set by the random flow spatial structure. Hence, both the fluid velocity and velocity gradient tend to follow a spatial Markov process along streamlines (Le Borgne et al. 2008a,b), as these properties decorrelate in space. By exploiting this property and the Protean coordinate frame, we derive a stochastic deformation model in steady random flows as a coupled continuous time random walk (CTRW). The CTRW framework allows accounting for the broad velocity distributions that characterise steady random flows and lead to anomalous transport behavior (Dentz et al. 2016a). A similar approach has been applied previously by Dentz et al. (2016b) to steady 2D flows, leading to significant insights into the mechanisms which control fluid deformation in steady flow, namely that fluid deformation evolves as a coupled Lévy process due to the coupling between shear deformation and velocity fluctuations along streamlines.

For 3D flows, the possibility of exponential fluid stretching augments the couplings between shear, stretching, velocity and deformation leading to wholly new deformation dynamics in steady 3D flows, which are captured by the presented CTRW model. By solving this CTRW model, we can then make quantitative predictions of longitudinal Λ l and transverse Λ t fluid deformation in steady 3D random flows from statistical characterisation of the fluid velocity and velocity gradients. This facilitates the quantification of fluid deformation (and hence mixing, reaction and dispersion) directly from e.g. experimental or computational flow datasets. We begin by consideration of fluid deformation in the Protean frame.

Fluid Deformation in the Protean Coordinate Frame

In this section we extend the Protean Coordinate frame of [START_REF] Adachi | Calculation of strain histories in Protean coordinate systems[END_REF][START_REF] Adachi | A note on the calculation of strain histories in orthogonal streamline coordinate systems[END_REF] to 3D steady flows and derive the evolution equations for the deformation gradient in streamline coordinates.

General Development

The evolution of the location x(t) of fluid element in the steady spatially heterogeneous flow field u(x) is given by the advection equation

dx(t; X) dt = v(t), v(t; X) ≡ u[x(t; X)], (3.1)
where we denote the reference material vectors in the Eulerian and Lagrangian frames by x(t) and X, respectively. The deformation gradient tensor F(t) quantifies how the infinitesimal vector dx(t; X) deforms from its reference state dx(t = 0; X) = dX as

dx = F(t) ⋅ dX (3.2)
and equivalently

F i j (t) ≡ ∂x i (t; X) ∂X j . (3.3)
Following this definition, the deformation gradient tensor F(t) evolves with travel time t along a Lagrangian trajectory (streamline) as

dF(t) dt = ϵ(t)F(t), F(0) = 1, (3.4)
where the velocity gradient tensor ϵ(t) = ∇v(t) ⊤ ≡ ∇u[x(t; X)] ⊤ and the superscript ⊤ denotes the transpose. Rather than use a curvilinear streamline coordinate system, we consider the transform of a fixed orthogonal Cartesian coordinate system into the Protean coordinate frame which consists of moving Cartesian coordinates. We denote spatial coordinates in the fixed Cartesian coordinate system as x = (x 1 , x 2 , x 3 ) ⊤ , and the reoriented Protean coordinate system by x′ = (x 1 ′ , x 2 ′ , x 3 ′ ) ⊤ , where the velocity vector v(x) = (v 1 , v 2 , v 3 ) ⊤ in the Cartesian frame transforms to v′(x′) = (v, 0, 0) ⊤ in the Protean frame, with v = |v|.

These two frames are related by the transform [START_REF] Truesdell | The non-linear field theories of mechanics[END_REF])

x′ = x 0 (t) + Q ⊤ (t)x, (3.5)
where x 0 (t) is an arbitrary translation vector, and Q(t) is a rotation matrix (or proper orthogonal transformation):

Q ⊤ (t) • Q(t) = 1, det[Q(t)] = 1.
From (3.5) the differential elements dx, dX transform respectively as dx′ = Q ⊤ (t)dx, dX′ = Q ⊤ (0)dX, and so from (3.2) the deformation gradient tensor then transforms as

F′(t) = Q ⊤ (t)F(t)Q(0) . (3.6)
Formally, F(t) is a psuedo-tensor [START_REF] Truesdell | The non-linear field theories of mechanics[END_REF]), as it is not objective as per (3.6), [START_REF] Ottino | The Kinematics of Mixing: Stretching, Chaos, and Transport[END_REF]. Differentiating (3.6) with respect to the Lagrangian travel time t yields

F′(t) ≠ Q ⊤ (t)F(t)Q(t)
dF′ dt = Q ˙⊤(t)F(t)Q(0) + Q ⊤ (t)ϵ(t)F(t)Q(0) = ϵ′(t)F′(t), (3.7)
where the transformed rate of strain tensor ϵ′(t) is defined as

ϵ′(t) ≡ Q ⊤ (t)ϵ(t)Q(t) + A(t), (3.8)
and the contribution due to a moving coordinate frame is

A(t) ≡ Q ˙⊤(t)Q(t) . As Q(t) is orthogonal, then Q ˙⊤(t)Q(t) + Q ⊤ (t)Q ˙(t) = 0
and so A(t) is skew-symmetric. Note the velocity vector v(t) = dx/dt also transforms as v′(t) = Q ⊤ (t)v(t). The basic idea of the Protean coordinate frame is to find an appropriate local reorientation Q(t) that renders the transformed velocity gradient ϵ′(t) upper triangular, simplifying solution of the deformation tensor F in (3.7).

Fluid Deformation in 3D Protean

Coordinates [START_REF] Adachi | Calculation of strain histories in Protean coordinate systems[END_REF][START_REF] Adachi | A note on the calculation of strain histories in orthogonal streamline coordinate systems[END_REF] shows that for For steady 2D flow, reorientation into streamline coordinates automatically yields an upper triangular velocity gradient tensor, hence the reorientation matrix Q(t) is simply given in terms of the velocity v = [v 1 , v 2 ] ⊤ in the fixed Cartesian frame as

Q(t) = 1 v v 1 -v 2 v 2 v 1 , (3.9) hence v′ = [v, 0] ⊤ , where v = v 1 2 + v 2 2 .
In contrast to steady 2D flow, the additional spatial dimension of steady 3D flows admits an additional degree of freedom for the specification of streamline coordinates. In 3D streamline coordinates, the base vector e 1 ′ (t) aligns with the velocity vector v(t), such that v′ (t) = [v(t), 0, 0] ⊤ , however the transverse vectors e 2 ′ (t) and e 3 ′ (t) are arbitrary up to a rotation about e 1 ′ (t) . As these base vectors are not material coordinates, this degree of freedom does not impact the governing kinematics.

For the Protean coordinate frame it is required to reorient the rate of strain tensor ϵ′(t) such that it is upper triangular, yielding explicit closed-form solution for the deformation gradient tensor F′(t) and elucidating the deformation dynamics. Whilst all square matrices are unitarily similar (mathematically equivalent) to an upper triangular matrix, it is unclear whether such a similarity transform is orthogonal, or corresponds to a reorientation of frame, or indeed corresponds to a reorientation into streamline coordinates. As demonstrated by [START_REF] Adachi | Calculation of strain histories in Protean coordinate systems[END_REF][START_REF] Adachi | A note on the calculation of strain histories in orthogonal streamline coordinate systems[END_REF], this condition is satisfied for all steady 2D flows, but this is still an open question for steady 3D flows.

For reasons that shall become apparent, we develop the 3D Protean coordinate frame by decomposing the reorientation matrix

Q(t) into two sequential reorientations Q 1 (t), Q 2 (t) as Q(t) ≡ Q 2 (t)Q 1 (t), (3.10)
where Q 1 (t) is a reorientation that aligns the velocity vector v with the Protean base vector e 1 ′ (t), and Q 2 (t) is a subsequent reorientation of the Protean frame about e 1 ′ (t) . As such, the reorientation Q 1 (t) renders the Protean frame a streamline coordinate frame, and Q 2 represents a (currently) arbitrary reorientation of this streamline frame about the velocity vector.

We begin by considering the 3D rotation matrix Q 1 (t) which is defined in terms of the unit rotation axis q(t) (which is orthogonal to both the velocity and e 1 vectors), and the associated reorientation angle θ(t), both of which are given in terms of the local velocity vector

v(t) = [v 1 , v 2 , v 3 ] ⊤ as q(t) = e 1 × v ∥ e 1 × v ∥ = 1 v 2 2 + v 3 2 {0, v 3 , -v 2 }, (3.11) cos θ(t) = e 1 ⋅ v ∥ e 1 ⋅ v ∥ = v 1 v . (3.12)
The rotation tensor Q 1 (t) then is given by

Q 1 (t) = cos θ(t)𝕀 + sin θ(t)(q) × ⊤ + [1 -cos θ(t)]q(t) ⊗ q(t), (3.13)
where 𝕀 is the identity matrix, (•) × denotes the cross product matrix and || • || is the ℓ 2 -norm.

The second reorientation matrix Q 2 (t) corresponds to an arbitrary reorientation (of angle α(t)) about the e 1 ′ axis, which is explictly

Q 2 (t) = 1 0 0 0 cos α(t) -sin α(t) 0 sin α(t) cos α(t) . (3.14)
Whilst the algebraic expressions (in terms of α, v 1 , v 2 , v 3 ) for the components of the composition rotation matrix

Q(t) = Q 2 (t)Q 1 (t)
are somewhat complicated, it is useful to note that Q(t) may be expressed in terms of the Protean frame basis vectors as

Q(t) = [e 1 ′ (t), e 2 ′ (t), e 3 ′ (t)], (3.15)
hence the moving coordinate contribution A(t) is given by

A(t) = Q ˙T(t)Q(t) = 0 e 2 ′ (t) ⋅ e ˙1 ′ (t) e 3 ′ (t) ⋅ e ˙1 ′ (t) -e 2 ′ (t) ⋅ e ˙1 ′ (t) 0 e 3 ′ (t) ⋅ e ˙2 ′ (t)
ϵ (t) ≡ Q T (t) ϵ (t)Q(t) = ϵ ′(t) -A(t), (3.17)
we find the components of A(t) may be related to ϵ t as follows. From (3.17) and (3.15), where from (3.13), (3.14)

ϵ i j (t) = e i ′(t) ϵ (t)e j ′(t), and since v ˙(t) = ϵ(t) ⋅ v(t), e 1 ′ (t) = v(t)/v(t), then e 2 ′ (t) ⋅ e ˙1 ′ (t) = e 2 ′ (t) ⋅ ϵ(t) ⋅ e 1 ′ (t) = ϵ 21 (t), (3.18) e 3 ′ (t) ⋅ e ˙1 ′ (t) = e 3 ′ (t) ⋅ ϵ(t) ⋅ e 1 ′ (t) = ϵ 31 (t)
ve 3 ′ ⋅ ∂e 2 ′ ∂v = 1 v + v 1 {0, v 3 , -v 2 } = v 3 cosα -v 2 sinα v + v 1 e 2 ′ - v 2 cosα + v 3 sinα v + v 1 e 3 ′ , (3.22) e 3 ′ ⋅ ∂e 2 ′ ∂α = -1, (3.23)
and so A 23 (t) is then

A 23 (t) = e 3 ′ ⋅ e ˙2 ′ = v 3 cos α -v 2 sin α v + v 1 ϵ 21 - v 2 cos α -v 3 sinα v + v 1 ϵ 31 - dα dt . (3.24)
Hence the equation A 23 t = ϵ 32 defines an ordinary differential equation (ODE) for the arbitrary orientation angle α(t) which renders ϵ′(t) upper triangular. Although the components of ϵ t vary with α, these may be expressed in terms of the components of ϵ (1) (which are independent of α), defined as

ϵ 1 t ≡ Q 1 T t ϵ t Q 1 t , (3.25)
where the components of ϵ t and ϵ (1) (t) are related as

ϵ 21 = ϵ 21 1 cos α + ϵ 31 1 sinα, (3.26) ϵ 31 = ϵ 31 1 cosα -ϵ 21 1 sinα, (3.27) ϵ 32 = ϵ 32 1 cos 2 α -ϵ 23 1 sin 2 α + ϵ 33 1 -ϵ 22 1 cos α sinα . (3.28)
Hence the condition A 23 = ϵ 32 is satisfied by the ODE

dα dt = g α, t = a t cos 2 α + b t sin 2 α + c t cosα sin α, (3.29)
where

a t = -ϵ 32 1 - v 2 v + v 1 ϵ 31 1 - v 3 v + v 1 ϵ 21 1 , (3.30) b t = ϵ 23 1 + v 2 v + v 1 ϵ 31 1 + v 3 v + v 1 ϵ 21 1 , (3.31) c t = ϵ 22 1 -ϵ 33 1 + 2v 2 v + v 1 ϵ 21 1 - 2v 3 v + v 1 ϵ 31 1 .
(3.32)

Evolution of Protean Orientation

Angle α(t)-Equation (3.29) describes a 1st-order ODE for the transverse orientation angle α(t) along a streamline which renders ϵ′ (t) upper triangular. Here the temporal derivative for α(t) is associated with both change in flow structure along a streamline and impact of a moving coordinate transform as encoded by A(t). This is analogous to the ODE system (A 5) for the continuous QR method (Appendix A), with the important difference that the Protean reorientation gives a closed form solution for e 1 = v/v which constrains two degrees of freedom (d.o.f) for the Protean frame and an ODE for the remaining d.o.f characterised by α(t), whereas the continuous QR method involves solution of the 3 degree of freedom ODE system (A 5), and requires unitary integrators to preserve orthogonality of 𝒬. These methods differ with respect to the initial

conditions 𝒬 0 = 𝕀, Q(0) = Q 1 (0) • Q 2 (0)
, where Q 1 (0) is given explicitly by (3.13), and Q 2 (0) is dependent upon the initial condition α(0) = α 0 as per (3.14), (3.29).

Whilst the ODE (3.29) may be satisfied for any arbitrary initial condition α(0) = α 0 , resulting in non-uniqueness of ϵ′(t), this non-autonomous ODE is locally dissipative as reflected by the divergence

∂g ∂α = b t -a t sin2α + c t cos2α, (3.33)
which admits maxima and minima of magnitude ±c t 1 + b t -a t /c t 2 respectively at

α t = 1 2 arctan b t -a t c t + π 2 sgn c t ∓ 1 2 . (3.34)
Hence for c(t) ≠ 0, the ODE (3.29) admits an attracting trajectory ℳ(t), which attracts solutions from all initial conditions α 0 as illustrated in Fig. 2(a). As this trajectory is also a solution of (3.29), it may be expressed as ℳ(t) = α(t; α 0,∞ ), where α 0,∞ is the initial condition at t = 0 which is already on ℳ(t). Whilst all solutions of (3.29) render ϵ′(t) upper triangular, only solutions along the attracting trajectory ℳ(t) represent asymptotic dynamics independent of the initial condition α 0 , hence we define the Protean frame as that which corresponds to ℳ(t) = α(t; α 0,∞ ).

As the attracting trajectory ℳ(t) maximizes dissipation -∂g/∂α over long times, the associated initial condition α 0,∞ is quantified by the limit α 0,∞ = lim τ→∞ α 0,τ , where α 0, τ t = argmin . This approach is particularly useful when Lagrangian data is only available over short times and explicitly identifies the inertial initial orientation angle α 0,∞ , allowing the Protean transform to be determined uniquely from t = 0.

Longitudinal and Transverse Deformation in 3D Steady

Flow-Given solution of (3.29) (such that A 23 = ϵ 32 ), reorientation into the Protean frame renders the streamline velocity gradient tensor ϵ′(t) upper triangular. From (3.7) it follows that F′(t) is also upper triangular, with diagonal components

F 11 ′ t = v t v 0 , (3.36a) F ii ′ t = exp ∫ 0 t dt′ ϵ ii ′ t′ , (3.36b)
for i = 2, 3. The off-diagonal elements are F ij (t) = 0 for i > j and else

F 12 ′ t = v t ∫ 0 t dt′ ϵ 12 ′ t′ F 22 ′ t′ v t′ , (3.36c) F 23 ′ t = F 22 ′ t ∫ 0 t dt′ ϵ 23 ′ t′ F 33 ′ t′ F 22 ′ t′ , (3.36d
)

F 13 ′ t = v t ∫ 0 t dt′ ϵ 12 ′ t′ F 23 ′ t′ + ϵ 13 ′ t′ F 33 ′ t′ v t′ . (3.36e)
Here the upper triangular form of ϵ′(t) simplifies solution of the deformation tensor F′(t)

such that the integrals (3.36) can be solved sequentially, rather than the coupled ODE (3.4).

It is interesting to note that the integrands of the shear deformations F 12 ′ t and F 13 ′ t are weighted by the inverse velocity. This implies that episodes of low velocity lead to an enhancement of shear induced deformation, see also the discussion in Dentz et al. (2016b) for 2D steady random flows. The impact of intermittent shear events in low velocity zones can be quantified using a stochastic deformation model based on continuous time random walks as outlined in the next section. For compactness of notation we henceforth omit these primes with the understanding that all quantities are in the Protean frame unless specified otherwise.

To illustrate how the deformation tensor controls longitudinal and transverse stretching of fluid elements, we decompose F into longitudinal and transverse components respectively as

F(t) = F l + F t ,
where

F l ≡ diag e 1 ⋅ F = F 11 F 12 F 13 0 0 0 0 0 0 , (3.37) F t ≡ diag e 2 + e 3 ⋅ F = 0 0 0 0 F 22 F 23 0 0 F 33 , (3.38)
where diag(a) is a diagonal matrix comprised of the vector a along the diagonal. From (3.2), a differential fluid line element δl(X, t) at Lagrangian position X then evolves with time t as

δl X, t = F X, t ⋅ δl X, 0 = F l X, t ⋅ δl X, 0 + F t X, t ⋅ δl X, 0 = δl l X, t + δl t X, t , (3.39)
where the line element may also be decomposed into the longitudinal and transverse components as δl l = δl l + δl t . Due to the upper triangular nature of F the length δl of these line elements can also be decomposed as

δl X, t ≡ δl X, t = δl X, 0 ⋅ F ⊤ X, t ⋅ F X, t ⋅ δl X, 0 = δl X, 0 ⋅ F l ⊤ X, t ⋅ F l X, t + F t ⊤ X, t ⋅ F t X, t ⋅ δl X, 0 = δl l X, t 2 + δl t X, t 2 .
(3.40)

Hence the total length of the line element is decomposed into longitudinal and transverse contributions as

δl 2 = δl l 2 + δl t 2 . (3.41)
We denote the length of the 1D lines ℋ 1D 2 , ℋ 1D

3 in Figure 1 respectively as

l 2 t , l 3 x 1 ,
where the respective arguments t, x 1 reflect the fact that 1D line ℋ 1D

(2) occurs at fixed time t since injection, whereas the 1D line ℋ 1D

(3) occurs at fixed distance x 1 downstream in the mean flow direction. Following the decomposition (3.41), l (2) (t) and l (3) x 1 are then given by the contour integrals along the injection line ℋ 1D (1) in Lagrangian space as (3.43) where t 1 X, x 1 is the time at which the fluid streamline at Lagrangian coordinate X reaches

l (2) (t) = ∫ ℋ 1D (1) δl l X, t 2 + δl t X, t 2 ds, (3.42) l (3) (x 1 ) = ∫ ℋ 1D (1) δl l X, t 1 X, x 1 ds,
x 1 , and ds is a differential increment associated with changes in X along ℋ 1D (1) . Hence stretching of the 1D material line ℋ 1D

(2) in Figure 1 is governed by both the transverse and longitudinal deformations F l , F t , whereas stretching of the 1D line ℋ 1D

(2) transverse to the mean flow direction is solely governed by the transverse deformation F t .

As discussed in §1, the deformation of these different lines is important for various applications. For example, Le [START_REF] Borgne | Stretching, coalescence, and mixing in porous media[END_REF][START_REF] Borgne | The lamellar description of mixing in porous media[END_REF] use the growth rate of l (2) (t) to predict the mixing of a pulsed tracer injection (illustrated as ℋ 1D

(2) in Figure 1 As such, we characterise longitudinal Λ l (t) and transverse Λ t (t) fluid deformation in terms of the corresponding deformation tensors as

Λ l t ≡ ∥ F l t ∥ = 〈 F 11 t 2 + F 12 t 2 + F 13 t 2 〉, (3.44) Λ t t ≡ ∥ F t t ∥ = 〈 F 22 t 2 + F 23 t 2 + F 33 t 2 〉, (3.45)
where the angled brackets denote an ensemble average over multiple realisations of a 3D steady random flow. We also define the ensemble-averaged length of an arbitrary fluid line as

ℓ t ≡ δl X, t , (3.46)
and in §4, we seek to derive stochastic models for the growth of Λ l (t), Λ t (t) and ℓ(t) in steady random 3D flows. In 2D steady flow, topological constraints (associated with conservation of area) prevent persistent growth of transverse deformation which simply fluctuates as

Λ t (t) = 〈F 22 (t)〉 = 1/〈F 11 (t)〉 = 〈v(0)/v(t)〉 = 1.
Conversely, longitudinal deformation can grow persistently as Λ l (t) ~ F 12 (t) 2 , and Dentz et al. (2016b) derive stochastic models for the evolution of Λ l (t) from velocity field statistical properties as a CTRW. In §4 we extend this approach to steady random 3D flows to yield stochastic models for Λ l (t), Λ t (t) and ℓ(t), which act as important inputs for models of fluid mixing and dispersion.

For illustration, we briefly consider the deformation along the axis of the streamline coordinate system. For a fluid element z(t) initially aligned with the 1-direction of the Protean coordinate system, i.e. z(0) = (z 0 , 0, 0) ⊤ we obtain the elongation

ℓ t = z t 2 ℓ t = z 0 v t v 0 . (3.47)
For ergodic flows (i.e. where any streamline eventually samples all of the flow structure), the average stretching is zero and elongation will tend asymptotically toward a constant. For a material element initially orientated along the 2-direction of the Protean coordinate system, i.e z(0) = (0, z 0 , 0) ⊤ we obtain the elongation

ℓ t = z 0 F 12 t 2 + F 22 t 2 . (3.48)
For an initial alignment with the 3-direction, z(0) = (0, 0, z 0 )

⊤ ℓ t = z 0 F 13 t 2 + F 23 t 2 + F 33 t 2 . (3.49)
In general, we have

ℓ t = z 0 F ⊤ t F t z 0 , (3.50)
where z(t = 0) ≡ z 0 . Thus, net elongation is only achieved for initial orientations away from the velocity tangent. For both chaotic and non-chaotic flows, the explicit structure of F′(t) facilitates the identification of the stochastic dynamics of fluid deformation in steady random flows and their quantification in terms of the fluid velocity and velocity gradient statistics.

Zero Helicity Density Flows

One important class of steady D=3 dimensional flows are zero helicity density (or complex lamellar [START_REF] Finnigan | A streamline coordinate system for distorted two-dimensional shear flows[END_REF])) flows such as the isotropic Darcy flow (Sposito 2001)

v x = -k x ∇ϕ x , (3.51)
which is commonly used to model Darcy-scale flow in heterogeneous porous media. Here v(x) is the specific discharge, the scalar k(x) represents hydraulic conductivity, and ϕ(x) is the flow potential (or velocity head). For such flows the helicity density [START_REF] Moffatt | The degree of knottedness of tangled vortex lines[END_REF] which measures the helical motion of the flow

h x ≡ v ⋅ ω = k x ∇ϕ x ⋅ ∇k x × ∇ϕ x , (3.52)
is identically zero (where ω is the vorticity vector). As shown by [START_REF] Kelvin | Reprint of Papers on Electrostatics and Magnetism[END_REF], all 3D zero helicity density flows are complex lamellar, and so may be posed in the form of an isotropic Darcy flow v(x) = -k(x)∇ϕ, where for porous media flow ϕ is the flow potential and k(x)

represents the (possibly heterogeneous) hydraulic conductivity. [START_REF] Arnol | On the topology of three-dimensional steady flows of an ideal fluid[END_REF]Arnol'd ( , 1965) ) shows that streamlines in complex lamellar steady zero helicity density flows are confined to a set of two orthogonal D=2 dimensional integral surfaces which can be interpreted as level sets of two stream functions ψ 1 , ψ 2 , with ∇ψ 1 • ∇ψ 2 = 0. As such, all zero helicity density flows may be represented as

v x = -k x ∇ϕ = ∇ψ 1 × ∇ψ 2 . (3.53)
As the gradients ∇ϕ, ∇ψ 1 , ∇ψ 2 are all mutually orthogonal, zero helicity density flows admit an orthogonal streamline coordinate system (ϕ, ψ 1 , ψ 2 ) with unit base vectors

e 1 = - ∇ϕ ∥ ∇ϕ ∥ = v ∥ v ∥ , e 2 = ∇ψ 1 ∥ ∇ψ 1 ∥ , e 3 = ∇ψ 2 ∥ ∇ψ 2 ∥ . (3.54)
If the streamfunctions ψ 1 , ψ 2 are known, there is no need to explicitly solve the orientation angle α(t) via (3.29), but rather the Protean frame can be determined from the coordinate directions given by (3.54). Note that the moving Protean coordinate frame differs from an orthogonal curvilinear streamline coordinate system based on (3.54), and we show in Appendix B that such a coordinate system does not yield an upper triangular velocity gradient.

The zero helicity density condition imposes important constraints upon the Lagrangian kinematics of these flows. First, as steady zero helicity-density flows must be non-chaotic due to integrability of the streamsurfaces ψ 1 , ψ 2 (Holm & Kimura 1991), hence fluid deformation must be sub-exponential (algebraic). From equation (3.36b), the ensemble average of the principal components ϵ ii of the velocity gradient deformation (which correspond to the Lyapunov exponents of the flow) must be zero. Secondly, zero helicity flows constrain the velocity gradient components such that ϵ 23 = ϵ 32 , due to the identity

h = v ⋅ ω = v ⋅ ε i jk : ϵ = v ⋅ ε i jk : ϵ = v ϵ 23 -ϵ 32 = 0, (3.55)
where ε ijk is the Levi-Civita tensor. Whilst equation (3.55) shows ϵ 23 = ϵ 32 in general for zero helicity flows, there exists a specific value of α, denoted α, which renders both ϵ 23 and ϵ 32 zero. This form of the velocity gradient (namely zero (2,3) and (3,2) components) is reflected by its representation in curvilinear streamline coordinate system shown in Appendix B, and can lead to a decoupling of fluid deformation between the (1,2) and (1,3) planes. However, as the Protean frame is a moving coordinate system, there is a non-zero contribution from A(t) to the (2,3) and (3,2) components of the velocity gradient. Typically, a different value of α is used than α to render ϵ 32 ′ = 0, and so for zero helicity flow

ϵ 23 ′ = 2 ϵ 23 , (3.56)
which is typically non-zero, hence F 23 ′ ≠ 0 in general. Note that as both the Protean coordinate frame and the curvilinear streamline coordinate system render the deformation tensor component F 23 non-zero (see Appendix B for details), the coupling between fluid deformation in the (1,2) and ( 1,3) planes is retained in both coordinate frames, hence this is not purely an artefact of the moving coordinate system.

The confinement of streamlines to integral streamsurfaces given by the level sets of ψ 1 , ψ 2 effectively means steady 3D zero helicity density flows can be considered as two superposed steady D=2 dimensional flows. This has several impacts upon the transport and deformation dynamics. First, as the 2D integral surfaces are topological cylinders or tori, streamlines confined within these surfaces cannot diverge exponentially in space. Thus, the principal transverse deformations F 22 , F 33 in the Protean frame may only fluctuate about the unit mean as where the indices r are governed by intermittency of the fluid shear events (Dentz et al. 2016b). We note that these constraints to not apply to anisotropic Darcy flows (i.e. v(x) = -K(x) • ∇ϕ(x), where K(x) is the tensorial hydraulic conductivity), as these no longer correspond to zero helicity density flows. Indeed, as demonstrated by Ye et al. ( 2015), these flows can give rise to chaotic advection and exponential stretching of material elements. We consider fluid deformation in such steady chaotic flows in the following Section.

F ii = 1, i = 1,

Fluid Deformation Along Stagnation

Lines-For streamlines which connect with a stagnation point (termed stagnation lines), the solution of the deformation tensor (3.36) diverges in the neighbourhood of the stagnation point as v → 0. To circumvent this issue, we consider fluid deformation along a stagnation line in the neighbourhood of a stagnation point x p , where the velocity gradient may be approximated as ϵ(t) ≈ ϵ 0 ≡ ϵ| x=x p . Note for streamlines the approaching a stagnation point from upstream, ϵ 11,0 must necessarily be negative. In the neighbourhood of the stagnation point, the velocity evolves as v(t)/v(t 0 ) ≈ exp(ϵ 11,0 (t -t 0 )), and the deformation gradient tensor can be solved directly from (3.4) as

F t ≈ F t 0 ⋅ exp ϵ 0 t -t 0 , (3.63)
where t 0 is the time at which the fluid element enters the neighbourhood of the stagnation point (defined as the region for which ϵ(t) ≈ ϵ 0 ). Note that (3.36) represents an explicit solution of the matrix exponential above, via the substitutions t ↦ tt 0 , v(t) ↦ v(t 0 ) exp(ϵ 11,0 (t -t 0 )), ϵ(t) ↦ ϵ 0 . Whilst the stochastic model for fluid deformation developed herein does apply to points of zero fluid velocity, the impact of deformation local to stagnation points can be included a posterori via (3.63).

Stochastic Fluid Deformation in Steady Random 3D Flow

In this Section, we develop a stochastic model for fluid deformation along streamlines in steady D = 3 dimensional flows that exhibit chaotic advection and exponential fluid stretching. This model is based on a continuous time random walk (CTRW) framework for particle transport in steady random flows, which is briefly summarized in the following.

Continous Time Random Walks for Fluid Velocities

The travel distance s(t) of a fluid particle along a streamline is given by

ds t dt = v s s t , dt s ds = 1 v s s , (4.1)
where v s (s) = v[t(s)]. We note that fluid velocities v s (s) can be represented as a Markov process that evolves with distance along fluid streamlines (Dentz et al. 2016a) for flow through heterogeneous porous and fractured media from pore to Darcy and regional scales [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Fiori | Ergodic transport trough aquifers of non-gaussian log conductivity distribution and occurence of anomalous behavior[END_REF]Le Borgne et al. 2008a,b;[START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF][START_REF] Holzner | Intermittent lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF][START_REF] Kang | Anomalous transport on regular fracture networks: Impact of conductivity heterogeneity and mixing at fracture intersections[END_REF]. This property stems from the fact that fluid velocities evolve on characteristics spatial scales imprinted in the spatial flow structure. Streamwise velocities v s (s) are correlated on the correlation scale ℓ c . For distances larger than ℓ c subsequent particle vleocities can be considered independent. Thus, particle motion along streamlines can be represented by the following recursion relations,

s n + 1 = s n + ℓ c , t n + 1 = t n + ℓ c v n , (4.2)
where we defined v n = v(s n ) for the n-th step. The position s(t) along a streamline is given in terms of (4.2) by s(t) = s n t , where n t = sup(n|t n ≤ t), the particle speed is given by v(t) = n n t . This coarse-grained picture is valid for times larger than the advection time scale τ u = ℓ c /ū, where ū is the mean flow velocity. We consider flow fields that are characterised by open ergodic streamlines. Thus, v n can be modeled as independent identically distributed random variables that are characterised by p s (v). The latter is related to the PDF p e (v) of Eulerian velocities v e through flux weighting as (Dentz et al. 2016a)

p s (v) = vp e (v) v e . (4.3)
The advective transition time over the characteristic distance ℓ c is defined by In order to illustrate this notion, we consider the joint evolution of a process A n and time t n along a streamline [START_REF] Scher | Stochastic transport in a disordered solid. I. Theory[END_REF],

τ n = ℓ c v n . (
A n + 1 = A n + ΔA n, t n + 1 = t n + τ n . (4.6)
The increments (ΔA, τ) are stepwise independent and in general coupled. They are characterized by the joint PDF ψ a (a, t), which can be written as

ψ a (a, t) = ψ a (a t)ψ(t), (4.7)
The PDF p a (a, t) of A(t) is then given by

p a (a, t) = ∫ 0 t dt′R(a, t′) ∫ t -t′ ∞ dt″ψ(t) (4.8)
where R(a, t) denotes the probability per time that A(t) has just assumed the value a. Thus, Eq. (4.8) can be read as follows. The probability p a (a, t) that A(t) has the value a at time t is equal to the proabality per time R(a, t′) that A(t) has arrived at a times the probability that the next transition takes longer than t -t′. The R(a, t) satisfies the Chapman-Kolmogorov type equation a′, t′), (4.9) where R 0 (a, t) is the joint distribution of A 0 and t 0 . Equations (4.8) and (4.9) ca be combined into the following generalized master equation [START_REF] Kenkre | Generalized master equations for continuous-time random walks[END_REF] for the evolution of p a (t)

R(a, t) = R 0 (a, t) + ∫ 0 t dt′ ∫ daψ(a -a′, t -t′)R(
∂ p a (a, t) ∂t = ∫ 0 t dt′ ∫ daΦ(a -a′, t -t′)[p a (a′, t′) -p a (a, t′)], (4.10)

where Φ(a, t) is defined in Laplace space as

Φ * (a, λ ) = λ ψ a *(a, λ ) 1 -ψ * ( λ ) . (4.11)
The general master equation (4.10) expresses the non-Markovian character of the evolution of the statistics of A(t).

In the following, we develop a stochastic model for deformation based on these relations. We anticipate that the velocity gradient statistics also follow a spatial Markov process with similar correlation structure. Spatial Markovianity of both the fluid velocity and velocity gradient then provide a basis for stochastic modelling of the deformation equations (3.36b)-(3.36e) as a coupled continuous time random walk (coupled CTRW) along streamlines, in a similar fashion as that developed by Dentz et al. (2016b) for deformation in steady 2D random flow. Based on the explicit expressions (3.36) of the deformation tensor, this stochastic approach is derived from first principles and so provides a link to Lagrangian velocity and deformation, which in turn may be linked to the Eulerian flow properties and medium characteristics (Dentz et al. 2016a;[START_REF] Fiori | Ergodic transport trough aquifers of non-gaussian log conductivity distribution and occurence of anomalous behavior[END_REF][START_REF] Edery | Origins of anomalous transport in heterogeneous media: Structural and dynamic controls[END_REF][START_REF] Tyukhova | Mechanisms of anomalous dispersion in flow through heterogeneous porous media[END_REF].

To this end, we recast (3.36) in terms of the distance along the streamline by using the transformation (4.1) from t to s. This gives for the diagonal components of

F(s) = F[t(s)] F 11 (s) = v s (s) v s (0) , (4.12a) F ii (s) = exp ∫ 0 s ds′ ϵ ii (s′) v s (s′) , (4.12b)
with i = 2, 3. For the off-diagonal elements, we obtain while F i j (s) = 0 for i > j. We denote ϵ i j (s) = ϵ i j [t(s)] . In the following, we investigate the evolution of stretching in different time regimes for chaotic steady random flow.

F 12 (s) = v s (s) ∫ 0 s ds′ ϵ 12 (s′)F 22 (s′) v s (s′) 2 , (4.12c)

Chaotic Steady Random Flow

For chaotic flows, the infinite-time Lyapunov exponent is defined by

λ ≡ lim t ∞ lim z 0 0 1 2t ln ℓ t 2 ℓ 0 2 . (4.13)
In the limit of t ≫ λ -1 , we may set

F 11 t = 1, F 22 t = exp ϵ t , F 33 t = exp -ϵ t , (4.14)
where we defined

ϵ = lim t ∞ 1 t ∫ 0 t dt′ϵ 22 t′ = -lim t ∞ 1 t ∫ 0 t dt′ϵ 33 t′ . (4.15)
The last equation is due to volume conservation. The long time behavior is fully dominated by exponential stretching. The characteristic stretching time scale for exponential stretching is given by τ ϵ = 1/ϵ, while the characteristic advection scale is τ u = ℓ c /ū with ū the mean and ℓ c the correlation scale of the steady random flow field u(x). The latter sets the end of the ballistic regimes, see below. If τ ϵ and τ u are well separated, τ ϵ ≫ τ u , we observe a subexponential stretching regime that is dominated by heterogeneous shear action. In the following, we analyse the behavior of deformation for time t ≪ τ ϵ . We first briefly discuss the ballistic regime for which t ≪ τ u . Then, we develop a CTRW based approach for deformation in the pre-exponential regime τ u ≪ t ≪ τ ϵ .

Ballistic

Regime: t ≪ τ u -We first consider the ballistic short time regime, for which t ≪ τ u . In this regime, the flow properties are essentially constant and we obtain the approximation for the deformation tensor (3.36)

F ii t = 1, (4.16a)
for i = 1, 2, 3 and for the off-diagonal elements with i > j

F 12 t = ϵ 12 t, F 23 t = ϵ 23 t, F 13 t = ϵ 12 ϵ 23 t/2 + ϵ 13 t . (4.16b)
As ϵ ij ~ 1/τ u , in this regime the transverse (3.45) and longitudinal (3.44) deformations as well the arbitrary fluid element ℓ(t) all evolve ballistically to leading order; Λ t (t), Λ l (t), ℓ(t) ∝ t.

4.2.2 Exponential stretching regime t ≫ τ ϵ -In the stretching regime, for times much larger than τ ϵ , the principal deformations scale as

F 11 t = 1, F 22 t = exp ϵt , F 33 t = exp -ϵt . (4.17)
In this regime the expressions for the off-diagonal elements of F(t) in (3.36) can be approximated by

F 12 t ≈ ϵ 12 ′ t exp ϵt ϵ , (4.18)
F 23 t ≈ exp ϵt ∫ 0 t dt′ϵ 23 ′ t′ exp -2ϵt′ , (4.19)

F 13 t ≈ F 12 t ϵ 12 t ∫ 0 t dt′ϵ 23 ′ t′ exp -2ϵt′ . (4.20)
Note that the integrals above are dominated through the exponential by the upper limit, which gives these approximations. In this regime the longitudinal (3.44) and transverse deformations (3.45), along with the elongation of a material element (3.50) all increase exponentially with time; Λ l (t), Λ t (t), ℓ(t) ∝ exp(ϵt). From (4.13) the Lyapunov exponent is then λ ≡ ϵ.

4.2.3

Intermediate regime t > τ u -For times t > τ u , we develop a CTRW approach to describe the impact of flow heterogeneity on the evolution of the deformation tensor, based on (4.2) and (4.12). In this regime, we approximate F 11 (s) ≈ 1, F 22 (s) ≈ exp λt s and F 33 ≈ exp [-λt(s)]. The off-diagonal components are then expressed as

F 12 s = v s s ℓ c r 2 s , F 23 s = F 22 s p s , (4.21a) F 13 s = v s s ∫ 0 s ds′ ϵ 12 s′ F 23 s′ v s s′ 2 + v s ℓ c r 3 s , (4.21b)
where we defined the auxiliary functions

dr i s ds = ℓ c ϵ 1i s F ii s v s s 2 , i = 2, 3, d p s ds = ϵ 23 s F 33 s v s s F 22 s , (4.22a)
with r i (s = 0) = 0 and p(s = 0) = 0. We approximated the off-diagonal components ϵ i j (s) for i > j as being δ-correlated on the scales of interest such that

ϵ i j s = σ i j 2 ℓ c ξ i j s , (4.23)
where ξ(s) denotes a white noise with zero mean and correlation 〈ξ ij (s)ξ ij (s′)〉 = δ(s -s′), the σ ij are characteristic deformation rates. We assume here that the deformation rates are independent of velocity magnitude, an assumption that can be easily relaxed. Discretising (4.22) on the scale ℓ c such that s n = nℓ c and using (4.23), we obtain

r i, n + 1 = r i, n + τ n 2 σ 1i exp[ -1 i λt n ]ξ 1i, n , i = 2, 3 (4.24a) p n + 1 = p n + τ n σ 23 exp -2λt n ξ 23, n . (4.24b)
These recursive relations describe a fully coupled continuous time random walk [START_REF] Zaburdaev | Lévy walks[END_REF] with non-stationary increments due to the explicit dependence on time t n . Note

that F i j (t) = F i j [s n t ] = F i j (n t ℓ c ) .
Thus, from (4.24) we obtain for the mean square deformations

〈F 12 t 2 〉 = σ 12 2 〈v 2 〉 ℓ c 2 B(t), (4.25) 〈F 23 t 2 〉 = σ 23 2 τ 2 4λ exp 2λt 1 -exp -4λt , (4.26) 〈F 13 t 2 〉 = v 2 σ 12 2 σ 23 2 τ 2 4λℓ c 2 A t + v 2 σ 13 2 ℓ c B t , (4.27)
where we define

A t = ∑ i = 0 n t -1 τ i 4 exp 2λi τ 1 -exp -4λi τ , (4.28) B t = ∑ i = 0 n t -1 τ i 4 exp 2λi τ . (4.29)
In order to further develop these expressions, we consider heavy-tailed velocity distributions that behave as

p s v ∝ v β -1 , (4.30)
for v < v 0 with v 0 a characteristic velocity and β > 1. Such power-law type behaviour of the velocity PDF for small velocities is characteristic of observed transport behaviour in heterogeneous porous media [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF]. From (4.30), the transition time PDF behaves as ψ(t) ∝ t -1-β for large t. Note that the moments 〈τ k 〉 of ψ(t) are finite for k < ⌊β⌋, where ⌊•⌋ denotes the floor function. For illustration, here we consider the case β > 2, which implies 〈τ 2 〉 < ∞.

Shear-dominated regime (τ u ≪ t ≪ τ ϵ ): We now determine the stretching behavior in the shear-dominated time regime τ u ≪ t ≪ τ ϵ . In this time regime, we approximate (4.28) as

A t = 4λ ∑ i = 0 n t -1 iτ i 4 , B t = ∑ i = 0 n t -1 τ i 4 . (4.31)
These expressions may be estimated as follows. We note that transitions that contribute to the elongation at time t must have durations shorter than t, i.e. τ n < t because

t n t = ∑ i = 0 n t -1 τ i < t .
Thus, we can approximate

A t ≈ 4λ ∫ 0 t dt′t′ 4 ψ t′ ∑ i = 0 n t -1 i ∝ t 6 -β , (4.32)
where we approximated 〈n t 〉 ≈ t/〈τ〉. Using the same approximations for B(t) in (4.31) gives B(t) ∝ t 5-β . These scalings are consistent with the ones known for coupled continuous time random walks [START_REF] Marcole | Scaling forms of particle densities for Lévy walks and strong anomalous diffusion[END_REF]. Hence in the regime τ u ≪ t ≪ τ ϵ the transverse (3.45) deformation evolves linearly with time, 〈Λ t (t) 2 〉 ∝ t because F 22 (t) and F 33 (t) are approximately constant whilst 〈F 23 (t) 2 〉 ∝ t according to (4.26). The longitudinal (3.44) and material (3.50) deformations evolve nonlinearly to leading order as 〈Λ l (t) 2 〉 ∝ 〈ℓ(t) 2 〉t 6-β according to (4.32).

Crossover (t ≳ τ ϵ ): We now focus on the crossover from the power-law to exponential stretching regimes. For t ≳ τ ϵ , we approximate A(t) ≈ B(t) in (4.28) because the second term is exponentially small, and B(t) is given by (4.29). Along the same lines as above, we derive for B(t)

B t ≈ ∫ 0 t dt′t′ 4 ψ t′ ∑ i = 0 n t -1 exp 2λi τ ∝ t 4 -β exp 2λt -1 . (4.33)
To test these expressions developed throughout this Section, we simulate the non-stationary coupled CTRW (4.24) for 10 5 fluid particles over 500 spatial steps over the range of Lyapunov exponents λ ∈ (10 -1 , 10 -2 , 10 -3 ). The velocity gradients ϵ 22 , ϵ 12 are modelled as Gaussian variables with mean λ and zero respectively, the mean travel time 〈τ〉 = 10 -1 , and

the variances σ 22 2 = σ 12 2 = 1.
Here the velocity PDF is given by the Nakagami distribution , where some deviation from the predicted exponential growth is observed for small β due to the heavy-tailed nature of ψ(τ).

p v (v; k; θ) = 2 Γ (k) k θ k v 2k -1 e - kv 2 θ , v > 0,

Fluid Deformation in a Model Random 3D Flow

To study fluid deformation in flow where the stretching regimes described in § §4.2.1-4. We consider 10 5 realisations of the VH flow for ψ = π/2, ψ = π/64, and for each realisation a particle trajectory is calculated over the time period t ∈ [0, 10 2 ] to precision 10 -14 via an implicit Gauss-Legendre method. The associated deformation tensor in the Eulerian frame is calculated via solution of (3.4) to the same precision via the discrete QR decomposition method (due to the large associated deformations, in excess of 10 100 ). The Protean velocity gradient tensor ϵ′(t) is then computed along these trajectories by solution of the attracting trajectory ℳ(t) (3.29) via an explicit Runge-Kutta method over fixed time-step (Δt = 10 -3 ) and the Protean deformation gradient tensor F′(t) is calculated by direct evaluation of the integrals (4.12). These results are compared with direct computation of the deformation tensor via (3.4) in the Cartesian frame, and the error between these methods is found to be of the same order of the integration precision.

As shown in Figure 5(a), the velocity v(s) for both values of ψ is well-described by the Nakagami distribution (4.34) which exhibits power-law scaling p v (v) ~ v -1+2k in the limit of small v. The long waiting times associated with this power-law scaling has significant implications for anomalous transport and stretching dynamics, as shall be shown. The autocorrelation structure of the velocity field along streamlines is shown in Figure 5 We observe similar behaviour for the velocity and velocity gradient distributions, crosscorrelation and autocorrelation functions for ψ = π/64, with the exception that the PDFs of ϵ ii ′ and ϵ i j ′ are no longer Gaussian, but rather exhibit exponentially decaying tails. As these distributions still posses finite moments, then via the central limit theorem the same approach can be employed as for ψ = π/2 to develop stochastic models of fluid deformation. As expected, fluid stretching is significantly weaker for the reduced helicity case ψ = π/64, where again ϵ ii ′ = 0, λ, -λ with λ ≈ 0.069. Whilst not shown, the autocorrelation structure of the velocity gradient components for ψ = π/2 and ψ = π/64 follows a similar first-order process to that of the fluid velocity, indicating that fluid deformation also follows a Markov process in space.

Along with spatial Markovianity of the Lagrangian velocity (and velocity gradient), the remarkably simple deformation structure (strongly decorrelated velocity gradient components) of the VH flow in the Protean frame facilitates explicit solution of the fluid deformation tensor as a CTRW using the methods outlined in §4. Although not shown here, we observe similar behaviour (spatial Markovianity, strongly decorrelated velocity gradient) for other steady random 3D flows, including the Kraichnan flow (Kraichnan 1970) and a steady random flow defined as v(x) ≡ ∇ g1 (x) × ∇ g2 (x), where g 1 (x), g 2 (x) are random steady fields.

The observation of largely uncorrelated velocity gradient components arises as the diagonal and off-diagonal components of the Protean velocity gradient characterize distinctly different deformations (fluid stretching and shear respectively). As the Protean frame separates transverse deformation into orthogonal contracting (ê 3 ) and expanding (ê 2 ) directions, the associated longitudinal ϵ 12 and ϵ 13 fluid shears are also independent.

Similarly, the transverse shear ϵ 23 , which characterises the helicity of the flow is also independent of these longitudinal shears. These observations suggest that the CTRW model of fluid deformation outlined in §3, §4 can be applied to a wide range of steady random 3D flows.

Given the difference in Lyapunov exponent for the strong (ψ = π/2) and weakly (ψ = π/64) VH flow (λ ≈ 0.679, λ ≈ 0.069 respectively), we expect a larger separation between the advective τ u and stretching τ ϵ = 1/ϵ timescales for the less chaotic VH flow, ψ = π/64.

Figure 7 shows the transition from ballistic (t ≪ τ u ) to shear-induced stretching (τ u ≪ t ≪ τ ϵ ) regimes, where the latter appears to persist beyond τ ϵ for the strongly helical (ψ = π/2) VH flow. To clearly observe the transition from power-law to exponential stretching as predicted by (4.33), we compare this prediction with numerical calculations of fluid deformation for both cases of the variable helicity flow in Figure 8 and the corresponding CTRW (4.24).

Evolution of Finite-Time Lyapunov Exponents

From §4 we may also approximate evolution of the finite-time Lyapunov exponent (FTLE) in chaotic flows as ν t ≡ 1 2t ln μ t , (5.5)

where µ(t) is the leading eigenvalue of the Cauchy-Green tensor C′(t) = F′(t) ⊤ • F′(t).

Under the approximations

F 11 ′ (t) ∼ 0, F 33 ′ (t) ∼ 0, F 22 ′ (t), F 23 ′ (t) ∼ exp(λt + σ 22 ξ (t)), F 12 ′ (t) , F 13 ′ (t) ∼ t 2 -β/2 exp(λt + σ 22 ξ (t)),
where ξ(t) is a Gaussian white noise with zero mean and unit variance, then to leading order

μ(t) ∼ exp(2λt + σ 22 ξ(t))(1 + t 4 -β ) . (5.6)
As such, the average FTLE then evolves as

v(t) ≈ λ + 2 - β 2 lnt t , (5.7)
and so as expected the FTLE converges to the infinite time Lyapunov exponent λ at rate given by the velocity PDF power-law exponent β. We may also compute the variance of ν(t)

as σ v 2 (t) = (v(t) -v(t) ) 2 = σ 22 2 t . (5.8)
As expected, both of these expressions provide accurate estimates of the FTLE mean and variance for all but short times, as reflected in Figure 9 and 10 respectively. Here the explicit role of algebraic stretching at times t < τ ϵ is shown to govern convergence of the ensemble average of the FTLE toward its infinite-time value λ.

Conclusions

We Under an appropriate streamline coordinate system both the velocity gradient and deformation gradient tensor are upper-triangular, where the diagonal components correspond to fluid stretching, and the off-diagonal terms correspond to fluid shear. This uppertriangular form of the velocity gradient tensor greatly simplifies solution of the deformation gradient tensor ODE, and facilitates derivation of an ab initio stochastic model of fluid deformation in random steady 3D flows. This approach differs significantly from conventional, empirical stochastic models of fluid deformation.

We apply this framework to several model flows and show that the statistics of the velocity gradient components are remarkably simple, with many components uncorrelated and Gaussian distributed. Similar to previous studies (Dentz et al. 2016b;Le Borgne et al. 2008a,b) we find that the fluid velocity and velocity gradient components follow a spatial Markov process, which forms the basis for a stochastic model for fluid deformation as a coupled CTRW along streamlines.

We develop analytic solutions to these CTRWs over various temporal regimes in both chaotic and non-chaotic flows, and these solutions show excellent agreement with direct numerical simulations. These solutions provide a concrete means of predicting transverse Λ t (t) and longitudinal Λ l (t) fluid deformation in random, steady flows from point-wise velocity statistics. We show that fluid deformation over all temporal regimes only depends upon three statistical parameters: the Lyapunov exponent λ, velocity PDF scaling β : lim v→0 p v (v) ~ v -1+β , and shear exponent α : ϵ 12 , ϵ 13 ~ v 5/2-α . For non-chaotic steady flows (such as isotropic Darcy flow), only β and α need be quantified. These parameters can be determined directly from high-resolution experimental (e.g. PIV Holzner et al. (2015)) or computational (e.g. CFD de Carvalho et al. ( 2017)) studies, or some (such as the velocity PDF scaling) may be estimated from e.g. breakthrough curve data [START_REF] Tyukhova | Mechanisms of anomalous dispersion in flow through heterogeneous porous media[END_REF]).

These statistical parameters completely quantify the growth rates of transverse Λ t (t) and longitudinal Λ l (t) fluid deformation. This quantitative link between the velocity field statistics and fluid deformation is especially useful for application porous media flows, where recent advances (i) quantitatively link fluid mixing to fluid deformation (de [START_REF] De Barros Felipe | Flow topology and scalar mixing in spatially heterogeneous flow fields[END_REF][START_REF] Borgne | The lamellar description of mixing in porous media[END_REF][START_REF] Lester | Chaotic mixing in three-dimensional porous media[END_REF][START_REF] Dentz | Evolution of a solute blob in heterogeneous porous media[END_REF], and (ii) quantitatively link medium properties (such as heterogeneity controls in Darcy scale flow) to velocity field statistics (Dentz et al. 2016a;[START_REF] Fiori | Ergodic transport trough aquifers of non-gaussian log conductivity distribution and occurence of anomalous behavior[END_REF][START_REF] Edery | Origins of anomalous transport in heterogeneous media: Structural and dynamic controls[END_REF][START_REF] Tyukhova | Mechanisms of anomalous dispersion in flow through heterogeneous porous media[END_REF]).

Solution of the governing CTRW uncovers a rich array of fluid deformation behaviour across the various regimes, ranging from linear and power-law to exponential and superexponential growth. In most cases there are marked differences between longitudinal Λ l (t) and transverse Λ t (t) deformation, reflecting the strong deformation anisotropy that results from the constraints inherent to steady 3D flows. For chaotic flows, Λ l (t) and Λ t (t) respectively grow quadratically and linearly with time in the ballistic regime t ≪ τ u , and both of these grow exponentially in the exponential regime t ≫ τ ϵ . In the shear-induced stretching regime τ u ≪ t ≪ τ ϵ , transverse deformation grows exponentially, whereas longitudinal deformation transitions from power-law to super-exponential growth due to coupling between shear and stretching. For non-chaotic, steady zero-helicity flows (such as isotropic Darcy flow), transverse and longitudinal fluid deformation both grows algebraically as power-laws via similar mechanisms to that of 2D steady flows (Dentz et al. 2016b), namely the intermittency of shear events.

These results uncover and quantify the evolution of transverse and longitudinal fluid deformation in steady random 3D flows, as illustrated in Figure 1. The ad innate stochastic model of fluid deformation identifies and quantifies the statistical parameters of the flow field which govern fluid deformation. The deformation measures Λ t (t), Λ l (t) provide critical inputs to models of fluid mixing and transport in heterogeneous steady flows. This approach provides a means to develop upscaled models of fluid mixing and transport in e.g. porous media flows which may be couched in terms of physical properties of the host medium.

  To address this shortcoming, in this study we derive an ab initio stochastic model of fluid deformation in steady, random 3D flows that fully resolve the longitudinal Λ l and transverse Λ t deformations. Rather than develop empirical models of fluid deformation in specific flows, we analyse the kinematics of deformation in steady 3D flows and derive a generic framework applicable to all such flows. This stochastic model takes the form of a coupled continuous time random walk (CTRW) along fluid streamlines, and represents an extension of previous work by Dentz et al. (2016b) on stochastic modelling of fluid deformation in steady 2D flows.We apply this model to several model flows and show it accurately resolves the deformation dynamics and provides a valuable link between the Eulerian flow properties and Lagrangian deformation structure. As this CTRW model may be readily applied to existent numerical or experimental datasets, it provides a means to fully characterise and predict Lagrangian fluid deformation in complex steady flows from the Eulerian velocity and velocity gradient statistics. In the context of porous media flow (e.g. pore-and Darcy-scale flow simulations or experiments) this link provides an important building block for the development of fluid transport, mixing and reaction models which may be couched directly in terms of medium properties, such as heterogeneity controls in Darcy scale flow (e.g.de Barros et al. (2012); Le[START_REF] Borgne | The lamellar description of mixing in porous media[END_REF]).

  -e 3 ′ (t) ⋅ e ˙1 ′ (t) -e 3 ′ (t) ⋅ e ˙2 ′ (t) 0 . (3.16) Defining the Protean velocity gradient in the absence of this contribution as ϵ t ,

α0 ∫ 0 τ

 0 Figure 2(b). This approach is particularly useful when Lagrangian data is only available over short times and explicitly identifies the inertial initial orientation angle α 0,∞ , allowing the Protean transform to be determined uniquely from t = 0.

  ) in a steady 2D Darcy flow. Conversely, Lester et al. (2016) use the growth rate of l (3) (t) to predict mixing of a continuously injected source in steady 3D pore-scale flow. These different deformation rates indicate that a single scalar cannot be used to characterise fluid deformation in steady random 3D flow. Instead it is necessary to characterise both longitudinal and transverse fluid deformation in such flows.

  2, 3, (3.57) and so the only persistent growth of fluid deformation in steady zero helicity flows only occurs via the longitudinal and transverse shear deformations F 12 , F 13 , F 23 . Secondly, as the streamlines of this steady flow are confined to 2D integral surfaces, exponential fluid stretching (such as occurs in chaotic advection) is not possible due to constraints associated with the Poincaré-Bendixson theorem. Hence the shear deformations may only grow algebraically in time, i.e. . (2016b) develop a CTRW model for shear deformation in steady random incompressible 2D flows and uncover algebraic stretching of F 12 ~ t r (again in Protean coordinates) and show that that the index r depends upon the coupling between shear deformation and velocity fluctuations along streamlines, ranging from diffusive r = 1/2 to superdiffusive r = 2 stretching. Hence fluid deformation is constrained to be algebraic in steady 3D zero helicity density flow, and evolves in a similar fashion to steady 2D flow Λ t t ∼ t r 23 , (3.61) Λ l t ∼ t max r 12 , r 13 , (3.62)

  4.4)Its distribution ψ(t) is obtained from (4.3) in terms of the Eulerian PDF p e (v) asψ(t) = ℓ c p e ( ℓ c /t) v e t 3. (4.5) Equations (4.2) constitute a continuous time random walk (CTRW) in that the time increment τ n ≡ ℓ c /v n at the nth random walk step is a continuous random variable unlike for discrete time random walks, which describe Markov processes in time. The CTRW is also called a semi-Markovian framework because the governing equations (4.2) are Markovian in step number, while any Markovian process A n when projected on time as A(t) = A n t is non-Markovian.

F

  23 (s) = F 22 (s) ∫ 0 s ds′ ϵ 23 (s′)F 33 (s′) v s (s′)F 22 (s′) , (4.12d) F 13 (s) = v s (s) ∫ 0 s ds′ ϵ 12 (s′)F 23 (s′) + ϵ 13 (s′)F 33 (s′) v s (s′) 2 , (4.12e)

  (4.34) with the parameters k = 1, θ = 10. The PDFs of the velocities in the numerical test cases used in §5 are very well characterised by this distribution. As such the PDF of τ = ℓ c /v is then ψ τ = ℓ c τ 2 p v ℓ c /τ , (4.35) hence ψ(t) ∝ t -1-β with β = 2k for t > τ u . The transition from algebraic to exponential growth is clearly shown in Figure 3(b)

1 N

 1 2.3 overlap, we consider a steady incompressible random flow with variable exponential stretching, the variable helicity (VH) flow. This flow is similar to that introduced by Dean et al.(2001) which involves the helicity parameter ψ ∈ [0, π/2] which controls the ensemble averaged helicity 〈v • ω〉 of the random flow. The velocity field of the VH flow is given byv x; ψ = ∑ n = χ n ψ × k n cos k ψ ⋅ x+ Ω n + χ n ψ × k n sin k ψ ⋅ x+ Ω n , (5.1)where the phase angle Ω is a random vector with independent, uniformly distributed components in [0, 2π], and χ(ψ), k(ψ) are the random vectors χ ψ = sin θ sin ψ sin ϕ, sin θ sin ψ cos ϕ, cos θ sin ψ , (5.2) k ψ = sin θ sin ψ sin ϕ, sin θ sin ψ cos ϕ, cos θ sin ψ , (5.3) where θ, ϕ are the independent, uniformly distributed variables θ ∈ [-π/2, π/2], ϕ ∈ [-π, π]. This flow is divergence free, and converges to a multi-Gaussian field with sufficiently large N; we use N=64 throughout. For ψ = 0, the VH flow is 2D and constrained to the x 1 , x 2 plane, and for ψ = π/2 the flow is a fully 3D flow. Hence fluid stretching increases monotonically from non-chaotic λ = 0 to globally chaotic λ ≈ ln 2 over the range ψ ∈ [0, π/2]. We consider fluid deformation for the two cases ψ = π/64, ψ = π/2, which corresponds to weak and strong exponential stretching respectively. Representative streamlines and velocity contours for both of these cases are shown in Figure 4, which clearly illustrates the quasi-2D and fully 3D flow structures for ψ = π/64, ψ = π/2 respectively.

  (b), indicating exponential decay which may be modelled as a first-order process for the correlation functionR s = v s′ v s′ -s ≈ exp -s/l c , (5.4)where the correlation length scale ℓ c is estimated as l c = ∫ 0 ∞ dsR s . This observation is consistent with the hypothesis used in (4.2) that steady random flows exhibit spatial Markovianity along streamlines (LeBorgne et al. 2008a,b).The distributions of the components of the velocity gradient tensor ϵ′ for both ψ = π/2 and ψ = π/64 are shown in Figure6, where for ψ = π/2 all of the components are Gaussian distributed. The diagonal components have approximately the same varianceσ i j 2 ≡ 〈( ϵ i j ′ -〈 ϵ i j ′ 〉) 2 〉and mean 〈 ϵ ii ′ = 0, λ, -λ for i = 1: 3, with the Lyapunov exponent λ ≈ 0.679, indicating exponential fluid stretching close to the theoretical limit ln 2 for autonomous and continuous 3 d.o.f. systems. All the components of ϵ′(t) are essentially uncorrelated (where the correlation r ϵ i j ′ , ϵ kl ′ < 10 -4 , except for the diagonal components which are negatively correlated due to the constraint of incompressibility:

  have studied the kinematics of fluid deformation in steady 3D random flows such as pore-and Darcy-scale flows in heterogeneous porous media. We have developed an ab initio stochastic model which predicts evolution of the fluid deformation gradient tensor F(t) as coupled continuous time random walk (CTRW) from the velocity and velocity gradient statistics. One important way in which steady flow differs from its unsteady counterpart is that steady flow imposes constraints upon the kinematics of fluid deformation, namely zero mean extensional deformation in the stream-wise direction. This constraint then renders fluid deformation in steady flows anisotropic with respect to the velocity direction, and so fluid deformation in steady 3D flow cannot be characterised in terms of a single scalar quantity. This anisotropy may be characterised in terms of the fluid deformation longitudinal Λ l (t) and transverse Λ t (t) to the mean flow direction, which respectively act as inputs to models of longitudinal (Le[START_REF] Borgne | The lamellar description of mixing in porous media[END_REF] and transverse[START_REF] Lester | Chaotic mixing in three-dimensional porous media[END_REF]) mixing in steady random 3D flows. As such, fluid deformation in such flows cannot be characterised in terms of a single scalar value.

Figure 1 .

 1 Figure 1. Evolution of a 2D material surface ℋ 2D (light grey sheet, blue online) with the mean flow direction 〈v〉 (black arrow) arising from a continuously injected 1D line source ℋ 1D 1 (black line). A 1D cross section ℋ 1D 3 (dark grey line, blue online) of this surface is only subject to transverse deformation, whereas a pulsed injection along ℋ 1D 1 subsequently forms the 1D line ℋ 1D (2) (light grey line, red online) which is subject to both transverse and longitudinal fluid deformation.

Figure 2 .

 2 Figure 2.(a) Convergence of the transverse orientation angle α(t) for different initial conditions α(0) = α 0 (thin black lines) toward the attracting trajectory ℳ(t) = α(t; α 0,∞ ) (thick gray line) for the ABC flow. Note the correspondence between ℳ(t) and the transverse angle α(t) associated with minimum (gray dashed line) divergence ∂g/∂α(α, t) responsible for creation of the attracting trajectory. The maximum divergence is also shown (black dashed line). (b) Convergence of the initial angle associated with the approximate attracting trajectory α 0,τ toward the infinite time limit α 0,∞ .

Figure 3 .

 3 Figure 3.Comparison between analytic predictions of B(t) (4.29) for t > τ u (gray dashed) with numerical simulations from the synthetic coupled CTRW (4.24) (solid dark grey, solid blue online) for β ∈ (1.5, 2.5, 3.5) (plots are offset for clarity with β increasing from bottom to top) and (a) λ = 10 -3 , σ 22 = 10 -3 , (b) λ = 10 -2 , σ 22 = 10 -2 , (c) λ = 10 -1 , σ 22 = 10 -1 .

Figure 4 .

 4 Figure 4. Typical particle trajectory in a realisation of the variable helicity flow for (a) ψ = π/2, (c) ψ = π/64, and contour plot of velocity magnitude distribution in x 3 = 0 plane for (b) ψ = π/2, (d) ψ = π/64.

Figure 5 .

 5 Figure 5.(a) Velocity probability distribution function (solid) with fitted Nakagami PDF (dashed), and (b) autocorrelation structure of fluid velocity v(s) for the variable helical flow with ψ = π/2 (light) and ψ = π/64 (dark).

Figure 6 .

 6 Figure 6. Distribution of the diagonal reoriented rate of strain components ϵ ii ′ t ϵ 11 ′ (dark grey), ϵ 22 ′ (medium grey), ϵ 33 ′ (light grey)) for the variable helicity flow with (a) ψ = π/2 and (b) ψ = π/64. Distribution of the off-diagonal reoriented rate of strain components ϵ i j ′ t ϵ 12 ′ (dark grey), ϵ 13 ′ (medium grey), ϵ 23 ′ (light grey)) for the variable helicity flow with (c) ψ = π/2 (d) ψ = π/64.

Figure 7 .

 7 Figure 7.Growth of the mean squares A(t) (black), B(t) (dark grey, blue online), 〈x 23 (t) 2 〉 (gray) with time t in the ballistic (t ≪ τ u ) and shear-induced stretching (τ u ≪ t ≪ τ ϵ ) regimes for the VH flow with (a) ψ = π/2, (b) ψ = π/64. The analytic scalings 〈x ij (t) 2 〉 ~ t 2 from (4.16b), and A(t) ~ t 6-β from (4.32) for these regimes are shown by the dashed grey lines.

Figure 8 .

 8 Figure 8. Comparison between analytic prediction of B(t) (4.33) (gray dashed line), numerical simulations from the synthetic CTRW (4.24) (black solid line) and deformation in the variable helicity flow (dark grey, blue online) with (a) ψ = π/2, (b) ψ = π/64. The dotted line in (b) also depicts the algebraic growth of 〈F 12 (t) 2 〉 ~ t 2 from (4.16b) in the ballistic regime (t ≪ τ u ).

Figure 9 .

 9 Figure 9.Comparison between analytic prediction of the mean FTLE 〈ν(t)〉 (5.7) (gray, dashed), the infinite-time Laypunov exponent λ (black, dashed) and 〈ν(t)〉 computed directly (dark grey, blue online) for the variable helicity flow with (a) ψ = π/2, (b) ψ = π/64.

Figure 10 .

 10 Figure 10. Comparison between analytic prediction of the FTLE variance σ v 2 (t) (5.8) (gray, dashed) and σ v 2 (t) computed directly (dark grey, blue online) for the variable helicity flow with (a) ψ = π/2, (b) ψ = π/64.

J Fluid Mech. Author manuscript; available in PMC 2019 May 01.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.