# Solute transport due to periodic loading in a soft porous material

Matilde Fiori, Satyajit Pramanik, Christopher W Macminn

► To cite this version:

Matilde Fiori, Satyajit Pramanik, Christopher W Macminn. Solute transport due to periodic loading in a soft porous material. Journal of Fluid Mechanics, 2025, 1009, pp.A15. 10.1017/jfm.2025.189 . hal-05050999

HAL Id: hal-05050999

https://hal.science/hal-05050999v1

Submitted on 29 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
J M

PAPERS

<figure><img src="image_2.png" /></figure>

Solute transport due to periodic loading in a soft
porous material

Matilde Fiori<fnref n="1" />,<fnref n="2" />, Satyajit Pramanik<fnref n="3" /> ID and Christopher W. MacMinn<fnref n="1" /> ID

**Corresponding author:** Christopher W. MacMinn, christopher.macminn@eng.ox.ac.uk

(Received 23 February 2024; revised 31 December 2024; accepted 1 February 2025)

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 including the geomechanics of contaminants in the subsurface and the biomechanics of nutrient transport in living tissues and scaffolds for tissue engineering. However, the basic features of this process have not previously been systematically investigated. Here, we fill this hole in the context of a one-dimensional 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 time scale – and amplitudes – from infinitesimal to large. We show that solute spreading over several loading cycles 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.

**Key words:** biomedical flows, dispersion, porous media

<footnote n="1">¹Department of Engineering Science, University of Oxford, Oxford OX1 3PJ, UK</footnote>

<footnote n="2">²Institut de Mécanique Des Fluides de Toulouse, IMFT, Université de Toulouse, CNRS, Toulouse 31400, France</footnote>

<footnote n="3">³Department of Mathematics, Indian Institute of Technology Guwahati, Guwahati, Assam 781039, India</footnote>
# 1. 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, Hou & Mow 1991; Preziosi, Joseph & Beavers 1996; Li, Borja & Regueiro 2004; Franceschini et al. 2006; Ehlers, Karajan & Markert 2009; Moeendarbary et al. 2013; Vuong, Yoshihara & Wall 2015; Borja & 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, Pramanik & MacMinn 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 (Saffman 1959; Scheidegger 1961; Whitaker 1967; Bear 1972; Brenner & Edwards 1993; Gelhar 1993; Whitaker 1998; 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 & 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) analogously to classical Taylor dispersion in a tube (Taylor 1953; Brenner & Stewartson 1980; Marbach & 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 (de Anna et al. 2013; Datta et al. 2013; Lester, Metcalfe & Trefry 2013, 2016a,b; Kree & Villermaux 2017; Gouze et al. 2021; Dentz, Hidalgo & Lester 2023; Souzy et al. 2020), thus inducing both longitudinal and transverse spreading (Scheidegger 1961; Gelhar & Axness 1983; Gelhar, Welty & Rehfeldt 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 in several different contexts. The impact of monotonic soil consolidation on transport has been studied extensively for its relevance to the management of landfills and other contaminated sediments, such as dredging or mining waste (e.g. Smith 2000; Peters & Smith 2002; Alshawabkeh & Rahbar 2006; Fox 2007a,b; Arega & Hayter 2008; Lewis, Pivonka & Smith 2009; Zhang et al. 2012, 2013; Xie et al. 2016; Pu, Fox & Shackelford 2018; Bonazzi, Jha & de Barros 2021). In that context, it is well known that consolidation enhances solute transport. Deformation has also been shown to increase mixing and reduce breakthrough time in the context of miscible viscous fingering (Tran & 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 one-dimensional (1-D) pore networks (Goldsztein & Santamarina 2004; Claria, Goldsztein & Santamarina 2012).
In a poroelastic material, solute transport due to small periodic deformations has been explored across a range of parameters, including compressibility and forcing frequency, for semi-infinite homogeneous systems (Pool, Dentz & Post 2016), finite homogeneous systems (Bonazzi et al. 2021) and finite heterogeneous systems (Trefry et al. 2019; Wu et al. 2020). The latter two studies focus in particular on the combined role of poroelasticity, heterogeneity and transient forcing in generating chaotic advection.

Periodic loading is also known to enhance the transport of nutrients in biological tissues (Ferguson, Ito & Pyrak-Nolte 2004; Gardiner et al. 2007; Zhang & Szeri 2008; Schmidt et al. 2010; Zhang 2011; Witt et al. 2014; DiDomenico et al. 2017). 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, Hung & Ateshian 2003; Cortez, Completo & Alves 2016; Fan et al. 2016; Kumar, Dey & Sekhar 2018), 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, Imparato & Netti 2008) and amplifies the role of hydrodynamic dispersion (Sengers, Oomens & Baaijens 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; Vaughan et al. 2013; DiDomenico et al. 2017).

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 Sengers et al. 2004). One context where dispersion is widely agreed to be important is in brain microcirculation (Kelley & 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, Borgne & Lorthois 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 Péclet number $Pe = VL/D_m$ is assumed to be small, where V is the characteristic fluid velocity, L the characteristic streamwise length scale and $D_m$ the molecular diffusivity. However, it is straightforward to show that Pe can be order 1 or larger in a tissue or gel subject to fast (0.1–1 Hz) and large (10–20 %) deformations (see table 2), 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; Vaughan et al. 2013; DiDomenico et al. 2017). 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 & Kinzelbach 2004; Kree & Villermaux 2017; Liang et al. 2018). Indeed, the typical pore size is ~ 10 nm in polymeric gels and in the extra-cellular matrix of tissues (e.g. around 6 nm in cartilage, Mow, Holmes & Lai (1984)) and can therefore be similar to (or smaller than) the size of large solute molecules (Maroudas 1970; DiDomenico, Lintz & Bonassar 2018), causing
solid-solute friction (Yao & Gu 2007; Ateshian et al. 2011). However, tissues and scaffolds are heterogeneous and multiscale materials; the presence of other components, such as collagen fibres, results in a ‘mesoscale’ of larger pores (e.g. 100–150 nm in cartilage, Maroudas (1975); Levick (1987); Federico & 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, Ritman & Dragomir-Daescu 2010; Mesallati et al. 2013; Lee et al. 2015). As a further counter-argument, we hypothesise that, even in pores that are small compared with the solute molecules, the irrelevance of pore-scale velocity gradients does not exclude velocity variations and streamline alterations in the overall network, which could cause longitudinal and transverse dispersion. This hypothesis is consistent with the quantification of tortuosity in several soft tissues (Maroudas 1970; Hrabe, Hrabětová & Segeth; 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 1-D, large-deformation poroelasticity model that includes rigorous nonlinear kinematics, deformation-dependent permeability and Hencky elasticity for the solid skeleton. In a companion study, we used this model to explore the poromechanics of large-amplitude periodic loading (Fiori et al. 2023). Here, we additionally consider solute transport due to advection, diffusion and dispersion. We first 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 numbers 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.

## 2. 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 caused by the former and responsible for the latter.

### 2.1. *Model problem*

We consider a 1-D sample of soft porous material of relaxed length *L* and relaxed porosity (fluid fraction) *φ*<sub>*f,0*</sub>. 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.
<figure><img src="image_3.png" /><figcaption>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).</figcaption></figure>

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], \qquad (2.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 ≤ A/L ≤ 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).

