# Solute Transport due to Periodic Loading in a Soft Porous Material

Matilde Fiori  Department of Engineering Science, University of Oxford,
Oxford, OX1 3PJ, UK  Satyajit Pramanik  Department of Mathematics, Indian
Institute of Technology Guwahati, Guwahati – 781039, Assam, India  Christopher
W. MacMinn  Department of Engineering Science, University of Oxford, Oxford,
OX1 3PJ, UK<EMAIL_ADDRESS>
###### Abstract

In soft porous media, deformation drives solute transport via the intrinsic
coupling between flow of the fluid and rearrangement of the pore structure.
Solute transport driven by periodic loading, in particular, can be of great
relevance in applications ranging from the geomechanics of contaminants in the
subsurface to the biomechanics of nutrient transport in living tissues,
scaffolds for tissue engineering, and biomedically employed hydrogels.
However, the basic features of this process have not previously been
systematically investigated. Here, we fill this hole in the context of a 1D
model problem. We do so by expanding the results from a companion study, in
which we explored the poromechanics of periodic deformations, by introducing
and analysing the impact of the resulting fluid and solid motion on solute
transport. We first characterise the independent roles of the three main
mechanisms of solute transport in porous media — advection, molecular
diffusion, and hydrodynamic dispersion — by examining their impacts on the
solute concentration profile during one loading cycle. We next explore the
impact of the transport parameters, showing how these alter the relative
importance of diffusion and dispersion. We then explore the loading parameters
by considering a range of loading periods — from slow to fast, relative to the
poroelastic timescale — and amplitudes — from infinitesimal to large. We show
that solute spreading over several loading cycle increases monotonically with
amplitude, but is maximised for intermediate periods because of the increasing
poromechanical localisation of the flow and deformation near the permeable
boundary as the period decreases.

##  I Introduction

Solutes spread and mix in deformable porous media in a variety of
geomechanical, industrial, and biological contexts. In general, the transport
of solutes in porous media is driven by molecular diffusion and by internal
fluid flow. In soft porous media, the latter is strongly coupled to external
mechanical loads through rearrangement of the pore structure (e.g., Mow  _et
al._ , 1980; Lai  _et al._ , 1991; Preziosi  _et al._ , 1996; Li  _et al._ ,
2004; Franceschini  _et al._ , 2006; Moeendarbary  _et al._ , 2013; Ehlers
_et al._ , 2009; Vuong  _et al._ , 2015; Borja and Choo, 2016). In many cases,
these loads are periodic; for example, compression due to surface loading can
induce the spreading of contaminants in soils, exacerbating environmental harm
and hindering remediation, while physiological loads can drive nutrient
transport and waste removal in biological tissues, thus potentially playing an
important role in cell growth and survival. In a companion study (Fiori  _et
al._ , 2023), we examined the poromechanics of periodic loading over a wide
range of loading periods and amplitudes. Here, we examine the implications of
those results for solute transport.

At the continuum (Darcy) scale, which is the framework of this study, solute
transport occurs through three fundamental mechanisms: advection, diffusion,
and hydrodynamic dispersion Bear (1972); Whitaker (1967); Saffman (1959);
Scheidegger (1961); Brenner and Edwards (1993); Whitaker (1998); Gelhar
(1993); Dentz  _et al._ (2011). Diffusion in a porous medium is weaker than in
a bulk fluid because of the tortuosity of the pore space Bear (1972);
Ghanbarian  _et al._ (2013); Tartakovsky and Dentz (2019). Both advection and
dispersion are driven by fluid flow. Advection is driven by the mean
interstitial fluid velocity. Dispersion results from the pore-scale deviations
from this Darcy-scale mean. In particular, dispersion is driven by two
phenomena: (i) analogous to classical Taylor dispersion in a tube Taylor
(1953); Brenner and Stewartson (1980); Marbach and Alim (2019), pore-scale
velocity gradients smear solute profiles along the flow direction, inducing
longitudinal spreading, and (ii) the morphology of the pore structure
introduces chaotic variability in the fluid streamlines Datta  _et al._
(2013); de Anna  _et al._ (2013); Lester  _et al._ (2013, 2016a, 2016b); Gouze
_et al._ (2021); Dentz  _et al._ (2023); Kree and Villermaux (2017); Souzy
_et al._ (2020), thus inducing both longitudinal and transverse spreading
Scheidegger (1961); Gelhar and Axness (1983); Gelhar  _et al._ (1992); Delgado
(2007). In soft porous media, therefore, deformation can enhance the transport
of solutes directly by driving fluid flow, thus leading to advection and
dispersion, and indirectly by distorting the pore space, and thus modifying
both dispersion and diffusion.

Solute transport in deformable porous media has been studied previously in
several different contexts. The impact of soil consolidation on transport has
been studied primarily for geomechanical applications (e.g., Fox, 2007a, b; Pu
_et al._ , 2018). In that context, it has been shown that solute transport is
enhanced by larger deformations Peters and Smith (2002); Lewis  _et al._
(2009); Alshawabkeh and Rahbar (2006); Zhang  _et al._ (2013); Bonazzi  _et
al._ (2021). Deformation has also been shown to increase mixing and reduce
breakthrough time in the context of miscible viscous fingering Tran and Jha
(2020). The key feature introduced by periodic loading is the continuously
fluctuating fluid flow, which can irreversibly modify diffusion and dispersion
even when the macroscopic advective component is perfectly reversible. The
role of periodic flow in enhancing solute transport and mixing has been
studied in rigid and compressible 1D pore networks Goldsztein and Santamarina
(2004); Claria  _et al._ (2012). In a poroelastic material, the impact of
small periodic deformations on solute transport has been analysed for a range
of frequencies Bonazzi  _et al._ (2021). In heterogeneous coastal aquifers,
periodic flow due to tidal loading has been shown to enhance solute transport
Wu  _et al._ (2020); Trefry  _et al._ (2019). However, all of these studies
examined a specific and relatively limited range of parameters.

In the context of biological tissues, periodic loading is known to enhance the
transport of nutrients (Ferguson  _et al._ , 2004; Gardiner  _et al._ , 2007;
Witt  _et al._ , 2014; Schmidt  _et al._ , 2010; Zhang, 2011; DiDomenico  _et
al._ , 2017; Zhang and Szeri, 2008). Similarly, periodic deformations are used
to enhance the infiltration of solutes into hydrogels Albro  _et al._ (2008);
Vaughan  _et al._ (2013) and other scaffolds for tissue engineering Mauck  _et
al._ (2003); Cortez  _et al._ (2016); Kumar  _et al._ (2018); Fan  _et al._
(2016), where the correlation between loading parameters, nutrient transport
and cell survival is of particular interest. Increasing the loading amplitude
and/or decreasing the loading period induces a transition from diffusion-
dominated to advection-dominated regimes Urciuolo  _et al._ (2008) and
amplifies the role of hydrodynamic dispersion Sengers  _et al._ (2004).
Decreasing the loading period also leads to localisation of flow and
deformation near permeable boundaries, resulting in larger velocities near the
surface that promote external solute infiltration Gardiner  _et al._ (2007);
Urciuolo  _et al._ (2008); DiDomenico  _et al._ (2017); Vaughan  _et al._
(2013).

In general, despite the established role of hydrodynamic dispersion in driving
the transport of solutes in porous media, dispersion is rarely included in
biomechanical models (with the notable exception of Ref. Sengers  _et al._
(2004)). One context where dispersion is widely agreed to be important is in
brain microcirculation Kelley and Thomas (2023). In the vascular network
within the brain, dispersion results from the shear-induced radial
concentration gradients in single vessels (e.g. Marbach and Alim (2019); Sharp
_et al._ (2019); Berg  _et al._ (2020); Troyetsky  _et al._ (2021);
Bojarskaite  _et al._ (2023)) and the progressive bifurcation of vessels into
smaller branches that can be modelled at the continuum scale as a porous
material (e.g. Zimmerman and Tartakovsky (2020); Goirand  _et al._ (2021)).

