3.1 General development

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

(3.1a,b ) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{x}(t;\boldsymbol{X})}{\text{d}t}=\boldsymbol{v}(t),\quad \boldsymbol{v}(t;\boldsymbol{X})\equiv \boldsymbol{u}[\boldsymbol{x}(t;\boldsymbol{X})],\end{eqnarray}$$

where we denote the reference material vectors in the Eulerian and Lagrangian frames by $\boldsymbol{x}(t)$ and $\boldsymbol{X}$ , respectively. The deformation gradient tensor $\unicode[STIX]{x1D641}(t)$ quantifies how the infinitesimal vector $\text{d}\boldsymbol{x}(t;\boldsymbol{X})$ deforms from its reference state $\text{d}\boldsymbol{x}(t=0;\boldsymbol{X})=\text{d}\boldsymbol{X}$ as

(3.2) $$\begin{eqnarray}\text{d}\boldsymbol{x}=\unicode[STIX]{x1D641}(t)\boldsymbol{\cdot }\text{d}\boldsymbol{X}\end{eqnarray}$$

and equivalently

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D60D}_{ij}(t)\equiv \frac{\unicode[STIX]{x2202}x_{i}(t;\boldsymbol{X})}{\unicode[STIX]{x2202}X_{j}}.\end{eqnarray}$$

Following this definition, the deformation gradient tensor $\unicode[STIX]{x1D641}(t)$ evolves with travel time $t$ along a Lagrangian trajectory (streamline) as

(3.4) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D641}(t)}{\text{d}t}=\unicode[STIX]{x1D750}(t)\unicode[STIX]{x1D641}(t),\quad \unicode[STIX]{x1D641}(0)=\mathbf{1},\end{eqnarray}$$

where the velocity gradient tensor $\unicode[STIX]{x1D750}(t)=\unicode[STIX]{x1D735}\boldsymbol{v}(t)^{\top }\equiv \unicode[STIX]{x1D735}\boldsymbol{u}[\boldsymbol{x}(t;\boldsymbol{X})]^{\top }$ and the superscript $\top$ 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 $\boldsymbol{x}=(x_{1},x_{2},x_{3})^{\top }$ , and the reoriented Protean coordinate system by $\boldsymbol{x}^{\prime }=(x_{1}^{\prime },x_{2}^{\prime },x_{3}^{\prime })^{\top }$ , where the velocity vector $\boldsymbol{v}(\boldsymbol{x})=(v_{1},v_{2},v_{3})^{\top }$ in the Cartesian frame transforms to $\boldsymbol{v}^{\prime }(\boldsymbol{x}^{\prime })=(v,0,0)^{\top }$ in the Protean frame, with $v=|\boldsymbol{v}|$ . These two frames are related by the transform (Truesdell & Noll Reference Truesdell and Noll1992)

(3.5) $$\begin{eqnarray}\boldsymbol{x}^{\prime }=\boldsymbol{x}_{0}(t)+\unicode[STIX]{x1D64C}^{\top }(t)\boldsymbol{x},\end{eqnarray}$$

where $\boldsymbol{x}_{0}(t)$ is an arbitrary translation vector, and $\unicode[STIX]{x1D64C}(t)$ is a rotation matrix (or proper orthogonal transformation): $\unicode[STIX]{x1D64C}^{\top }(t)\boldsymbol{\cdot }\unicode[STIX]{x1D64C}(t)=\mathbf{1}$ , $\det [\unicode[STIX]{x1D64C}(t)]=1$ . From (3.5) the differential elements $\text{d}\boldsymbol{x}$ , $\text{d}\boldsymbol{X}$ transform respectively as $\text{d}\boldsymbol{x}^{\prime }=\unicode[STIX]{x1D64C}^{\top }(t)\,\text{d}\boldsymbol{x}$ , $\text{d}\boldsymbol{X}^{\prime }=\unicode[STIX]{x1D64C}^{\top }(0)\,\text{d}\boldsymbol{X}$ , and so from (3.2) the deformation gradient tensor then transforms as

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D641}^{\prime }(t)=\unicode[STIX]{x1D64C}^{\top }(t)\unicode[STIX]{x1D641}(t)\unicode[STIX]{x1D64C}(0).\end{eqnarray}$$

Formally, $\unicode[STIX]{x1D641}(t)$ is a pseudo-tensor (Truesdell & Noll Reference Truesdell and Noll1992), as it is not objective as per (3.6), $\unicode[STIX]{x1D641}^{\prime }(t)\neq \unicode[STIX]{x1D64C}^{\top }(t)\unicode[STIX]{x1D641}(t)\unicode[STIX]{x1D64C}(t)$ (Ottino Reference Ottino1989). Differentiating (3.6) with respect to the Lagrangian travel time $t$ yields

(3.7) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D641}^{\prime }}{\text{d}t}=\dot{\unicode[STIX]{x1D64C}}^{\top }(t)\unicode[STIX]{x1D641}(t)\unicode[STIX]{x1D64C}(0)+\unicode[STIX]{x1D64C}^{\top }(t)\unicode[STIX]{x1D750}(t)\unicode[STIX]{x1D641}(t)\unicode[STIX]{x1D64C}(0)=\unicode[STIX]{x1D750}^{\prime }(t)\unicode[STIX]{x1D641}^{\prime }(t),\end{eqnarray}$$

where the transformed rate of strain tensor $\unicode[STIX]{x1D750}^{\prime }(t)$ is defined as

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D750}^{\prime }(t)\equiv \unicode[STIX]{x1D64C}^{\top }(t)\unicode[STIX]{x1D750}(t)\unicode[STIX]{x1D64C}(t)+\unicode[STIX]{x1D63C}(t),\end{eqnarray}$$

and the contribution due to a moving coordinate frame is $\unicode[STIX]{x1D63C}(t)\equiv \dot{\unicode[STIX]{x1D64C}}^{\top }(t)\unicode[STIX]{x1D64C}(t)$ . As $\unicode[STIX]{x1D64C}(t)$ is orthogonal, then $\dot{\unicode[STIX]{x1D64C}}^{\top }(t)\unicode[STIX]{x1D64C}(t)+\unicode[STIX]{x1D64C}^{\top }(t)\dot{\unicode[STIX]{x1D64C}}(t)=0$ and so $\unicode[STIX]{x1D63C}(t)$ is skew-symmetric. Note that the velocity vector $\boldsymbol{v}(t)=\text{d}\boldsymbol{x}/\text{d}t$ also transforms as $\boldsymbol{v}^{\prime }(t)=\unicode[STIX]{x1D64C}^{\top }(t)\boldsymbol{v}(t)$ . The basic idea of the Protean coordinate frame is to find an appropriate local reorientation $\unicode[STIX]{x1D64C}(t)$ that renders the transformed velocity gradient $\unicode[STIX]{x1D750}^{\prime }(t)$ upper triangular, simplifying solution of the deformation tensor $\unicode[STIX]{x1D641}$ in (3.7).

3.2 Fluid deformation in 3-D Protean coordinates

Adachi (Reference Adachi1983, Reference Adachi1986) shows that for steady 2-D flow, reorientation into streamline coordinates automatically yields an upper triangular velocity gradient tensor, hence the reorientation matrix $\unicode[STIX]{x1D64C}(t)$ is simply given in terms of the velocity $\boldsymbol{v}=[v_{1},v_{2}]^{\top }$ in the fixed Cartesian frame as

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D64C}(t)=\frac{1}{v}\left(\begin{array}{@{}cc@{}}v_{1} & -v_{2}\\ v_{2} & v_{1}\end{array}\right),\end{eqnarray}$$

hence $\boldsymbol{v}^{\prime }=[v,0]^{\top }$ , where $v=\sqrt{v_{1}^{2}+v_{2}^{2}}$ .

In contrast to steady 2-D flow, the additional spatial dimension of steady 3-D flows admits an additional degree of freedom for the specification of streamline coordinates. In 3-D streamline coordinates, the base vector $\boldsymbol{e}_{1}^{\prime }(t)$ aligns with the velocity vector $\boldsymbol{v}(t)$ , such that $\boldsymbol{v}^{\prime }(t)=[v(t),0,0]^{\top }$ , but the transverse vectors $\boldsymbol{e}_{2}^{\prime }(t)$ and $\boldsymbol{e}_{3}^{\prime }(t)$ are arbitrary up to a rotation about $\boldsymbol{e}_{1}^{\prime }(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 $\unicode[STIX]{x1D750}^{\prime }(t)$ such that it is upper triangular, yielding explicit closed-form solution for the deformation gradient tensor $\unicode[STIX]{x1D641}^{\prime }(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 Adachi (Reference Adachi1983, Reference Adachi1986), this condition is satisfied for all steady 2-D flows, but this is still an open question for steady 3-D flows.

For reasons that will become apparent, we develop the 3-D Protean coordinate frame by decomposing the reorientation matrix $\unicode[STIX]{x1D64C}(t)$ into two sequential reorientations $\unicode[STIX]{x1D64C}_{1}(t)$ , $\unicode[STIX]{x1D64C}_{2}(t)$ as

(3.10) $$\begin{eqnarray}\unicode[STIX]{x1D64C}(t)\equiv \unicode[STIX]{x1D64C}_{2}(t)\unicode[STIX]{x1D64C}_{1}(t),\end{eqnarray}$$

where $\unicode[STIX]{x1D64C}_{1}(t)$ is a reorientation that aligns the velocity vector $\boldsymbol{v}$ with the Protean base vector $\boldsymbol{e}_{1}^{\prime }(t)$ , and $\unicode[STIX]{x1D64C}_{2}(t)$ is a subsequent reorientation of the Protean frame about $\boldsymbol{e}_{1}^{\prime }(t)$ . As such, the reorientation $\unicode[STIX]{x1D64C}_{1}(t)$ renders the Protean frame a streamline coordinate frame, and $\unicode[STIX]{x1D64C}_{2}$ represents a (currently) arbitrary reorientation of this streamline frame about the velocity vector.

We begin by considering the 3-D rotation matrix $\unicode[STIX]{x1D64C}_{1}(t)$ which is defined in terms of the unit rotation axis $\boldsymbol{q}(t)$ (which is orthogonal to both the velocity and $\boldsymbol{e}_{1}$ vectors), and the associated reorientation angle $\unicode[STIX]{x1D703}(t)$ , both of which are given in terms of the local velocity vector $\boldsymbol{v}(t)=[v_{1},v_{2},v_{3}]^{\top }$ as

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}(t)=\frac{\boldsymbol{e}_{1}\times \boldsymbol{v}}{\Vert \boldsymbol{e}_{1}\times \boldsymbol{v}\Vert }=\frac{1}{\sqrt{v_{2}^{2}+v_{3}^{2}}}\{0,v_{3},-v_{2}\}, & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle \cos \unicode[STIX]{x1D703}(t)=\frac{\boldsymbol{e}_{1}\boldsymbol{\cdot }\boldsymbol{v}}{\Vert \boldsymbol{e}_{1}\boldsymbol{\cdot }\boldsymbol{v}\Vert }=\frac{v_{1}}{v}. & \displaystyle\end{eqnarray}$$

The rotation tensor $\unicode[STIX]{x1D64C}_{1}(t)$ then is given by

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D64C}_{1}(t)=\cos \unicode[STIX]{x1D703}(t)\unicode[STIX]{x1D644}+\sin \unicode[STIX]{x1D703}(t)(\boldsymbol{q})_{\times }^{\top }+[1-\cos \unicode[STIX]{x1D703}(t)]\boldsymbol{q}(t)\otimes \boldsymbol{q}(t),\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the identity matrix, $(\cdot )_{\times }$ denotes the cross-product matrix and $\Vert \cdot \Vert$ is the $\ell ^{2}$ -norm. The second reorientation matrix $\unicode[STIX]{x1D64C}_{2}(t)$ corresponds to an arbitrary reorientation (of angle $\unicode[STIX]{x1D6FC}(t)$ ) about the $\boldsymbol{e}_{1}^{\prime }$ axis, which is explicitly

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D64C}_{2}(t)=\left(\begin{array}{@{}ccc@{}}1 & 0 & 0\\ 0 & \cos \unicode[STIX]{x1D6FC}(t) & -\text{sin}\,\unicode[STIX]{x1D6FC}(t)\\ 0 & \sin \unicode[STIX]{x1D6FC}(t) & \cos \unicode[STIX]{x1D6FC}(t)\end{array}\right).\end{eqnarray}$$

