2.1. Model problem

We consider a 1-D sample of soft porous material of relaxed length $L$ and relaxed porosity (fluid fraction) $\phi _{f,0}$ . The left boundary of the material (at $x=a(t))$ is moving and permeable, whereas the right boundary (at $x=L$ ) is fixed and impermeable. The position of the left boundary, $a(t)$ , is imposed to be

(2.1) \begin{equation} a(t)= \frac {A}{2} \left [1-\cos \left (\frac {2\pi t}{T}\right )\right ], \end{equation}

where $A$ and $T$ are the amplitude and period of loading, respectively. We consider imposed deformations ranging from small to large macroscopic strains ( $-0.4\,\%$ to $-20\,\%$ or $0.004\leq {}A/L\leq {}0.2$ ). We take the fluid and solid to be individually incompressible, such that changes in bulk volume correspond directly to the movement of fluid into and out of the pore space. We presented and analysed the poromechanics of this scenario in detail in a companion study (Fiori et al. Reference Fiori, Pramanik and MacMinn2023). We now introduce a strip of passive solute of initial width $l$ located at the right boundary and we study the impact of this periodic, displacement-driven deformation on the evolution of the solute distribution (figure 1).

Figure 1. We consider a 1-D sample of soft porous material of relaxed length $L$ , subject to a periodic, displacement-driven loading at its left boundary (white arrows). The left boundary is permeable, thus allowing fluid flow in or out (pale blue squiggles) to accommodate the loading. The right boundary is fixed and impermeable. The solute is initially localised against the right boundary in a strip of width $l$ (dark blue).

2.2. Kinematics

We consider an Eulerian reference frame, in which the solid displacement is $\mathbf {u}_{\mathbf {s}}=\mathbf {x}-\mathbf {X}(\mathbf {x},t)$ , with $\mathbf {X}(\mathbf {x},t)$ the reference position of the solid material point that at time $t$ occupies position $\mathbf {x}$ . We choose our reference configuration to be the relaxed configuration, such that $\mathbf {X}(\mathbf {x},0)=\mathbf {x}$ and $\mathbf {u}_{\mathbf {s}}(\mathbf {x},0)=0$ . The true volume fractions of fluid and solid are $\phi _f$ and $\phi _{s}$ , respectively, where $\phi _{f}+\phi _{s}=1$ . In this uniaxial setting, the solid displacement and the solid and fluid velocities are one-dimensional and given by

(2.2) \begin{equation} \mathbf{u}_{\mathbf{s}}=u_{s}(x,t) {\hat {\mathbf{e}}_{\mathbf {x}}}, \; \mathbf {v}_{\mathbf {s}} = v_{s}(x,t) {\hat {\mathbf{e}}_{\mathbf {x}}}, \; \mathbf {v}_{\mathbf {f}} = v_{f}(x,t) {\hat{\mathbf{e}}_{\mathbf {x}}}, \end{equation}

where $\mathbf {v}_{\mathbf {f}}$ and $\mathbf {v}_{\mathbf {s}}$ are the fluid and solid velocities, respectively, $u_s$ , $v_s$ and $v_f$ are the $x$ -components of these fields and ${\hat {\mathbf e}_{\mathbf {x}}}$ is the unit vector in the $x$ -direction. The local current volume per unit reference volume is measured by the Jacobian determinant, which in this uniaxial setting is given by $J= (1-\partial {u_s}/\partial {x})^{-1}$ . For incompressible constituents and uniform initial porosity $\phi _{f,0}$ , the local change in volume relates to the change in porosity as

(2.3) \begin{equation} J(x,t) = \frac {1-\phi _{f,0}}{1-\phi _f} \quad \to \quad \frac {\partial {u_s}}{\partial {x}} =\frac {\phi _f-\phi _{f,0}}{1-\phi _{f,0}}. \end{equation}

Continuity can be written

(2.4) \begin{equation} \frac {\partial {\phi _f}}{\partial {t}} +\frac {\partial }{\partial {x}}(\phi _f v_f) = 0 \quad \mathrm {and}\quad \frac {\partial {\phi _f}}{\partial {t}} -\frac {\partial }{\partial {x}}{[(1-\phi _f)v_s]} = 0, \end{equation}

which together imply that the total flux $q=\phi _fv_f+(1-\phi _f)v_s$ is uniform in space, $\partial {q}/\partial {x}=0$ .

2.3. Fluid flow

We assume that the fluid flows relative to the solid according to Darcy’s law

(2.5) \begin{equation} \phi _ f (v_{f}-v_{s}) = - \frac {k(\phi _f)}{\mu }\frac {\partial p}{\partial x}, \end{equation}

where $k(\phi _f)$ is the permeability of the solid skeleton, $\mu$ is the dynamic viscosity of the fluid and $p$ is the fluid (pore) pressure, and where we have neglected gravity. As in Fiori et al. (Reference Fiori, Pramanik and MacMinn2023), we take the permeability to be deformation-dependent according to a normalised Kozeny–Carman relation, $k(\phi _{f}) = k_0 ({(1-\phi _{f,0})^2}/{\phi _{f,0}^3}) ({\phi _{f}^3}/{(1-\phi _{f})^2})$ , where $k_0\equiv {}k(\phi _{f,0})$ is the permeability of the initial state. We discuss this choice in detail in Fiori et al. (Reference Fiori, Pramanik and MacMinn2023).

Combining (2.4) and (2.5), we arrive at the nonlinear flow equations

(2.6) \begin{equation} \frac {\partial {\phi _f}}{\partial {t}} +\frac {\partial }{\partial {x}}\bigg [{\phi _f q}-(1-\phi _f) \frac {k(\phi _f)}{\mu }\frac {\partial {p}}{\partial {x}}\bigg ]=0 \quad \mathrm {and}\quad \frac {\partial {q}}{\partial {x}}=0, \end{equation}

where the total flux $q$ is again

(2.7) \begin{equation} q\equiv \phi _f v_f + (1-\phi _f)v_s, \end{equation}

and the fluid and solid velocities are given by

(2.8) \begin{equation} v_f=q-\frac {(1-\phi _f)}{\phi _f}\frac {k(\phi _f)}{\mu }\frac {\partial p}{\partial x} \quad \mathrm {and}\quad v_s=q+\frac {k(\phi _f)}{\mu }\frac {\partial p}{\partial x}. \end{equation}

Note that the fluid flux is

(2.9) \begin{equation} q_f=\phi _f v_f. \end{equation}

2.4. Mechanical equilibrium and elasticity law