Dispersion is typically neglected in the context of tissues and gels for two
main reasons. First, fluid flow is often assumed to be slow, implying that
transport is dominated by diffusion. In other words, the Peclet number
$\mathrm{Pe}=VL/D$ is assumed to be small, where $V$ is the fluid velocity,
$L$ the streamwise length scale, and $D$ the diffusivity. However, it is
straightforward to show that $\mathrm{Pe}$ can be of order 1 or larger in a
tissue or gel subject to fast ($0.1-1$ Hz) and large ($10-20\%$) deformations
(see table II), suggesting that dispersion may be important or even dominant
in some scenarios Delgado (2007). Indeed, many studies highlight a transition
from diffusion-dominated to advection-dominated transport without
acknowledging the potential role of dispersion Gardiner  _et al._ (2007);
Urciuolo  _et al._ (2008); DiDomenico  _et al._ (2017); Vaughan  _et al._
(2013). With an analogous argument, Davit  _et al._ (2013) illustrated the
importance of including dispersion in models for solute transport in biofilms.
The second typical reason for neglecting dispersion in tissues and gels is the
assumption that the longitudinal and transverse dispersivities themselves are
negligible. This expectation is a result of physical insight derived from
transport in granular materials, where the dispersivity is typically taken to
be proportional to the pore size Saffman (1959); Oswald and Kinzelbach (2004);
Kree and Villermaux (2017); Liang  _et al._ (2018). Indeed, the typical pore
size is $\sim 10$ nm in polymeric gels and in the extra-cellular matrix of
tissues (e.g., around 6 nm in cartilage Mow  _et al._ (1984)) and can
therefore be similar to (or smaller than) the size of large solute molecules
Maroudas (1970); DiDomenico  _et al._ (2018), originating solid-solute
friction Ateshian and Weiss (2013); Yao and Gu (2007). However, tissues and
scaffolds are heterogeneous and multiscale materials; the presence of other
components, such as collagen fibres, originates a “mesoscale” of larger pores
(e.g., 100–150 nm in cartilage Maroudas (1975); Levick (1987); Federico and
Herzog (2008)), where even larger solute molecules can pass DiDomenico  _et
al._ (2017, 2018) and where dispersion is likely to play a much larger role.
The same is true for double-porosity scaffolds and gels, where additional
channels and/or macroscopic pores are included to enhance fluid flow
throughout the scaffold depth Buijs  _et al._ (2010); Lee  _et al._ (2015);
Mesallati  _et al._ (2013). As a further counter-argument, we hypothesise
that, even in pores that are small compared to the solute molecules, the
irrelevance of pore-scale velocity gradients does not exclude velocity
variations and streamline alterations in the overall network, which could
originate longitudinal and transverse dispersion. This hypothesis is
consistent with the quantification of tortuosity in several soft tissues
Maroudas (1970); Hrabe  _et al._ (2004); Zhang and Szeri (2005).

Thus, the impact of periodic loading on solute transport in soft porous media
has been addressed with various approaches and assumptions across a variety of
specific applications in soils, tissues, hydrogels, and scaffolds. However, no
single study has yet provided a comprehensive understanding across a wide
range of loading frequencies and amplitudes. Moreover, the impact of
hydrodynamic dispersion remains relatively unexplored and therefore poorly
understood, particularly in the context of biological and biomedical
applications. Here, we study the transport and mixing of solutes due to
arbitrarily large, periodic deformations of a soft porous material. For the
flow and deformation, we adopt a one-dimensional, large-deformation
poroelasticity model that includes rigorous nonlinear kinematics, deformation-
dependent permeability, and Hencky elasticity for the solid skeleton. We used
this model to explore the poromechanics of large-amplitude periodic loading in
a companion study (Fiori  _et al._ , 2023). Here, we introduce a model for
solute transport due to advection, diffusion, and dispersion. We study the
separate roles of advection, diffusion, and dispersion during one loading
cycle. We then consider the impact of the transport and loading parameters on
transport and mixing over longer time periods and/or larger number of loading
cycles. We report the impact of a wide range of loading amplitudes and periods
on each transport mechanism and observe how transport depends on the
poromechanical response through its impact on local fluid flow. When
dispersion is negligible, we show that diffusion is insensitive to loading
period but slightly suppressed by increased loading amplitude. With
dispersion, larger amplitudes always boost solute spreading; however,
progressively shorter periods impact transport and mixing in more complex
ways: fast loading promotes spreading by inducing large fluid velocities, but
very fast loading hinders spreading by progressively localising the flow and
deformation. We show that the competition between these two effects results in
maximum solute transport and mixing for intermediate loading periods.

##  II Theoretical Model

Our model combines large-deformation poroelasticity with solute transport. The
coupling between periodic deformations and solute movement occurs primarily
via the fluid flow, which is originated by the former and responsible for the
latter.

###  II.1 Model Problem

We consider a one-dimensional 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:

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

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._ (2023). 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 1D 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).

###  II.2 Kinematics

We consider an Eulerian reference frame, in which the solid displacement is
$\mathbf{u_{s}}=\mathbf{x}-\mathbf{X}(\mathbf{x},t)$, with
$\mathbf{X}(\mathbf{x},t)$ the reference position of the 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_{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

$\mathbf{u_{s}}=u_{s}(x,t)\mathbf{\hat{e}_{x}},\;\mathbf{v_{s}}=v_{s}(x,t)\mathbf{\hat{e}_{x}},\;\mathbf{v_{f}}=v_{f}(x,t)\mathbf{\hat{e}_{x}},$
(2)

where $\mathbf{v_{f}}$ and $\mathbf{v_{s}}$ are the fluid and solid
velocities, respectively, $u_{s}$, $v_{s}$, and $v_{f}$ are the $x$-components
of these fields, and $\mathbf{\hat{e}_{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

$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}}.$
(3)

Continuity can be written

$\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,$
(4)

which together imply that the total flux $q=\phi_{f}v_{f}+(1-\phi_{f})v_{s}$
is uniform in space, $\partial{q}/\partial{x}=0$.

###  II.3 Fluid flow

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

$\phi_{f}(v_{f}-v_{s})=-\frac{k(\phi_{f})}{\mu}\frac{\partial p}{\partial x},$
(5)

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 Ref. Fiori  _et al._ (2023), we take
the permeability to be deformation-dependent according to a normalised Kozeny-
Carman relation,
$k(\phi_{f})=k_{0}\frac{(1-\phi_{f,0})^{2}}{\phi_{f,0}^{3}}\frac{\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 Ref. Fiori  _et al._ (2023).

Combining equations (4) and (5), we arrive at the nonlinear flow equations:

$\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,$
(6)

where the total flux $q$ is again

$q\equiv\phi_{f}v_{f}+(1-\phi_{f})v_{s}$ (7)

and the fluid and solid velocities are given by

$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}.$ (8)

Note that the fluid flux is

$q_{f}=\phi_{f}v_{f}.$ (9)

###  II.4 Mechanical equilibrium and elasticity law

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

$\frac{\partial\sigma^{\prime}}{\partial x}=\frac{\partial p}{\partial x},$
(10)

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

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:

$\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,$
(11)

where the nonlinear composite constitutive function

$D_{f}(\phi_{f})=(1-\phi_{f})\frac{k(\phi_{f})}{\mu}\frac{\mathrm{d}\sigma^{\prime}}{\mathrm{d}\phi_{f}}$
(12)

is the poroelastic diffusivity.

We use Hencky hyperelasticity Hencky (1931) to capture kinematic nonlinearity,
as discussed in detail in Ref. Fiori  _et al._ (2023). For a uniaxial
deformation, the relevant component of the effective stress is then MacMinn
_et al._ (2016); Auton and MacMinn (2018)

$\sigma^{\prime}=\mathcal{M}\frac{\ln(J)}{J}=\mathcal{M}\bigg{(}\frac{1-\phi_{f}}{1-\phi_{f,0}}\bigg{)}\ln\bigg{(}\frac{1-\phi_{f,0}}{1-\phi_{f}}\bigg{)},$
(13)

where $\mathcal{M}$ is the $p$-wave or oedometric modulus. With appropriate
initial and boundary conditions and the normalised Kozeny-Carman permeability
law, equations (11), (12), and (13) comprise a closed model for the evolution
of the porosity.

###  II.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 1D, conservation of mass for the solute can be written

$\frac{\partial}{\partial t}(\phi_{f}c)+\frac{\partial}{\partial
x}\bigg{[}{\phi_{f}cv_{f}}-\phi_{f}\mathcal{D}\frac{\partial c}{\partial
x}\bigg{]}=0.$ (14)

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 (1961)). 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