Whilst the algebraic expressions (in terms of $\unicode[STIX]{x1D6FC}$ , $v_{1}$ , $v_{2}$ , $v_{3}$ ) for the components of the composition rotation matrix $\unicode[STIX]{x1D64C}(t)=\unicode[STIX]{x1D64C}_{2}(t)\unicode[STIX]{x1D64C}_{1}(t)$ are somewhat complicated, it is useful to note that $\unicode[STIX]{x1D64C}(t)$ may be expressed in terms of the Protean frame basis vectors as

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D64C}(t)=[\boldsymbol{e}_{1}^{\prime }(t),\boldsymbol{e}_{2}^{\prime }(t),\boldsymbol{e}_{3}^{\prime }(t)],\end{eqnarray}$$

hence the moving coordinate contribution $\unicode[STIX]{x1D63C}(t)$ is given by

(3.16) $$\begin{eqnarray}\unicode[STIX]{x1D63C}(t)=\dot{\unicode[STIX]{x1D64C}}^{\top }(t)\unicode[STIX]{x1D64C}(t)=\left(\begin{array}{@{}ccc@{}}0 & \boldsymbol{e}_{2}^{\prime }(t)\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{1}^{\prime }(t) & \boldsymbol{e}_{3}^{\prime }(t)\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{1}^{\prime }(t)\\ -\boldsymbol{e}_{2}^{\prime }(t)\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{1}^{\prime }(t) & 0 & \boldsymbol{e}_{3}^{\prime }(t)\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{2}^{\prime }(t)\\ -\boldsymbol{e}_{3}^{\prime }(t)\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{1}^{\prime }(t) & -\boldsymbol{e}_{3}^{\prime }(t)\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{2}^{\prime }(t) & 0\end{array}\right).\end{eqnarray}$$

Defining the Protean velocity gradient in the absence of this contribution as $\tilde{\unicode[STIX]{x1D750}}(t)$ ,

(3.17) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D750}}(t)\equiv \unicode[STIX]{x1D64C}^{\top }(t)\unicode[STIX]{x1D750}(t)\unicode[STIX]{x1D64C}(t)=\unicode[STIX]{x1D750}^{\prime }(t)-\unicode[STIX]{x1D63C}(t),\end{eqnarray}$$

we find the components of $\unicode[STIX]{x1D63C}(t)$ may be related to $\tilde{\unicode[STIX]{x1D750}}(t)$ as follows. From (3.17) and (3.15), $\tilde{\unicode[STIX]{x1D716}}_{ij}(t)=\boldsymbol{e}_{i}^{\prime }(t)\unicode[STIX]{x1D750}(t)\boldsymbol{e}_{j}^{\prime }(t)$ , and since $\dot{\boldsymbol{v}}(t)=\unicode[STIX]{x1D750}(t)\boldsymbol{\cdot }\boldsymbol{v}(t)$ , $\boldsymbol{e}_{1}^{\prime }(t)=\boldsymbol{v}(t)/v(t)$ , then

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{e}_{2}^{\prime }(t)\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{1}^{\prime }(t)=\boldsymbol{e}_{2}^{\prime }(t)\boldsymbol{\cdot }\unicode[STIX]{x1D750}(t)\boldsymbol{\cdot }\boldsymbol{e}_{1}^{\prime }(t)=\tilde{\unicode[STIX]{x1D716}}_{21}(t), & \displaystyle\end{eqnarray}$$

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{e}_{3}^{\prime }(t)\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{1}^{\prime }(t)=\boldsymbol{e}_{3}^{\prime }(t)\boldsymbol{\cdot }\unicode[STIX]{x1D716}(t)\boldsymbol{\cdot }\boldsymbol{e}_{1}^{\prime }(t)=\tilde{\unicode[STIX]{x1D716}}_{31}(t). & \displaystyle\end{eqnarray}$$

From these relationships the Protean velocity gradient tensor $\unicode[STIX]{x1D750}^{\prime }(t)=\tilde{\unicode[STIX]{x1D750}}+\unicode[STIX]{x1D63C}(t)$ is then

(3.20) $$\begin{eqnarray}\unicode[STIX]{x1D750}^{\prime }(t)=\left(\begin{array}{@{}ccc@{}}\tilde{\unicode[STIX]{x1D716}}_{11} & \tilde{\unicode[STIX]{x1D716}}_{12}+\tilde{\unicode[STIX]{x1D716}}_{21} & \tilde{\unicode[STIX]{x1D716}}_{13}+\tilde{\unicode[STIX]{x1D716}}_{31}\\ 0 & \tilde{\unicode[STIX]{x1D716}}_{22} & \tilde{\unicode[STIX]{x1D716}}_{23}+\unicode[STIX]{x1D608}_{23}(t)\\ 0 & \tilde{\unicode[STIX]{x1D716}}_{32}-\unicode[STIX]{x1D608}_{23}(t) & -\tilde{\unicode[STIX]{x1D716}}_{22}-\tilde{\unicode[STIX]{x1D716}}_{11}\end{array}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D608}_{23}(t)=\boldsymbol{e}_{3}^{\prime }(t)\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{2}^{\prime }(t)$ . As the reorientation $\unicode[STIX]{x1D64C}_{2}$ (3.14) only impacts the $\unicode[STIX]{x1D716}_{22}^{\prime }$ , $\unicode[STIX]{x1D716}_{23}^{\prime }$ , $\unicode[STIX]{x1D716}_{32}^{\prime }$ , $\unicode[STIX]{x1D716}_{33}^{\prime }$ components, reorientation (via $\unicode[STIX]{x1D64C}_{1}(t)$ ) into streamline coordinates automatically renders the $\unicode[STIX]{x1D716}_{21}^{\prime }$ , $\unicode[STIX]{x1D716}_{31}^{\prime }$ components of the velocity gradient tensor zero. It is this transform which renders the Protean velocity gradient upper triangular for 2-D steady flows.

Whilst $\unicode[STIX]{x1D716}_{32}^{\prime }=\tilde{\unicode[STIX]{x1D716}}_{32}-\unicode[STIX]{x1D608}_{23}$ is non-zero in general, $\unicode[STIX]{x1D750}^{\prime }(t)$ may be rendered upper triangular by appropriate manipulation of the arbitrary reorientation angle $\unicode[STIX]{x1D6FC}(t)$ such that $\unicode[STIX]{x1D608}_{23}=\tilde{\unicode[STIX]{x1D716}}_{32}$ . In a similar manner to (3.18), (3.19), we may express $\unicode[STIX]{x1D608}_{23}$ in terms of the components of $\unicode[STIX]{x1D750}(t)$ as

(3.21) $$\begin{eqnarray}\boldsymbol{e}_{3}^{\prime }(t)\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{2}^{\prime }(t)=v(t)\boldsymbol{e}_{3}^{\prime }(t)\left[\frac{\unicode[STIX]{x2202}\boldsymbol{e}_{2}^{\prime }(t)}{\unicode[STIX]{x2202}\boldsymbol{v}}\unicode[STIX]{x1D750}(t)\right]\boldsymbol{e}_{1}^{\prime }(t)+\boldsymbol{e}_{3}^{\prime }(t)\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{e}_{2}^{\prime }(t)}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}(t)}\frac{\text{d}\unicode[STIX]{x1D6FC}}{\text{d}t},\end{eqnarray}$$

where from (3.13), (3.14)

(3.22) $$\begin{eqnarray}\displaystyle & \displaystyle v\boldsymbol{e}_{3}^{\prime }\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{e}_{2}^{\prime }}{\unicode[STIX]{x2202}\boldsymbol{v}}=\frac{1}{v+v_{1}}\{0,v_{3},-v_{2}\}=\frac{v_{3}\cos \unicode[STIX]{x1D6FC}-v_{2}\sin \unicode[STIX]{x1D6FC}}{v+v_{1}}\boldsymbol{e}_{2}^{\prime }-\frac{v_{2}\cos \unicode[STIX]{x1D6FC}+v_{3}\sin \unicode[STIX]{x1D6FC}}{v+v_{1}}\boldsymbol{e}_{3}^{\prime },\qquad & \displaystyle\end{eqnarray}$$

(3.23) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{e}_{3}^{\prime }\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{e}_{2}^{\prime }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}=-1, & \displaystyle\end{eqnarray}$$

and so $\unicode[STIX]{x1D608}_{23}(t)$ is then

(3.24) $$\begin{eqnarray}\unicode[STIX]{x1D608}_{23}(t)=\boldsymbol{e}_{3}^{\prime }\boldsymbol{\cdot }\dot{\boldsymbol{e}}_{2}^{\prime }=\frac{v_{3}\cos \unicode[STIX]{x1D6FC}-v_{2}\sin \unicode[STIX]{x1D6FC}}{v+v_{1}}\tilde{\unicode[STIX]{x1D716}}_{21}-\frac{v_{2}\cos \unicode[STIX]{x1D6FC}+v_{3}\sin \unicode[STIX]{x1D6FC}}{v+v_{1}}\tilde{\unicode[STIX]{x1D716}}_{31}-\frac{\text{d}\unicode[STIX]{x1D6FC}}{\text{d}t}.\end{eqnarray}$$

Hence the equation $\unicode[STIX]{x1D608}_{23}(t)=\tilde{\unicode[STIX]{x1D716}}_{32}$ defines an ordinary differential equation (ODE) for the arbitrary orientation angle $\unicode[STIX]{x1D6FC}(t)$ which renders $\unicode[STIX]{x1D750}^{\prime }(t)$ upper triangular. Although the components of $\tilde{\unicode[STIX]{x1D750}}(t)$ vary with $\unicode[STIX]{x1D6FC}$ , these may be expressed in terms of the components of $\unicode[STIX]{x1D750}^{(1)}$ (which are independent of $\unicode[STIX]{x1D6FC}$ ), defined as

(3.25) $$\begin{eqnarray}\unicode[STIX]{x1D750}^{(1)}(t)\equiv \unicode[STIX]{x1D64C}_{1}^{\top }(t)\unicode[STIX]{x1D750}(t)\unicode[STIX]{x1D64C}_{1}(t),\end{eqnarray}$$

where the components of $\tilde{\unicode[STIX]{x1D750}}(t)$ and $\unicode[STIX]{x1D750}^{(1)}(t)$ are related as