## 2.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 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}_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)\hat{\mathbf{e}}_x, \quad \mathbf{v}_s = v_s(x, t)\hat{\mathbf{e}}_x, \quad \mathbf{v}_f = v_f(x, t)\hat{\mathbf{e}}_x, \qquad (2.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 $\hat{\mathbf{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} \rightarrow \frac{\partial u_s}{\partial x} = \frac{\phi_f - \phi_{f,0}}{1 - \phi_{f,0}}. \qquad (2.3)$$

Continuity can be written

$$\frac{\partial \phi_f}{\partial t} + \frac{\partial}{\partial x}(\phi_f v_f) = 0 \quad \text{and} \quad \frac{\partial \phi_f}{\partial t} - \frac{\partial}{\partial x}[(1 - \phi_f) v_s] = 0, \qquad (2.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$.

## 2.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}, \qquad (2.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 Fiori et al. (2023), 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. (2023).

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

$$ \frac{\partial \phi_f}{\partial t} + \frac{\partial}{\partial x} \left[ \phi_f q - (1 - \phi_f) \frac{k(\phi_f)}{\mu} \frac{\partial p}{\partial x} \right] = 0 \quad \text{and} \quad \frac{\partial q}{\partial x} = 0, \qquad (2.6) $$

where the total flux $q$ is again

$$ q = \phi_f v_f + (1 - \phi_f) v_s, \qquad (2.7) $$

and the fluid and solid velocities are given by

$$ v_f = q - \frac{(1 - \phi_f) k(\phi_f)}{\mu} \frac{\partial p}{\partial x} \quad \text{and} \quad v_s = q + \frac{k(\phi_f)}{\mu} \frac{\partial p}{\partial x}. \qquad (2.8) $$

Note that the fluid flux is

$$ q_f = \phi_f v_f. \qquad (2.9) $$

## 2.4. Mechanical equilibrium and elasticity law

Neglecting inertia, gravity and other body forces, mechanical equilibrium can be expressed as $\nabla \cdot \boldsymbol{\sigma} = \nabla \cdot \boldsymbol{\sigma}' - \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

$$ \frac{\partial \boldsymbol{\sigma}'}{\partial x} = \frac{\partial p}{\partial x}, \qquad (2.10) $$

where $\boldsymbol{\sigma}'$ is the $xx$ component of $\boldsymbol{\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 $\boldsymbol{\sigma}' = \boldsymbol{\sigma}'(\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} \left[ \phi_f q - D_f(\phi_f) \frac{\partial \phi_f}{\partial x} \right] = 0 \quad \text{and} \quad \frac{\partial q}{\partial x} = 0, \qquad (2.11) $$

where the nonlinear composite constitutive function

$$ D_f(\phi_f) = (1 - \phi_f) \frac{k(\phi_f)}{\mu} \frac{d\boldsymbol{\sigma}'}{d\phi_f} \qquad (2.12) $$

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. Davis and Russel 1989; Peppin, Elliott & Worster 2006; Style & Peppin 2011; Bouchaudy and Salmon 2019; Worster, Peppin & Wettlaufer 2021).

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

$$ \sigma' = M \frac{\ln(J)}{J} = M \left( \frac{1 - \phi_f}{1 - \phi_{f,0}} \right) \ln \left( \frac{1 - \phi_{f,0}}{1 - \phi_f} \right), \qquad (2.13) $$

where $M$ is the p-wave or oedometric modulus (MacMinn, Dufresne & Wettlaufer 2016). 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. 2016). We motivate and discuss our constitutive choices and explore the poroelastic diffusivity in more detail in Fiori et al. (2023).

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

$$ \frac{\partial}{\partial t}(\phi_f c) + \frac{\partial}{\partial x} \left[ \phi_f c v_f - \phi_f D \frac{\partial c}{\partial x} \right] = 0. \qquad (2.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 $D$ can be written

$$ D = D_m + D_h, \qquad (2.15) $$

where $D_m$ and $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 one dimension by

$$ D_h = \alpha |v_f - v_s|, \qquad (2.16) $$

where $\alpha$ is the longitudinal dispersivity (Scheidegger 1961; Brenner & Edwards 1993; Gelhar 1993; Whitaker 1998). 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. 2016, 2017) and compaction has been shown to have a non-trivial impact on dispersion in bead packs and packed beds (Charlaix, Hulin & Plona 1987; Östergren & Trägårdh 2000; Liu et al. 2024). 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

$$ \phi_f(x, 0) = \phi_{f,0} \text{ and } u_s(x, 0) = 0. \qquad (2.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} \{\tanh[s(x - L + l)] + 1\}, \qquad (2.18) $$

where $s$ is a steepness parameter.

### 2.6.2. Left boundary

For $t > 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

$$ u_s(a, t) = a(t), \quad v_s(a, t) = \frac{da}{dt} \quad \text{and} \quad p(a, t) = 0. \qquad (2.19) $$

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

$$ c(a, t) = 0. \qquad (2.20) $$

### 2.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 \text{and} \quad \left. \frac{\partial c}{\partial x} \right|_{x=L} = 0. \qquad (2.21) $$

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

$$ v_f = -\frac{(1 - \phi_f)}{\phi_f} v_s. \qquad (2.22) $$
## 2.7. Scaling and summary

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

$$ \tilde{x} = \frac{x}{L}, \quad \tilde{u}_s = \frac{u_s}{L}, \quad \tilde{t} = \frac{t}{T_{pe}}, \quad \tilde{\sigma}' = \frac{\sigma'}{M}, \quad \tilde{p} = \frac{p}{M}, \quad \tilde{k} = \frac{k(\phi)}{k_0}, \quad \tilde{v}_f = \frac{v_f}{L/T_{pe}}, \quad \tilde{v}_s = \frac{v_s}{L/T_{pe}}, \quad (2.23) $$

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

We then scale quantities related to solute transport as

$$ \tilde{c} = \frac{c}{c_0}, \quad \tilde{l} = \frac{l}{L}, \quad \tilde{\alpha} = \frac{\alpha}{L}. \quad (2.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}} \left[ \tilde{D}_f(\phi_f) \frac{\partial \phi_f}{\partial \tilde{x}} \right] = 0, \quad (2.25) $$

where

$$ \tilde{D}_f = \frac{D_f}{D_{f,0}} = (1 - \phi_f) \tilde{k}(\phi_f) \frac{d\tilde{\sigma}'}{d\phi_f}, \quad (2.26) $$

and

$$ \frac{\partial}{\partial \tilde{t}}(\phi_f \tilde{c}) + \frac{\partial}{\partial \tilde{x}} \left[ \phi_f \tilde{c} \tilde{v}_f - \phi_f \tilde{D} \frac{\partial \tilde{c}}{\partial \tilde{x}} \right] = 0. \quad (2.27) $$

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

$$ \tilde{D} = \frac{\mathcal{D}}{D_m} = \text{Pe}^{-1} + \tilde{\alpha} |\tilde{v}_f - \tilde{v}_s|, \quad (2.28) $$

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

The initial conditions are

$$ \tilde{a}(0) = 0, \quad \phi_f(\tilde{x}, 0) = \phi_{f,0}, \quad (2.29) $$

and

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

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

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

and

$$ \tilde{u}_s(1, \tilde{t}) = \tilde{v}_s(1, \tilde{t}) = \tilde{v}_f(1, \tilde{t}) = 0 \quad \text{and} \quad \left. \frac{\partial \tilde{c}}{\partial \tilde{x}} \right|_{\tilde{x}=1} = 0, \quad (2.32) $$

where $\tilde{A} = A/L$ and $\tilde{T} = T/T_{pe}$.
M. Fiori, S. Pramanik and C.W. MacMinn

As shown in Fiori et al. (2023), 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

$$T_{adv} = \frac{L}{v_f^*} = \frac{LT\phi_{f,0}}{2A(1 - \phi_{f,0})} \quad \rightarrow \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}}, \quad (2.33)$$