$\mathcal{D}=\mathcal{D}_{m}+\mathcal{D}_{h},$ (15)

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._ (2018, 2023). The most widely
used model for the macroscopic dispersive flux is Fickian, as above, with a
velocity-dependent dispersion coefficient given in 1D by Scheidegger (1961);
Brenner and Edwards (1993); Whitaker (1998); Gelhar (1993)

$\mathcal{D}_{h}=\alpha|v_{f}-v_{s}|,$ (16)

where $\alpha$ is the longitudinal dispersivity. Note that the dispersive flux
is therefore 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 the 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. This coupling has not yet been
studied and will be the subject of future work. Here, we take $\alpha$ to be
constant for simplicity.

###  II.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.

####  II.6.1 Initial conditions

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

$\phi_{f}(x,0)=\phi_{f,0}\textrm{ and }u_{s}(x,0)=0.$ (17)

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

$c(x,0)=\frac{c_{0}}{2}\left\\{\tanh{[s(x-L+l)}]+1\right\\},$ (18)

where $s$ is a steepness parameter.

####  II.6.2 Left boundary

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

$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.$
(19)

We take the fluid outside the domain to be “clean”, such that

$c(a,t)=0.$ (20)

####  II.6.3 Right boundary

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

$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.$ (21)

Equation (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 (7) then implies that the fluid and the solid always locally move in
opposite directions,

$v_{f}=-\frac{(1-\phi_{f})}{\phi_{f}}v_{s}.$ (22)

###  II.7 Scaling and summary

As in in Ref. Fiori  _et al._ (2023), we apply the following non-
dimensionalization to the poromechanical model

$\begin{split}\tilde{x}=\frac{x}{L},\;\tilde{u}_{s}=\frac{u_{s}}{L},\;\tilde{t}=\frac{t}{T_{\mathrm{p}e}},\;\tilde{\sigma}^{\prime}_{x}=\frac{\sigma^{\prime}_{x}}{\mathcal{M}},\;\tilde{p}=\frac{p}{\mathcal{M}},\;\tilde{k}=\frac{k(\phi)}{k_{0}},\;\tilde{v}_{f}=\frac{v_{f}}{L/T_{\mathrm{p}e}},\;\tilde{v}_{s}=\frac{v_{s}}{L/T_{\mathrm{p}e}},\end{split}$
(23)

where $T_{\mathrm{pe}}=L^{2}/D_{f,0}=\mu{}L^{2}/(k_{0}\mathcal{M})$ is the
classical poroelastic timescale 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

$\tilde{c}=\frac{c}{c_{0}},\;\tilde{l}=\frac{l}{L},\;\tilde{\alpha}=\frac{\alpha}{L}.$
(24)

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

$\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,$
(25)

where

$\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}},$
(26)

and

$\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.$
(27)

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

$\tilde{\mathcal{D}}=\frac{\mathcal{D}}{\mathcal{D}_{m}}=\mathrm{Pe}^{-1}+\tilde{\alpha}|\tilde{v}_{f}-\tilde{v}_{s}|,$
(28)

where
$\mathrm{Pe}=\frac{L^{2}/T_{\mathrm{p}e}}{\mathcal{D}_{m}}=\frac{k_{0}\mathcal{M}}{\mu\mathcal{D}_{m}}$
is the Péclet number, which measures the importance of deformation-driven
advection relative to molecular diffusion. The initial conditions are

$\tilde{a}(0)=0,\;\phi_{f}(\tilde{x},0)=\phi_{f,0},$ (29)

and

$\tilde{c}(\tilde{x},0)=\frac{1}{2}\\{\tanh{[\tilde{s}(\tilde{x}-1+\tilde{l})}]+1\\},$
(30)

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

$\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,$
(31)

and

$\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,$
(32)

where $\tilde{A}=A/L$ and $\tilde{T}=T/T_{\mathrm{pe}}$. Below, we consider
only dimensionless quantities, dropping the tildes for convenience.

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. Below,
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.

##  III Solute transport and mixing

###  III.1 Quantification of solute transport and mixing

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. 1.

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 interaction of the concentration profile with the
left boundary increases strongly.

  2. 2.

The degree of mixing $\chi$ measures the degree to which the initial
concentration profile has homogenised, and is related to the spatial variance
of the concentration profile about its center of mass. We express the degree
of mixing in terms of the variance of the concentration profile by
generalising the definitions of Jha  _et al._ (2011) to account for a porosity
field that varies in space. Considering the fluid-volume-weighted average
$\langle\ast\rangle_{f}$, defined as

$\displaystyle\langle\ast\rangle_{f}=\frac{\int_{a}^{1}\phi_{f}\ast\mathrm{d}x}{\int_{a}^{1}\phi_{f}\;\mathrm{d}x},$
(33)

the variance of the concentration profile is then

$\sigma^{2}(t)=\langle{}c^{2}\rangle{}_{f}-\langle{}c\rangle{}_{f}^{2}$ (34)

and the degree of mixing is

$\chi(t)=1-\frac{\sigma^{2}(t)}{\sigma_{\rm max}^{2}},$ (35)

where $\sigma_{\rm 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).

Figure 2: Schematic representation of the travel distance or mixing length
$\delta$, which measures the distance traveled 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, for example, Tan and Homsy
(1988) or Mishra  _et al._ (2008)). 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$.

###  III.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
LABEL:tab:quantitative). 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 section
III.4, $A$ and $T$ in section III.5, and $\phi_{f,0}$ and $l$ in Appendix D,
respectively.

We choose a baseline amplitude $A=0.1$, corresponding to moderate-large
deformations. We choose a baseline period $T=6\pi$, which, following our
companion study (Ref. Fiori  _et al._ (2023)), ensures that the poromechanics
are quasi-static (i.e., “slow-loading”, see the first part of section III.5).
The baseline values of $\mathrm{Pe}^{-1}$ and $\alpha$ are in the range of
those proposed by Sengers  _et al._ (2004) 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.

[ caption=Baseline parameter values. , label = tab:quantitative, entry = none,
] width = colspec = Q[185]Q[77]Q[83]Q[235]Q[121]Q[142]Q[90], cells = c,
vline2-7 = -, hline2 = -, Parameter & $A$ $T$ $\mathrm{Pe}^{-1}$ $\alpha$
$\phi_{f,0}$ $l$
Baseline value $0.1$ $6\pi$ $3\times 10^{-5}$ $0.01$ $0.75$ $0.05$

###  III.3 Qualitative impacts of periodic loading on solute transport

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_{\mathrm{adv}}=\phi_{f}v_{f}c$,
$q_{\mathrm{diff}}=\phi_{f}\mathrm{Pe}^{-1}(\partial c/\partial x)$, and
$q_{\mathrm{disp}}=\phi_{f}\alpha|v_{f}-v_{s}|(\partial c/\partial x)$,
respectively.

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 curve is proportional to
$t^{-\frac{1}{2}}$.