(3.26) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D716}}_{21}=\unicode[STIX]{x1D716}_{21}^{(1)}\cos \unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D716}_{31}^{(1)}\sin \unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

(3.27) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D716}}_{31}=\unicode[STIX]{x1D716}_{31}^{(1)}\cos \unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D716}_{21}^{(1)}\sin \unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

(3.28) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D716}}_{32}=\unicode[STIX]{x1D716}_{32}^{(1)}\cos ^{2}\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D716}_{23}^{(1)}\sin ^{2}\unicode[STIX]{x1D6FC}+(\unicode[STIX]{x1D716}_{33}^{(1)}-\unicode[STIX]{x1D716}_{22}^{(1)})\cos \unicode[STIX]{x1D6FC}\sin \unicode[STIX]{x1D6FC}. & \displaystyle\end{eqnarray}$$

Hence the condition $\unicode[STIX]{x1D608}_{23}=\tilde{\unicode[STIX]{x1D716}}_{32}$ is satisfied by the ODE

(3.29) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D6FC}}{\text{d}t} & = & \displaystyle g(\unicode[STIX]{x1D6FC},t)\nonumber\\ \displaystyle & = & \displaystyle a(t)\cos ^{2}\unicode[STIX]{x1D6FC}+b(t)\sin ^{2}\unicode[STIX]{x1D6FC}+c(t)\cos \unicode[STIX]{x1D6FC}\sin \unicode[STIX]{x1D6FC},\end{eqnarray}$$

where

(3.30) $$\begin{eqnarray}\displaystyle & \displaystyle a(t)=-\unicode[STIX]{x1D716}_{32}^{(1)}-\frac{v_{2}}{v+v_{1}}\unicode[STIX]{x1D716}_{31}^{(1)}-\frac{v_{3}}{v+v_{1}}\unicode[STIX]{x1D716}_{21}^{(1)}, & \displaystyle\end{eqnarray}$$

(3.31) $$\begin{eqnarray}\displaystyle & \displaystyle b(t)=\unicode[STIX]{x1D716}_{23}^{(1)}+\frac{v_{2}}{v+v_{1}}\unicode[STIX]{x1D716}_{31}^{(1)}+\frac{v_{3}}{v+v_{1}}\unicode[STIX]{x1D716}_{21}^{(1)}, & \displaystyle\end{eqnarray}$$

(3.32) $$\begin{eqnarray}\displaystyle & \displaystyle c(t)=\unicode[STIX]{x1D716}_{22}^{(1)}-\unicode[STIX]{x1D716}_{33}^{(1)}+\frac{2v_{2}}{v+v_{1}}\unicode[STIX]{x1D716}_{21}^{(1)}-\frac{2v_{3}}{v+v_{1}}\unicode[STIX]{x1D716}_{31}^{(1)}. & \displaystyle\end{eqnarray}$$

3.2.1 Evolution of Protean orientation angle $\unicode[STIX]{x1D6FC}(t)$

Equation (3.29) describes a first-order ODE for the transverse orientation angle $\unicode[STIX]{x1D6FC}(t)$ along a streamline which renders $\unicode[STIX]{x1D750}^{\prime }(t)$ upper triangular. Here the temporal derivative for $\unicode[STIX]{x1D6FC}(t)$ is associated with both change in flow structure along a streamline and impact of a moving coordinate transform as encoded by $\unicode[STIX]{x1D63C}(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 $\boldsymbol{e}_{1}=\boldsymbol{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 $\unicode[STIX]{x1D6FC}(t)$ , whereas the continuous QR method involves solution of the three-degree-of-freedom (3-d.o.f.) ODE system (A 5), and requires unitary integrators to preserve orthogonality of ${\mathcal{Q}}$ . These methods differ with respect to the initial conditions ${\mathcal{Q}}(0)=\unicode[STIX]{x1D644}$ , $\unicode[STIX]{x1D64C}(0)=\unicode[STIX]{x1D64C}_{1}(0)\boldsymbol{\cdot }\unicode[STIX]{x1D64C}_{2}(0)$ , where $\unicode[STIX]{x1D64C}_{1}(0)$ is given explicitly by (3.13), and $\unicode[STIX]{x1D64C}_{2}(0)$ is dependent upon the initial condition $\unicode[STIX]{x1D6FC}(0)=\unicode[STIX]{x1D6FC}_{0}$ as per (3.14), (3.29).

Figure 2. (a) Convergence of the transverse orientation angle $\unicode[STIX]{x1D6FC}(t)$ for different initial conditions $\unicode[STIX]{x1D6FC}(0)=\unicode[STIX]{x1D6FC}_{0}$ (thin black lines) toward the attracting trajectory ${\mathcal{M}}(t)=\unicode[STIX]{x1D6FC}(t;\unicode[STIX]{x1D6FC}_{0,\infty })$ (thick grey line) for the ABC flow. Note the correspondence between ${\mathcal{M}}(t)$ and the transverse angle $\unicode[STIX]{x1D6FC}(t)$ associated with minimum (grey dashed line) divergence $\unicode[STIX]{x2202}g/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FC},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 $\unicode[STIX]{x1D6FC}_{0,\unicode[STIX]{x1D70F}}$ toward the infinite-time limit $\unicode[STIX]{x1D6FC}_{0,\infty }$ .

Whilst the ODE (3.29) may be satisfied for any arbitrary initial condition $\unicode[STIX]{x1D6FC}(0)=\unicode[STIX]{x1D6FC}_{0}$ , resulting in non-uniqueness of $\unicode[STIX]{x1D750}^{\prime }(t)$ , this non-autonomous ODE is locally dissipative as reflected by the divergence

(3.33) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}=[b(t)-a(t)]\sin 2\unicode[STIX]{x1D6FC}+c(t)\cos 2\unicode[STIX]{x1D6FC},\end{eqnarray}$$

which admits maxima and minima of magnitude $\pm c(t)\sqrt{1+[b(t)-a(t)/c(t)]^{2}}$ respectively at

(3.34) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}(t)=\frac{1}{2}\arctan \left[\frac{b(t)-a(t)}{c(t)}\right]+\frac{\unicode[STIX]{x03C0}}{2}\frac{\text{sgn}\,c(t)\mp 1}{2}.\end{eqnarray}$$

Hence for $c(t)\neq 0$ , the ODE (3.29) admits an attracting trajectory ${\mathcal{M}}(t)$ , which attracts solutions from all initial conditions $\unicode[STIX]{x1D6FC}_{0}$ as illustrated in figure 2(a). As this trajectory is also a solution of (3.29), it may be expressed as ${\mathcal{M}}(t)=\unicode[STIX]{x1D6FC}(t;\unicode[STIX]{x1D6FC}_{0,\infty })$ , where $\unicode[STIX]{x1D6FC}_{0,\infty }$ is the initial condition at $t=0$ which is already on ${\mathcal{M}}(t)$ . Whilst all solutions of (3.29) render $\unicode[STIX]{x1D750}^{\prime }(t)$ upper triangular, only solutions along the attracting trajectory ${\mathcal{M}}(t)$ represent asymptotic dynamics independent of the initial condition $\unicode[STIX]{x1D6FC}_{0}$ , hence we define the Protean frame as that which corresponds to ${\mathcal{M}}(t)=\unicode[STIX]{x1D6FC}(t;\unicode[STIX]{x1D6FC}_{0,\infty })$ .

As the attracting trajectory ${\mathcal{M}}(t)$ maximises dissipation $-\unicode[STIX]{x2202}g/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}$ over long times, the associated initial condition $\unicode[STIX]{x1D6FC}_{0,\infty }$ is quantified by the limit $\unicode[STIX]{x1D6FC}_{0,\infty }=\lim _{\unicode[STIX]{x1D70F}\rightarrow \infty }\unicode[STIX]{x1D6FC}_{0,\unicode[STIX]{x1D70F}}$ , where

(3.35) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{0,\unicode[STIX]{x1D70F}}(t)=\arg \min _{\unicode[STIX]{x1D6FC}_{0}}\int _{0}^{\unicode[STIX]{x1D70F}}\frac{\unicode[STIX]{x2202}g(\unicode[STIX]{x1D6FC}(t^{\prime };\unicode[STIX]{x1D6FC}_{0}),t^{\prime })}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\,\text{d}t^{\prime }.\end{eqnarray}$$

Whilst ${\mathcal{M}}(t)$ can be identified by evolving (3.29) until acceptable convergence is obtained, ${\mathcal{M}}(t)$ may be identified at shorter times $\unicode[STIX]{x1D70F}\sim 1/|c|$ via the approximation (3.35), as shown in figure 2(b). This approach is particularly useful when Lagrangian data are only available over short times and explicitly identifies the inertial initial orientation angle $\unicode[STIX]{x1D6FC}_{0,\infty }$ , allowing the Protean transform to be determined uniquely from $t=0$ .

3.2.2 Longitudinal and transverse deformation in 3-D steady flow

Given solution of (3.29) (such that $\unicode[STIX]{x1D608}_{23}=\tilde{\unicode[STIX]{x1D716}}_{32}$ ), reorientation into the Protean frame renders the streamline velocity gradient tensor $\unicode[STIX]{x1D750}^{\prime }(t)$ upper triangular. From (3.7) it follows that $\unicode[STIX]{x1D641}^{\prime }(t)$ is also upper triangular, with diagonal components

(3.36a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D60D}_{11}^{\prime }(t)=\frac{v(t)}{v(0)}, & \displaystyle\end{eqnarray}$$

(3.36b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D60D}_{ii}^{\prime }(t)=\exp \left[\int _{0}^{t}\text{d}t^{\prime }\unicode[STIX]{x1D716}_{ii}^{\prime }(t^{\prime })\right], & \displaystyle\end{eqnarray}$$

for $i=2,3$ . The off-diagonal elements are $\unicode[STIX]{x1D60D}_{ij}(t)=0$ for $i>j$ and else

(3.36c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D60D}_{12}^{\prime }(t)=v(t)\int _{0}^{t}\text{d}t^{\prime }\frac{\unicode[STIX]{x1D716}_{12}^{\prime }(t^{\prime })\unicode[STIX]{x1D60D}_{22}^{\prime }(t^{\prime })}{v(t^{\prime })}, & \displaystyle\end{eqnarray}$$

(3.36d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D60D}_{23}^{\prime }(t)=\unicode[STIX]{x1D60D}_{22}^{\prime }(t)\int _{0}^{t}\text{d}t^{\prime }\frac{\unicode[STIX]{x1D716}_{23}^{\prime }(t^{\prime })\unicode[STIX]{x1D60D}_{33}^{\prime }(t^{\prime })}{\unicode[STIX]{x1D60D}_{22}^{\prime }(t^{\prime })}, & \displaystyle\end{eqnarray}$$

(3.36e ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D60D}_{13}^{\prime }(t)=v(t)\int _{0}^{t}\text{d}t^{\prime }\frac{\unicode[STIX]{x1D716}_{12}^{\prime }(t^{\prime })\unicode[STIX]{x1D60D}_{23}^{\prime }(t^{\prime })+\unicode[STIX]{x1D716}_{13}^{\prime }(t^{\prime })\unicode[STIX]{x1D60D}_{33}^{\prime }(t^{\prime })}{v(t^{\prime })}. & \displaystyle\end{eqnarray}$$