Neglecting inertia, gravity and other body forces, mechanical equilibrium can be expressed as ${\boldsymbol {\nabla }\cdot {\boldsymbol {\sigma }}}= {\boldsymbol {\nabla }\cdot {\boldsymbol {\sigma '}}}-\boldsymbol {\nabla } p= 0$ , where $\boldsymbol {\sigma }$ is the true Cauchy total stress, decomposed into contributions from the fluid pressure $p$ and from Terzaghi’s effective stress $\boldsymbol {\sigma '}$ . In one dimension, mechanical equilibrium reads

(2.10) \begin{equation} \frac {\partial \sigma '}{\partial x}=\frac {\partial p}{\partial x}, \end{equation}

where $\sigma ^\prime$ is the $xx$ component of $\boldsymbol {\sigma ^\prime}$ .

We take the solid skeleton to be elastic, with no viscous or dissipative behaviours. Since any elasticity law can be written in the form $\sigma ^\prime =\sigma ^\prime (\phi _f)$ for a uniaxial deformation, this problem can be described by a nonlinear advection-diffusion equation

(2.11) \begin{equation} \frac {\partial {\phi _f}}{\partial {t}} +\frac {\partial }{\partial {x}}\bigg [{\phi _f q}-D_f(\phi _f)\frac {\partial {\phi _f}}{\partial {x}}\bigg ]=0 \quad \mathrm {and}\quad \frac {\partial {q}}{\partial {x}}=0, \end{equation}

where the nonlinear composite constitutive function

(2.12) \begin{equation} D_f(\phi _f)=(1-\phi _f)\frac {k(\phi _f)}{\mu }\frac {\mathrm {d}\sigma ^\prime }{\mathrm {d}\phi _f} \end{equation}

is the poroelastic diffusivity. Note that a very similar model is used for the solidification of colloidal suspensions in applications such as filtration and sedimentation, for which the poroelastic diffusivity $D_f(\phi _f)$ (i.e. the ‘solids diffusivity’) is characterised as a composite material property (e.g. Reference Davis and RusselDavis and Russel 1989; Peppin, Elliott & Worster Reference Peppin, Elliott and Worster2006; Style & Peppin Reference Style and Peppin2011; Bouchaudy and Salmon Reference Bouchaudy and Salmon2019; Worster, Peppin & Wettlaufer Reference Worster, Peppin and Wettlaufer2021).

We use Hencky hyperelasticity (Hencky Reference Hencky1931) as a simple, large-deformation model that captures kinematic nonlinearity. For a uniaxial deformation, the relevant component of the effective stress is then

(2.13) \begin{equation} \sigma ^\prime = \mathcal {M}\frac {\ln (J)}{J}=\mathcal {M} \left (\frac {1-\phi _f}{1-\phi _{f,0}}\right ) \ln \left (\frac {1-\phi _{f,0}}{1-\phi _f}\right ), \end{equation}

where $\mathcal {M}$ is the $p$ -wave or oedometric modulus (MacMinn, Dufresne & Wettlaufer Reference MacMinn, Dufresne and Wettlaufer2016). Note that, for these constitutive choices of Kozeny–Carman permeability and Hencky elasticity, the permeability at the left boundary can vanish for sufficiently large $A$ and/or small $T$ , because the poroelastic diffusivity remains finite rather than diverging as $\phi _f\to 0$ (Hewitt et al. Reference Hewitt, Paterson, Balmforth and Martinez2016). We motivate and discuss our constitutive choices and explore the poroelastic diffusivity in more detail in Fiori et al. (Reference Fiori, Pramanik and MacMinn2023).

With appropriate initial conditions, boundary conditions, and the normalised Kozeny–Carman permeability law, equations (2.11), (2.12) and (2.13) comprise a closed model for the evolution of the porosity.

2.5. Solute transport

We now consider the transport of solute. We denote the true local solute concentration in the fluid phase by $c$ (amount of solute per unit current fluid volume). We take the solute to be passive and charge neutral, with no chemical or other interaction with the solid or fluid phases, so that neither the fluid properties nor the solid properties depend on $c$ . The flow and mechanics above are then independent of the transport problem.

In one dimension, it is well known that conservation of mass for a passive solute can be written

(2.14) \begin{equation} \frac {\partial }{\partial t} (\phi _f c) +\frac {\partial }{\partial x}\left [{\phi _f c v_f }- \phi _f \mathcal {D} \frac {\partial c}{\partial x}\right ]=0. \end{equation}

The first term in the square brackets is the Darcy-scale solute flux due to advection, which occurs here entirely in response to the deformation. The second term in the square brackets combines molecular diffusion and hydrodynamic dispersion, thus taking the latter to be a Fickian process (e.g. Scheidegger Reference Scheidegger1961). The latter term is multiplied by the porosity $\phi _f$ since solute movements only occur in the fluid phase. The coefficient $\mathcal {D}$ can be written

(2.15) \begin{equation} \mathcal {D} =\mathcal {D}_m+\mathcal {D}_h, \end{equation}

where $\mathcal {D}_m$ and $\mathcal {D}_h$ are the coefficients of molecular diffusion and hydrodynamic dispersion, respectively. Dispersion, in which pore-scale velocity gradients and the tortuosity of the pore space lead to macroscopic spreading of solute, depends sensitively on flow conditions and the details of the pore structure in ways that are not yet fully understood, even for rigid porous materials (Dentz et al. Reference Urciuolo, Imparato and Netti2018, Reference DiDomenico, Wang and Bonassar2023). The most widely used model for the macroscopic dispersive flux is Fickian, as above, with a velocity-dependent dispersion coefficient given in one dimension by

(2.16) \begin{equation} \mathcal {D}_h=\alpha |v_f-v_s|, \end{equation}

where $\alpha$ is the longitudinal dispersivity (Scheidegger Reference Scheidegger1961; Brenner & Edwards (Reference Brenner and Edwards1993; Gelhar Reference Gelhar1993; Whitaker Reference Whitaker1998). Note that the dispersive flux is proportional to $|v_f-v_s|$ , unlike the advective flux, because dispersion is driven by flow of fluid through the pore structure (i.e. $v_f=v_s\neq 0$ would lead to advection but no dispersion). Note also that, unlike the advective flux, the diffusive and dispersive fluxes are independent of the direction of the fluid flow.

The dispersivity $\alpha$ is typically taken to be a constant material property for a given pore structure. In a deforming porous material, and particularly for moderate to large deformations, it is likely that $\alpha$ should be deformation-dependent to account for the evolving pore structure. For example, particle–particle interactions and rearrangements are known to drive enhanced dispersion in dense suspensions (Souzy et al. Reference Vaughan, Galie, Stegemann and Grotberg2016, Reference Vuong, Yoshihara and Wall2017) and compaction has been shown to have a non-trivial impact on dispersion in bead packs and packed beds (Charlaix, Hulin & Plona Reference Charlaix, Hulin and Plona1987; Östergren & Trägårdh Reference Östergren and Trägårdh2000; Liu et al. Reference Liu, Gong, Xiao and Wang2024). We expect similar but even larger effects in poroelastic materials under large deformations, which may ultimately require novel dispersion models, but these phenomena are beyond the scope of the present study. Here, we take $\alpha$ to be a constant for simplicity.

2.6. Initial and boundary conditions

We next specify initial and boundary conditions for the solid, the fluid and the solute. Recall that the left and right boundaries of the solid are at $x=a(t)$ and $x=L$ , respectively.

2.6.1. Initial conditions

Equation (2.1) implies that $a(0)=0$ , and thus that the initial porosity is uniform and equal to the relaxed porosity

(2.17) \begin{equation} \phi _f(x,0)= \phi _{f,0} \textrm { and } u_s(x,0)=0. \end{equation}

We take the solute to be initially localised against the right boundary in a strip of width $l$ and concentration $c_0$ , such that

(2.18) \begin{equation} c(x,0)=\frac {c_0}{2} \left \{ \tanh {[s(x-L+l)}]+1 \right \}, \end{equation}

where $s$ is a steepness parameter.

2.6.2. Left boundary

For $t\gt 0$ , we apply a displacement-controlled loading at the left boundary according to equation (2.1). We take this moving boundary to be fluid and solute permeable. The associated boundary conditions are

(2.19) \begin{equation} u_s(a,t)=a(t), \quad v_s(a,t)=\frac {\mathrm {d}a}{\mathrm {d}t}\quad \mathrm {and}\quad p(a,t)=0. \end{equation}

We take the fluid outside the domain to be ‘clean’, such that

(2.20) \begin{equation} c(a,t)=0. \end{equation}

2.6.3. Right boundary

We take the right boundary to be fixed and impermeable, such that

(2.21) \begin{equation} u_s(L,t)=v_s(L,t)=v_f(L,t)=0 \quad \mathrm {and} \quad \frac {\partial c}{\partial x}\bigg |_{x=L}=0. \end{equation}

Equation (2.21) and the requirement that $q$ be uniform in space imply that there can be no net flow from left to right in our problem, $q\equiv 0$ . Equation (2.7) then implies that the fluid and the solid always locally move in opposite directions

(2.22) \begin{equation} v_f= -\frac {(1-\phi _f)}{\phi _f} v_s. \end{equation}

2.7. Scaling and summary

As in Fiori et al. (Reference Fiori, Pramanik and MacMinn2023), we apply the following non-dimensionalisation to the poromechanical model:

(2.23) \begin{equation} \begin{aligned} \tilde {x}=\frac {x}{L},\; \tilde {u}_s=\frac {u_s}{L},\; \tilde {t}=\frac {t}{T_{pe}},\; \tilde {\sigma }^\prime=\frac {\sigma '}{\mathcal {M}},\; \tilde {p}=\frac {p}{\mathcal {M}},\; \tilde {k}=\frac {k(\phi )}{k_0},\; \tilde {v}_f=\frac {v_f}{L/T_{pe}}, \; \tilde {v}_s=\frac {v_s}{L/T_{pe}}, \end{aligned} \end{equation}

where $T_{{pe}}=L^2/D_{f,0}=\mu {}L^2/(k_0\mathcal {M})$ is the classical poroelastic time scale for the relaxation of pressure over a distance $L$ and $D_{f,0}=k_0\mathcal {M}/\mu$ is the constant linear-poroelastic diffusivity.

We then scale quantities related to solute transport as

(2.24) \begin{equation} \tilde {c}=\frac {c}{c_0}, \; \tilde {l}=\frac {l}{L}, \; \tilde {\alpha }=\frac {\alpha }{L}. \end{equation}

Taking $q\equiv {}0$ , as noted above, the full problem can then be rewritten in dimensionless form as

(2.25) \begin{equation} \frac {\partial {\phi _f}}{\partial {\tilde {t}}} -\frac {\partial }{\partial {\tilde {x}}}\bigg [\tilde {D}_f(\phi _f)\frac {\partial {\phi _f}}{\partial {\tilde {x}}}\bigg ]=0, \end{equation}

where

(2.26) \begin{equation} \tilde {D}_f=\frac {D_f}{D_{f,0}}=(1-\phi _f)\tilde {k}(\phi _f)\frac {\mathrm {d}\tilde {\sigma }^\prime }{\mathrm {d}\phi _f}, \end{equation}

and

(2.27) \begin{equation} \frac {\partial }{\partial \tilde {t}}(\phi _f \tilde {c}) +\frac {\partial }{\partial \tilde {x}}\bigg [{\phi _f \tilde {c} \tilde {v}_f }- \phi _f \tilde {\mathcal {D}} \frac {\partial \tilde {c}}{\partial \tilde {x}}\bigg ]=0. \end{equation}

The dimensionless coefficient of diffusion/dispersion $\tilde {\mathcal {D}}$ is

(2.28) \begin{equation} \tilde {\mathcal {D}}=\frac {\mathcal {D}}{\mathcal {D}_m}=\mathrm {Pe}^{-1}+\tilde {\alpha } | \tilde {v}_f-\tilde {v}_s|, \end{equation}

where $\mathrm {Pe}=({L^2/T_{pe}})/{\mathcal {D}_m}= {k_0 \mathcal {M}}/({\mu \mathcal {D}_m})$ is the Péclet number, which measures the importance of poroelastic-relaxation-driven advection relative to molecular diffusion.

The initial conditions are

(2.29) \begin{equation} \tilde {a}(0)=0, \; \phi _f(\tilde {x},0)= \phi _{f,0}, \end{equation}

and

(2.30) \begin{equation} \tilde {c}(\tilde {x},0)=\frac {1}{2} \{ \tanh {[\tilde {s}(\tilde {x}-1+\tilde {l})}]+1\}, \end{equation}

where we take $\tilde {s}=sL=60$ . The boundary conditions are

(2.31) \begin{equation} \tilde {u}_s(\tilde {a},\tilde {t})=\tilde {a}(\tilde {t})= \frac {\tilde {A}}{2} \bigg [1-\cos \left (\frac {2\pi \tilde {t}}{\tilde {T}}\right )\bigg ] \,,\,\,\tilde {v}_s(\tilde {a},\tilde {t})=\frac {\mathrm {d} \tilde {a}}{\mathrm {d}\tilde {t}}\,,\,\ \quad \tilde {p}(\tilde {a},\tilde {t})=0 \; \mathrm {and} \quad \tilde {c}(\tilde {a},\tilde {t})=0 ,\end{equation}

and

(2.32) \begin{equation} \tilde {u}_s(1,\tilde {t}) =\tilde {v}_s(1,\tilde {t}) =\tilde {v}_f(1,\tilde {t})=0 \quad \mathrm {and}\quad \frac {\partial \tilde {c}}{\partial \tilde {x}}\bigg |_{\tilde {x}=1}=0, \end{equation}

where $\tilde {A}=A/L$ and $\tilde {T}=T/T_{{pe}}$ .

As shown in Fiori et al. (Reference Fiori, Pramanik and MacMinn2023), the forcing considered here will drive a typical solid velocity of size $v_s^*=2A/T$ and thus a typical fluid velocity of size $v_f^*= ( ({1-\phi _{f,0}})/{\phi _{f,0}} )(2A/T)$ . The characteristic advection time $T_{ {adv}}$ and diffusion time $T_{{diff}}$ are then

(2.33) \begin{equation} T_{ {adv}}=\frac {L}{v_f^*} =\frac {LT\phi _{f,0}}{2A(1-\phi _{f,0})} \quad \to \quad \tilde {T}_{ {adv}} =\frac {T_{ {adv}}}{ T_{{pe}}} =\frac {\tilde {T}\phi _{f,0}}{2\tilde {A}(1-\phi _{f,0})} \propto \frac {\tilde {T}}{\tilde {A}}, \end{equation}

and

(2.34) \begin{equation} T_{{diff}}=\frac {L^2}{\mathcal {D}_m} \quad \to \quad \tilde {T}_{{diff}}=\frac {{T}_{{diff}}}{T_{{pe}}} = \frac {D_{f,0}}{\mathcal {D}_m}=\mathrm {Pe} .\end{equation}

Recall that the Péclet number – as defined above – quantifies the rate of advection due to poromechanical relaxation relative to the rate of molecular diffusion. The characteristic times above suggest that the balance between loading-driven advection and molecular diffusion is better measured by an effective Péclet number $\mathrm {Pe}_{{eff}}$ ,

(2.35) \begin{equation} \mathrm {Pe}_{{eff}} =\mathrm {Pe}\,\frac {\tilde {A}}{\tilde {T}} \propto \frac {\tilde {T}_{{diff}}}{\tilde {T}_{ {adv}}}. \end{equation}

In our results below, we explore a wide range of $\mathrm {Pe}_{{\it eff}}$ : from $\sim$ 1 to $\sim 10^5$ . We show in table 2 that this range is biologically relevant.

The above model describes uniaxial flow, mechanics and solute transport in a poroelastic material subject to periodic deformations. The kinematics are rigorous and thus nonlinear, the elasticity law is Hencky elasticity and the permeability law is the normalised Kozeny–Carman formula. Solute transport occurs via advection, molecular diffusion and hydrodynamic dispersion. We solve this system numerically in MATLAB using compact finite differences in space and an implicit Runge–Kutta method in time, as described in more detail in Appendix A. We provide an example code in Fiori, Pramanik & MacMinn (Reference Fiori, Pramanik and MacMinn2025). Below, we consider only dimensionless quantities, dropping the tildes for convenience.

3.1. Quantification of solute transport and mixing

Figure 2. Schematic representation of the travel distance or mixing length $\delta$ , which measures the distance travelled by the left edge of the concentration profile during the time $t$ . For solute initially localised in a finite strip at the right, we calculate $\delta (t)$ by choosing a small threshold concentration $c_{\delta }$ and then finding the leftmost position $x_{\delta }(t)$ where that concentration occurs. Then, $\delta (t) = |x_{\delta }(t) - x_{\delta }(0) |$ (see, e.g. Tan & Homsy Reference Tan and Homsy1988; Mishra, Martin & De Wit Reference Mishra, Martin and De Wit2008). Here, we show $c(x,0)$ (dashed curve), $c(x,t)$ (solid curve) and the corresponding $\delta (t)$ . The value of $c_\delta$ is arbitrary and should have no qualitative impact on the results. In the results shown below, we take $c_{\delta }= 0.01$ .

We begin with some qualitative examples that illustrate the impact of deformation on each transport mechanism individually. We also assess solute transport and mixing quantitatively via two metrics:

1. (i) The travel distance or mixing length $\delta$ measures the distance travelled by the left edge of the concentration profile (figure 2). The travel distance can range from $0$ to $1-l$ , but it becomes less meaningful as it approaches $1-l-A$ , by which point the concentration profile interacts strongly with the left boundary.

2. (ii) The degree of mixing $\chi$ measures the degree to which the initial concentration profile has homogenised, and is closely related to the variance of the concentration distribution. We express the degree of mixing in terms of the variance of the concentration distribution by generalising the standard definition (see, e.g. Danckwerts Reference Danckwerts1952; Jha, Cueto-Felgueroso & Juanes Reference Jha, Cueto-Felgueroso and Juanes2011) to account for a porosity field that varies in space. Considering the fluid-volume-weighted average $\langle \ast \rangle _{f}$ , defined as

   (3.1) \begin{align} \langle \ast \rangle _f=\frac {\int _a^1\,\phi _f\ast \,\mathrm {d}x}{\int _a^1\,\phi _f\,\mathrm {d}x}, \end{align}
   the variance of the concentration distribution is then
   (3.2) \begin{equation} \sigma ^2(t) = \langle {}c^2\rangle {}_{f} - \langle {}c\rangle {}_{f}^2 ,\end{equation}
   and the degree of mixing is
   (3.3) \begin{equation} \chi (t) = 1 - \frac {\sigma ^2(t)}{\sigma _{\textrm { max}}^2}, \end{equation}
   where $\sigma _ {{max}}^2=\sigma ^2(t=0)$ in this case. Note that $\chi$ can range from $0$ to $1$ , where the former corresponds to no mixing (i.e. the initial state by definition) and the latter is characteristic of a completely mixed configuration (i.e. spatially uniform concentration).

3.2. Baseline values

For a given total loading time, $\delta$ and $\chi$ depend on the transport parameters $\mathrm {Pe}^{-1}$ and $\alpha$ ; the loading parameters $A$ and $T$ ; the initial porosity $\phi _{f,0}$ ; and the initial width of the solute strip $l$ . We choose a baseline value for each parameter (table 1). We use these baseline values in all of the results presented below, except where explicitly noted otherwise. We explore the impact of individually changing $\mathrm {Pe}^{-1}$ and $\alpha$ in § 3.4, $A$ and $T$ in § 3.5 and $\phi _{f,0}$ and $l$ in Appendix D.

Table 1. Baseline parameter values.

We choose a baseline amplitude $A=0.1$ , corresponding to moderately large deformations. We choose a baseline period $T=6\pi$ , which, following our companion study (Fiori et al. Reference Fiori, Pramanik and MacMinn2023), ensures that the poromechanics are quasi-static (i.e. ‘slow loading’ see the first part of § 3.5). The fluid flux $q_f$ and the relative velocity $|v_f-v_s|$ for this baseline case are shown in figures 7(d) and 7(h), respectively. The baseline values of $\mathrm {Pe}^{-1}$ and $\alpha$ are in the range of those proposed by Sengers et al. (Reference Sengers, Oomens and Baaijens2004) for cartilage constructs, with the specific values chosen to ensure that diffusion dominates over dispersion for the slowest period considered in this study (see Appendix C). The baseline value for $\phi _{f,0}$ is representative of hydrogels or soft biological tissues, whereas $l$ is arbitrarily chosen to be a small fraction of the domain length.

3.3. Qualitative impacts of periodic loading on solute transport

Figure 3. Evolution of the solute flux across $x=1-l$ during 5 loading cycles. We show the total flux of solute (solid black) and the separate contributions of advection (dotted blue), molecular diffusion (dash-dotted green) and hydrodynamic dispersion (dashed red) for $A=0.4, \alpha =0.025$ . Note that $A$ and $\alpha$ are higher than the baseline values to better illustrate the roles of advection and dispersion. The solid grey envelope is proportional to $t^{-\frac {1}{2}}$ .

We begin by isolating and comparing the solute transport mechanisms. To illustrate the contribution of each mechanism, we consider the time evolution of their separate contributions to the total solute flux at $x=1-l$ , which is the initial left edge of the concentration profile, during five loading cycles (figure 3). The individual contribution of the advective, diffusive and dispersive solute fluxes are $q_{ {adv}}=\phi _f v_f c$ , $q_{{diff}}=\phi _f \mathrm {Pe}^{-1} (\partial c/\partial x)$ and $q_{ {disp}}=\phi _f \alpha |v_f-v_s| (\partial c/\partial x)$ , respectively. During the loading half of each cycle ( $\dot {a}\gt 0$ ), all three fluxes are negative, implying that all three mechanisms drive solute to the left. During the unloading half of each cycle ( $\dot {a}\lt 0$ ), however, the flow changes direction and the advective flux changes sign (now positive, meaning to the right), whereas the diffusive and dispersive fluxes remain negative (still to the left). The flow and deformation are periodic after an initial transient that decays exponentially (see Fiori et al. Reference Fiori, Pramanik and MacMinn2023), in which case the net contribution of advection over one full cycle is zero (see figure 4 b). Thus, net transport at the end of each cycle depends on the cumulative amount of diffusion and dispersion. Diffusion and dispersion are strongest at early times, when the concentration gradient is largest, and decay over time as $t^{-1/2}$ . The strengths of diffusion and dispersion are proportional to $\mathrm {Pe}^{-1}$ and $\alpha$ , respectively.

Figure 4. Evolution of the concentration profile during one cycle (red to blue through white) for four cases: (a) diffusion only ( $A=\alpha =0, \mathrm {Pe}^{-1}=3\times 10^{-5}$ ); (b) advection only ( $A=0.4, \mathrm {Pe}^{-1}=\alpha =0$ ); (c) advection and diffusion ( $A=0.4, \mathrm {Pe}^{-1}=3\times 10^{-5}, \alpha =0$ ); (d) advection, diffusion and dispersion ( $A=0.4, \mathrm {Pe}^{-1}=3\times 10^{-5}, \alpha =0.025$ ). We plot concentration against the spatial coordinate $x$ and split the evolution into two phases, loading ( $\dot {a}\gt 0$ , first half of the cycle, dark to light red) and unloading ( $\dot {a}\lt 0$ , second half, light to dark blue). In panel (b), the unloading curves (dashed) overlap with the loading curves (solid). The initial profile is shown in black. For each case, we also show the evolution of $\delta$ throughout the loading cycle (insets); in all cases, the dotted curves are for diffusion without loading (with the dashed reference line showing linearity with $\sqrt {t/T}$ ), the dash-dot curves are for advection only, the thin solid curves are for advection and diffusion and the thick solid curve is for advection, diffusion and dispersion. Note that $A$ and $\alpha$ are higher than the baseline values to better illustrate the roles of advection and dispersion.

We next plot the evolution of the concentration profile during the first cycle (figure 4). We consider four cases: molecular diffusion only, in which $A=0$ (no loading); advection only, in which $\mathrm {Pe}^{-1}=\alpha =0$ ; advection and molecular diffusion only, in which $\alpha =0$ ; and the general case, including all three mechanisms. For diffusion only (figure 4 a) solute spreading is driven exclusively by concentration gradients and the travel distance $\delta$ grows as $\delta \propto \sqrt {t}$ after an initial (slower) phase in which the profile adjusts from its initial condition toward classical self similarity (see Appendix B). When a deformation is applied (figure 4b–d), four main factors impact the movement of the solute: (i) the motion of the fluid drives advection; (ii) the motion of the fluid through the pore space drives dispersion; (iii) the decrease in porosity weakly hinders diffusion and dispersion since $\phi _f$ is a prefactor in both of those fluxes; and (iv) the stretched solute profile, with the same quantity of fluid (and solute) now occupying a larger spatial extent, weakens concentration gradients. The latter two effects become increasingly strong during loading and then decreasingly strong during unloading. The fourth mechanism is most obvious in the case of advection only (figure 4 b), where the motion of the solute is perfectly reversible and $\delta =0$ at the end of the loading cycle. The fact that loading weakly suppresses molecular diffusion via the third and fourth mechanisms is apparent in figure 4 c, where the final value of $\delta$ is lower for diffusion with loading than for diffusion without loading. When dispersion is included (figure 4 d), transport is greatly amplified despite the weak suppression of diffusion.

We next consider several quantitative measures of transport and mixing.

3.4. Impact of diffusion and dispersion coefficients on transport and mixing

Figure 5. Impact of $\mathrm {Pe}^{-1}$ on the evolution of $\delta$ over 5 loading cycles. (a) We plot the evolution of $\delta$ with $\sqrt {t}$ for nine different values of $\mathrm {Pe}^{-1}$ $\in [3 \times 10^{-8},3 \times 10^{-4}]$ (dark to light). Note that the curves for the two smallest values of $\mathrm {Pe}^{-1}$ overlap. In each case, delta is roughly linear in $\sqrt {t}$ with a slope that increases monotonically with $\mathrm {Pe}^{-1}$ . The dashed curve indicates the baseline value of $\mathrm {Pe}^{-1}$ . (b) We plot the final value of $\delta$ at $t=5T$ as function of $\mathrm {Pe}^{-1}$ . The open circle indicates the baseline value of $\mathrm {Pe}^{-1}$ . The black dashed curve is our estimate $\delta _{{est, SL}}$ from equation (3.4).

Figure 6. Impact of $\alpha$ on the evolution of $\delta$ over 5 loading cycles. (a) We plot the evolution of $\delta$ with $\sqrt {t}$ for nine different values of $\alpha$ $\in [ 10^{-5},10^{-1}]$ (dark to light). Note that the curves for the two smallest values of $\alpha$ overlap. In each case, delta is roughly linear in $\sqrt {t}$ with a slope that increases monotonically with $\alpha$ . The dashed curve indicates the baseline value of $\alpha$ . (b) We plot the final value of $\delta$ at $t=5T$ as function of $\alpha$ . The open circle indicates the baseline value of $\alpha$ . The black dashed curve is our estimate $\delta _{{est, SL}}$ from equation (3.4).

We next isolate the roles of $\mathrm {Pe}^{-1}$ and $\alpha$ . For that purpose, we focus on five cycles and consider a wide range of $\mathrm {Pe}^{-1}$ and $\alpha$ . We consider the roles of $\phi _{f,0}$ and $l$ in Appendix D. The impact of changing $\mathrm {Pe}^{-1}$ and $\alpha$ on $\delta$ are shown in figures 5 and 6, respectively. All of the curves are roughly linear in $\sqrt {t}$ after an initial transient (i.e. spreading is Fickian on average), exhibiting fluctuations with a period $T$ because of the phenomena described in § 3.3: loading decreases the porosity, forcing the solute to spread (advection to the left), and unloading reverses this process. Larger values of $\mathrm {Pe}^{-1}$ and $\alpha$ enhance the diffusive and dispersive fluxes, respectively, and hence drive faster spreading, as should be expected. For sufficiently small $\mathrm {Pe}^{-1}$ , dispersion dominates diffusion and the rate of spreading becomes independent of $\mathrm {Pe}^{-1}$ (figure 5). Similarly, diffusion dominates dispersion for sufficiently small $\alpha$ and the rate of spreading becomes independent of $\alpha$ (figure 6).

We next introduce a naive estimate for the travel distance during slow loading, $\delta _{{est, SL}}$ , based on the assumption that net transport is Fickian on average. That is, we assume that $\delta _{ {est,SL}} \propto \sqrt {\mathcal {D}_{ {eff}}t}$ for some effective diffusion/dispersion coefficient $\mathcal {D}_{{eff}}\approx {}\mathrm {Pe}^{-1}+\alpha |v_f-v_s|$ . During slow loading, the quantity $|v_f-v_s|$ is proportional to $A/T$ and decreases linearly from left to right (see Fiori et al. Reference Fiori, Pramanik and MacMinn2023). Thus, we assume that, on average, $|v_f-v_s|\sim {}C_2A/T$ for some constant $C_2$ . The resulting estimate is then

(3.4) \begin{equation} \delta _{{est, SL}} =C_1 f\left (\frac {\overline {\langle {\phi _f}\rangle }}{\phi _{f,0}},\mathrm {Pe}^{-1}\right ) \sqrt { 4 \left (\mathrm {Pe}^{-1} + C_2 \alpha \frac {A}{T}\right ) t}, \end{equation}

where $C_1$ is a constant and the function $f(\overline {\langle {\phi _f}\rangle }/\phi _{f,0},\mathrm {Pe}^{-1})$ is an empirical prefactor to capture the impact of the average deformation on diffusive spreading, as discussed in more detail below (see figure 10), with $\overline {\langle {\phi _f}\rangle }$ the overall average porosity (see equation 3.6). By fitting equation (3.4) to the results shown in figures 5, 6 and 10, we find $C_1=1.29$ , $C_2=0.25$ and

(3.5) \begin{equation} f\left (\frac {\overline {\langle {\phi _f}\rangle }}{\phi _{f,0}},\mathrm {Pe}^{-1}\right )=1+C_3\ln (\mathrm {Pe}^{-1})\left (1-\frac {\overline {\langle {\phi _f}\rangle }}{\phi _{f,0}}\right ) ,\end{equation}

with $C_3=0.123$ . We discuss this functional form in detail below, around figure 10. With $\delta _{{est,SL}}$ fully specified, we compare this prediction with the results in figures 5 and 6 (dashed black curves), where it provides a reasonable qualitative and quantitative estimate of $\delta$ across four orders of magnitude in both $\alpha$ and $\mathrm {Pe}^{-1}$ . Note that the function $f$ is a constant of order 1 in figure 6 because it does not depend on $\alpha$ , whereas $f$ varies by a few per cent across the full range of $\mathrm {Pe}^{-1}$ in figure 5 because it is a weak function of $\mathrm {Pe}^{-1}$ .

This estimate ignores the periodic velocity field by assuming that dispersion occurs according to the time-averaged magnitude of $|v_f-v_s|$ , our estimate for which is based on slow loading. During faster loading, $|v_f-v_s|$ will be increasingly localised near the left boundary and suppressed in the interior of the material. We consider the impact of localisation in the next section.

3.5. Quantitative impacts of periodic loading on solute transport

We next consider the effects of the loading parameters $A$ and $T$ . To help interpret these results, we first consider the poromechanical response.

3.5.1. Poromechanical response to periodic loading

In our companion study (Fiori et al. Reference Fiori, Pramanik and MacMinn2023), we explored the impact of $A$ and $T$ on the poromechanical response. During slow loading ( $T\gg 1$ ), the time scale of the loading is much slower than the poroelastic response of the material and the response is quasi-static for any amplitude. The porosity is uniform in space throughout the cycle, returning to its undeformed value at the end of each cycle. The displacement, the fluid velocity, the solid velocity and the Darcy flux all decrease linearly from a spatial extremum at the piston to zero at the right boundary. During fast loading ( $T\ll 1$ ), the time scale of the loading is much faster than the poroelastic response of the material. As a result, the deformation is non-uniform in space and increasingly localised near the left (permeable) boundary as the period decreases. The left portion of the domain experiences both compression and tension, whereas the right portion is in compression at all times. For very fast loading ( $T\lll 1$ ), the deformation is entirely localised near the left boundary and decays exponentially with $x$ , such that the right portion of the material is in static compression. As the amplitude of the deformation increases, the change in porosity at the left boundary with respect to the relaxed state becomes increasingly asymmetric between loading and unloading, with a larger decrease (compression) during loading than the respective increase (tension) during unloading.

Figure 7. Evolution of (a)–(d) fluid flux $q_f$ and (e)–(h) $|v_f-v_s|$ at ten different values of $X=x-u_s(X,t)$ from $0$ to $1$ (dark to light) during one cycle for four different values of $T$ (columns). We distinguish between the loading half of the cycle ( $\dot {a}\gt 0$ ; solid curves) and the unloading half of the cycle ( $\dot {a}\lt 0$ ; dashed curves). (i) Normalised time average of $|v_f-v_s|$ as a function of $X$ for the same four values of $T$ (increasing dark to light). (j) Maximum in time of $|v_f - v_s|$ at $x=1-l$ as a function of $T/\pi$ for $A=0.05$ . The dashed black curve shows the slow-loading prediction of $\pi {}Al/(\phi _{f,0}T)$ and the dotted black curve shows the very-fast-loading prediction of $[\pi {}A/(\phi _{f,0}T)]\exp [-(1-l)\sqrt {\pi /T}]$ (see equation (3.7)).

Figure 8. Evolution of (a)–(d) fluid flux $q_f$ and (e)–(h) relative velocity $|v_f-v_s|$ at ten different values of $X=x-u_s(X,t)$ from $0$ to $1$ (dark to light) during one cycle for $T=0.1\pi$ , and for four different values of $A$ (columns). We distinguish between the loading half of the cycle ( $\dot {a}\gt 0$ ; solid curves) and the unloading half of the cycle ( $\dot {a}\lt 0$ ; dashed curves). (i) Normalised time average of $|v_f-v_s|$ as function of $X$ for the same four values of $A$ (increasing dark to light). (j) Maximum in time of $|v_f-v_s|$ at $x=1-l$ as a function of $A$ .

As noted above, the poromechanical response impacts solute transport through the motion of the fluid and through the changes in porosity. In figure 7 and figure 8, respectively, we show the impact of $T$ and $A$ on the fluid flux $q_f$ (driving advection) and on the relative velocity $|v_f-v_s|$ (driving dispersion). As should be expected, the fluid flux exhibits localisation for fast loading and asymmetry in loading and unloading for large amplitudes. In figure 7, we fix $A$ to the baseline value and consider four values of $T$ . For slow loading, $|v_f - v_s|\sim 2A(1-x)/T$ (see Fiori et al. (Reference Fiori, Pramanik and MacMinn2023)). For very large values of $T$ (e.g. figure 7 d, h), the deformation is uniform and very slow, and $|v_f - v_s|$ is low, especially toward the right (lighter shades) where the solute is positioned. Diffusion dominates over dispersion, even with a large amplitude. As $T$ decreases from $6\pi$ to $0.1\pi$ (e.g. figure 7 b,c,f,g), $|v_f - v_s|$ increases throughout the domain. As $T$ decreases further, however, the deformation is increasingly localised near the left boundary: figures 7(a) and 7(e) show that both $q_f$ and $|v_f-v_s|$ are orders of magnitude larger near the left boundary (darkest curves) than near the right boundary (lightest curves). This localisation is highlighted in figure 7(i), where we plot the time-averaged profile of $|v_f - v_s|$ for each $T$ , and in figure 7(j), where we plot the maximum value of $|v_f - v_s|$ at $x=1-l$ against $T/\pi$ . Near the right boundary, where the solute is located, the relative velocity increases and then decreases with $T$ , exhibiting a maximum around $T=0.1\pi$ .

In figure 8, we show the impact of $A$ on the same quantities for a fixed period, $T=0.1\pi$ . The magnitudes of $q_f$ and $|v_f-v_s|$ increase monotonically with $A$ (figure 8 a–h) and the normalised time-averaged value of $|v_f-v_s|$ is relatively insensitive to $A$ (figure 8 i), suggesting that $|v_f-v_s|$ is essentially proportional $A$ . This suggestion is confirmed in figure 8(j).

3.5.2. Solute transport for different loading amplitudes and periods

Figure 9. (a) Travel distance $\delta _{\mathcal {D}}$ and (b) degree of mixing $\chi _{\mathcal {D}}$ as a function of $T$ and for a wide range of $A$ (evenly spaced from 0.004 to 0.2, increasing dark to light) after a loading time of $12\pi$ followed by a relaxation time of 1 (total time $12\pi +1$ ) for advection and molecular diffusion but no dispersion ( $\alpha =0$ ). Note that the values of $T$ are selected to provide an integer number of loading cycles in a total time of $12\pi$ , but the results are plotted as continuous curves for visual clarity. This constraint leads to periods ranging from $T=0.01\pi$ applied for 1200 cycles to $T=12\pi$ applied for 1 cycle. Dashed lines (darkest colour) correspond to diffusion with no loading ( $A=0$ ). Note that the minimum porosity in the domain occurs near the piston and decreases monotonically with increasing $A$ and with decreasing $T$ . Each curve ends on the left at the value of $T$ for which the minimum porosity vanishes and the simulations fail (see § 2.4).

We showed in § 3.3 how the three transport mechanisms act individually on the concentration profile. We now extend this analysis to examine the roles of $A$ and $T$ . In figure 9, we show the travel distance $\delta _{\mathcal {D}}$ and the degree of mixing $\chi _{\mathcal {D}}$ after a fixed total loading time of $12\pi$ followed by a relaxation time of 1 (total time $12\pi +1$ ), for advection and molecular diffusion but no dispersion ( $\alpha =0$ ). We include results over a wide range of $T$ – from very fast to slow loading – and $A$ – from small to large deformations.

We illustrated in § 3.3 that the contribution of advection is reversible over one loading cycle, independent of $A$ and $T$ . However, as noted above, both the porosity field and the concentration gradients do depend on $A$ and $T$ . Hence, molecular diffusion is expected to vary weakly with $A$ and $T$ . The porosity $\phi _f$ is on average lower than the initial value $\phi _{f,0}$ , because the loading has a non-zero mean – the material is on average compressed. As noted in Fiori et al. (Reference Fiori, Pramanik and MacMinn2023), the overall average porosity $\overline {\langle {\phi _f}\rangle }$ over any integer number of cycles is given by

(3.6) \begin{equation} \overline {\langle {\phi _f}\rangle } = \frac {1}{mT}\int _{nT}^{(n+m)T} \langle {\phi _f}\rangle \,\mathrm {d}t = 1 - \frac {1-\phi _{f,0}}{\sqrt {1-A}}, \end{equation}

where $m$ is an arbitrary positive integer and $n$ is an arbitrary non-negative integer. Note that $\overline {\langle {\phi _f}\rangle }$ decreases with $A$ but is independent of $T$ , as is also true of $\delta _{\mathcal {D}}$ and $\chi _{\mathcal {D}}$ in figure 9, except for very fast loading. For slow loading, it is straightforward to show that the deformation stretches concentration gradients by, on average, a factor of $\phi _{f,0}/\overline {\langle {\phi _f}\rangle }$ .

Figure 10. (a) Variation of $\delta _{\mathcal {D}}$ (red), $\chi _{\mathcal {D}}$ (blue) and $\overline {\langle {\phi _f}\rangle }$ (black) with $A$ for $T=12\pi$ , and where all three quantities are normalised by their values at $A=0$ . We also plot the variation of $\delta _{\mathcal {D}}$ with $A$ , again normalised by its value at $A=0$ , for (b) five values of $\phi _{f,0}\in [0.5,0.95]$ (light to dark) and (c) five values of $\mathrm {Pe}^{-1}\in [3 \times 10^{-6},3 \times 10^{-4}]$ (light to dark). Black dashed lines are the empirical function $f(\overline {\langle {\phi _f}\rangle }/\phi _{f,0},\mathrm {Pe}^{-1})$ from $\delta _{{est,SL}}$ (see (3.4) and (3.5)).

In figure 10(a), we fix $T=12\pi$ and plot $\delta _{\mathcal {D}}$ , $\chi _{\mathcal {D}}$ and $\overline {\langle {\phi _f}\rangle }$ against $A$ ; all three quantities are normalised by their values at $A=0$ , demonstrating that they exhibit a qualitatively similar decay with $A$ . We then plot $\delta _{\mathcal {D}}/\delta _{\mathcal {D}, A=0}$ against $A$ for several values of $\phi _{f,0}$ (figure 10 b) and $\mathrm {Pe}^{-1}$ (figure 10 c). Panel (b) shows that increasing $\phi _{f,0}$ strongly mitigates the decay with $A$ , while panel (c) shows that increasing $\mathrm {Pe}^{-1}$ has a similar but much weaker effect.

With no dispersion ( $\alpha =0$ ), the ratio of $\delta _{{est,SL}}$ to its value for $A=0$ is precisely $f (\overline {\langle {\phi _f}\rangle }/\phi _{f,0},\mathrm {Pe}^{-1} )$ . Thus, the curves in figures 10(b) and 10(c) correspond to profiles of $f$ against $A$ for different values of $\phi _{f,0}$ and $\mathrm {Pe}^{-1}$ . As indicated by the functional dependence of $f$ , we hypothesise that $f$ depends on both $A$ and $\phi _{f,0}$ exclusively through the quantity $\overline {\langle {\phi _f}\rangle }/\phi _{f,0}$ , which, as noted above, measures the (inverse of the) average stretching of concentration gradients due to the deformation (see § 3.3). As indicated in equation (3.5), we find that $f$ is linear in $\overline {\langle {\phi _f}\rangle }/\phi _{f,0}$ to a very close approximation; figure 10(b) demonstrates excellent agreement between this expression for $f$ and our numerical results at fixed $\mathrm {Pe}^{-1}$ .

The dependence of $f$ on $\mathrm {Pe}^{-1}$ is weaker. As indicated in equation (3.5) and illustrated in figure 10(c), we find that a logarithmic variation with $\mathrm {Pe}^{-1}$ provides reasonable qualitative agreement between our expression for $f$ and our numerical results, with good quantitative agreement for larger values of $\mathrm {Pe}^{-1}$ . Note that our expression for $f$ has just one fitting parameter, $C_3$ , the value of which is strongly constrained by the functional form of $f$ and by the clear linear trend with $\overline {\langle {\phi _f}\rangle }/\phi _{f,0}$ ; of course, a better fit is possible with a different functional form and more fitting parameters.

Although the precise functional form chosen here for $f$ is ultimately ad hoc, it is clear that stronger stretching of the concentration profile (i.e. increasing $A$ or $\phi _{f,0}$ , decreasing $\overline {\langle {\phi _f}\rangle }/\phi _{f,0}$ ) will more strongly hinder diffusion (decreasing $f$ ). The role of $\mathrm {Pe}^{-1}$ is more difficult to rationalise. We find that stretching of concentration gradients more strongly hinders diffusion at smaller values of $\mathrm {Pe}^{-1}$ , corresponding to weaker diffusion. At smaller $\mathrm {Pe}^{-1}$ , the initially steep concentration gradient will decay more slowly and it may be the case that diffusion is more sensitive to the stretching of these steeper gradients than shallower ones.

Figure 11. Time evolution of $\delta _{\mathcal {D}}$ at the end of a long series of cycles for six of the smallest values of $T$ considered here, $T=0.015\pi$ to $0.8\pi$ (dark to light) with advection and molecular diffusion but no dispersion ( $\alpha =0$ ). The inset focuses on the very last portion of the main plot. Black stars mark the end of the last cycle of periodic loading and the beginning of the relaxation phase (total time of 1), during which the material returns to its undeformed state. Note that the horizontal axis is on a log scale.

Figures 9 and 10 confirm that $\delta _{\mathcal {D}}$ and $\chi _{\mathcal {D}}$ decrease weakly but monotonically with $A$ , and are essentially independent of $T$ for all but the smallest values of $T$ . For those smallest values, both $\delta _{\mathcal {D}}$ and $\chi _{\mathcal {D}}$ increase with $T$ . This effect is not related to $\overline {\langle {\phi _f}\rangle }$ , which is independent of $T$ (see (3.6)). We explore this behaviour in figure 11 by plotting the evolution of $\delta _{\mathcal {D}}$ over time for a fixed amplitude $A=0.06$ and for several values of $T$ , from $T=0.015\pi$ to $T=0.8\pi$ . For the slowest case ( $T=0.8\pi$ ), $\delta _{\mathcal {D}}$ is minimum at the end of the loading time ( $12\pi$ , marked by a black star) and then increases due to diffusion during the relaxation time (a further time of 1). For smaller values of $T$ , $\delta _{\mathcal {D}}$ instead decreases during the initial stages of relaxation, immediately after the end of loading. This decrease is due to the relaxation of the static far-field compression induced for very fast loading (e.g.  $T \lesssim 0.1 \pi$ ). Since the deformation is much faster than the material response, the domain is never fully relaxed during periodic loading and the right portion, in particular, is in a state of static compression that contributes to an additional stretch of the concentration gradients that hinders diffusion. Once the loading stops, this compressed material relaxes by drawing fluid in, thereby further retracting the solute profile to a degree that is stronger and longer as $T$ decreases, such that $\delta _{\mathcal {D}}$ decreases as $T$ decreases in very fast loading (figures 9 and 11).

Figure 12. (a) Travel distance $\delta$ and (b) degree of mixing $\chi$ as functions of $T$ and for a wide range of $A$ (increasing dark to light) after a loading time of $12\pi$ followed by a relaxation time of 1 (total time $12\pi +1$ ), as in figure 9, but now with all three transport mechanisms simultaneously active. The ranges of amplitudes and periods are the same as in figure 9. Portions shown in grey scale indicate simulations where the solute reaches the left boundary and begins to leave the domain.

We next repeat the preceding analysis, but now including hydrodynamic dispersion (figure 12). Both $\delta$ and $\chi$ are greatly enhanced by dispersion relative to the results in figure 9 across much of the range of $T$ . Recall that the strength of dispersion is expected to scale with $A/T$ . For very slow loading ( $T\gtrsim {}12\pi$ ), the contribution of dispersion is negligible relative to that of diffusion and the values of $\delta$ and $\chi$ converge toward their values without dispersion. As $T$ decreases, the contribution of dispersion grows and increasingly dominates over diffusion, reaching a peak around $T\approx 0.1\pi$ . For these parameters, the values of $\delta$ in the peak region are almost one order of magnitude larger than without dispersion (cf. Figure 9). Across the full range of $T$ where dispersion dominates, $\delta$ and $\chi$ increase with $A$ , as should be expected from figure 8(i): the relative velocity in the interior increases monotonically with $A$ . As $T$ further decreases, $\delta$ and $\chi$ instead begin to decrease and for very fast loading ( $T\lesssim {}0.01\pi$ ) approach their values without dispersion. This trend can be linked to the impact of $T$ on the ratio of $|v_f - v_s|$ at $x=1-l$ , shown in figure 7(i): for faster loading, the deformation is increasingly localised near the piston and the region occupied by the solutes is increasingly not engaged. We examine these observations in more detail in figure 13, by plotting the evolution of $\delta$ throughout the loading time for a fixed amplitude and for several periods between $T=0.03\pi$ and $T=0.8\pi$ , thus spanning the peak in figure 12. In all cases, $\delta$ exhibits oscillations with period $T$ on top of a roughly Fickian growth. As $T$ decreases, these oscillations decrease in magnitude as they increase in frequency, consistent with the deformation being increasingly localised at the left. The overall rate of spreading increases as $T$ decreases from $0.8\pi$ to $0.1\pi$ , for which the increase in frequency leads to a net increase in dispersive flux despite the decrease in magnitude (see figure 7). As $T$ decreases further, the decrease in local magnitude begins to dominate the increase in frequency and the rate of spreading instead slows (see again figure 7).

Figure 13. Evolution of $\delta$ over the entire loading time for $A=0.06$ and for five values of $T$ (dark to light, values as indicated) when advection, molecular diffusion and hydrodynamic dispersion are simultaneously active. We also show the case of diffusion only (no loading, lightest curve).

We next compare the numerical values of $\delta$ shown in figure 12 with a revised estimate $\delta _{{est}}$ that accounts for diffusion, dispersion and localisation. An analytical expression for the Darcy flux during very fast loading at very low amplitude ( $T\ll 1$ and $A\ll 1$ ; from linear poroelasticity and analogous to Stokes’s classical ‘second problem’) can be derived from the results of Fiori et al. (Reference Fiori, Pramanik and MacMinn2023) by calculating the solid velocity as $v_{s, {vf}}\approx {}\mathrm {d}u_{s, {vf}}/\mathrm {d}t$ , where the solid displacement $u_{s, {vf}}$ is given in their equation (3.12) and the subscript (‘ $\mathrm {vf}$ ’) stands for ‘very fast’. Equation (2.22) (also their equation 2.20) implies that the Darcy flux is then $\phi _{f, {vf}}(v_{f, {vf}}-v_{s, {vf}})=-v_{s, {vf}}$ , such that

(3.7) \begin{equation} \begin{aligned} \phi _{f, {vf}} (v_{f, {vf}} -v_{s, {vf}}) = - &\frac {\pi A}{T}\left [ \exp \left (-x\sqrt \frac {\pi }{T}\right )\sin \left (\frac {2\pi {}t}{T}-x\sqrt \frac {\pi }{T}\right )\right ]. \end{aligned} \end{equation}

Figure 14. Estimated travel distance $\delta _{{est}}$ as a function of $T$ and for a wide range of $A$ (increasing dark to light) after a loading time of $12\pi$ followed by a relaxation time of 1 (total time $12\pi +1$ ), as in figure 12, but here calculated via equation (3.8). The dashed black line corresponds to the maxima for $\delta$ from figure 12, for comparison.

Equation (3.7) suggests that, in this regime, the amplitude of the Darcy flux decays exponentially with distance from the left boundary and with $T^{-1/2}$ . Based on this solution, we reformulate equation (3.4) to include this exponential localisation in the magnitude of the dispersive flux

(3.8) \begin{equation} \delta _{ {est}}=C_1 f\left (\frac {\overline {\langle {\phi _f}\rangle }}{\phi _{f,0}} ,\mathrm {Pe}^{-1}\right ) \sqrt { 4 \left [\mathrm {Pe}^{-1} + C_2 \alpha \frac {A}{T} \exp \left ( {- (1-l) \sqrt {\frac {\pi }{T}}}\right )\right ] t}, \end{equation}

where $f (({\overline {\langle {\phi _f}\rangle }}/{\phi _{f,0}}) ,\mathrm {Pe}^{-1} )$ is as defined in (3.5). Figure 14 shows that $\delta _{{est}}$ captures the qualitative trends observed for $\delta$ in figure 12 (for which we report the maxima for a direct comparison), with the central peak that is, however, lower and slightly shifted compared with the one observed for $\delta$ . We do not expect strong quantitative agreement because this estimate ignores the details of the spatial variation in the dispersive flux, which decreases linearly from left to right even for slow loading, as well as the oscillatory nature of the flow and, most importantly, is based on an analytical solution that is only valid for small deformations and in the very-fast-loading region of the parameter space.

Hence, figures 12–14 highlight two competing mechanisms: progressively faster loading enhances dispersion by promoting large dispersive fluxes in general, but also progressively localises the flow and deformation near the left boundary, suppressing the flux (and thus dispersion) in the interior. The competition between these two behaviours is the origin of the local maximum in $\delta$ and $\chi$ with $T$ and is qualitatively captured by the estimated travel distance $\delta _{{est}}$ .

In summary, this analysis reveals the link between transport fluxes and loading parameters and, consequently, the poromechanics of periodic loading. Through the travel distance $\delta _{{est}}$ , we are able to capture the three main transport regimes, corresponding to analogous regions of the $A-T$ domain shown in figure 12:

1. (i) $T\gtrsim 0.3\pi$ : the loading is slow and the strength of dispersion is roughly proportional to $A/T$ . Thus, faster and larger deformations progressively increase the strength of dispersion relative to diffusion. In particular, $|v_f-v_s|$ reaches higher peaks as $T$ decreases, promoting hydrodynamic dispersion, which is the dominant mechanism for this region.

2. (ii) $0.1 \pi \lesssim T \lesssim 0.3 \pi$ : the loading is fast, with an optimal balance between magnitude and depth of the fluid fluxes; dispersion is at its peak.

3. (iii) $T \lesssim 0.1 \pi$ : the loading is fast and the strong localisation of the flow and deformation near the left boundary dominates, increasingly suppressing dispersion in the interior as $T$ further decreases.

Note that the quantitative values of $\delta$ and $\chi$ also depend on the specific values of $\mathrm {Pe}^{-1}$ and $\alpha$ , but the evolution of the porosity and fluid velocity with $A$ and $T$ does not. As a result, varying $\mathrm {Pe}^{-1}$ and $\alpha$ would change the width and height of the curves in figure 12, but would not change the qualitative features of these plots or the position of the central peak.

3.6. Péclet number in biological examples

We conclude by reporting typical dimensionless loading parameters and diffusion coefficients for several examples of soft biological tissues (table 2). Based on these values, we calculate the Péclet number and the effective Péclet number, as defined in § 2.7, to understand the range of loading scenarios where advection and dispersion are likely to be non-negligible. For simplicity, we use a moderate dimensionless amplitude $\tilde {A}=0.1$ in all cases. The detailed parameters to calculate $\tilde {T}$ and $\mathrm {Pe}$ are reported, for the same biological tissues, in Fiori et al. (Reference Fiori, Pramanik and MacMinn2023). In all cases, the values of $\mathrm {Pe}_{\it eff}$ range from moderate to high, indicating that advection and potentially dispersion are unlikely to be negligible in these systems.

Table 2. Loading and transport parameters for some examples of biological materials.