During the loading half of each cycle ($\dot{a}>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}<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 Ref Fiori  _et al._ (2023)), in which case the
net contribution of advection over one full cycle is zero (see figure 4b).
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 strength of diffusion and dispersion are proportional to
$\mathrm{Pe}^{-1}$ and $\alpha$, respectively.

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.

Figure 4: Evolution of the concentration profile during one cycle (red to blue
through white) for four cases: (a) diffusion with no loading
($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 Lagrangian spatial coordinate $X=x-u_{s}$ and
split the evolution into two phases, loading ($\dot{a}>0$, first half of the
cycle, dark to light red) and unloading ($\dot{a}<0$, second half, light to
dark blue). The initial profile is shown in black. For each case, we also show
the associated evolution of $\delta$ (insets); in all cases, the dotted curves
are for diffusion without loading (with the dashed reference line to show
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.

For diffusion only (figure 4a) 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 fluid moves, driving advection;
(ii) the fluid moves through the pore space driving dispersion; (iii) the
porosity decreases, weakly hindering diffusion and dispersion since $\phi_{f}$
is a prefactor in both of those fluxes; and (iv) the solute spreads since a
lower porosity means that the same quantity of fluid (and solute) must now
occupy a larger spatial extent, thus weakening 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 4b), 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 4c, where the final value of $\delta$ is
lower for diffusion with loading than for diffusion without loading. When
dispersion is included (figure 4d), transport is greatly amplified despite the
weak suppression of diffusion.

We next consider several quantitative measures of transport and mixing.

###  III.4 Impact of diffusion and dispersion coefficients on transport and
mixing

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$. In Appendix D, we consider the roles of
$\phi_{f,0}$ and $l$. The impact of changing $\mathrm{Pe}^{-1}$ and $\alpha$
on $\delta$ are shown in figure 5 and figure 6, respectively.

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 line shows the approximate theoretical prediction
$\delta_{\mathrm{est,SL}}$ from equation (36).  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 line shows the approximate
theoretical prediction $\delta_{\mathrm{est,SL}}$ from equation (36).

All of the curves are roughly linear in $\sqrt{t}$ (i.e., spreading is Fickian
on average), exhibiting fluctuations with a period $T$ because of the
phenomenon described in section III.3: loading decreases the porosity, forcing
the solute to spread (advection to the left), and unloading increases the
porosity, forcing the solute to contract (advection to the right). Larger
values of $\mathrm{Pe}^{-1}$ and $\alpha$ enhance the respective fluxes 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). In both figure 5 and figure 6, we
compare the value of $\delta$ in each case with a naive Fickian estimate
$\delta_{\mathrm{est,SL}}$, in which we neglect advection but include
dispersion:

$\delta_{\mathrm{est,SL}}=C_{1}\sqrt{4\left(\mathrm{Pe}^{-1}+C_{2}\alpha\frac{A}{T}\right)t}.$
(36)

The constants $C_{1}$ and $C_{2}$ are derived fitting the curves corresponding
to $\delta_{\mathrm{est,SL}}$ to the data in figure 6 and figure 5,
respectively, with $C_{1}=1.28$ and $C_{2}=0.25$. Note that this estimate
ignores the periodic nature of the flow and assumes that the time-averaged
magnitude of $|v_{f}-v_{s}|$ decreases linearly from left to right during slow
loading (see Ref. Fiori  _et al._ (2023)) and therefore is proportional to
$A/T$ throughout the material. During faster loading, however, $|v_{f}-v_{s}|$
is increasingly localised near the left boundary and suppressed in the
interior of the material. This estimate is therefore only appropriate for slow
loading (SL). We consider the impact of localisation in detail in the next
section.

###  III.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 provide a brief summary of the key findings of our
companion study before considering solute transport and mixing.

##### Poromechanical response to periodic loading

In our companion study (Fiori  _et al._ , 2023), we explored the impact of $A$
and $T$ on the poromechanical response. During slow loading ($T\gg 1$), the
timescale 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
timescale 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.

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 $|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 details, in figure
7, we fix $A$ to the baseline value and consider four values of $T$. For slow
loading, $|v_{f}-v_{s}|\sim A/T$ (see equation (47) and Ref. Fiori  _et al._
(2023)). For very large values of $T$ (e.g. figure 7d and h), the deformation
is uniform and very slow, and $|v_{f}-v_{s}|$ is low, especially towards the
right (lighter lines) where the solute is positioned. Diffusion dominates over
dispersion, even with a large amplitude. As $T$ decreases from $10\pi$ to
$0.1\pi$ (e.g., figure 7b,c,f,g), the entire domain is characterised by an
increasingly higher $|v_{f}-v_{s}|$. As $T$ decreases further, however, the
deformation is increasingly localised at the left boundary: figure 7a and e,
show that both $q_{f}$ and $|v_{f}-v_{s}|$ are orders of magnitude larger at
the left boundary (darkest curves) than near the right boundary (lightest
curves). This localisation is highlighted in figure 7i, where we plot the
maximum value of $|v_{f}-v_{s}|$ at $x=1-l$ against $T$: towards the right
boundary, where the solute is located, the strength of dispersion 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$. As $A$ increases, $q_{f}$ and $|v_{f}-v_{s}|$ progressively
increase, with increasingly larger and lower magnitudes during loading and
unloading, respectively. For fixed $T$, the strength of dispersion increases
monotonically with $A$ (figure 8i).

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}>0$; solid
curves) and the unloading half of the cycle ($\dot{a}<0$; dashed curves). (i)
Maximum in time of $|v_{f}-v_{s}|$ at $x=1-l$ as a function of $T$.  Figure 8:
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 $T=0.1\pi$, and for four different values of $A$
(columns). We distinguish between the loading half of the cycle ($\dot{a}>0$;
solid curves) and the unloading half of the cycle ($\dot{a}<0$; dashed
curves). (i) Maximum in time of $|v_{f}-v_{s}|$ at $x=1-l$ as a function of
$A$.

##### Solute transport for different loading amplitude and period

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 color) 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.

We showed in section III.3 how the three transport mechanisms act individually
on the concentration profile. We now extend this analysis to examine the role
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 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 section III.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 Ref. Fiori  _et al._ (2023), the overall
average porosity $\langle\overline{\phi}_{f}\rangle$ over any integer number
of cycles is given by:

$\langle\overline{\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}},$
(37)

where $\langle{\phi_{f}}\rangle$ indicates the spatial average. Note that
$\langle\overline{\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. In figure 10, we fix $T=12\pi$ and
plot the normalised values of $\delta_{\mathcal{D}}$, $\chi_{\mathcal{D}}$,
and $\langle\overline{\phi}_{f}\rangle$ against $A$, demonstrating that all
three quantities exhibit a qualitatively similar decay with $A$.

Figure 10: Normalised $\delta_{\mathcal{D}}$ (dashed red curve),
$\chi_{\mathcal{D}}$ (dotted blue curve), and
$\langle\overline{\phi}_{f}\rangle$ (solid black curve) against $A$ for
$T=12\pi$.

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$ (see (37)).
This effect is not related to $\langle\overline{\phi}_{f}\rangle$, which is
independent of $T$. 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 of the smallest values of $T$, from $T=0.015\pi$ to
$T=0.8\pi$.

Figure 11: Time evolution of $\delta_{\mathcal{D}}$ at the end of a long
series of cycles for six 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.

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 very weakly as $T$ decreases in very fast
loading (figure 10).

Figure 12: (a) Travel distance $\delta$ and (b) degree of mixing $\chi$
(right) 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 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 greyscale 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>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). As $T$ further decreases, $\delta$ and
$\chi$ instead begin to decrease and again approach their values without
dispersion for very fast loading ($T<~{}0.01\pi$). Across the full range of
$T$ where dispersion dominates, $\delta$ and $\chi$ increase with $A$. 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 peaks 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
naive analytical estimate $\delta_{\mathrm{est}}$ that includes both
dispersion and localisation.