Here the upper triangular form of $\unicode[STIX]{x1D750}^{\prime }(t)$ simplifies solution of the deformation tensor $\unicode[STIX]{x1D641}^{\prime }(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 $\unicode[STIX]{x1D60D}_{12}^{\prime }(t)$ and $\unicode[STIX]{x1D60D}_{13}^{\prime }(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. (Reference Dentz, Lester, Borgne and de Barros2016b ) for 2-D 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 $\unicode[STIX]{x1D641}$ into longitudinal and transverse components respectively as $\unicode[STIX]{x1D641}(t)=\unicode[STIX]{x1D641}_{l}+\unicode[STIX]{x1D641}_{t}$ , where

(3.37) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D641}_{l}\equiv \text{diag}(\boldsymbol{e}_{1})\boldsymbol{\cdot }\unicode[STIX]{x1D641}=\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D60D}_{11} & \unicode[STIX]{x1D60D}_{12} & \unicode[STIX]{x1D60D}_{13}\\ 0 & 0 & 0\\ 0 & 0 & 0\\ \end{array}\right), & \displaystyle\end{eqnarray}$$

(3.38) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D641}_{t}\equiv \text{diag}(\boldsymbol{e}_{2}+\boldsymbol{e}_{3})\boldsymbol{\cdot }\unicode[STIX]{x1D641}=\left(\begin{array}{@{}ccc@{}}0 & 0 & 0\\ 0 & \unicode[STIX]{x1D60D}_{22} & \unicode[STIX]{x1D60D}_{23}\\ 0 & 0 & \unicode[STIX]{x1D60D}_{33}\end{array}\right), & \displaystyle\end{eqnarray}$$

where $\text{diag}(\boldsymbol{a})$ is a diagonal matrix comprising the vector $\boldsymbol{a}$ along the diagonal. From (3.2), a differential fluid line element $\unicode[STIX]{x1D6FF}\boldsymbol{l}(\boldsymbol{X},t)$ at Lagrangian position $\boldsymbol{X}$ then evolves with time $t$ as

(3.39) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{l}(\boldsymbol{X},t) & = & \displaystyle \unicode[STIX]{x1D641}(\boldsymbol{X},t)\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{l}(\boldsymbol{X},0)\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D641}_{l}(\boldsymbol{X},t)\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{l}(\boldsymbol{X},0)+\unicode[STIX]{x1D641}_{t}(\boldsymbol{X},t)\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{l}(\boldsymbol{X},0)=\unicode[STIX]{x1D6FF}\boldsymbol{l}_{l}(\boldsymbol{X},t)+\unicode[STIX]{x1D6FF}\boldsymbol{l}_{t}(\boldsymbol{X},t),\end{eqnarray}$$

where the line element may also be decomposed into the longitudinal and transverse components as $\unicode[STIX]{x1D6FF}\boldsymbol{l}_{l}=\unicode[STIX]{x1D6FF}\boldsymbol{l}_{l}+\unicode[STIX]{x1D6FF}\boldsymbol{l}_{t}$ . Due to the upper triangular nature of $\unicode[STIX]{x1D641}$ the length $\unicode[STIX]{x1D6FF}l$ of these line elements can also be decomposed as

(3.40) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}l(\boldsymbol{X},t)\equiv |\unicode[STIX]{x1D6FF}\boldsymbol{l}(\boldsymbol{X},t)| & = & \displaystyle \sqrt{\unicode[STIX]{x1D6FF}\boldsymbol{l}(\boldsymbol{X},0)\boldsymbol{\cdot }\unicode[STIX]{x1D641}^{\top }(\boldsymbol{X},t)\boldsymbol{\cdot }\unicode[STIX]{x1D641}(\boldsymbol{X},t)\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{l}(\boldsymbol{X},0)}\nonumber\\ \displaystyle & = & \displaystyle \sqrt{\unicode[STIX]{x1D6FF}\boldsymbol{l}(\boldsymbol{X},0)\boldsymbol{\cdot }(\unicode[STIX]{x1D641}_{l}^{\top }(\boldsymbol{X},t)\boldsymbol{\cdot }\unicode[STIX]{x1D641}_{l}(\boldsymbol{X},t)+\unicode[STIX]{x1D641}_{t}^{\top }(\boldsymbol{X},t)\boldsymbol{\cdot }\unicode[STIX]{x1D641}_{t}(\boldsymbol{X},t))\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{l}(\boldsymbol{X},0)}\nonumber\\ \displaystyle & = & \displaystyle \sqrt{\unicode[STIX]{x1D6FF}l_{l}(\boldsymbol{X},t)^{2}+\unicode[STIX]{x1D6FF}l_{t}(\boldsymbol{X},t)^{2}}.\end{eqnarray}$$

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

(3.41) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}l^{2}=\unicode[STIX]{x1D6FF}l_{l}^{2}+\unicode[STIX]{x1D6FF}l_{t}^{2}.\end{eqnarray}$$

We denote the length of the 1-D lines ${\mathcal{H}}_{1\text{D}}^{(2)}$ , ${\mathcal{H}}_{1\text{D}}^{(3)}$ in figure 1 respectively as $l^{(2)}(t)$ , $l^{(3)}(\bar{x}_{1})$ , where the respective arguments $t$ , $\bar{x}_{1}$ reflect the fact that the 1-D line ${\mathcal{H}}_{1\text{D}}^{(2)}$ occurs at fixed time $t$ since injection, whereas the 1-D line ${\mathcal{H}}_{1\text{D}}^{(3)}$ occurs at fixed distance $\bar{x}_{1}$ downstream in the mean flow direction. Following the decomposition (3.41), $l^{(2)}(t)$ and $l^{(3)}(\bar{x}_{1})$ are then given by the contour integrals along the injection line ${\mathcal{H}}_{1\text{D}}^{(1)}$ in Lagrangian space as

(3.42) $$\begin{eqnarray}\displaystyle & \displaystyle l^{(2)}(t)=\int _{{\mathcal{H}}_{1\text{D}}^{(1)}}\sqrt{\unicode[STIX]{x1D6FF}l_{l}(\boldsymbol{X},t)^{2}+\unicode[STIX]{x1D6FF}l_{t}(\boldsymbol{X},t)^{2}}\,\text{d}s, & \displaystyle\end{eqnarray}$$

(3.43) $$\begin{eqnarray}\displaystyle & \displaystyle l^{(3)}(\bar{x}_{1})=\int _{{\mathcal{H}}_{1\text{D}}^{(1)}}\unicode[STIX]{x1D6FF}l_{l}(\boldsymbol{X},t_{1}(\boldsymbol{X},\bar{x}_{1}))\,\text{d}s, & \displaystyle\end{eqnarray}$$

where $\bar{t}_{1}(\boldsymbol{X},\bar{x}_{1})$ is the time at which the fluid streamline at Lagrangian coordinate $\boldsymbol{X}$ reaches $\bar{x}_{1}$ , and $\text{d}s$ is a differential increment associated with changes in $\boldsymbol{X}$ along ${\mathcal{H}}_{1\text{D}}^{(1)}$ . Hence stretching of the 1-D material line ${\mathcal{H}}_{1\text{D}}^{(2)}$ in figure 1 is governed by both the transverse and longitudinal deformations $\unicode[STIX]{x1D641}_{l}$ , $\unicode[STIX]{x1D641}_{t}$ , whereas stretching of the 1-D line ${\mathcal{H}}_{1\text{D}}^{(2)}$ transverse to the mean flow direction is solely governed by the transverse deformation $\unicode[STIX]{x1D641}_{t}$ .

As discussed in § 1, the deformation of these different lines is important for various applications. For example, Le Borgne et al. (Reference Le Borgne, Dentz and Villermaux2013, Reference Le Borgne, Dentz and Villermaux2015) use the growth rate of $l^{(2)}(t)$ to predict the mixing of a pulsed tracer injection (illustrated as ${\mathcal{H}}_{1\text{D}}^{(2)}$ in figure 1) in a steady 2-D Darcy flow. Conversely, Lester et al. (Reference Lester, Dentz and Le Borgne2016) use the growth rate of $l^{(3)}(t)$ to predict mixing of a continuously injected source in steady 3-D pore-scale flow. These different deformation rates indicate that a single scalar cannot be used to characterise fluid deformation in steady random 3-D flow. Instead it is necessary to characterise both longitudinal and transverse fluid deformation in such flows.

As such, we characterise longitudinal $\unicode[STIX]{x1D6EC}_{l}(t)$ and transverse $\unicode[STIX]{x1D6EC}_{t}(t)$ fluid deformation in terms of the corresponding deformation tensors as

(3.44) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6EC}_{l}(t)\equiv \langle \Vert \unicode[STIX]{x1D641}_{l}(t)\Vert \rangle =\langle \sqrt{\unicode[STIX]{x1D60D}_{11}(t)^{2}+\unicode[STIX]{x1D60D}_{12}(t)^{2}+\unicode[STIX]{x1D60D}_{13}(t)^{2}}\rangle , & \displaystyle\end{eqnarray}$$

(3.45) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6EC}_{t}(t)\equiv \langle \Vert \unicode[STIX]{x1D641}_{t}(t)\Vert \rangle =\langle \sqrt{\unicode[STIX]{x1D60D}_{22}(t)^{2}+\unicode[STIX]{x1D60D}_{23}(t)^{2}+\unicode[STIX]{x1D60D}_{33}(t)^{2}}\rangle , & \displaystyle\end{eqnarray}$$

where the angled brackets denote an ensemble average over multiple realisations of a 3-D steady random flow. We also define the ensemble-averaged length of an arbitrary fluid line as

(3.46) $$\begin{eqnarray}\ell (t)\equiv \langle \unicode[STIX]{x1D6FF}l(\boldsymbol{X},t)\rangle ,\end{eqnarray}$$