and

$$T_{diff} = \frac{L^2}{D_m} \quad \rightarrow \quad \tilde{T}_{diff} = \frac{T_{diff}}{T_{pe}} = \frac{D_{f,0}}{D_m} = Pe. \quad (2.34)$$

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 $Pe_{eff}$,

$$Pe_{eff} = Pe \frac{\tilde{A}}{\tilde{T}} \propto \frac{\tilde{T}_{diff}}{\tilde{T}_{adv}}. \quad (2.35)$$

In our results below, we explore a wide range of $Pe_{eff}$: from ~1 to ~$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 (2025). Below, we consider only dimensionless quantities, dropping the tildes for convenience.

## 3. Solute transport and mixing

### 3.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:

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

(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 1952; Jha, Cueto-Felgueroso & Juanes 2011) to account for a porosity field that varies in space. Considering the fluid-volume-weighted average $\langle \ast \rangle_f$,
<figure><img src="image_4.png" /><figcaption>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 1988; Mishra, Martin & De Wit 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$.</figcaption></figure>

defined as

$$ \langle * \rangle_f = \frac{\int_a^1 \phi_f * dx}{\int_a^1 \phi_f dx}, \qquad (3.1) $$

the variance of the concentration distribution is then

$$ \sigma^2(t) = \langle c^2 \rangle_f - \langle c \rangle_f^2, \qquad (3.2) $$

and the degree of mixing is

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

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 $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 $Pe^{-1}$ and $\alpha$ in § 3.4, $A$ and $T$ in § 3.5 and $\phi_{f,0}$ and $l$ in Appendix D.

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. 2023), 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 $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
<table><caption>Table 1. Baseline parameter values.</caption><thead><tr><th>Parameter</th><th>A</th><th>T</th><th>Pe<sup>-1</sup></th><th>α</th><th>φ<sub>f,0</sub></th><th>l</th></tr></thead><tbody><tr><td>Baseline value</td><td>0.1</td><td>6π</td><td>3 × 10<sup>-5</sup></td><td>0.01</td><td>0.75</td><td>0.05</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>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$, $α = 0.025$. Note that A and α are higher than the baseline values to better illustrate the roles of advection and dispersion. The solid grey envelope is proportional to $t^{-1/2}$.</figcaption></figure>

slowest period considered in this study (see Appendix C). The baseline value for $φ_{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

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} = φ_f v_f c$, $q_{diff} = φ_f Pe^{-1} (\partial c / \partial x)$ and $q_{disp} = φ_f α |v_f - v_s| (\partial c / \partial x)$, respectively. 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 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 strengths of diffusion and dispersion are proportional to $Pe^{-1}$ and $α$, respectively.
<figure><img src="image_6.png" /></figure>

Figure 4. Evolution of the concentration profile during one cycle (red to blue through white) for four cases: (a) diffusion only ($A = \alpha = 0, Pe^{-1} = 3 \times 10^{-5}$); (b) advection only ($A = 0.4, Pe^{-1} = \alpha = 0$); (c) advection and diffusion ($A = 0.4, 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} > 0$, first half of the cycle, dark to light red) and unloading ($\dot{a} < 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 $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 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 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 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
<figure><img src="image_7.png" /></figure>

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

<figure><img src="image_8.png" /></figure>

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

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.

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

We next isolate the roles of $Pe^{-1}$ and $\alpha$. For that purpose, we focus on five cycles and consider a wide range of $Pe^{-1}$ and $\alpha$. We consider the roles of $\phi_{f,0}$ and $l$ in Appendix D. The impact of changing $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 $Pe^{-1}$ and $\alpha$ enhance the diffusive and dispersive fluxes, respectively, and hence drive faster spreading, as should be expected. For sufficiently small $Pe^{-1}$, dispersion dominates diffusion and the rate of spreading becomes independent of $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{D_{eff}}$ for some effective diffusion/dispersion coefficient $D_{eff} \approx 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. 2023). Thus, we assume that, on average, $|v_f - v_s| \sim C_2 A/T$ for some constant $C_2$. The resulting estimate is then

$$ \delta_{est,SL} = C_1 f \left( \frac{\langle \phi_f \rangle}{\phi_{f,0}}, \text{Pe}^{-1} \right) \sqrt{4 \left( \text{Pe}^{-1} + C_2 \alpha \frac{A}{T} \right) t}, \quad (3.4) $$

where $C_1$ is a constant and the function $f(\langle \phi_f \rangle/\phi_{f,0}, \text{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 $\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

$$ f \left( \frac{\langle \phi_f \rangle}{\phi_{f,0}}, \text{Pe}^{-1} \right) = 1 + C_3 \ln(\text{Pe}^{-1}) \left( 1 - \frac{\langle \phi_f \rangle}{\phi_{f,0}} \right), \quad (3.5) $$

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 $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 $Pe^{-1}$ in figure 5 because it is a weak function of $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. 2023), 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 \ll 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 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. (2023)). For very large values of $T$ (e.g. figure 7d, $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 7b, 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 8a–h) and the normalised time-averaged value of $|v_f - v_s|$ is relatively insensitive to $A$ (figure 8i), suggesting that $|v_f - v_s|$ is essentially proportional to $A$. This suggestion is confirmed in figure 8(j).

### 3.5.2. Solute transport for different loading amplitudes and periods

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_D$ and the degree of mixing $\chi_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. (2023), the overall average porosity $\langle \phi_f \rangle$ over any integer number
<figure><img src="image_9.png" /><figcaption>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) 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 A/(\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)).</figcaption></figure>

of cycles is given by

$$ \overline{\langle \phi_f \rangle} = \frac{1}{mT} \int_{nT}^{(n+m)T} \langle \phi_f \rangle \, dt = 1 - \frac{1 - \phi_{f,0}}{\sqrt{1 - A}}, \quad (3.6) $$

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_D$ and $\chi_D$ in
<figure><img src="image_10.png" /><figcaption>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} > 0$; solid curves) and the unloading half of the cycle ($\dot{a} < 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$.</figcaption></figure>

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}/\langle\phi_f\rangle$.