Figure 14: Estimated travel distance $\delta_{\mathrm{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 (39).

Following Ref. Fiori  _et al._ (2023), the analytical solution 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”) is

$\phi_{f}(v_{f,\mathrm{vf}}-v_{s,\mathrm{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],$
(38)

where the subscript (“$\mathrm{vf}$”) stands for “very fast”. Equation (38)
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 (36) to include localisation:

$\delta_{\mathrm{est}}=C_{1}\sqrt{4\left[\mathrm{Pe}^{-1}+C_{2}\alpha\frac{A}{T}\exp\left({-(1-l)\sqrt{\frac{\pi}{T}}}\right)\right]t}.$
(39)

Figure 14 shows that $\delta_{\mathrm{est}}$ captures the qualitative trends
observed for $\delta$ in figure 12, with the central peak that is however
larger and slightly shifted compared to 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, and also the oscillatory
nature of the flow. 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$.

In summary, this analysis reveals the link between transport fluxes and
loading parameters and, consequently, with the poromechanics of periodic
loading. We identify two principal regions of the $A-T$ domain shown in Figure
12:

  * •

$T\gtrsim 0.1\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, as $T$ decreases $|v_{f}-v_{s}|$ reaches higher peaks, promoting
hydrodynamic dispersion, which is the dominant mechanism for this region.

  * •

$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.

###  III.6 Péclet number in biological examples

We conclude by reporting typical loading parameters, material properties and
diffusion coefficients for several examples of soft biological tissues (table
III.6). Based on these values, we calculate the Péclet number as defined in
section II.7, which quantifies the importance of poromechanical relaxation
relative to molecular diffusion, and an effective Péclet number
$\mathrm{Pe_{eff}}=\mathrm{Pe}A/T$, which additionally accounts for the impact
of loading on advection to understand the range of loading amplitudes and
frequencies where advection and dispersion are likely to be non-negligible.
For simplicity, we use a moderate dimensionless amplitude $A=0.1$ in all
cases. In all cases, the values of $\mathrm{Pe_{eff}}$ range from moderate to
high, indicating that advection and potentially dispersion are probably not
negligible in these systems.

Tissue |  $L$ [m] | $\tilde{A}$ |  $k_{0}$ [m2] |  $\mathcal{M}$ [Pa] |  $T_{\mathrm{pe}}$ [s] |  Typical loading frequency [Hz] | $\tilde{T}$ |  $\mathcal{D}_{m}$ [$\frac{\mathrm{m}^{2}}{\mathrm{s}}$] | $\mathrm{Pe}$ | $\mathrm{Pe_{eff}}=\mathrm{Pe}\frac{\tilde{A}}{\tilde{T}}$
---|---|---|---|---|---|---|---|---|---|---
Brain ECM Kedarasetti  _et al._ (2020) | $2\times 10^{-4}$ | $0.1$ | $2\times 10^{-15}$ | $2\times 10^{3}$ | $10$ |  $0.3$ – $10$ |  $0.003\pi$ – $0.1\pi$ | $1.4\times 10^{-10}$ | $30$ |  $3\times 10^{1}$ – $10^{3}$
Cartilage Ferguson  _et al._ (2004); Gardiner  _et al._ (2007); DiDomenico  _et al._ (2017) | $2\times 10^{-3}$ | $0.1$ | $7.5\times 10^{-18}$ | $5\times 10^{6}$ | $10^{2}$ |  $0.001$ – $0.1$ (sitting) $0.1$ – $1$ (running) |  $0.03\pi$ – $3\pi$ $0.003\pi$ – $0.03\pi$ | $5\times 10^{-12}$ | $7\times 10^{3}$ |  $10^{2}$ – $10^{3}$ $10^{3}$ – $10^{5}$
Intervertebral disk (Anulus F.) Ferguson  _et al._ (2004) | $10^{-2}$ | $0.1$ | $7.5\times 10^{-19}$ | $2.5\times 10^{6}$ | $5\times 10^{4}$ |  $2\times 10^{-5}$ (wake cycle) $0.001$ – $0.1$ (sitting) $0.1$ – $1$ (running) |  $0.3\pi$ $6\times 10^{-5}\pi$ – $6\times 10^{-3}\pi$ $6\times 10^{-6}\pi$ – $6\times 10^{-5}\pi$ | $5\times 10^{-12}$ | $4\times 10^{2}$ |  $4\times 10^{1}$ – $2\times 10^{3}$ $2\times 10^{3}$ – $2\times 10^{5}$ $2\times 10^{5}$ – $2\times 10^{6}$
Cartilage scaffold (bioreactor) Sengers  _et al._ (2004) | $2\times 10^{-3}$ | $0.1$ | $10^{-17}$ | $10^{5}$ | $4\times 10^{3}$ |  $0.001$ – $1$ |  $0.0001\pi$ – $0.1\pi$ | $10^{-9}$ | $10^{4}$ |  $4\times 10^{1}$ – $4\times 10^{4}$

##  IV Conclusions

We have derived physical insight into solute transport and mixing in a
periodically compressed soft porous material under large deformations. To do
so, we used a 1D continuum model, formulated following a rigorous nonlinear
kinematic approach and considering Hencky elasticity as the constitutive law
for the solid skeleton. Overall, we demonstrated that the characteristic
rearranging of the porous structure — resulting in a strong coupling between
mechanical loading and fluid flow — originates phenomena that are not
reversible, despite the cyclic nature of the applied load.

Our analysis is linked to a companion study, where we characterise the
mechanical response of a soft porous medium under the same loading scenarios
Fiori  _et al._ (2023). That study showed that, depending on the loading
period, the deformation can belong to either a slow-loading regime – in which
the deformation is uniform throughout the domain – or a fast-loading regime –
in which the deformation is increasingly localised near the permeable boundary
(i.e., the left boundary). Fiori  _et al._ (2023). Here, we analysed how these
different mechanical behaviours related to the loading parameters influence
solute transport and mixing.

We first studied the nature of the transport phenomena, which act in the same
direction during loading — from right to left — whereas only advection changes
direction during unloading (figure 3). Next, we focused on a single loading
cycle and compared the evolution of the solute concentration profile for four
cases: (1) molecular diffusion only (no loading); (2) advection only (no
diffusion or dispersion); (3) advection and diffusion only (no dispersion);
and (4) a general case where all the transport mechanisms are active. We found
that advection is reversible at the end of each loading cycle, and that
deformation weakly suppresses diffusion by decreasing the porosity and by
stretching the concentration gradients.

We then introduced two variables — travel distance $\delta$ and degree of
mixing $\chi$ — and compared them for different transport coefficients and
different loading parameters over several loading cycles. We showed that
diffusion and dispersion are roughly Fickian on average. For advection and
diffusion only (no dispersion), larger amplitudes increasingly reduce the
average porosity and hence increasingly suppress diffusion. This case is
essentially independent of $T$ except for very small values of $T$, where the
material is additionally compressed and $\delta$ and $\chi$ decrease as $T$
decreases. When including dispersion, solute transport increases monotonically
with $A$ but varies non-monotonically with $T$. The latter is due to the
progressive localisation of the deformation at the left boundary, which
reduces the intensity of the dispersive flux at the right boundary, where the
solutes are positioned.

Our results have two important applications: the prediction of concentration
profiles and the control of concentration profiles. The former is useful in
cases where the conditions are fixed (e.g., to predict nutrient distributions
in biological tissues in vivo), whereas the latter could enable the
optimization of the stimulation and material features to enhance or suppress
spreading and mixing, or to reach a desired final solute profile (e.g., when
designing scaffolds and bioreactors for tissue engineering).

We have taken the diffusion and dispersion coefficients themselves to be
constants. As noted above, future work should consider the impact of the
changing pore structure on the values of these material properties. It is
likely that, much like $k$, both $D_{m}$ and $\alpha$ should be deformation
dependent for moderate to large deformations.

Declaration of interests: The authors report no conflicts of interest.

###### Acknowledgements.

This work was supported by the European Research Council (ERC) under the
European Union’s Horizon 2020 Programme [Grant No. 805469]. S.P. was supported
by Start-Up Research Grant (SRG/2021/001269) by the Science and Engineering
Research Board, Department of Science and Technology, Government of India. For
the purpose of Open Access, the authors have applied a CC BY public copyright
licence to any Author Accepted Manuscript (AAM) version arising from this
submission. We thank Luis Cueto-Felgueroso for the helpful discussion on
compact finite differences.

##  Appendix A Numerical method

Our model equations are solved using finite difference methods. We use a
sixth-order compact finite difference approximation for spatial derivatives
and a third-order compact finite difference approximation for derivative
boundary conditions Lele (1992). For time integration, an implicit Runge-Kutta
method is used via MATLAB’s built-in solver ODE15s Shampine and Reichelt
(1997).

To account for the moving boundary, we work in a scaled coordinate system. We
start from a general conservation law of the form:

$\frac{\partial\Phi}{\partial t}+\frac{\partial}{\partial x}[F(\Phi)]=0.$ (40)

We then rescale the spatial coordinate $x$, introducing the variable

$\xi=\frac{x-a}{1-a},$ (41)

which maps $a\leq{}x\leq{}1$ to $0\leq{}\xi\leq{}1$. Equation 40 can then be
rewritten in terms of $\xi$ as

$\frac{\partial\Phi}{\partial
t}-\frac{(1-\xi)}{(1-a)}\dot{a}\frac{\partial\Phi}{\partial\xi}+\bigg{(}\frac{1}{1-a}\bigg{)}\frac{\partial
F(\Phi)}{\partial\xi}=0.$ (42)

When solving equation (25), we then take $\Phi=\phi_{f}$ and

${F(\phi_{f})}=-\tilde{D}_{f}(\phi_{f})\frac{\partial{\phi_{f}}}{\partial{\tilde{x}}},$
(43)

while for equation 27 we take ${\Phi=\phi_{f}\tilde{c}}$ and

$F(\phi_{f}\tilde{c})=(\phi_{f}\tilde{c})\tilde{v}_{f}-\phi_{f}\tilde{\mathcal{D}}\frac{\partial\tilde{c}}{\partial\tilde{x}}.$
(44)

For our spatial discretisation, we perform a convergence analysis in the
number of grid points $N_{x}$ (see figure 15) by calculating the Root Mean
Square (RMS) relative eror in $c(t,x=1)$ for each solution with respect the
solution for $N_{x}=1000$.

Figure 15: Convergence analysis: RMS relative error in $c(t,x=1)$ relative to
the solution for $N_{x}=1000$. (a) We fix $A=0.02$ and consider different
values of $T$, from very fast to slow. (b) We fix $T=0.1\pi$ and consider
different values of $A$, from small to large.

To balance between accuracy and computational time, we choose for all our
simulations $N_{x}=300$, with a RMS error of less than $10^{-3}$. We fix our
absolute and relative error tolerances for time integration to be $10^{-10}$.

As a reference, the cases of pure molecular diffusion and pure advection are
compared with analytical solutions for this problem (see section B), resulting
in good agreement between numerical and analytical results.

##  Appendix B Analytical solution for molecular diffusion and advection

The analytical solution for the molecular diffusion of a delta function in a
semi-infinite material, as formulated by Crank (1979), is:

$c(x,t)=c_{0}\bigg{\\{}1-\frac{1}{2}\left[\mathrm{erfc}\left(\frac{l+x-1}{2\sqrt{\mathrm{Pe}^{-1}t}}\right)+\mathrm{erfc}\left(\frac{l-x+1}{2\sqrt{\mathrm{Pe}^{-1}t}}\right)\right]\bigg{\\}}.$
(45)

In case of slow-loading, the analytical solution for advection
($\mathrm{Pe}^{-1}=\alpha=0$) is

$\tilde{c}(\tilde{x},t)=\frac{1}{2}\bigg{\\{}\tanh{\bigg{[}\tilde{s}\left((\tilde{x}-1)\frac{\phi_{f,0}}{\phi_{f}}+\tilde{l}\right)}\bigg{]}+1\bigg{\\}}.$
(46)

In figure 16, we compare these analytical solutions to the numerical solutions
obtained with the method described in Appendix A. Note that, for the case of
diffusion, we consider an initial solute profile that is smoother than the one
for which the analytical solution is formulated. Hence, the two solutions show
an initial discrepancy that vanishes with time, as the profile adjusts towards
classical self similarity.

Figure 16: Qualitative comparison between the analytical solution (solid red
curves) and numerical solution (dotted black curves) for diffusion only (a-b)
and for advection (c-d). We show (a - c) the evolution of the concentration
profiles $c(x,t)$ in time (dark to light) and (b-d) the evolution of $\delta$
over time.

##  Appendix C Dispersive flux

In this Appendix, we justify the choice of fixing $\mathrm{Pe}^{-1}$ and
$\alpha$ to the specific baseline values reported in III.2. We do so by
quantifying the strength of dispersion relative to diffusion using the ratio
of the dispersive solute flux $q_{\mathrm{disp}}$ to the diffusive solute flux
$q_{\mathrm{diff}}$ which we define as

$\frac{q_{\mathrm{disp}}}{q_{\mathrm{diff}}}=\alpha\mathrm{Pe}|{v}_{f}-{v}_{s}|,$
(47)

and which measures the relative importance of these two mechanisms. In figure
17, we compare $q_{\mathrm{disp}}/q_{\mathrm{diff}}$ for the slow-loading and
low-amplitude case and for the fast-loading and high-amplitude case. In the
former, diffusion prevails over dispersion throughout the domain, while in the
latter dispersion is the dominant mechanism. Hence, the specific baseline
values chosen for $\mathrm{Pe}^{-1}$ and $\alpha$ are such that our results
span the range from diffusion-dominated to dispersion-dominated transport
across the range of $A$ and $T$ considered.

Figure 17: (a) Lowest and (b) largest values of
$q_{\mathrm{disp}}/q_{\mathrm{diff}}$ considered in this study. The flux is
plot at ten different values of $X=x-u_{s}(X,t)$ from $0$ to $1$ (dark to
light) during one cycle. We distinguish between the loading half of the cycle
($\dot{a}>0$; solid curves) and the unloading half of the cycle ($\dot{a}<0$;
dashed curves).

##  Appendix D Initial porosity and solute strip width

Results of varying $\phi_{f,0}$ and $l$ are shown in figure 18.

Figure 18: Evolution of $\delta$ and $\chi$ over 5 cycles for different values
for (a,b) the initial porosity $\phi_{f,0}$ and (c,d) the initial solute
stripe width $l$.

Decreasing $\phi_{f,0}$ leads to larger oscillations in fluid flux for a given
solid velocity (see eq 22) In this way, both advection and hydrodynamic
dispersion are enhanced. However, the resulting impact on travel distance in
figure 18 is relatively small because dispersion is relatively weak for these
parameters.

Increasing the initial solute amount leads to a sharp increase in the travel
distance: this is generally expected because the total concentration of solute
is higher and therefore all the transport mechanisms are amplified. In
addition, the flux fluid increases monotonically in magnitude from right to
left, so larger values of l expose the solute front to stronger advection and
dispersion. However, note that a wider solute strip is associated with a lower
degree of mixing because the variance of the solute compared to the initial
variance is lower when there is more solute (see equation 35).

## References

  * Mow  _et al._ (1980) V. C. Mow, S. C. Kuei, W. M. Lai, and C. G. Armstrong, Journal of Biomechanical Engineering 102, 73 (1980).
  * Lai  _et al._ (1991) W. M. Lai, J. S. Hou, and V. C. Mow, Journal of Biomechanical Engineering 113, 245 (1991).
  * Preziosi  _et al._ (1996) L. Preziosi, D. Joseph, and G. Beavers, International Journal of Multiphase Flow 22, 1205 (1996).
  * Li  _et al._ (2004) C. Li, R. I. Borja, and R. A. Regueiro, Computer Methods in Applied Mechanics and Engineering 193, 3837 (2004).
  * Franceschini  _et al._ (2006) G. Franceschini, D. Bigoni, P. Regitnig, and G. Holzapfel, Journal of the Mechanics and Physics of Solids 54, 2592 (2006).
  * Moeendarbary  _et al._ (2013) E. Moeendarbary, L. Valon, M. Fritzsche, A. R. Harris, D. A. Moulding, A. J. Thrasher, E. Stride, L. Mahadevan, and G. T. Charras, Nature Materials 12, 253 (2013).
  * Ehlers  _et al._ (2009) W. Ehlers, N. Karajan, and B. Markert, Biomechanics and Modeling in Mechanobiology 8, 233 (2009).
  * Vuong  _et al._ (2015) A. T. Vuong, L. Yoshihara, and W. A. Wall, Computer Methods in Applied Mechanics and Engineering 283, 1240 (2015).
  * Borja and Choo (2016) R. I. Borja and J. Choo, Computer Methods in Applied Mechanics and Engineering 309, 653 (2016).
  * Fiori  _et al._ (2023) M. Fiori, S. Pramanik, and C. W. MacMinn, Journal of Fluid Mechanics 974, A2 (2023).
  * Bear (1972) J. Bear,  _Dynamics of fluids in porous media_ (American Elsevier Publishing Company, 1972).
  * Whitaker (1967) S. Whitaker, AIChE Journal 13, 420 (1967).
  * Saffman (1959) P. G. Saffman, Journal of Fluid Mechanics 6, 321 (1959).
  * Scheidegger (1961) A. E. Scheidegger, Journal of Geophysical Research 66, 3273 (1961).
  * Brenner and Edwards (1993) H. Brenner and D. Edwards,  _Macrotransport Processes_ , Butterworth-Heinemann series in chemical engineering (Elsevier Science, 1993).
  * Whitaker (1998) S. Whitaker,  _The Method of Volume Averaging_ , Theory and Applications of Transport in Porous Media (Springer Netherlands, 1998).
  * Gelhar (1993) L. Gelhar,  _Stochastic Subsurface Hydrology_ (Prentice-Hall, 1993).
  * Dentz  _et al._ (2011) M. Dentz, T. Le Borgne, A. Englert, and B. Bijeljic, Journal of Contaminant Hydrology 120-121, 1 (2011).
  * Ghanbarian  _et al._ (2013) B. Ghanbarian, A. G. Hunt, R. P. Ewing, and M. Sahimi, Soil science society of America journal 77, 1461 (2013).
  * Tartakovsky and Dentz (2019) D. M. Tartakovsky and M. Dentz, Transport in Porous Media 130, 105 (2019).
  * Taylor (1953) G. I. Taylor, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 219, 186 (1953).
  * Brenner and Stewartson (1980) H. Brenner and K. Stewartson, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 297, 81 (1980).
  * Marbach and Alim (2019) S. Marbach and K. Alim, Physical Review Fluids 4, 114202 (2019).
  * Datta  _et al._ (2013) S. S. Datta, H. Chiang, T. S. Ramakrishnan, and D. A. Weitz, Physical Review Letters 111, 064501 (2013).
  * de Anna  _et al._ (2013) P. de Anna, T. Le Borgne, M. Dentz, A. M. Tartakovsky, D. Bolster, and P. Davy, Physical Review Letters 110, 184502 (2013).
  * Lester  _et al._ (2013) D. R. Lester, G. Metcalfe, and M. G. Trefry, Physical Review Letters 111, 1 (2013).
  * Lester  _et al._ (2016a) D. R. Lester, M. Dentz, and T. Le Borgne, Journal of Fluid Mechanics 803, 144–174 (2016a).
  * Lester  _et al._ (2016b) D. R. Lester, M. G. Trefry, and G. Metcalfe, Advances in Water Resources 97, 175 (2016b).
  * Gouze  _et al._ (2021) P. Gouze, A. Puyguiraud, T. Porcher, and M. Dentz, Frontiers in Water 3 (2021).
  * Dentz  _et al._ (2023) M. Dentz, J. J. Hidalgo, and D. R. Lester, Transport in Porous Media 146, 5 – 53 (2023).
  * Kree and Villermaux (2017) M. Kree and E. Villermaux, Physical Review Fluids 2, 104502 (2017).
  * Souzy  _et al._ (2020) M. Souzy, H. Lhuissier, Y. Méheust, T. Le borgne, and B. Metzger, Journal of Fluid Mechanics 891 (2020).
  * Gelhar and Axness (1983) L. W. Gelhar and C. L. Axness, Water Resources Research 19, 161 (1983).
  * Gelhar  _et al._ (1992) L. W. Gelhar, C. Welty, and K. R. Rehfeldt, Water Resources Research 28, 1955 (1992).
  * Delgado (2007) J. Delgado, Chemical Engineering Research and Design 85, 1245 (2007).
  * Fox (2007a) P. J. Fox, Journal of Geotechnical and Geoenvironmental Engineering 133, 3 (2007a).
  * Fox (2007b) P. J. Fox, Journal of Geotechnical and Geoenvironmental Engineering 133, 16 (2007b).
  * Pu  _et al._ (2018) H. Pu, P. J. Fox, and C. D. Shackelford, Journal of Geotechnical and Geoenvironmental Engineering 144, 1 (2018).
  * Peters and Smith (2002) G. P. Peters and D. W. Smith, International Journal for Numerical and Analytical Methods in Geomechanics 26, 683 (2002).
  * Lewis  _et al._ (2009) T. W. Lewis, P. Pivonka, and D. W. Smith, International Journal for Numerical and Analytical Methods in Geomechanics 33, 95 (2009).
  * Alshawabkeh and Rahbar (2006) A. N. Alshawabkeh and N. Rahbar, Journal of Geotechnical and Geoenvironmental Engineering 132, 1001 (2006).
  * Zhang  _et al._ (2013) H. J. Zhang, D. S. Jeng, D. A. Barry, B. R. Seymour, and L. Li, Journal of Hydrology 479, 189 (2013).
  * Bonazzi  _et al._ (2021) A. Bonazzi, B. Jha, and F. P. J. de Barros, International Journal for Numerical and Analytical Methods in Geomechanics 45, 307 (2021).
  * Tran and Jha (2020) M. Tran and B. Jha, Advances in Water Resources 136, 103485 (2020).
  * Goldsztein and Santamarina (2004) G. H. Goldsztein and J. C. Santamarina, Applied Physics Letters 85, 2432 (2004).
  * Claria  _et al._ (2012) J. J. Claria, G. H. Goldsztein, and J. C. Santamarina, Transport in Porous Media 93, 51 (2012).
  * Wu  _et al._ (2020) J. Wu, D. R. Lester, M. G. Trefry, and G. Metcalfe, Water Resources Research 56, 1 (2020).
  * Trefry  _et al._ (2019) M. G. Trefry, D. R. Lester, G. Metcalfe, and J. Wu, Water Resources Research 55, 3347 (2019).
  * Ferguson  _et al._ (2004) S. J. Ferguson, K. Ito, and L. J. Pyrak-Nolte, Journal of Biomechanics 37, 213 (2004).
  * Gardiner  _et al._ (2007) B. Gardiner, D. Smith, P. Pivonka, A. Grodzinsky, E. Frank, and L. Zhang, Computer Methods in Biomechanics and Biomedical Engineering 10, 265 (2007).
  * Witt  _et al._ (2014) F. Witt, G. N. Duda, C. Bergmann, and A. Petersen, Tissue Engineering - Part A 20, 486 (2014).
  * Schmidt  _et al._ (2010) H. Schmidt, A. Shirazi-Adl, F. Galbusera, and H.-J. Wilke, Journal of Biomechanics 43, 1849 (2010).
  * Zhang (2011) L. Zhang, International Journal of Applied Mechanics 3, 507 (2011).
  * DiDomenico  _et al._ (2017) C. D. DiDomenico, Z. X. Wang, and L. J. Bonassar, Journal of Biomechanical Engineering 139, 1 (2017).
  * Zhang and Szeri (2008) L. Zhang and A. Z. Szeri, Journal of Biomechanics 41, 430 (2008).
  * Albro  _et al._ (2008) M. B. Albro, N. O. Chahine, R. Li, K. Yeager, C. T. Hung, and G. A. Ateshian, Journal of Biomechanics 41, 3152 (2008).
  * Vaughan  _et al._ (2013) B. L. Vaughan, P. A. Galie, J. P. Stegemann, and J. B. Grotberg, Biophysical Journal 105, 2188 (2013).
  * Mauck  _et al._ (2003) R. L. Mauck, C. T. Hung, and G. A. Ateshian, Journal of Biomechanical Engineering 125, 602 (2003).
  * Cortez  _et al._ (2016) S. Cortez, A. Completo, and J. L. Alves, ECCOMAS Congress 2016 - Proceedings of the 7th European Congress on Computational Methods in Applied Sciences and Engineering 1, 27 (2016).
  * Kumar  _et al._ (2018) P. Kumar, B. Dey, and G. P. Raja Sekhar, International Journal of Engineering Science 127, 201 (2018).
  * Fan  _et al._ (2016) L. Fan, S. Pei, X. Lucas Lu, and L. Wang, Bone Research 4, 16032 (2016).
  * Urciuolo  _et al._ (2008) F. Urciuolo, G. Imparato, and P. A. Netti, AIChE Journal 54, 824 (2008).
  * Sengers  _et al._ (2004) B. G. Sengers, C. W. J. Oomens, and F. P. T. Baaijens, Journal of Biomechanical Engineering 126, 82 (2004).
  * Kelley and Thomas (2023) D. H. Kelley and J. H. Thomas, Annual Review of Fluid Mechanics 55, 237 (2023).
  * Sharp  _et al._ (2019) K. M. Sharp, R. o. Carare, and A. Martin, Bryn, Fluids and Barriers of the CNS 16 (2019).
  * Berg  _et al._ (2020) M. Berg, Y. Davit, M. Quintard, and S. Lorthois, Journal of Fluid Mechanics 884, A39 (2020).
  * Troyetsky  _et al._ (2021) D. E. Troyetsky, J. Tithof, J. H. Thomas, and D. H. Kelley, Scientific Reports 11 (2021).
  * Bojarskaite  _et al._ (2023) L. Bojarskaite, D. M. Bjørnstad, A. Vallet, K. M. Gullestad Binder, C. Cunen, K. Heuser, M. Kuchta, K.-A. Mardal, and R. Enger, Nature Communications 14, 953 (2023).
  * Zimmerman and Tartakovsky (2020) R. A. Zimmerman and D. M. Tartakovsky, Journal of Fluid Mechanics 901, A24 (2020).
  * Goirand  _et al._ (2021) F. Goirand, T. Le Borgne, and S. Lorthois, Nature Communications 12 (2021).
  * Davit  _et al._ (2013) Y. Davit, H. Byrne, J. Osborne, J. Pitt-Francis, D. Gavaghan, and M. Quintard, Physical Review E 87, 012718 (2013).
  * Oswald and Kinzelbach (2004) S. Oswald and W. Kinzelbach, Journal of Hydrology 290, 22 (2004).
  * Liang  _et al._ (2018) Y. Liang, B. Wen, M. A. Hesse, and D. DiCarlo, Geophysical Research Letters 45, 9690 (2018).
  * Mow  _et al._ (1984) V. C. Mow, M. H. Holmes, and W. M. Lai, Journal of Biomechanics 17, 377 (1984).
  * Maroudas (1970) A. Maroudas, Biophysical Journal 10, 365 (1970).
  * DiDomenico  _et al._ (2018) C. D. DiDomenico, M. Lintz, and L. Bonassar, Nature Reviews Rheumatology 14, 393 (2018).
  * Ateshian and Weiss (2013) G. A. Ateshian and J. A. Weiss, in  _Computer Models in Biomechanics: From Nano to Macro_ (Springer, 2013) pp. 231–249.
  * Yao and Gu (2007) H. Yao and W. Y. Gu, Biomechanics and Modeling in Mechanobiology 6, 63 (2007).
  * Maroudas (1975) A. Maroudas, Biorheology 12, 233 (1975).
  * Levick (1987) J. Levick, Quarterly Journal of Experimental Physiology: Translation and Integration 72, 409 (1987).
  * Federico and Herzog (2008) S. Federico and W. Herzog, International Journal of Solids and Structures 45, 2160 (2008).
  * DiDomenico  _et al._ (2017) C. D. DiDomenico, A. Goodearl, A. Yarilina, V. Sun, S. Mitra, A. S. Sterman, and L. J. Bonassar, Journal of Biomechanical Engineering 139, 091005 (2017).
  * Buijs  _et al._ (2010) J. O. D. Buijs, E. L. Ritman, and D. Dragomir-Daescu, Tissue Engineering - Part C: Methods 16, 1145 (2010).
  * Lee  _et al._ (2015) M. K. Lee, M. H. Rich, K. Baek, J. Lee, and H. Kong, Scientific Reports 5, 1 (2015).
  * Mesallati  _et al._ (2013) T. Mesallati, C. T. Buckley, T. Nagel, and D. J. Kelly, Biomechanics and Modeling in Mechanobiology 12, 889 (2013).
  * Hrabe  _et al._ (2004) J. Hrabe, S. Hrabĕtová, and K. Segeth, Biophysical journal 87, 1606 (2004).
  * Zhang and Szeri (2005) L. Zhang and A. Szeri, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 461, 2021 (2005).
  * Hencky (1931) H. Hencky, Journal of Rheology 2, 169 (1931).
  * MacMinn  _et al._ (2016) C. W. MacMinn, E. R. Dufresne, and J. S. Wettlaufer, Physical Review Applied 5, 044020 (2016).
  * Auton and MacMinn (2018) L. C. Auton and C. W. MacMinn, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 474, 20180284 (2018).
  * Dentz  _et al._ (2018) M. Dentz, F. P. De Barros, T. Le Borgne, and D. R. Lester, Journal of Fluid Mechanics 853, 621 (2018).
  * Jha  _et al._ (2011) B. Jha, L. Cueto-Felgueroso, and R. Juanes, Physical Review E 84, 066312 (2011).
  * Tan and Homsy (1988) C. T. Tan and G. M. Homsy, The Physics of Fluids 31, 1330 (1988).
  * Mishra  _et al._ (2008) M. Mishra, M. Martin, and A. De Wit, Physical Review E 78, 066306 (2008).
  * Kedarasetti  _et al._ (2020) R. T. Kedarasetti, P. J. Drew, and F. Costanzo, Fluids and Barriers of the CNS 19, 34 (2020).
  * Lele (1992) S. K. Lele, Journal of computational physics 103, 16 (1992).
  * Shampine and Reichelt (1997) L. F. Shampine and M. W. Reichelt, SIAM Journal on Scientific Computing 18, 1 (1997).
  * Crank (1979) J. Crank,  _The mathematics of diffusion_ (Oxford university press, 1979).