and in § 4, we seek to derive stochastic models for the growth of $\unicode[STIX]{x1D6EC}_{l}(t)$ , $\unicode[STIX]{x1D6EC}_{t}(t)$ and $\ell (t)$ in steady random 3-D flows. In 2-D steady flow, topological constraints (associated with conservation of area) prevent persistent growth of transverse deformation which simply fluctuates as $\unicode[STIX]{x1D6EC}_{t}(t)=\langle \unicode[STIX]{x1D60D}_{22}(t)\rangle =1/\langle \unicode[STIX]{x1D60D}_{11}(t)\rangle =\langle v(0)/v(t)\rangle =1$ . Conversely, longitudinal deformation can grow persistently as $\unicode[STIX]{x1D6EC}_{l}(t)\sim \unicode[STIX]{x1D60D}_{12}(t)^{2}$ , and Dentz et al. (Reference Dentz, Lester, Borgne and de Barros2016b ) derive stochastic models for the evolution of $\unicode[STIX]{x1D6EC}_{l}(t)$ from velocity field statistical properties as a CTRW. In § 4 we extend this approach to steady random 3-D flows to yield stochastic models for $\unicode[STIX]{x1D6EC}_{l}(t)$ , $\unicode[STIX]{x1D6EC}_{t}(t)$ and $\ell (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 $\boldsymbol{z}(t)$ initially aligned with the $1$ -direction of the Protean coordinate system, i.e.  $\boldsymbol{z}(0)=(z_{0},0,0)^{\top }$ we obtain the elongation $\ell (t)=\sqrt{\boldsymbol{z}(t)^{2}}$

(3.47) $$\begin{eqnarray}\ell (t)=z_{0}\frac{v(t)}{v(0)}.\end{eqnarray}$$

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.  $\boldsymbol{z}(0)=(0,z_{0},0)^{\top }$ , we obtain the elongation

(3.48) $$\begin{eqnarray}\ell (t)=z_{0}\sqrt{\unicode[STIX]{x1D60D}_{12}(t)^{2}+\unicode[STIX]{x1D60D}_{22}(t)^{2}}.\end{eqnarray}$$

For an initial alignment with the $3$ -direction, $\boldsymbol{z}(0)=(0,0,z_{0})^{\top }$

(3.49) $$\begin{eqnarray}\ell (t)=z_{0}\sqrt{\unicode[STIX]{x1D60D}_{13}(t)^{2}+\unicode[STIX]{x1D60D}_{23}(t)^{2}+\unicode[STIX]{x1D60D}_{33}(t)^{2}}.\end{eqnarray}$$

In general, we have

(3.50) $$\begin{eqnarray}\ell (t)=\sqrt{\boldsymbol{z}_{0}\unicode[STIX]{x1D641}^{\top }(t)\unicode[STIX]{x1D641}(t)\boldsymbol{z}_{0}},\end{eqnarray}$$

where $\boldsymbol{z}(t=0)\equiv \boldsymbol{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 $\unicode[STIX]{x1D641}^{\prime }(t)$ facilitates the identification of the stochastic dynamics of fluid deformation in steady random flows and its quantification in terms of the fluid velocity and velocity gradient statistics.

3.3 Zero-helicity-density flows

One important class of steady $D=3$ dimensional flows are zero helicity density (or complex lamellar (Finnigan Reference Finnigan1983)) flows such as the isotropic Darcy flow (Sposito Reference Sposito2001)

(3.51) $$\begin{eqnarray}\boldsymbol{v}(\boldsymbol{x})=-k(\boldsymbol{x})\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}(\boldsymbol{x}),\end{eqnarray}$$

which is commonly used to model Darcy-scale flow in heterogeneous porous media. Here $\boldsymbol{v}(\boldsymbol{x})$ is the specific discharge, the scalar $k(\boldsymbol{x})$ represents hydraulic conductivity and $\unicode[STIX]{x1D719}(\boldsymbol{x})$ is the flow potential (or velocity head). For such flows the helicity density (Moffatt Reference Moffatt1969) which measures the helical motion of the flow,

(3.52) $$\begin{eqnarray}h(\boldsymbol{x})\equiv \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D74E}=k(\boldsymbol{x})\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}(\boldsymbol{x})\boldsymbol{\cdot }(\unicode[STIX]{x1D735}k(\boldsymbol{x})\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D719}(\boldsymbol{x})),\end{eqnarray}$$

is identically zero (where $\unicode[STIX]{x1D74E}$ is the vorticity vector). As shown by Kelvin (Reference Kelvin1884), all 3-D zero-helicity-density flows are complex lamellar, and so may be posed in the form of an isotropic Darcy flow $\boldsymbol{v}(\boldsymbol{x})=-k(\boldsymbol{x})\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ , where for porous media flow $\unicode[STIX]{x1D719}$ is the flow potential and $k(\boldsymbol{x})$ represents the (possibly heterogeneous) hydraulic conductivity. Arnol’d (Reference Arnol’d1965, Reference Arnol’d1966) 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 $\unicode[STIX]{x1D713}_{1}$ , $\unicode[STIX]{x1D713}_{2}$ , with $\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{2}=0$ . As such, all zero helicity density flows may be represented as

(3.53) $$\begin{eqnarray}\boldsymbol{v}(\boldsymbol{x})=-k(\boldsymbol{x})\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{1}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{2}.\end{eqnarray}$$

As the gradients $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ , $\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{1}$ , $\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{2}$ are all mutually orthogonal, zero-helicity-density flows admit an orthogonal streamline coordinate system ( $\unicode[STIX]{x1D719}$ , $\unicode[STIX]{x1D713}_{1}$ , $\unicode[STIX]{x1D713}_{2}$ ) with unit base vectors

(3.54a-c ) $$\begin{eqnarray}\hat{\boldsymbol{e}}_{1}=-\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}}{\Vert \unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\Vert }=\frac{\boldsymbol{v}}{\Vert \boldsymbol{v}\Vert },\quad \hat{\boldsymbol{e}}_{2}=\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{1}}{\Vert \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{1}\Vert },\quad \hat{\boldsymbol{e}}_{3}=\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{2}}{\Vert \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{2}\Vert }.\end{eqnarray}$$

If the streamfunctions $\unicode[STIX]{x1D713}_{1}$ , $\unicode[STIX]{x1D713}_{2}$ are known, there is no need to explicitly solve the orientation angle $\unicode[STIX]{x1D6FC}(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 $\unicode[STIX]{x1D713}_{1}$ , $\unicode[STIX]{x1D713}_{2}$ (Holm & Kimura Reference Holm and Kimura1991), hence fluid deformation must be sub-exponential (algebraic). From (3.36b ), the ensemble average of the principal components $\unicode[STIX]{x1D716}_{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 $\tilde{\unicode[STIX]{x1D716}}_{23}=\tilde{\unicode[STIX]{x1D716}}_{32}$ , due to the identity

(3.55) $$\begin{eqnarray}h=\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D74E}=\boldsymbol{v}\boldsymbol{\cdot }(\unicode[STIX]{x1D700}_{ijk}:\unicode[STIX]{x1D750})=\tilde{\boldsymbol{v}}\boldsymbol{\cdot }(\unicode[STIX]{x1D700}_{ijk}:\tilde{\unicode[STIX]{x1D750}})=v(\tilde{\unicode[STIX]{x1D716}}_{23}-\tilde{\unicode[STIX]{x1D716}}_{32})=0,\end{eqnarray}$$

where $\unicode[STIX]{x1D700}_{ijk}$ is the Levi-Civita tensor. Whilst (3.55) shows $\tilde{\unicode[STIX]{x1D716}}_{23}=\tilde{\unicode[STIX]{x1D716}}_{32}$ in general for zero helicity flows, there exists a specific value of $\unicode[STIX]{x1D6FC}$ , denoted $\tilde{\unicode[STIX]{x1D6FC}}$ , which renders both $\tilde{\unicode[STIX]{x1D716}}_{23}$ and $\tilde{\unicode[STIX]{x1D716}}_{32}$ zero. This form of the velocity gradient (namely zero (2,3) and (3,2) components) is reflected by its representation in the 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 $\unicode[STIX]{x1D63C}(t)$ to the (2,3) and (3,2) components of the velocity gradient. Typically, a different value of $\unicode[STIX]{x1D6FC}$ is used than $\tilde{\unicode[STIX]{x1D6FC}}$ to render $\unicode[STIX]{x1D716}_{32}^{\prime }=0$ , and so for zero helicity flow

(3.56) $$\begin{eqnarray}\unicode[STIX]{x1D716}_{23}^{\prime }=2\tilde{\unicode[STIX]{x1D716}}_{23},\end{eqnarray}$$

which is typically non-zero, hence $\unicode[STIX]{x1D60D}_{23}^{\prime }\neq 0$ in general. Note that as both the Protean coordinate frame and the curvilinear streamline coordinate system render the deformation tensor component $\unicode[STIX]{x1D60D}_{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 $\unicode[STIX]{x1D713}_{1}$ , $\unicode[STIX]{x1D713}_{2}$ effectively means that steady 3-D 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 2-D integral surfaces are topological cylinders or tori, streamlines confined within these surfaces cannot diverge exponentially in space. Thus, the principal transverse deformations $\unicode[STIX]{x1D60D}_{22}$ , $\unicode[STIX]{x1D60D}_{33}$ in the Protean frame may only fluctuate about the unit mean as

(3.57) $$\begin{eqnarray}\langle \unicode[STIX]{x1D60D}_{ii}\rangle =1,\quad i=1,2,3,\end{eqnarray}$$

and so the only persistent growth of fluid deformation in steady zero helicity flows only occurs via the longitudinal and transverse shear deformations $\unicode[STIX]{x1D60D}_{12}$ , $\unicode[STIX]{x1D60D}_{13}$ , $\unicode[STIX]{x1D60D}_{23}$ . Secondly, as the streamlines of this steady flow are confined to 2-D 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.

(3.58) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{t\rightarrow \infty }\langle \unicode[STIX]{x1D60D}_{12}(s)\rangle \sim t^{r_{12}}, & \displaystyle\end{eqnarray}$$

(3.59) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{t\rightarrow \infty }\langle \unicode[STIX]{x1D60D}_{13}(s)\rangle \sim t^{r_{13}}, & \displaystyle\end{eqnarray}$$

(3.60) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{t\rightarrow \infty }\langle \unicode[STIX]{x1D60D}_{23}(s)\rangle \sim t^{r_{23}}. & \displaystyle\end{eqnarray}$$

Dentz et al. (Reference Dentz, Lester, Borgne and de Barros2016b ) develop a CTRW model for shear deformation in steady random incompressible 2-D flows and uncover algebraic stretching of $\unicode[STIX]{x1D60D}_{12}\sim t^{r}$ (again in Protean coordinates) and show that the index $r$ depends upon the coupling between shear deformation and velocity fluctuations long streamlines, ranging from diffusive $r=1/2$ to superdiffusive $r=2$ stretching.

Hence fluid deformation is constrained to be algebraic in steady 3-D zero-helicity-density flow, and evolves in a similar fashion to steady 2-D flow,

(3.61) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6EC}_{t}(t)\sim t^{r_{23}}, & \displaystyle\end{eqnarray}$$

(3.62) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6EC}_{l}(t)\sim t^{max(r_{12},r_{13})}, & \displaystyle\end{eqnarray}$$

where the indices $r$ are governed by intermittency of the fluid shear events (Dentz et al. Reference Dentz, Lester, Borgne and de Barros2016b ). We note that these constraints do not apply to anisotropic Darcy flows (i.e.  $\boldsymbol{v}(\boldsymbol{x})=-\boldsymbol{K}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}(\boldsymbol{x})$ , where $\boldsymbol{K}(\boldsymbol{x})$ is the tensorial hydraulic conductivity), as these no longer correspond to zero-helicity-density flows. Indeed, as demonstrated by Ye et al. (Reference Ye, Chiogna, Cirpka, Grathwohl and Rolle2015), 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.

3.3.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\rightarrow 0$ . To circumvent this issue, we consider fluid deformation along a stagnation line in the neighbourhood of a stagnation point $\boldsymbol{x}_{p}$ , where the velocity gradient may be approximated as $\unicode[STIX]{x1D750}(t)\approx \unicode[STIX]{x1D750}_{0}\equiv \unicode[STIX]{x1D750}|_{\boldsymbol{x}=\boldsymbol{x}_{p}}$ . Note that for streamlines approaching a stagnation point from upstream, $\unicode[STIX]{x1D716}_{11,0}$ must necessarily be negative. In the neighbourhood of the stagnation point, the velocity evolves as $v(t)/v(t_{0})\approx \exp (\unicode[STIX]{x1D716}_{11,0}(t-t_{0}))$ , and the deformation gradient tensor can be solved directly from (3.4) as

(3.63) $$\begin{eqnarray}\unicode[STIX]{x1D641}(t)\approx \unicode[STIX]{x1D641}(t_{0})\boldsymbol{\cdot }\exp (\unicode[STIX]{x1D750}_{0}(t-t_{0})),\end{eqnarray}$$