In figure 10(a), we fix $T = 12\pi$ and plot $\delta_D$, $\chi_D$ and $\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_D/\delta_{D,A=0}$ against $A$ for several values of $\phi_{f,0}$ (figure 10b) and $\text{Pe}^{-1}$ (figure 10c). Panel (b) shows that increasing $\phi_{f,0}$ strongly mitigates the decay with $A$, while panel (c) shows that increasing $\text{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(\langle\phi_f\rangle/\phi_{f,0}, \text{Pe}^{-1})$. Thus, the curves in figures 10(b) and 10(c) correspond to profiles
<figure><img src="image_11.png" /><figcaption>Figure 9. (a) Travel distance $\delta_D$ and (b) degree of mixing $\chi_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).</figcaption></figure>

of $f$ against $A$ for different values of $\phi_{f,0}$ and $\text{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 $\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 $\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 $\text{Pe}^{-1}$.

The dependence of $f$ on $\text{Pe}^{-1}$ is weaker. As indicated in equation (3.5) and illustrated in figure 10(c), we find that a logarithmic variation with $\text{Pe}^{-1}$ provides reasonable qualitative agreement between our expression for $f$ and our numerical results, with good quantitative agreement for larger values of $\text{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 $\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 $\langle\phi_f\rangle/\phi_{f,0}$) will more strongly hinder diffusion (decreasing $f$). The role of $\text{Pe}^{-1}$ is more difficult to rationalise. We find that stretching of concentration gradients more strongly hinders diffusion at smaller values of $\text{Pe}^{-1}$, corresponding to weaker diffusion. At smaller $\text{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.

Figures 9 and 10 confirm that $\delta_D$ and $\chi_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_D$ and $\chi_D$ increase with $T$. This effect is not related to $\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_D$ over time for a fixed amplitude $A = 0.06$ and for several values of $T$, from
<figure><img src="image_12.png" /><figcaption>Figure 10. (a) Variation of $\delta_D$ (red), $\chi_D$ (blue) and $\overline{(\phi_f)}$ (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_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 $\text{Pe}^{-1} \in [3 \times 10^{-6}, 3 \times 10^{-4}]$ (light to dark). Black dashed lines are the empirical function $f(\overline{(\phi_f)}/\phi_{f,0}, \text{Pe}^{-1})$ from $\delta_{est,SL}$ (see (3.4) and (3.5)).</figcaption></figure>

$T = 0.015\pi$ to $T = 0.8\pi$. For the slowest case ($T = 0.8\pi$), $\delta_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_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_D$ decreases as $T$ decreases in very fast loading (figures 9 and 11).

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
<figure><img src="image_13.png" /><figcaption>Figure 11. Time evolution of $\delta_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.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>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.</figcaption></figure>

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
<figure><img src="image_15.png" /><figcaption>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).</figcaption></figure>

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

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. (2023) by calculating the solid velocity as $v_{s,vf} \approx du_{s,vf}/dt$, where the solid displacement $u_{s,vf}$ is given in their equation (3.12) and the subscript (‘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

$$ \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]. \quad (3.7) $$

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

$$ \delta_{est} = C_1 f \left( \frac{\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, \quad (3.8) $$
<figure><img src="image_16.png" /></figure>

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.

where $f((\overline{\langle\phi_f\rangle}/\phi_{f,0}), \text{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:

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

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

(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 $\text{Pe}^{-1}$ and $\alpha$, but the evolution of the porosity and fluid velocity with $A$ and $T$ does not. As a
<table><caption>Table 2. Loading and transport parameters for some examples of biological materials.</caption><thead><tr><th>Tissue</th><th>D<sub>m</sub> [m<sup>2</sup>s<sup>-1</sup>]</th><th>Pe</th><th>T&#775;/&pi;</th><th>Pe<sub>eff</sub> = Pe(A&#775;/T&#775;)</th></tr></thead><tbody><tr><td>Brain ECM (Kedarasetti et al. 2020)</td><td>1.4 &times; 10<sup>-10</sup></td><td>30</td><td>0.003 &ndash; 0.1</td><td>3 &times; 10<sup>1</sup> &ndash; 1 &times; 10<sup>3</sup></td></tr><tr><td>Cartilage (Ferguson et al. 2004; Gardiner et al. 2007; DiDomenico et al. 2017)</td><td>5 &times; 10<sup>-12</sup></td><td>7 &times; 10<sup>3</sup></td><td>0.03 &ndash; 3 (sitting)<br>0.003 &ndash; 0.03 (running)</td><td>10<sup>2</sup> &ndash; 10<sup>3</sup><br>10<sup>3</sup> &ndash; 10<sup>5</sup></td></tr><tr><td>Intervertebral disc (Annulus F.) (Ferguson et al. 2004)</td><td>5 &times; 10<sup>-12</sup></td><td>4 &times; 10<sup>2</sup></td><td>0.3 (wake cycle)<br>6 &times; 10<sup>-5</sup> &ndash; 6 &times; 10<sup>-3</sup> (sitting)<br>6 &times; 10<sup>-6</sup> &ndash; 6 &times; 10<sup>-5</sup> (running)</td><td>4 &times; 10<sup>1</sup> &ndash; 2 &times; 10<sup>3</sup><br>2 &times; 10<sup>3</sup> &ndash; 2 &times; 10<sup>5</sup><br>2 &times; 10<sup>5</sup> &ndash; 2 &times; 10<sup>6</sup></td></tr><tr><td>Cartilage scaffold (bioreactor) (Sengers et al. 2004)</td><td>10<sup>-9</sup></td><td>10<sup>4</sup></td><td>10<sup>-4</sup> &ndash; 10<sup>-1</sup></td><td>4 &times; 10<sup>1</sup> &ndash; 4 &times; 10<sup>4</sup></td></tr></tbody></table>

result, varying Pe⁻¹ and α 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 Pe are reported, for the same biological tissues, in Fiori et al. (2023). In all cases, the values of Pe$_{eff}$ range from moderate to high, indicating that advection and potentially dispersion are unlikely to be negligible in these systems.

## 4. 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 1-D 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 – causing 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); see 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: (i) molecular diffusion only (no loading); (ii) advection only (no diffusion or dispersion); (iii) advection and diffusion only (no dispersion); and (iv) 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 considered 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. We formalise these trends in the expression of the characteristic transport length $\delta_{est}$, which summarises the different regimes according to the transport and loading parameters.

For small deformations, it may be possible to derive an effective equation for net transport in this system. For example, Pool et al. (2016) used a multiple-scale technique to develop an effective transport equation for periodic loading in a semi-infinite domain that is valid for times much larger than the forcing period. The same approach could potentially be adapted to the present scenario, where the small-deformation analytical solution for the flow field is a full Fourier series rather than a single exponentially damped mode because the domain is bounded. This analysis is beyond the scope of the present work, but it would provide valuable qualitative insight into deformation-driven transport.

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 optimisation 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 and, ultimately, on the underlying model for hydrodynamic dispersion. At the least, it is likely that, much like $k$, both $D_m$ and $\alpha$ should be deformation-dependent for moderate to large deformations.

Future work should also include the expansion of this study to 2-D and 3-D media with heterogeneous properties. Already for small deformations, periodic deformations have been shown to cause complex transport and mixing dynamics in heterogeneous media (Trefry et al. 2019; Wu et al. 2020, 2024) and in response to multi-modal forcing (Trefry et al. 2020). Complex media and forcings are particularly relevant to biological systems; examples include heterogeneous structures originating in the multi-scale nature of tissues and/or multi-modal loading originating from the co-existence of different vital cyclic pulsations, such as the cardiac and the respiratory cycles. Finally, future work should also focus on an experimental validation of the results presented here.

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

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

## 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, we use an implicit Runge-Kutta method via MATLAB's built-in solver ODE15s (Shampine & 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(x, t)}{\partial t} + \frac{\partial}{\partial x} [F(\Phi(x, t))] = 0. \quad (A1) $$

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

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

which maps $a \le x \le 1$ to $0 \le \xi \le 1$. Performing a standard change of variables from $(x, t)$ to $(\xi, t)$, partial derivatives transform according to

$$ \frac{\partial}{\partial t} \rightarrow \frac{\partial}{\partial t} + \frac{\partial \xi}{\partial t} \frac{\partial}{\partial \xi} \quad \text{and} \quad \frac{\partial}{\partial x} \rightarrow \frac{\partial \xi}{\partial x} \frac{\partial}{\partial \xi} \quad (A3) $$

and equation (A1) can be rewritten as

$$ \frac{\partial \Phi}{\partial t} - \frac{(1-\xi)}{(1-a)} \frac{\partial \Phi}{\partial \xi} + \left( \frac{1}{1-a} \right) \frac{\partial F(\Phi)}{\partial \xi} = 0. \quad (A4) $$

<figure><img src="image_17.png" /><figcaption>Figure 15. Convergence analysis: r.m.s. 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.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>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 δ over time.</figcaption></figure>

When solving equation (2.25), we then take $\Phi = \phi_f$ and

$$F(\phi_f) = -\tilde{D}_f(\phi_f) \frac{\partial \phi_f}{\partial \tilde{x}}, \qquad (A5)$$

while for equation (2.27) we take $\Phi = \phi_f \tilde{c}$ and

$$F(\phi_f \tilde{c}) = (\phi_f \tilde{c}) \tilde{v}_f - \phi_f \tilde{D} \frac{\partial \tilde{c}}{\partial \tilde{x}}. \qquad (A6)$$

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 (r.m.s.) relative error in $c(t, x = 1)$ for each solution with respect to the solution for $N_x = 1000$. To balance between accuracy and computational time, we choose for all our simulations $N_x = 300$, with a r.m.s. 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 Appendix B), resulting in good agreement between numerical and analytical results.
<figure><img src="image_19.png" /><figcaption>Figure 17. (a) Smallest and (b) largest values of $q_{disp}/q_{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).</figcaption></figure>

## Appendix B. Analytical solution for molecular diffusion and advection

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

$$c(x, t) = c_0 \left\{ 1 - \frac{1}{2} \left[ \operatorname{erfc} \left( \frac{l+x-1}{2\sqrt{\mathrm{Pe}^{-1}t}} \right) + \operatorname{erfc} \left( \frac{l-x+1}{2\sqrt{\mathrm{Pe}^{-1}t}} \right) \right] \right\}. \quad (\mathrm{B}1)$$

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

$$c(x, t) = \frac{1}{2} \left\{ \tanh \left[ s \left( (x-1) \frac{\phi_{f,0}}{\phi_f} + l \right) \right] + 1 \right\}. \quad (\mathrm{B}2)$$

In figure 16, we compare these analytical solutions with 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.

## 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 table 1. We do so by quantifying the strength of dispersion relative to diffusion using the ratio of the dispersive solute flux $q_{disp}$ to the diffusive solute flux $q_{diff}$ which we define as

$$\frac{q_{disp}}{q_{diff}} = \alpha \mathrm{Pe} |v_f - v_s|, \quad (\mathrm{C}1)$$

and which measures the relative importance of these two mechanisms. In figure 17, we compare $q_{disp}/q_{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
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>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$.</figcaption></figure>

domain, while in the latter dispersion is the dominant mechanism. Hence, the specific baseline values chosen for $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.

## Appendix D. Initial porosity and solute strip width

Results of varying $\phi_{f,0}$ and $l$ are shown in figure 18. Decreasing $\phi_{f,0}$ leads to larger oscillations in fluid flux for a given solid velocity (see equation (2.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 fluid flux 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 with the initial variance is lower when there is more solute (see equation 3.3).
REFERENCES

ALBRO, M.B., CHAHINE, N.O., LI, R., YEAGER, K., HUNG, C.T. & ATESHIAN, G.A. 2008 Dynamic loading of deformable porous media can induce active solute transport. *J. Biomech.* **41** (15), 3152–3157.

ALSHAWABKEH, A.N. & RAHBAR, N. 2006 Parametric study of one-dimensional solute transport in deformable porous media. *J. Geotech. Geoenviron. Engng* **132** (8), 1001–1010.

DE ANNA, P., LE BORGNE, T., DENTZ, M., TARTAKOVSKY, A.M., BOLSTER, D. & DAVY, P. 2013 Flow intermittency, dispersion, and correlated continuous time random walks in porous media. *Phys. Rev. Lett.* **110** (18), 184502.

AREGA, F. & HAYTER, E. 2008 Coupled consolidation and contaminant transport model for simulating migration of contaminants through the sediment and a cap. *Appl. Math. Model.* **32** (11), 2413–2428.

ATESHIAN, G.A., ALBRO, M.B., MAAS, S. & WEISS, J.A. 2011 Finite element implementation of mechanochemical phenomena in neutral deformable porous media under finite deformation. *J. Biomech.* *Eng.* **133** (8), 081005.

AUTON, L.C. & MACMIN, C.W. 2018 From arteries to boreholes: transient response of a poroelastic cylinder to fluid injection. *Proc. R. Soc. Lond. A: Math. Phys. Engng Sci.* **474**, 20180284.

BEAR, J. 1972 *Dynamics of Fluids in Porous Media*. American Elsevier Publishing Company.

BERG, M., DAVIT, Y., QUINTARD, M. & LORTHOIS, S. 2020 Modelling solute transport in the brain microcirculation: is it really well mixed inside the blood vessels? *J. Fluid Mech.* **884**, A39.

BOJARSKAITE, L., VALLET, A., BJØRNSTAD, D.M., GULLESTAD BINDER, K.M., CUNEN, CÉLINE, HEUSER, K., KUCHTA, M., MARDAL, K.-A. & ENGER, R. 2023 Sleep cycle-dependent vascular dynamics in male mice and the predicted effects on perivascular cerebrospinal fluid flow and solute transport. *Nat. Commun.* **14** (1), 953.

BONAZZI, A., JHA, B. & DE BARROS, F.P.J. 2021 Transport analysis in deformable porous media through integral transforms. *Intl J. Numer. Anal. Meth. Geomech.* **45** (3), 307–324.

BORJA, R.I. & CHOO, J. 2016 Cam-clay plasticity, part VIII: a constitutive framework for porous materials with evolving internal structure. *Comput. Meth. Appl. Mech. Engng* **309**, 653–679.

BOUCHAUDY, A. & SALMON, J.-B. 2019 Drying-induced stresses before solidification in colloidal dispersions: *in situ* measurements. *Soft Matter* **15**, 2768–2781.

BRENNER, H. & EDWARDS, D. 1993 *Macrotransport Processes, Butterworth-Heinemann series in chemical engineering*. Elsevier.

BRENNER, H. & STEWARTSON, K. 1980 Dispersion resulting from flow through spatially periodic porous media. *Phil. Trans. R. Soc. Lond. A Math Phys. Sci.* **297** (1430), 81–133.

BUIJS, J.O.D., RITMAN, E.L. & DRAGOMIR-DAESCU, D. 2010 Validation of a fluid-structure interaction model of solute transport in pores of cyclically deformed tissue scaffolds. *Tissue Engng C: Methods* **16**, 1145–1156.

CHARLAIX, E., HULIN, J.P. & PLONA, T.J. 1987 Experimental study of tracer dispersion in sintered glass porous materials of variable compaction. *Phys. Fluids* **30** (6), 1690–1698.

CLARIA, J.J., GOLDSZTEIN, G.H. & SANTAMARINA, J.C. 2012 AC diffusion: transport in porous networks subjected to zero-time-average advective flow. *Transp. Porous Med.* **93** (1), 51–61.

CORTEZ, S., COMPLETO, A. & ALVES, J.L. 2016 The influence of mechanical stimulus on nutrient transport and cell growth in engineered cartilage: A finite element approach. In *Proceedings of the VII European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS)* (ed. M. Papadrakakis, V. Papadopoulos, G. Stefanou & V. Plevris), 5–10 June, Crete Island, Greece.

CRANK, J. 1979 *The Mathematics of Diffusion*. Oxford University Press.

DANCKWERTS, P.V. 1952 The definition and measurement of some characteristics of mixtures. *Flow Turbul. Combust.* **3**, 279–296.

DATTA, S.S., CHIANG, H., RAMAKRISHNAN, T.S. & WEITZ, D.A. 2013 Spatial fluctuations of fluid velocities in flow through a three-dimensional porous medium. *Phys. Rev. Lett.* **111** (6), 064501.

DAVIS, K.E. & RUSSEL, W.B. 1989 An asymptotic description of transient settling and ultrafiltration of colloidal dispersions. *Phys. Fluid A: Fluid Dyn.* **1** (1), 82–100.

DAVIT, Y., BYRNE, H., OSBORNE, J., PITT-FRANCIS, J., GAVAGHAN, D. & QUINTARD, M. 2013 Hydrodynamic dispersion within porous biofilms. *Phys. Rev. E* **87** (1), 012718.

DELGADO, J. 2007 *Chem. Engng Res. Design* **85**, 1245.

DENTZ, M., DE BARROS, F.P., LE BORGNE, T. & LESTER, D.R. 2018 Evolution of solute blobs in heterogeneous porous media. *J. Fluid Mech.* **853**, 621–646.

DENTZ, M., HIDALGO, J.J. & LESTER, D.R. 2023 Mixing in porous media: concepts and approaches across scales. *Transp. Porous Med.* **146** (1-2), 5–53.

DENTZ, M., LE BORGNE, T., ENGLERT, A. & BIJELJIC, B. 2011 Mixing, spreading and reaction in heterogeneous media: a brief review. *J. Contam. Hydrol.* **120-121**, 1–17.
DIDOMENICO, C.D., GOODEARL, A., YARILINA, A., SUN, V., MITRA, S., STERMAN, A.S. & BONASSAR, L.J. 2017 The effect of antibody size and mechanical loading on solute diffusion through the articular surface of cartilage. *J. Biomech. Engng* **139** (9), 091005.

DIDOMENICO, C.D., LINTZ, M. & BONASSAR, L. 2018 Molecular transport in articular cartilage – what have we learned from the past 50 years? *Nat. Rev. Rheumatol.* **14** (7), 393–403.

DIDOMENICO, C.D., WANG, Z.X. & BONASSAR, L.J. 2017 Cyclic mechanical loading enhances transport of antibodies into articular cartilage. *J. Biomech. Engng* **139** (1), 1–7.

EHLERS, W., KARAJAN, N. & MARKERT, B. 2009 An extended biphasic model for charged hydrated tissues with application to the intervertebral disc. *Biomech. Model. Mechanobiol.* **8** (3), 233–251.

FAN, L., PEI, S., LUCAS LU, X. & WANG, L. 2016 A multiscale 3D finite element analysis of fluid/solute transport in mechanically loaded bone. *Bone Res.* **4** (1), 16032.

FEDERICO, S. & HERZOG, W. 2008 On the permeability of fibre-reinforced porous materials. *Intl J. Solids Struct.* **45** (7–8), 2160–2172.

FERGUSON, S.J., ITO, K. & PYRAK-NOLTE, L.J. 2004 Fluid flow and convective transport of solutes within the intervertebral disc. *J. Biomech.* **37** (2), 213–221.

FIORI, M., PRAMANIK, S. & MACMINN, C.W. 2023 Flow and deformation due to periodic loading in a soft porous material. *J. Fluid Mech.* **974**, A2.

FIORI, M., PRAMANIK, S. & MACMINN, C.W. 2025 Solute transport due to periodic loading in a soft porous material: example code. Zenodo, http://dx.doi.org/10.5281/zenodo.14853162.

FOX, P.J. 2007a Coupled large strain consolidation and solute transport. I: model development. *J. Geotech. Geoenviron. Engng* **133** (1), 3–2007.

FOX, P.J. 2007b Coupled large strain consolidation and solute transport. II: model verification and simulation results. *J. Geotech. Geoenviron. Engng* **133** (1), 16–2007.

FRANCESCHINI, G., BIGONI, D., REGITNIG, P. & HOLZAPFEL, G. 2006 Brain tissue deforms similarly to filled elastomers and follows consolidation theory. *J. Mech. Phys. Solids* **54** (12), 2592–2620.

GARDINER, B., SMITH, D., PIVONKA, P., GRODZINSKY, A., FRANK, E. & ZHANG, L. 2007 Solute transport in cartilage undergoing cyclic deformation. *Comput. Method. Biomed. Engng* **10** (4), 265–278.

GELHAR, L. 1993 *Stochastic Subsurface Hydrology*. Prentice Hall.

GELHAR, L.W. & AXNESS, C.L. 1983 Three-dimensional stochastic analysis of macrodispersion in aquifers. *Water Resour. Res.* **19** (1), 161–180.

GELHAR, L.W., WELTY, C. & REHFEEDT, K.R. 1992 A critical review of data on field-scale dispersion in aquifers. *Water Resour. Res.* **28** (7), 1955–1974.

GHANBARIAN, B., HUNT, A.G., EWING, R.P. & SAHIMI, M. 2013 Tortuosity in porous media: a critical review. *Soil Sci. Soc. Am. J.* **77** (5), 1461–1477.

GOIRAND, F., BORGNE, T.LE & LORTHOIS, S. 2021 Network-driven anomalous transport is a fundamental component of brain microvascular dysfunction. *Nat. Commun.* **12** (7295).

GOLDSZTEIN, G.H. & SANTAMARINA, J.C. 2004 Solute transport during cyclic flow in saturated porous media. *Appl. Phys. Lett.* **85** (12), 2432–2434.

GOUZE, P., PUYGUIRAUD, A., PORCHER, T. & DENTZ, M. 2021 Modeling longitudinal dispersion in variable porosity porous media: control of velocity distribution and microstructures. *Front. Water* **3**, 766338.

HENCKY, H. 1931 The law of elasticity for isotropic and quasi-isotropic substances by finite deformations. *J. Rheol.* **2** (2), 169–176.

HEWITT, D.R., PATERSON, D.T., BALMFORTH, N.J. & MARTINEZ, D.M. 2016 Dewatering of fibre suspensions by pressure filtration. *Phys. Fluids* **28** (6), 063304.

HRABE, J., HRABĚTOVÁ, S. & SEGETH, K. 2004 A model of effective diffusion and tortuosity in the extracellular space of the brain. *Biophys. J.* **87** (3), 1606–1617.

JHA, B., CUETO-FELGUEROSO, L. & JUANES, R. 2011 Quantifying mixing in viscously unstable porous media flows. *Phys. Rev. E* **84** (6), 066312.

KEDARASETTI, R.T., DREW, P.J. & COSTANZO, F. 2020 Arterial vasodilation drives convective fluid flow in the brain: a poroelastic model. *Fluids Barriers CNS* **19** (1), 34.

KELLEY, D.H. & THOMAS, J.H. 2023 Cerebrospinal fluid flow. *Annu. Rev. Fluid Mech.* **55** (1), 237–264.

KREE, M. & VILLERMAUX, E. 2017 Scalar mixtures in porous media. *Phys. Rev. Fluids* **2** (10), 104502.

KUMAR, P., DEY, B. & RAJA SEKHAR, G.P. 2018 Nutrient transport through deformable cylindrical scaffold inside a bioreactor: an application to tissue engineering. *Intl J. Engng Sci.* **127**, 201–216.

LAI, W.M., HOU, J.S. & MOW, V.C. 1991 A triphasic theory for the swelling and deformation behaviors of articular cartilage. *J. Biomech. Engng* **113** (3), 245–258.

LEE, M.K., RICH, M.H., BAEK, K., LEE, J. & KONG, H. 2015 Bioinspired tuning of hydrogel permeability-rigidity dependency for 3D cell culture. *Sci. Rep. UK* **5**, 1–7.
LELE, S.K. 1992 Compact finite difference schemes with spectral-like resolution. *J. Comput. Phys.* **103** (1), 16–42.

LESTER, D.R., DENTZ, M. & LE BORGNE, T. 2016a Chaotic mixing in three-dimensional porous media. *J. Fluid Mech.* **803**, 144–174.

LESTER, D.R., METCALFE, G. & TREFRY, M.G. 2013 Is chaotic advection inherent to porous media flow? *Phys. Rev. Lett.* **111** (17), 1–5.

LESTER, D.R., TREFRY, M.G. & METCALFE, G. 2016b Chaotic advection at the pore scale: mechanisms, upscaling and implications for macroscopic transport. *Adv. Water Resour.* **97**, 175–2016.

LEVICK, J. 1987 Flow through interstitium and other fibrous matrices. *Q. J. Expl Physiol. Translation and Integration* **72** (4), 409–437.

LEWIS, T.W., PIVONKA, P. & SMITH, D.W. 2009 Theoretical investigation of the effects of consolidation on contaminant transport through clay barriers. *Intl J. Numer. Anal. Meth. Geomech.* **33** (1), 95–116.

LI, C., BORJA, R.I. & REGUEIRO, R.A. 2004 Dynamics of porous media at finite strain. *Comput. Meth. Appl. Mech. Engng* **193** (36–38), 3837–3870.

LIANG, Y., WEN, B., HESSE, M.A. & DICARLO, D. 2018 Effect of dispersion on solutal convection in porous media. *Geophys. Res. Lett.* **45** (18), 9690–9698.

LIU, Y., GONG, W., XIAO, H. & WANG, M. 2024 Non-monotonic effect of compaction on longitudinal dispersion coefficient of porous media. *J. Fluid Mech.* **988**, R2.

MACMINN, C.W., DUFRESNE, E.R. & WETTLAUFER, J.S. 2016 Large deformations of a soft porous material. *Phys. Rev. Appl.* **5** (4), 044020.

MARBACH, S. & ALIM, K. 2019 Active control of dispersion within a channel with flow and pulsating walls. *Phys. Rev. Fluids* **4** (11), 114202.

MAROUDAS, A. 1970 Distribution and diffusion of solutes in articular cartilage. *Biophys. J.* **10** (5), 365–379.

MAROUDAS, A. 1975 Biophysical chemistry of cartilaginous tissues with special reference to solute and fluid transport. *Biorheology* **12** (3–4), 233–248.

MAUCK, R.L., HUNG, C.T. & ATESHIAN, G.A. 2003 Modeling of neutral solute transport in a dynamically loaded porous permeable gel: implications for articular cartilage biosynthesis and tissue engineering. *J. Biomech. Engng* **125** (5), 602–614.

MESALLATI, T., BUCKLEY, C.T., NAGEL, T. & KELLY, D.J. 2013 Scaffold architecture determines chondrocyte response to externally applied dynamic compression. *Biomech. Model. Mechanbiol.* **12** (5), 889–899.

MISHRA, M., MARTIN, M. & DE WIT, A. 2008 Differences in miscible viscous fingering of finite width slices with positive or negative log-mobility ratio. *Phys. Rev. E* **78** (6), 066306.

MOEENDARBARY, E., VALON, L., FRITZSCHE, M., HARRIS, A.R., MOULDING, D.A., THRASHER, A.J., STRIDE, E., MAHADEVAN, L. & CHARRAS, G.T. 2013 The cytoplasm of living cells behaves as a poroelastic material. *Nat. Mater.* **12** (3), 253–261.

MOW, V.C., HOLMES, M.H. & LAI, W.M. 1984 Fluid transport and mechanical properties of articular cartilage: a review. *J. Biomech.* **17** (5), 377–394.

MOW, V.C., KUEI, S.C., LAI, W.M. & ARMSTRONG, C.G. 1980 Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiments. *J. Biomech. Engng* **102** (1), 73–84.

ÖSTERGREN, K.C.E. & TRÄGÅRDH, C. 2000 Characterization of hydrodynamic dispersion in a chromatographic column under compression. *Chem. Engng J.* **79** (2), 103–111.

OSWALD, S. & KINZELBACH, W. 2004 Three-dimensional physical benchmark experiments to test variable-density flow models. *J. Hydrol* **290** (1-2), 22–42.

PEPPIN, S., ELLIOTT, J. & WORSTER, M.G. 2006 Solidification of colloidal suspensions. *J. Fluid Mech.* **554**, 147.

PETERS, G.P. & SMITH, D.W. 2002 Solute transport through a deforming porous medium. *Intl J. Numer. Anal. Meth. Geomech.* **26** (7), 683–717.

POOL, M. & DENTZ, M. 2018 Effects of heterogeneity, connectivity, and density variations on mixing and chemical reactions under temporally fluctuating flow conditions and the formation of reaction patterns. *Water Resour. Res.* **54** (1), 186–204.

POOL, M., DENTZ, M. & POST, V.E.A. 2016 Transient forcing effects on mixing of two fluids for a stable stratification. *Water Resour. Res.* **52** (9), 7178–7197.

PREZIOSI, L., JOSEPH, D. & BEAVERS, G. 1996 Infiltration of initially dry, deformable porous media. *Intl J. Multiphase Flow* **22** (6), 1205–1222.

PU, H., FOX, P.J. & SHACKELFORD, C.D. 2018 *J. Geotech. Geoenviron. Engng* **144**, 1.

SAFFMAN, P.G. 1959 A theory of dispersion in a porous medium. *J. Fluid Mech.* **6** (03), 321.

SCHEIDEGGER, A.E. 1961 General theory of dispersion in porous media. *J. Geophys. Res.* **66** (10), 3273–3278.
SCHMIDT, H., SHIRAZI-ADL, A., GALBUSERA, F. & WILKE, H.-J. 2010 Response analysis of the lumbar spine during regular daily activities: a finite element analysis. *J. Biomech.* **43** (10), 1849–1856.

SENGERS, B.G., OOMENS, C.W.J. & BAAIJENS, F.P.T. 2004 An integrated finite-element approach to mechanics, transport and biosynthesis in tissue engineering. *J. Biomech. Engng* **126** (1), 82–91.

SHAMPINE, L.F. & REICHELT, M.W. 1997 The MATLAB ODE suite. *SIAM J. Sci. Comput.* **18** (1), 1–22.

SHARP, K.M., CARARE, R.O. & BRYN, A. M, 2019 Dispersion in porous media in oscillatory flow between flat plates: applications to intrathecal, periarterial and paraarterial solute transport in the central nervous system. *Fluids and Barriers of the CNS* **16** (13).

SMITH, D.W. 2000 One-dimensional contaminant transport through a deforming porous medium: theory and a solution for a quasi-steady-state problem. *Intl J. Numer. Anal. Meth Geomech.* **24** (8), 693–722.

SOUZY, M., LHUISSIER, H., MÉHEUST, Y., LE BORGNE, T. & METZGER, B. 2020 Velocity distributions, dispersion and stretching in three-dimensional porous media. *J. Fluid Mech.* **891**, A16.

SOUZY, M., LHUISSIER, H., VILLERMAUX, E. & METZGER, B. 2017 Stretching and mixing in sheared particulate suspensions. *J. Fluid Mech.* **812**, 611–635.

SOUZY, M., PHAM, P. & METZGER, B. 2016 Taylor’s experiment in a periodically sheared particulate suspension. *Phys. Rev. Fluids* **1** (4).

STYLE, R.W. & PEPPIN, S.S.L. 2011 *Proc. R. Soc. Lond. A: Math. Phys. Engng Sci.* **467**, 174.

TAN, C.T. & HOMSY, G.M. 1988 Simulation of nonlinear viscous fingering in miscible displacement. *Phys. Fluids* **31** (6), 1330–1338.

TARTAKOVSKY, D.M. & DENTZ, M. 2019 Diffusion in porous media: phenomena and mechanisms. *Transp. Porous Med.* **130** (1), 105–127.

TAYLOR, G.I. 1953 Dispersion of soluble matter in solvent flowing slowly through a tube. *Proc. R. Soc. Lond. A. Math Phys. Sci.* **219** (1137), 186–203.

TRAN, M. & JHA, B. 2020 Coupling between transport and geomechanics affects spreading and mixing during viscous fingering in deformable aquifers. *Adv. Water Resour.* **136**, 103485.

TREFRY, M.G., LESTER, D.R., METCALFE, G. & WU, J. 2019 Temporal fluctuations and poroelasticity can generate chaotic advection in natural groundwater systems. *Water Resour. Res.* **55** (4), 3347–3374.

TREFRY, M.G., LESTER, D.R., METCALFE, G. & WU, J. 2020 Lagrangian complexity persists with multimodal flow forcing in compressible porous systems. *Transp. Porous Med.* **135** (3), 555–586.

TROYETSKY, D.E., TITHOF, J., THOMAS, J.H. & KELLEY, D.H. 2021 Dispersion as a waste-clearance mechanism in flow through penetrating perivascular spaces in the brain. *Sci. Rep. UK* **11** (4595).

URCIUOLO, F., IMPARATO, G. & NETTI, P.A. 2008 Effect of dynamic loading on solute transport in soft gels implication for drug delivery. *AIChE J.* **54** (3), 824–834.

VAUGHAN, B.L., GALIE, P.A., STEGEMANN, J.P. & GROTBERG, J.B. 2013 A poroelastic model describing nutrient transport and cell stresses within a cyclically strained collagen hydrogel. *Biophys. J.* **105** (9), 2188–2198.

VUONG, A.T., YOSHIHARA, L. & WALL, W.A. 2015 A general approach for modeling interacting flow through porous media under finite deformations. *Comput. Meth. Appl. Mech. Engng* **283**, 1240–1259.

WHITAKER, S. 1967 Diffusion and dispersion in porous media. *AIChE J.* **13** (3), 420–427.

WHITAKER, S. 1998 *The Method of Volume Averaging, Theory and Applications of Transport in Porous Media*. Springer.

WITT, F., DUDA, G.N., BERGMANN, C. & PETERSEN, A. 2014 *Tissue Engng A* **20**, 486.

WORSTER, M.G., PEPPIN, S. & WETTLAUFER, J. 2021 Colloidal mushy layers. *J. Fluid Mech.* **914**, A28.

WU, J., LESTER, D., TREFRY, M.G. & METCALFE, G. 2024 Lagrangian coherent structures control solute dispersion in heterogeneous poroelastic media. *Phys. Rev. Fluids* **9** (4), 044501.

WU, J., LESTER, D.R., TREFRY, M.G. & METCALFE, G. 2020 When do complex transport dynamics arise in natural groundwater systems? *Water Resour. Res.* **56** (2), 1–17.

XIE, H., YAN, H., FENG, S., WANG, Q. & CHEN, P. 2016 An analytical model for contaminant transport in landfill composite liners considering coupled effect of consolidation, diffusion, and degradation. *Environ. Sci. Pollut. Res.* **23** (19), 19362–19375.

YAO, H. & GU, W.Y. 2007 Convection and diffusion in charged hydrated soft tissues: a mixture theory approach. *Biomech. Model. Mechanbiol.* **6** (1–2), 63–72.

ZHANG, H.J., JENG, D.S., BARRY, D.A., SEYMOUR, B.R. & LI, L. 2013 Solute transport in nearly saturated porous media under landfill clay liners: a finite deformation approach. *J. Hydrol.* **479**, 189–199.

ZHANG, H.J., JENG, D.-S., SEYMOUR, B.R., BARRY, D.A. & LI, L. 2012 Solute transport in partially-saturated deformable porous media: application to a landfill clay liner. *Adv. Water Resour.* **40**, 1–10.

ZHANG, L. 2011 Solute transport in cyclic deformed heterogeneous articular cartilage. *Intl J. Appl. Mech.* **3** (3), 507–524.
ZHANG, L. & SZERI, A.Z. 2005 Transport of neutral solute in articular cartilage: effects of loading and particle size. *Proc. R. Soc. Lond. A* **461** (2059), 2021–2042.

ZHANG, L. & SZERI, A.Z. 2008 Transport of neutral solute in articular cartilage: effect of microstructure anisotropy. *J. Biomech.* **41** (2), 430–437.

ZIMMERMAN, R.A. & TARTAKOVSKY, D.M. 2020 Solute dispersion in bifurcating networks. *J. Fluid Mech.* **901**, A24.