where $t_{0}$ is the time at which the fluid element enters the neighbourhood of the stagnation point (defined as the region for which $\unicode[STIX]{x1D750}(t)\approx \unicode[STIX]{x1D750}_{0}$ ). Note that (3.36) represents an explicit solution of the matrix exponential above, via the substitutions $t\mapsto t-t_{0}$ , $v(t)\mapsto v(t_{0})\exp (\unicode[STIX]{x1D716}_{11,0}(t-t_{0}))$ , $\unicode[STIX]{x1D750}(t)\mapsto \unicode[STIX]{x1D750}_{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 posteriori via (3.63).

4.1 Continous-time random walks for fluid velocities

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

(4.1a,b ) $$\begin{eqnarray}\frac{\text{d}s(t)}{\text{d}t}=v_{s}[s(t)],\quad \frac{\text{d}t(s)}{\text{d}s}=\frac{1}{v_{s}(s)},\end{eqnarray}$$

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. Reference Dentz, Kang, Comolli, Le Borgne and Lester2016a ) for flow through heterogeneous porous and fractured media from pore to Darcy and regional scales (Berkowitz et al. Reference Berkowitz, Cortis, Dentz and Scher2006; Fiori et al. Reference Fiori, Jankovic, Dagan and Cvetkovic2007; Le Borgne et al. Reference Le Borgne, Dentz and Carrera2008a ,Reference Le Borgne, Dentz and Carrera b ; Bijeljic, Mostaghimi & Blunt Reference Bijeljic, Mostaghimi and Blunt2011; De Anna et al. Reference De Anna, Le Borgne, Dentz, Tartakovsky, Bolster and Davy2013; Holzner et al. Reference Holzner, Morales, Willmann and Dentz2015; Kang et al. Reference Kang, Dentz, Le Borgne and Juanes2015). This property stems from the fact that fluid velocities evolve on characteristic spatial scales imprinted in the spatial flow structure. Streamwise velocities $v_{s}(s)$ are correlated on the correlation scale $\ell _{c}$ . For distances larger than $\ell _{c}$ subsequent particle velocities can be considered independent. Thus, particle motion along streamlines can be represented by the recursion relations

(4.2a,b ) $$\begin{eqnarray}s_{n+1}=s_{n}+\ell _{c},\quad t_{n+1}=t_{n}+\frac{\ell _{c}}{v_{n}},\end{eqnarray}$$

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}\leqslant 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 $\unicode[STIX]{x1D70F}_{u}=\ell _{c}/\overline{u}$ , where $\overline{u}$ is the mean flow velocity. We consider flow fields that are characterised by open ergodic streamlines. Thus, $v_{n}$ can be modelled as independent identically distributed random variables that are characterised by $p_{s}(v)$ . The latter is related to the probability density function (PDF) $p_{e}(v)$ of Eulerian velocities $v_{e}$ through flux weighting as (Dentz et al. Reference Dentz, Kang, Comolli, Le Borgne and Lester2016a )

(4.3) $$\begin{eqnarray}p_{s}(v)=\frac{vp_{e}(v)}{\langle v_{e}\rangle }.\end{eqnarray}$$

The advective transition time over the characteristic distance $\ell _{c}$ is defined by

(4.4) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{n}=\frac{\ell _{c}}{v_{n}}.\end{eqnarray}$$

Its distribution $\unicode[STIX]{x1D713}(t)$ is obtained from (4.3) in terms of the Eulerian PDF $p_{e}(v)$ as

(4.5) $$\begin{eqnarray}\unicode[STIX]{x1D713}(t)=\frac{\ell _{c}p_{e}(\ell _{c}/t)}{\langle v_{e}\rangle t^{3}}.\end{eqnarray}$$

Equations (4.2) constitute a continuous-time random walk (CTRW) in that the time increment $\unicode[STIX]{x1D70F}_{n}\equiv \ell _{c}/v_{n}$ at the $n$ th 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.

In order to illustrate this notion, we consider the joint evolution of a process $A_{n}$ and time $t_{n}$ along a streamline (Scher & Lax Reference Scher and Lax1973),

(4.6a,b ) $$\begin{eqnarray}A_{n+1}=A_{n}+\unicode[STIX]{x0394}A_{n},\quad t_{n+1}=t_{n}+\unicode[STIX]{x1D70F}_{n}.\end{eqnarray}$$

The increments $(\unicode[STIX]{x0394}A,\unicode[STIX]{x1D70F})$ are stepwise independent and in general coupled. They are characterised by the joint PDF $\unicode[STIX]{x1D713}_{a}(a,t)$ , which can be written as

(4.7) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{a}(a,t)=\unicode[STIX]{x1D713}_{a}(a|t)\unicode[STIX]{x1D713}(t),\end{eqnarray}$$

The PDF $p_{a}(a,t)$ of $A(t)$ is then given by

(4.8) $$\begin{eqnarray}p_{a}(a,t)=\int _{0}^{t}\text{d}t^{\prime }R(a,t^{\prime })\int _{t-t^{\prime }}^{\infty }\text{d}t^{\prime \prime }\unicode[STIX]{x1D713}(t),\end{eqnarray}$$

where $R(a,t)$ denotes the probability per time that $A(t)$ has just assumed the value $a$ . Thus, equation (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 probability per time $R(a,t^{\prime })$ that $A(t)$ has arrived at $a$ times the probability that the next transition takes longer than $t-t^{\prime }$ . The $R(a,t)$ satisfies the Chapman–Kolmogorov-type equation

(4.9) $$\begin{eqnarray}R(a,t)=R_{0}(a,t)+\int _{0}^{t}\text{d}t^{\prime }\int \text{d}a\unicode[STIX]{x1D713}(a-a^{\prime },t-t^{\prime })R(a^{\prime },t^{\prime }),\end{eqnarray}$$

where $R_{0}(a,t)$ is the joint distribution of $A_{0}$ and $t_{0}$ . Equations (4.8) and (4.9) can be combined into the following generalised master equation (Kenkre, Montroll & Shlesinger Reference Kenkre, Montroll and Shlesinger1973) for the evolution of $p_{a}(t)$ :

(4.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p_{a}(a,t)}{\unicode[STIX]{x2202}t}=\int _{0}^{t}\text{d}t^{\prime }\int \text{d}a\unicode[STIX]{x1D6F7}(a-a^{\prime },t-t^{\prime })[p_{a}(a^{\prime },t^{\prime })-p_{a}(a,t^{\prime })],\end{eqnarray}$$

where $\unicode[STIX]{x1D6F7}(a,t)$ is defined in Laplace space as

(4.11) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}^{\ast }(a,\unicode[STIX]{x1D706})=\frac{\unicode[STIX]{x1D706}\unicode[STIX]{x1D713}_{a}^{\ast }(a,\unicode[STIX]{x1D706})}{1-\unicode[STIX]{x1D713}^{\ast }(\unicode[STIX]{x1D706})}.\end{eqnarray}$$

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 provides a basis for stochastic modelling of the deformation equations (3.36b,e ) as a coupled continuous-time random walk (coupled CTRW) along streamlines, in a similar fashion to that developed by Dentz et al. (Reference Dentz, Lester, Borgne and de Barros2016b ) for deformation in steady 2-D 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 (Fiori et al. Reference Fiori, Jankovic, Dagan and Cvetkovic2007; Edery et al. Reference Edery, Guadagnini, Scher and Berkowitz2014; Dentz et al. Reference Dentz, Kang, Comolli, Le Borgne and Lester2016a ; Tyukhova et al. Reference Tyukhova, Dentz, Kinzelbach and Willmann2016).

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 $\hat{\unicode[STIX]{x1D641}}(s)=\unicode[STIX]{x1D641}[t(s)]$

(4.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D60D}}_{11}(s)=\frac{v_{s}(s)}{v_{s}(0)}, & \displaystyle\end{eqnarray}$$

(4.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D60D}}_{ii}(s)=\exp \left[\int _{0}^{s}\text{d}s^{\prime }\frac{\hat{\unicode[STIX]{x1D716}}_{ii}(s^{\prime })}{v_{s}(s^{\prime })}\right], & \displaystyle\end{eqnarray}$$

with $i=2,3$ . For the off-diagonal elements, we obtain

(4.12c ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D60D}}_{12}(s)=v_{s}(s)\int _{0}^{s}\text{d}s^{\prime }\frac{\hat{\unicode[STIX]{x1D716}}_{12}(s^{\prime })\hat{\unicode[STIX]{x1D60D}}_{22}(s^{\prime })}{v_{s}(s^{\prime })^{2}}, & \displaystyle\end{eqnarray}$$

(4.12d ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D60D}}_{23}(s)=\hat{\unicode[STIX]{x1D60D}}_{22}(s)\int _{0}^{s}\text{d}s^{\prime }\frac{\hat{\unicode[STIX]{x1D716}}_{23}(s^{\prime })\hat{\unicode[STIX]{x1D60D}}_{33}(s^{\prime })}{v_{s}(s^{\prime })\unicode[STIX]{x1D60D}_{22}(s^{\prime })}, & \displaystyle\end{eqnarray}$$

(4.12e ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D60D}}_{13}(s)=v_{s}(s)\int _{0}^{s}\text{d}s^{\prime }\frac{\hat{\unicode[STIX]{x1D716}}_{12}(s^{\prime })\hat{\unicode[STIX]{x1D60D}}_{23}(s^{\prime })+\hat{\unicode[STIX]{x1D716}}_{13}(s^{\prime })\hat{\unicode[STIX]{x1D60D}}_{33}(s^{\prime })}{v_{s}(s^{\prime })^{2}}, & \displaystyle\end{eqnarray}$$

while $\hat{\unicode[STIX]{x1D60D}}_{ij}(s)=0$ for $i>j$ . We denote $\hat{\unicode[STIX]{x1D716}}_{ij}(s)=\unicode[STIX]{x1D716}_{ij}[t(s)]$ . In the following, we investigate the evolution of stretching in different time regimes for chaotic steady random flow.

4.2 Chaotic steady random flow

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

(4.13) $$\begin{eqnarray}\unicode[STIX]{x1D706}\equiv \lim _{t\rightarrow \infty }\lim _{z_{0}\rightarrow 0}\frac{1}{2t}\ln \left[\frac{\ell (t)^{2}}{\ell (0)^{2}}\right].\end{eqnarray}$$

In the limit of $t\gg \unicode[STIX]{x1D706}^{-1}$ , we may set

(4.14a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D60D}_{11}(t)=1,\quad \unicode[STIX]{x1D60D}_{22}(t)=\exp (\unicode[STIX]{x1D716}t),\quad \unicode[STIX]{x1D60D}_{33}(t)=\exp (-\unicode[STIX]{x1D716}t),\end{eqnarray}$$

where we defined

(4.15) $$\begin{eqnarray}\unicode[STIX]{x1D716}=\lim _{t\rightarrow \infty }\frac{1}{t}\int _{0}^{t}\text{d}t^{\prime }\unicode[STIX]{x1D716}_{22}(t^{\prime })=-\!\lim _{t\rightarrow \infty }\frac{1}{t}\int _{0}^{t}\text{d}t^{\prime }\unicode[STIX]{x1D716}_{33}(t^{\prime }).\end{eqnarray}$$

The last equation is due to volume conservation. The long-time behaviour is fully dominated by exponential stretching. The characteristic stretching time scale for exponential stretching is given by $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}=1/\unicode[STIX]{x1D716}$ , while the characteristic advection scale is $\unicode[STIX]{x1D70F}_{u}=\ell _{c}/\overline{u}$ with $\overline{u}$ the mean and $\ell _{c}$ the correlation scale of the steady random flow field $\boldsymbol{u}(\boldsymbol{x})$ . The latter sets the end of the ballistic regimes: see below. If $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}$ and $\unicode[STIX]{x1D70F}_{u}$ are well separated, $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}\gg \unicode[STIX]{x1D70F}_{u}$ , we observe a subexponential stretching regime that is dominated by heterogeneous shear action. In the following, we analyse the behaviour of deformation for time $t\ll \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}$ . We first briefly discuss the ballistic regime for which $t\ll \unicode[STIX]{x1D70F}_{u}$ . Then, we develop a CTRW based approach for deformation in the pre-exponential regime $\unicode[STIX]{x1D70F}_{u}\ll t\ll \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}$ .

4.2.1 Ballistic regime: $t\ll \unicode[STIX]{x1D70F}_{u}$

We first consider the ballistic short-time regime, for which $t\ll \unicode[STIX]{x1D70F}_{u}$ . In this regime, the flow properties are essentially constant and we obtain the approximation for the deformation tensor (3.36),

(4.16a ) $$\begin{eqnarray}\unicode[STIX]{x1D60D}_{ii}(t)=1,\end{eqnarray}$$

for $i=1,2,3$ and for the off-diagonal elements with $i>j$ ,

(4.16b-d ) $$\begin{eqnarray}\unicode[STIX]{x1D60D}_{12}(t)=\unicode[STIX]{x1D716}_{12}t,\quad \unicode[STIX]{x1D60D}_{23}(t)=\unicode[STIX]{x1D716}_{23}t,\quad \unicode[STIX]{x1D60D}_{13}(t)=(\unicode[STIX]{x1D716}_{12}\unicode[STIX]{x1D716}_{23}t/2+\unicode[STIX]{x1D716}_{13})t.\end{eqnarray}$$

As $\unicode[STIX]{x1D716}_{ij}\sim 1/\unicode[STIX]{x1D70F}_{u}$ , in this regime the transverse (3.45) and longitudinal (3.44) deformations as well as the arbitrary fluid element $\ell (t)$ all evolve ballistically to leading order: $\unicode[STIX]{x1D6EC}_{t}(t)$ , $\unicode[STIX]{x1D6EC}_{l}(t)$ , $\ell (t)\propto t$ .

4.2.2 Exponential stretching regime $t\gg \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}$

In the stretching regime, for times much larger than $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}$ , the principal deformations scale as

(4.17a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D60D}_{11}(t)=1,\quad \unicode[STIX]{x1D60D}_{22}(t)=\exp (\unicode[STIX]{x1D716}t),\quad \unicode[STIX]{x1D60D}_{33}(t)=\exp (-\unicode[STIX]{x1D716}t).\end{eqnarray}$$

In this regime the expressions for the off-diagonal elements of $\unicode[STIX]{x1D641}(t)$ in (3.36) can be approximated by

(4.18) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D60D}_{12}(t)\approx \frac{\unicode[STIX]{x1D716}_{12}^{\prime }(t)\exp (\unicode[STIX]{x1D716}t)}{\unicode[STIX]{x1D716}}, & \displaystyle\end{eqnarray}$$

(4.19) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D60D}_{23}(t)\approx \exp (\unicode[STIX]{x1D716}t)\int _{0}^{t}\text{d}t^{\prime }\unicode[STIX]{x1D716}_{23}^{\prime }(t^{\prime })\exp (-2\unicode[STIX]{x1D716}t^{\prime }), & \displaystyle\end{eqnarray}$$

(4.20) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D60D}_{13}(t)\approx \unicode[STIX]{x1D60D}_{12}(t)\unicode[STIX]{x1D716}_{12}(t)\int _{0}^{t}\text{d}t^{\prime }\unicode[STIX]{x1D716}_{23}^{\prime }(t^{\prime })\exp (-2\unicode[STIX]{x1D716}t^{\prime }). & \displaystyle\end{eqnarray}$$

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: $\unicode[STIX]{x1D6EC}_{l}(t)$ , $\unicode[STIX]{x1D6EC}_{t}(t)$ , $\ell (t)\propto \exp (\unicode[STIX]{x1D716}t)$ . From (4.13) the Lyapunov exponent is then $\unicode[STIX]{x1D706}\equiv \unicode[STIX]{x1D716}$ .

4.2.3 Intermediate regime $t>\unicode[STIX]{x1D70F}_{u}$

For times $t>\unicode[STIX]{x1D70F}_{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 $\hat{\unicode[STIX]{x1D60D}}_{11}(s)\approx 1$ , $\hat{\unicode[STIX]{x1D60D}}_{22}(s)\approx \exp [\unicode[STIX]{x1D706}t(s)]$ and $\unicode[STIX]{x1D60D}_{33}\approx \exp [-\unicode[STIX]{x1D706}t(s)]$ . The off-diagonal components are then expressed as

(4.21a,b ) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D60D}}_{12}(s)=\frac{v_{s}(s)}{\ell _{c}}r_{2}(s),\quad \hat{\unicode[STIX]{x1D60D}}_{23}(s)=\hat{\unicode[STIX]{x1D60D}}_{22}(s)p(s), & & \displaystyle\end{eqnarray}$$

(4.21c ) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D60D}}_{13}(s)=v_{s}(s)\int _{0}^{s}\text{d}s^{\prime }\frac{\hat{\unicode[STIX]{x1D716}}_{12}(s^{\prime })\hat{\unicode[STIX]{x1D60D}}_{23}(s^{\prime })}{v_{s}(s^{\prime })^{2}}+\frac{v(s)}{\ell _{c}}r_{3}(s), & & \displaystyle\end{eqnarray}$$

where we defined the auxiliary functions

(4.22a,b ) $$\begin{eqnarray}\frac{\text{d}r_{i}(s)}{\text{d}s}=\ell _{c}\frac{\hat{\unicode[STIX]{x1D716}}_{1i}(s)\hat{\unicode[STIX]{x1D60D}}_{ii}(s)}{v_{s}(s)^{2}},\quad i=2,3,\quad \frac{\text{d}p(s)}{\text{d}s}=\frac{\hat{\unicode[STIX]{x1D716}}_{23}(s)\hat{\unicode[STIX]{x1D60D}}_{33}(s)}{v_{s}(s)\unicode[STIX]{x1D60D}_{22}(s)},\end{eqnarray}$$

with $r_{i}(s=0)=0$ and $p(s=0)=0$ . We approximated the off-diagonal components $\hat{\unicode[STIX]{x1D716}}_{ij}(s)$ for $i>j$ as being $\unicode[STIX]{x1D6FF}$ -correlated on the scales of interest such that

(4.23) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D716}}_{ij}(s)=\sqrt{\frac{\unicode[STIX]{x1D70E}_{ij}^{2}}{\ell _{c}}}\unicode[STIX]{x1D709}_{ij}(s),\end{eqnarray}$$

where $\unicode[STIX]{x1D709}(s)$ denotes white noise with zero mean and correlation $\langle \unicode[STIX]{x1D709}_{ij}(s)\unicode[STIX]{x1D709}_{ij}(s^{\prime })\rangle =\unicode[STIX]{x1D6FF}(s-s^{\prime })$ , and the $\unicode[STIX]{x1D70E}_{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 $\ell _{c}$ such that $s_{n}=n\ell _{c}$ and using (4.23), we obtain

(4.24a ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{i,n+1}=r_{i,n}+\unicode[STIX]{x1D70F}_{n}^{2}\unicode[STIX]{x1D70E}_{1i}\exp [(-1)^{i}\unicode[STIX]{x1D706}t_{n}]\unicode[STIX]{x1D709}_{1i,n},\quad i=2,3, & \displaystyle\end{eqnarray}$$

(4.24b ) $$\begin{eqnarray}\displaystyle & \displaystyle p_{n+1}=p_{n}+\unicode[STIX]{x1D70F}_{n}\unicode[STIX]{x1D70E}_{23}\exp (-2\unicode[STIX]{x1D706}t_{n})\unicode[STIX]{x1D709}_{23,n}. & \displaystyle\end{eqnarray}$$

These recursive relations describe a fully coupled continuous time random walk (Zaburdaev, Denisov & Klafter Reference Zaburdaev, Denisov and Klafter2015) with non-stationary increments due to the explicit dependence on time $t_{n}$ . Note that $\unicode[STIX]{x1D60D}_{ij}(t)=\hat{\unicode[STIX]{x1D60D}}_{ij}[s_{n_{t}}]=\hat{\unicode[STIX]{x1D60D}}_{ij}(n_{t}\ell _{c})$ . Thus, from (4.24) we obtain for the mean square deformations

(4.25) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D60D}_{12}(t)^{2}\rangle =\frac{\unicode[STIX]{x1D70E}_{12}^{2}\langle v^{2}\rangle }{\ell _{c}^{2}}B(t), & \displaystyle\end{eqnarray}$$

(4.26) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D60D}_{23}(t)^{2}\rangle =\frac{\unicode[STIX]{x1D70E}_{23}^{2}\langle \unicode[STIX]{x1D70F}^{2}\rangle }{4\unicode[STIX]{x1D706}}\exp (2\unicode[STIX]{x1D706}t)[1-\exp (-4\unicode[STIX]{x1D706}t)], & \displaystyle\end{eqnarray}$$

(4.27) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D60D}_{13}(t)^{2}\rangle =\frac{\langle v^{2}\rangle \unicode[STIX]{x1D70E}_{12}^{2}\unicode[STIX]{x1D70E}_{23}^{2}\langle \unicode[STIX]{x1D70F}^{2}\rangle }{4\unicode[STIX]{x1D706}\ell _{c}^{2}}A(t)+\frac{\langle v^{2}\rangle \unicode[STIX]{x1D70E}_{13}^{2}}{\ell _{c}}B(t), & \displaystyle\end{eqnarray}$$

where we define

(4.28) $$\begin{eqnarray}\displaystyle & \displaystyle A(t)=\left\langle \mathop{\sum }_{i=0}^{n_{t}-1}\unicode[STIX]{x1D70F}_{i}^{4}\exp (2\unicode[STIX]{x1D706}\text{i}\langle \unicode[STIX]{x1D70F}\rangle )[1-\exp (-4\unicode[STIX]{x1D706}\text{i}\langle \unicode[STIX]{x1D70F}\rangle )]\right\rangle , & \displaystyle\end{eqnarray}$$

(4.29) $$\begin{eqnarray}\displaystyle & \displaystyle B(t)=\left\langle \mathop{\sum }_{i=0}^{n_{t}-1}\unicode[STIX]{x1D70F}_{i}^{4}\exp (2\unicode[STIX]{x1D706}\text{i}\langle \unicode[STIX]{x1D70F}\rangle )\right\rangle . & \displaystyle\end{eqnarray}$$

In order to further develop these expressions, we consider heavy-tailed velocity distributions that behave as

(4.30) $$\begin{eqnarray}p_{s}(v)\propto v^{\unicode[STIX]{x1D6FD}-1},\end{eqnarray}$$

for $v<v_{0}$ with $v_{0}$ a characteristic velocity and $\unicode[STIX]{x1D6FD}>1$ . Such power-law type behaviour of the velocity PDF for small velocities is characteristic of observed transport behaviour in heterogeneous porous media (Berkowitz et al. Reference Berkowitz, Cortis, Dentz and Scher2006). From (4.30), the transition time PDF behaves as $\unicode[STIX]{x1D713}(t)\propto t^{-1-\unicode[STIX]{x1D6FD}}$ for large $t$ . Note that the moments $\langle \unicode[STIX]{x1D70F}^{k}\rangle$ of $\unicode[STIX]{x1D713}(t)$ are finite for $k<\lfloor \unicode[STIX]{x1D6FD}\rfloor$ , where $\lfloor \cdot \rfloor$ denotes the floor function. For illustration, here we consider the case $\unicode[STIX]{x1D6FD}>2$ , which implies $\langle \unicode[STIX]{x1D70F}^{2}\rangle <\infty$ .

Shear-dominated regime $(\unicode[STIX]{x1D70F}_{u}\ll t\ll \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}})$ . We now determine the stretching behaviour in the shear-dominated time regime $\unicode[STIX]{x1D70F}_{u}\ll t\ll \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}$ . In this time regime, we approximate (4.28) as

(4.31a,b ) $$\begin{eqnarray}A(t)=4\unicode[STIX]{x1D706}\left\langle \mathop{\sum }_{i=0}^{n_{t}-1}\text{i}\unicode[STIX]{x1D70F}_{i}^{4}\right\rangle ,\quad B(t)=\left\langle \mathop{\sum }_{i=0}^{n_{t}-1}\unicode[STIX]{x1D70F}_{i}^{4}\right\rangle .\end{eqnarray}$$

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.  $\unicode[STIX]{x1D70F}_{n}<t$ because $t_{n_{t}}=\sum _{i=0}^{n_{t}-1}\unicode[STIX]{x1D70F}_{i}<t$ . Thus, we can approximate

(4.32) $$\begin{eqnarray}A(t)\approx 4\unicode[STIX]{x1D706}\int _{0}^{t}\text{d}t^{\prime }{t^{\prime }}^{4}\unicode[STIX]{x1D713}(t^{\prime })\mathop{\sum }_{i=0}^{\langle n_{t}\rangle -1}\text{i}\propto t^{6-\unicode[STIX]{x1D6FD}},\end{eqnarray}$$

where we approximated $\langle n_{t}\rangle \approx t/\langle \unicode[STIX]{x1D70F}\rangle$ . Using the same approximations for $B(t)$ in (4.31) gives $B(t)\propto t^{5-\unicode[STIX]{x1D6FD}}$ . These scalings are consistent with the ones known for coupled continuous-time random walks (Dentz et al. Reference Dentz, Le Borgne, Lester and de Barros2015). Hence, in the regime $\unicode[STIX]{x1D70F}_{u}\ll t\ll \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}$ the transverse (3.45) deformation evolves linearly with time, $\langle \unicode[STIX]{x1D6EC}_{t}(t)^{2}\rangle \propto t$ , because $\unicode[STIX]{x1D60D}_{22}(t)$ and $\unicode[STIX]{x1D60D}_{33}(t)$ are approximately constant whilst $\langle \unicode[STIX]{x1D60D}_{23}(t)^{2}\rangle \propto t$ according to (4.26). The longitudinal (3.44) and material (3.50) deformations evolve nonlinearly to leading order as $\langle \unicode[STIX]{x1D6EC}_{l}(t)^{2}\rangle \propto \langle \ell (t)^{2}\rangle t^{6-\unicode[STIX]{x1D6FD}}$ according to (4.32).

Cross-over $(t\gtrsim \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}})$ . We now focus on the cross-over from the power-law to exponential stretching regimes. For $t\gtrsim \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}$ , we approximate $A(t)\approx 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)$

(4.33) $$\begin{eqnarray}B(t)=\int _{0}^{t}\text{d}t^{\prime }{t^{\prime }}^{4}\unicode[STIX]{x1D713}(t^{\prime })\mathop{\sum }_{i=0}^{\langle n_{t}\rangle -1}\exp (2\unicode[STIX]{x1D706}\text{i}\langle \unicode[STIX]{x1D70F}\rangle )\propto t^{4-\unicode[STIX]{x1D6FD}}[\exp (2\unicode[STIX]{x1D706}t)-1].\end{eqnarray}$$

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 $\unicode[STIX]{x1D706}\in (10^{-1},10^{-2},10^{-3})$ . The velocity gradients $\hat{\unicode[STIX]{x1D716}}_{22}$ , $\hat{\unicode[STIX]{x1D716}}_{12}$ are modelled as Gaussian variables with mean $\unicode[STIX]{x1D706}$ and zero respectively, the mean travel time $\langle \unicode[STIX]{x1D70F}\rangle =10^{-1}$ , and the variances $\unicode[STIX]{x1D70E}_{22}^{2}=\unicode[STIX]{x1D70E}_{12}^{2}=1$ . Here the velocity PDF is given by the Nakagami distribution

(4.34) $$\begin{eqnarray}p_{v}(v;k,\unicode[STIX]{x1D703})=\frac{2}{\unicode[STIX]{x1D6E4}(k)}\left(\frac{k}{\unicode[STIX]{x1D703}}\right)^{k}v^{2k-1}\text{e}^{-kv^{2}/\unicode[STIX]{x1D703}},\quad v>0,\end{eqnarray}$$

with the parameters $k=1$ , $\unicode[STIX]{x1D703}=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 $\unicode[STIX]{x1D70F}=\ell _{c}/v$ is then

(4.35) $$\begin{eqnarray}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D70F})=\frac{\ell _{c}}{\unicode[STIX]{x1D70F}^{2}}p_{v}(\ell _{c}/\unicode[STIX]{x1D70F}),\end{eqnarray}$$

hence $\unicode[STIX]{x1D713}(t)\propto t^{-1-\unicode[STIX]{x1D6FD}}$ with $\unicode[STIX]{x1D6FD}=2k$ for $t>\unicode[STIX]{x1D70F}_{u}$ . The transition from algebraic to exponential growth is clearly shown in figure 3(b), where some deviation from the predicted exponential growth is observed for small $\unicode[STIX]{x1D6FD}$ due to the heavy-tailed nature of $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D70F})$ .

Figure 3. Comparison between analytic predictions of $B(t)$ (4.29) for $t>\unicode[STIX]{x1D70F}_{u}$ (grey dashed) with numerical simulations from the synthetic coupled CTRW (4.24) (solid dark grey, solid blue online) for $\unicode[STIX]{x1D6FD}\in (1.5,2.5,3.5)$ (plots are offset for clarity with $\unicode[STIX]{x1D6FD}$ increasing from bottom to top) and (a) $\unicode[STIX]{x1D706}=10^{-3}$ , $\unicode[STIX]{x1D70E}_{22}=10^{-3}$ , (b) $\unicode[STIX]{x1D706}=10^{-2}$ , $\unicode[STIX]{x1D70E}_{22}=10^{-2}$ , (c) $\unicode[STIX]{x1D706}=10^{-1}$ , $\unicode[STIX]{x1D70E}_{22}=10^{-1}$ .

5.1 Evolution of finite-time Lyapunov exponents

From § 4 we may also approximate evolution of the finite-time Lyapunov exponent (FTLE) in chaotic flows as

(5.5) $$\begin{eqnarray}\unicode[STIX]{x1D708}(t)\equiv \frac{1}{2t}\ln \unicode[STIX]{x1D707}(t),\end{eqnarray}$$

where $\unicode[STIX]{x1D707}(t)$ is the leading eigenvalue of the Cauchy–Green tensor $\unicode[STIX]{x1D63E}^{\prime }(t)=\unicode[STIX]{x1D641}^{\prime }(t)^{\top }\boldsymbol{\cdot }\unicode[STIX]{x1D641}^{\prime }(t)$ . Under the approximations $\unicode[STIX]{x1D60D}_{11}^{\prime }(t)\sim 0$ , $\unicode[STIX]{x1D60D}_{33}^{\prime }(t)\sim 0$ , $\unicode[STIX]{x1D60D}_{22}^{\prime }(t)$ , $|\unicode[STIX]{x1D60D}_{23}^{\prime }(t)|\sim \exp (\unicode[STIX]{x1D706}t+\unicode[STIX]{x1D70E}_{22}\unicode[STIX]{x1D709}(t))$ , $|\unicode[STIX]{x1D60D}_{12}^{\prime }(t)|$ , $|\unicode[STIX]{x1D60D}_{13}^{\prime }(t)|\sim t^{2-\unicode[STIX]{x1D6FD}/2}\exp (\unicode[STIX]{x1D706}t+\unicode[STIX]{x1D70E}_{22}\unicode[STIX]{x1D709}(t))$ , where $\unicode[STIX]{x1D709}(t)$ is a Gaussian white noise with zero mean and unit variance, then to leading order

(5.6) $$\begin{eqnarray}\unicode[STIX]{x1D707}(t)\sim \exp (2\unicode[STIX]{x1D706}t+\unicode[STIX]{x1D70E}_{22}\unicode[STIX]{x1D709}(t))(1+t^{4-\unicode[STIX]{x1D6FD}}).\end{eqnarray}$$

As such, the average FTLE then evolves as

(5.7) $$\begin{eqnarray}\langle \unicode[STIX]{x1D708}(t)\rangle \approx \unicode[STIX]{x1D706}+\left(2-\frac{\unicode[STIX]{x1D6FD}}{2}\right)\frac{\ln t}{t},\end{eqnarray}$$

and so as expected the FTLE converges to the infinite-time Lyapunov exponent $\unicode[STIX]{x1D706}$ at rate given by the velocity PDF power-law exponent $\unicode[STIX]{x1D6FD}$ . We may also compute the variance of $\unicode[STIX]{x1D708}(t)$ as

(5.8) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D708}}^{2}(t)=\langle (\unicode[STIX]{x1D708}(t)-\langle \unicode[STIX]{x1D708}(t)\rangle )^{2}\rangle =\frac{\unicode[STIX]{x1D70E}_{22}^{2}}{t}.\end{eqnarray}$$

As expected, both of these expressions provide accurate estimates of the FTLE mean and variance for all but short times, as reflected in figures 9 and 10 respectively. Here the explicit role of algebraic stretching at times $t<\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D716}}$ is shown to govern convergence of the ensemble average of the FTLE toward its infinite-time value $\unicode[STIX]{x1D706}$ .

Figure 9. Comparison between analytic prediction of the mean FTLE $\langle \unicode[STIX]{x1D708}(t)\rangle$ (5.7) (grey, dashed), the infinite-time Lyapunov exponent $\unicode[STIX]{x1D706}$ (black, dashed) and $\langle \unicode[STIX]{x1D708}(t)\rangle$ computed directly (dark grey, blue online) for the variable-helicity flow with (a) $\unicode[STIX]{x1D713}=\unicode[STIX]{x03C0}/2$ , (b) $\unicode[STIX]{x1D713}=\unicode[STIX]{x03C0}/64$ .

Figure 10. Comparison between analytic prediction of the FTLE variance $\langle \unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D708}}^{2}(t)\rangle$ (5.8) (grey, dashed) and $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D708}}^{2}(t)$ computed directly (dark grey, blue online) for the variable-helicity flow with (a) $\unicode[STIX]{x1D713}=\unicode[STIX]{x03C0}/2$ , (b) $\unicode[STIX]{x1D713}=\unicode[STIX]{x03C0}/64$ .