# On the fractional transversely isotropic functionally graded nature of soft
biological tissues

G. Sachin<EMAIL_ADDRESS>S. Natarajan<EMAIL_ADDRESS>O. Barrera<EMAIL_ADDRESS>+44 7585041147. Department of
Mechanical Engineering, Indian Institute of Technology Madras, Chennai,
600036, India.  School of Engineering, Computing and Mathematics, Oxford
Brookes University, Headington, Oxford OX3 0BP, United Kingdom.  Department of
Engineering Science, University of Oxford, Parks Road, OX1 3PJ, Oxford, United
Kingdom.

###### Abstract

This paper focuses on the origin of the poroelastic anisotropic behaviour of
the meniscal tissue and its spatially varying properties. We present confined
compression creep test results on samples extracted from three parts of the
tissue (Central body, Anterior horn and Posterior horn) in three orientations
(Circumferential, Radial and Vertical). We show that a poroelastic model in
which the fluid flow evolution is ruled by non-integer order operators
(fractional Darcy’s law) provides accurate agreement with the experimental
creep data. The model is validated against two additional sets of experimental
data: stress relaxation and fluid loss during the consolidation process
measured as weight reduction. Results show that the meniscus can be considered
as a transversely isotropic poroelastic material. This behaviour is due to the
fluid flow rate being about three times higher in the circumferential
direction than in the radial and vertical directions in the body region of the
meniscus. In the anterior horn, the elastic properties are transversely
isotropic, with the aggregate modulus higher in the radial direction than in
the circumferential and vertical directions. The 3D fractional poroelastic
model is implemented in finite element software and quantities such as flux of
interstitial fluid during the consolidation process, a non-trivial
experimental measure, are determined.

###### keywords:

Meniscus , Fractional poroelasticity , Confined compression tests , Anomalous
diffusion , Finite element simulation

##  1 Introduction

Biological soft tissues, such as the meniscal tissue (Figure 1), exhibit a
hierarchical porous solid matrix and an interstitial fluid flowing into the
pores [1, 2, 3, 4, 5]. The overall mechanical behaviour depends not only on
the solid matrix deformation but also on the movement of the fluid in and out
of the pores during the deformation. Furthermore, the internal architecture of
the structure constituting the solid matrix and collagen determines the
anisotropic behaviour of both the elastic and the transport properties [6].
The poroelastic theory is employed to model the response of biphasic tissue
and better understand the anisotropic nature of the meniscus on the
functioning of the knee.

(a)

(b)

(c)

(d)

Figure 1: (a) Lateral and Medial menisci in the knee joint, (b) Regions and
the directions in the meniscus, (c) 3D porous microstructure of meniscus
cylindrical sample (d)Micro CT scan and Multi-photon microscopy of meniscus
tissue.

The menisci situated between the femoral and tibial cartilage of the knee
joint (Figure 1(a)) have three main regions, viz., anterior horn, central body
and posterior horn (Figure 1(b)). They perform several functions like load
bearing, lubrication, energy absorption and stability [7, 8]. In light of our
recent experimental, theoretical and numerical studies [1, 2, 9, 10, 11], we
understand that each region has different functions. In [2, 1], the unique
structure of the body region that fulfils load-bearing and energy absorption
capabilities is discussed. This is achieved through a sandwich-like structure
with two thin, stiff outside layers which take the load-bearing role and a
thick, softer internal layer which acts as an energy absorption element, an
effective natural damper [2, 11]. The architectural arrangement of the solid
collagen-based matrix appears to be different and gradually changing in the
outside and the internal layers. The outside layers show a dense distribution
of collagen fibres with random orientation and low mean diameter. The elastic
modulus is one order of magnitude higher than the internal layer. Whereas the
internal layer of the meniscus relies on a hierarchical anisotropic network of
collagen channels mainly aligned along the circumferential direction in which
fluid flows during the physiological deformation process (Figure 1). In [5],
we have shown that the permeability changes from the anterior/posterior horns
to the body region, with the body region being more permeable. Here, we show
that the permeability tensor in this region is also transversely isotropic and
that the permeability in the circumferential direction (the preferential
direction of the collagen channel) is higher than the other two directions
(Figure 1). Fluid flowing inside these channels determines the time-dependent
behaviour of the tissue. As the tissue deforms under the action of loads,
these channels change morphology, which results in change in permeability as a
function of channels dimension, porosity, and tortuosity. Thus, making the
permeability tensor not only depending on the pressure but also varying
spatially and temporally because of the local variation of applied pressure
gradients. Recent finite element studies of the human knee have shown that
anterior and posterior horns might have a stability function [12, 10] as they
deform to allow the meniscus to be extruded so that the body region has a
higher contact area with the cartilage and can then carry most of the load.

Appropriate modelling of the fluid flow behaviour in the different portions of
the meniscal tissue when subject to a range of loading conditions is essential
to gain insight into the biomechanical function of the tissue and design
appropriate artificial tissue substitutes. Current barriers to the clinical
and functional outcomes of the to-date devices [13] are that the transport and
structure-properties relationship of the native tissues is still not well
understood. In particular, the relationship between the “micro/nano”-scopic
composition and structure of the meniscus and its coarse-scale anisotropic
behaviour structure is still an open question [14].

Experimental tests, such as confined compression tests (relaxation and/or
creep) are currently used to characterize the material parameters such as
elastic modulus at equilibrium and permeability. The main models used to
identify these two parameters are based on biphasic theory. In [15], the
authors examined, through a confined compression test, the effect of fluid in
cartilage is bearing up to 90% of the applied load. They used a linear
biphasic model to find out the material properties in terms of the aggregate
modulus and permeability. In [16], authors performed a creep indentation test
on different regions of the meniscus of different animal species. The
different material properties such as the aggregate modulus, Poisson’s ratio,
shear modulus and permeability were estimated by using the biphasic model to
fit the experimental data.

Our recent experimental findings on human meniscal samples show that the
biphasic model was not sufficient to reproduce the observed experimental
behaviour [5]. To address this, [17] proposed a poroelastic model wherein the
pore pressure diffusion equation is derived by adopting a modified version of
Darcy’s law involving fractional derivatives. It follows that the hydraulic
permeability becomes “anomalous” and the rate of fluid flow is governed by the
order of the derivative. It is thus clear that the anomalous hydraulic
permeability, the order of derivative tensors and their anisotropic nature,
which govern the transport of fluid within the complex porous structure of the
meniscus, are key properties of the meniscal tissue. Though a finite element
implementation of the fractional poroelastic model is reported in [17],
however, the model is restricted to estimating only the pore pressure field.

The purpose of this paper is to present for the first time creep data of
confined compression tests performed in different regions of the meniscal
tissue (anterior-horn, central-body and posterior-horn) and in different
directions (radial, circumferential and vertical). We adopt the fractional
poroelastic model presented in [17, 5] to fit the experimental data to recover
the aggregate modulus, anomalous permeability and order of derivative and
their dependence on the directions (Aggregate modulus, anomalous permeability
and order of derivative tensors). For the first time, we validate the model in
terms of Stress, weight loss and displacement fields by using the parameters
recovered from fitting the creep curve to match two additional sets of data
for the same samples, i.e., relaxation and weight loss data. The fractional
poroelastic model is implemented in commercial software Abaqus using UMATHT
and UMAT. The numerical simulation results show a good agreement with
analytical solutions regarding both displacements and pore pressure fields
which were lacking in the literature. This enabled us to run consolidation
tests numerically and calculate fields such as the flux of interstitial fluid
during the consolidation process, which is a difficult experimental measure to
obtain. Additionally, we perform statistical analysis of the parameters using
ANOVA method to compare the variation of the material parameters in different
directions.

##  2 Material and Methods

###  2.1 Confined consolidation creep tests

Menisci were harvested from patients undergoing total knee arthroplasty under
ethically approved protocol as reported in [5]. Samples labelled as “degraded”
by the gross investigation of the surgeon were discarded. Menisci were thawed
in Phosphate-buffered solution (PBS) for thirty minutes to recover their
hydration. A total of 29 meniscal test samples were harvested. The samples
were cylindrical in shape with a 3mm diameter and 3-4 mm in height. Samples
were extracted from the central body, anterior horn and posterior horn regions
along each principal direction, i.e. circumferential, radial and vertical
(Figure 1).

(a)

(b)

(c)

Figure 2: (a) Schematic representation of confined compression test, (b)
applied load for confined compression creep test and (c) measured
displacement.

The confined compression test set-up consisted of a chamber with transparent
and insulated walls and a porous plate at the bottom(Figure 2(a)). The sample
to be tested is confined inside the chamber. A load of 0.5 N (corresponding to
a stress of about 0.07 MPa which is a physiological value for menisci [18]) at
the top of the chamber is applied. It is applied through a piston with a stage
velocity of 0.3 % height of the sample till the load 0f 0.5 N is reached and
then the load is kept constant for 400s as shown in Figure 2. The applied
force and the obtained displacement that is measured at the top surface are
shown in Figures 2(b)-2(c), respectively.

###  2.2 Fractional consolidation - governing equations & finite element
implementation

#### Governing equations

Consider a poroelastic material occupying $\Omega\subset\mathbb{R}^{d}$,
bounded by $\partial\Omega\subset\mathbb{R}^{d-1}$, where $d=1,2,3$. The
boundary accommodates the following decomposition:
$\partial\Omega=\partial\Omega_{N}\cup\partial\Omega_{D}$ and
$\partial\Omega_{N}\cap\partial\Omega_{D}=\emptyset$, where Neumann and
Dirichlet boundary conditions are applied on $\partial\Omega_{N}$ and
$\partial\Omega_{D}$, respectively. In the absence of body force and inertia,
the governing differential for a poroelastic material is given by:

$\displaystyle\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}$
$\displaystyle=\mathbf{0}$ (1a) $\displaystyle\dfrac{\partial\zeta}{\partial
t}+\boldsymbol{\nabla}\cdot\mathbf{J}_{p}$ $\displaystyle=0$ (1b)

supplemented with the following boundary conditions:
$\mathbf{u}=\hat{\mathbf{u}},~{}{\rm and}~{}p=\hat{p}$ on $\partial\Omega_{D}$
and $\boldsymbol{\sigma}\cdot\mathbf{n}=F,~{}{\rm
and}~{}\mathbf{J}_{p}\cdot\mathbf{n}=0$ on $\partial\Omega_{N}$. In Equation
(1), $\zeta$ is the variation of the fluid content, $p$ is the pore pressure,
$\mathbf{J}_{p}$ is the fluid flux and $\boldsymbol{\sigma}$ is the stress.
The relation between the stresses, displacements and pressure; and between the
flux and the pressure is given by:

$\displaystyle\boldsymbol{\sigma}$
$\displaystyle=2G\boldsymbol{\epsilon}+\lambda\mathrm{Tr}(\boldsymbol{\epsilon})\textbf{I}-\alpha
p\textbf{I}$ (2a) $\displaystyle\mathbf{J}_{p}$
$\displaystyle=-\lambda_{\beta}\prescript{}{0}{D}_{t}^{\beta}\boldsymbol{\nabla}p$
(2b)

where $\lambda=\dfrac{(3K-2G)}{3}$ is Lamé’s constant, $G,K$ are the shear and
the bulk modulus, respectively, $\alpha$ is the Biot coefficient and
$\prescript{}{0}{D}_{t}^{\beta}$ is the fractional derivative with respect to
time of order $\beta$. It is important to note that in case of $\beta=$ 0,
$\lambda_{\beta}$ is with dimension of $[L]^{4}[F]^{-1}[T]^{-1}$. In the case
of $\beta\neq 0$, $\lambda_{\beta}$ has dimension of
$[L]^{4}[F]^{-1}[T]^{\beta-1}$. The pore pressure diffusion equation (c.f.
Equation (2b)), written in terms of dilatational strain, $\epsilon_{d}$, is
given by:

$\dfrac{\partial p}{\partial
t}=\dfrac{K_{u}-K}{\alpha^{2}}\lambda_{\beta}\prescript{}{0}{D}_{t}^{\beta}\boldsymbol{\nabla}^{2}p-\dfrac{K_{u}-K}{\alpha}\dfrac{\partial\epsilon_{d}}{\partial
t}$ (3)

where $K_{u}$ is the undrained bulk modulus.

During confined compression tests, a cylinder of meniscal tissue is inserted
inside a cylindrical chamber with a porous base, see Figure 2(a) for geometry
and boundary conditions. Upon applying pressure at the top of the sample, the
fluid contained inside the sample will flow out of the base. This can be
modelled as a one-dimensional version of Equation (2) as the only non-zero
components of strain will be $\epsilon_{zz}$. Equation (2) will be modified to
[17]:

$\displaystyle\left(K+\frac{4G}{3}\right)\dfrac{\partial\epsilon_{zz}}{\partial
z}-\alpha\dfrac{\partial p}{\partial z}=0$ (4a) $\displaystyle\dfrac{\partial
p}{\partial
t}=\dfrac{K_{u}-K}{\alpha^{2}}\lambda_{\beta}\prescript{}{0}{D}_{t}^{\beta}\dfrac{\partial^{2}p}{\partial
z^{2}}-\dfrac{K_{u}-K}{\alpha}\dfrac{\partial\epsilon}{\partial t}$ (4b)

and the corresponding boundary conditions are, $\forall~{}0\leq z\leq h$,
$p(z,t=0)=\gamma P_{A};~{}P(z=0,t)=P_{A},~{}p(z=h,t)=0,~{}\dfrac{\partial
p}{\partial z}\Big{|}_{z=0}=0$ and $u(z=h,t)=0.$ where $h$ is the initial
height of the sample. A constant compressive stress,
$\sigma_{zz}(z=0,t)=-P_{A}$ in the $z-$ direction is applied to the cylinder
at $z=0$. The initial pore pressure is derived assuming that the material is
under undrained conditions at the initial time of loading, given by:

$p(z,t=0)=P_{A}\frac{3(K_{u}-K)}{\alpha\left(4G+3K_{u}\right)}$ (5)

It is shown that the analytical solution of the pore pressure is given by
[17]:

$p(z,t)=P_{A}\gamma\sum\limits_{n=1,3}^{\infty}E_{1-\beta,1}\left(-{\frac{n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}}\right)c_{n}\cos\frac{n\pi
z}{2h}$ (6)

where $E_{1-\beta,1}$ is the Mittag-Leffler function and

$\displaystyle\gamma$
$\displaystyle=\dfrac{3\left(K_{u}-K\right)}{\alpha\left(4G+3K_{u}\right)};\qquad
c_{n}=\dfrac{4}{n\pi}\big{(}-1\big{)}^{\frac{n-1}{2}}$ (7)
$\displaystyle\bar{\lambda}$
$\displaystyle=\lambda_{\beta}\dfrac{\left(4G+3K\right)\left(K_{u}-K\right)}{\alpha^{2}\left(4G+3K_{u}\right)}$

The analytical solution for the displacement is given by:

$u(z,t)=\frac{3P_{A}}{3K+4G}\Bigg{[}(h-z)+\gamma\alpha\sum\limits_{n=1,3}^{\infty}E_{1-\beta,1}\left(-{\dfrac{n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}}\right)\dfrac{8h}{(n\pi)^{2}}\bigg{\\{}(-1)^{\dfrac{n-1}{2}}\sin\frac{n\pi
z}{2h}-1\bigg{\\}}\Bigg{]}$ (8)

The detailed derivation is given in Appendix A.
The expression for the flux is found by substituting the solution for the pore
pressure (c.f. Equation (6)) in the flux definition (c.f. Equation (2b)) and
solving the fractional derivative as follows:

$J_{p}=\lambda_{\beta}P_{A}\gamma
t^{-\beta}\sum\limits_{n=1,3}^{\infty}E_{1-\beta,1-\beta}\Bigg{(}\dfrac{-n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}\Bigg{)}\dfrac{2}{h}\sin\left(\dfrac{n\pi
z}{2h}\right)$ (9)

The weight of the sample over time $W(t)$ is related to the initial weight of
the sample $W_{0}$ and the fluid flux $J_{p}$ as follows:

$W(t)=W_{0}-\int\limits_{0}^{t}J_{p}\cdot
w_{s}~{}\boldsymbol{A}~{}\mathrm{d}t\\\ $ (10)

where $w_{s}$ is the specific weight of the fluid, and $\boldsymbol{A}$ is the
cross-section over which the fluid flows out. Substituting Equation (9) in
Equation (10) and evaluating the integral analytically, we obtain the
following relation:

$W(t)=W_{0}-P_{A}Aw_{s}\gamma\frac{2}{h}\sum\limits_{n=1,3}^{\infty}\left[\lambda_{\beta}t^{1-\beta}E_{1-\beta,2-\beta}\left(-{\frac{\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}}\right)\right]\\\
$ (11)

For detailed derivation, refer to Appendix B.

#### Finite Element implementation and verification

Biot’s model with fractional Darcy’s law is implemented in Abaqus using UMATHT
and UMAT subroutines using the method proposed by Barrera [17].In [17] only
the pore pressure field is discussed, here for the first time we show that
both the pore pressure and the displacement fields agree with the analytical
solution. Using the similarity between the governing equations of
thermoelasticity and poroelasticity [17], UMATHT of Abaqus, which is generally
meant for thermoelasticity, is used for the case of poroelasticity. In Abaqus,
Temperature can be treated as Pore pressure, and the heat flux can be treated
as heat flux. UMAT is used for transferring the stress data to the UMATHT for
modelling the coupling behaviour. UVARM can be used for transferring stress
values. However, due to the order of implementation in Abaqus, which is
elastic model-UMATHT-UVARM, the stress values from UVARM reach the UMATHT in
the next iteration, which introduces errors in the computation. To avoid this
error, UMAT is used, which allows for sharing stress values in the same
increment. UMAT is written for a linear elastic model. The coefficient
coupling the strain and pore pressure $\alpha_{T}=\alpha/3K$ is given as the
thermal expansion coefficient.

Table 1: Material parameters Bulk modulus, $K$ | 1.67$\times$105Pa
---|---
Shear modulus, $G$ | 7.69$\times$ 104 Pa ($E=$ 0.2$\times$106Pa, $\nu=$ 0.3)
Skempton coefficient, $B$ | 0.88
Biot coefficient, $\alpha$ | 0.65
Permeability, $\lambda_{\beta}$ | 8.33$\times$10-8m4/Ns1-β
Fractional Order, $\beta$ | 0, 0.1, 0.2, 0.4

Abaqus model is validated by comparing with theoretical solutions of a
confined compression test with material properties given in Table 1 [17]. For
this, a cylindrical computational model with height $h=$ 3mm and diameter $d=$
3mm is considered. The domain is discretized with 8-noded trilinear hexahedral
elements with 4 degrees of freedom per node (C3D8T elements). The
displacements in all three directions are restrained on the bottom face and
the displacements on the side faces are restrained laterally. The pore
pressure is zero at the bottom and an instantaneous load, $P_{A}=$ 0.07 MPa is
applied on the top surface. Based on a systematic study, a total of 28,110
8-noded hexahedral elements, were found to be adequate to model the
poroelastic response. The displacements and the pore pressure obtained from
the numerical simulation are compared against the theoretical results. Figures
3(a)-3(b) shows the comparison of the pore pressure and displacement as a
function for depth for different time steps for $\beta=$ 0 between theoretical
solution and numerical computation. The influence of $\beta$ on the pore
pressure and the displacement is shown in Figures 3(c)-3(d). It is seen from
Figure 3 that there is a very good agreement with the numerical result and
theoretical solution and that the fractional order has a strong influence on
the pore pressure and the displacement. Large jumps in the pore pressure and
displacement in the initial time in Figures 3(c), 3(d) are due to the
numerical approximation of step function with ramp in the first time
increment.

(a)

(b)

(c)

(d)

Figure 3: Comparison of numerical results with theoretical results: (a-b)
variation of pore pressure and displacement with depth at different time steps
and $\beta=$ 0 and (c-d) influence of $\beta$ on the pore pressure and
displacement as a function of time for a point on the top surface.

##  3 Results

###  3.1 Material parameters: anomalous permeability, order of derivative and
aggregate modulus

The proposed fractional poroelastic model is used to characterize the response
of the meniscus. To this, we use the expression in Equation (8) for fitting
the experimental data. The meniscus is assumed to be incompressible material,
due to which $B=1,\alpha=1,$ and $K_{u}\rightarrow\infty$ reducing Equation
(8) to:

$u(z,t)=\dfrac{P_{A}}{M}\Bigg{[}(h-z)+\sum\limits_{n=1,3}^{\infty}E_{1-\beta,1}\left(-{\dfrac{n^{2}\pi^{2}\lambda_{\beta}Mt^{1-\beta}}{4h^{2}}}\right)\dfrac{8h}{(n\pi)^{2}}\bigg{(}(-1)^{\dfrac{n-1}{2}}\sin\dfrac{n\pi
z}{2h}-1\bigg{)}\Bigg{]}$ (12)

where, $M=(3K+4G)/3$, is the aggregate modulus. Equation (12) is used to find
the material properties of meniscus tissue by fitting with the creep
displacement vs time data. It is noticed that this model requires only three
material properties, viz., the aggregate modulus ($M$), the fractional order
($\beta$) and the permeability ($\lambda_{\beta}$). On using $\beta=$ 0, this
theory converts into the classical Biot’s theory. The experiments were
displacement controlled, i.e., a constant piston velocity was used until the
load reached the required value. So, ramp loading is neglected for fitting and
approximated as step loading. Displacements measured in the experiments were
fitted with the analytical solution Equation (12) at $z=0$, i.e.,

$u=\frac{P_{A}}{M}\bigg{[}h-\sum_{n=1,3}^{\infty}E_{1-\beta,1}\bigg{(}\frac{-n^{2}\pi^{2}\lambda_{\beta}Mt^{1-\beta}}{4h^{2}}\bigg{)}\frac{8h}{(n\pi)^{2}}\bigg{]}$
(13)

Fitting is done in MATLAB using the inbuilt function ‘fminsearch’ with
material properties as variables and the RMS error between the analytical and
the experimental data as the function to be minimized. Fitting plots for four
samples are shown in Figure 5. Table 2 shows the material properties of the
meniscus obtained from the experimental results for samples in both body and
anterior horn for three different directions, viz., circumferential, vertical
and radial. The last two characters of the sample name indicate the location
and orientation information, for example, BC represents Body Circular. The
other notations employed are A: Anterior horn, P: Posterior horn, R: Radial,
and V: vertical directions.

(a)

(b)

(c)

Figure 4: Comparison between Biot’s theory with classical and fractional
Darcy’s law to fit creep displacement experiment data.

To understand the fractional order’s influence, we compared the classical
Biot’s poroelastic model and the proposed one (c.f. Section 2.2). The material
properties obtained for one sample TK11BC through a fitting for classical
Biot’s model111Refer Appendix C for material properties corresponding to the
classical Biot’s theory for all the samples. are $M=$ 3.56$\times$10-1 MPa &
$\lambda=$ 6.95$\times$10-13 m4/Ns and whilst that with the fractional
framework are: $M=$ 1.27$\times$10-1 MPa, $\beta=$ 0.73 & $\lambda_{\beta}=$
2.95$\times$10-12 m4/Ns1-β. The RMS error using Classical Biot’s theory is
6.53$\times$10-5, and the RMS error while using Biot’s theory with fractional
Darcy’s law is 1.42$\times$10-5. Figure 4 shows a comparison of the
displacement, flux out and pore pressure as a function of time for classical
Biot’s theory and fractional model with experiments. It is opined that the
classical Biot’s model does not fit the experimental results well and
reemphasizes a need for a fractional poroelastic framework.

(a) TK11 Body Circumferential

(b) TK11 Body Radial

(c) TK11 Body Vertical

(d) TK16 Anteriorhorn Radial

Figure 5: Displacement confined creep data fitting with fractional poroelasticity and modelling using FEM, (a) Sample taken from the body region of TK11 sample in the circumferential direction, (b) Sample taken from the body region of TK11 sample in the radial direction, (c) Sample taken from the body region of TK11 sample in the vertical direction, (d) Sample taken from the anterior horn region of TK16 sample in the radial direction.  Table 2: Meniscus material parameters obtained from fitting with creep data S No. | Sample | h |  M$\times$105 | $\boldsymbol{\beta}$ |  $\boldsymbol{\lambda_{\beta}}$$\times$10-12 | RMS error
---|---|---|---|---|---|---
|  | (mm) | (Pa) |  |  (m2/Pa.s${}^{1-\beta})$ |  ($\times$10-5)
1 | TK11BC | 3.7 | 1.27 | 0.73 | 2.95 | 1.42
2 | TK11BR | 3.3 | 0.71 | 0.59 | 1.48 | 1.61
3 | TK11BV | 4.1 | 0.83 | 0.53 | 1.56 | 1.74
4 | TK16BC1 | 3.2 | 0.76 | 0.7 | 2.13 | 1.39
5 | TK16BR2 | 2.9 | 0.82 | 0.59 | 0.89 | 1.30
6 | TK16BV | 3.2 | 0.62 | 0.59 | 0.85 | 1.48
7 | TK16AC | 3.2 | 0.60 | 0.63 | 0.97 | 1.97
8 | TK16AR1 | 3.33 | 0.85 | 0.58 | 1.01 | 2.02
9 | TK16AV | 3.33 | 0.57 | 0.5 | 0.37 | 1.63
10 | TK16BC2 | 3.33 | 0.90 | 0.76 | 4.28 | 1.53
11 | TK16BR1 | 2.9 | 0.55 | 0.64 | 0.91 | 1.83
12 | TK16BV2 | 3.13 | 0.61 | 0.63 | 0.93 | 1.56
13 | TK16PV | 2.67 | 1.41 | 0.75 | 1.99 | 1.40
14 | TK17BC | 3.2 | 0.10 | 0.73 | 1.03 | 2.58
15 | TK17BR | 2.2 | 0.41 | 0.72 | 1.36 | 1.55
16 | TK17BV | 2.3 | 0.86 | 0.54 | 0.90 | 1.31
17 | TK17AC | 2.73 | 1.87 | 0.79 | 4.96 | 1.49
18 | TK17AR | 2.7 | 1.94 | 0.76 | 3.05 | 1.44
19 | TK18BC | 4.6 | 1.34 | 0.74 | 141.03 | 3.81
20 | TK18BR | 2.9 | 0.66 | 0.72 | 1.90 | 1.78
21 | TK18BV | 2.8 | 0.68 | 0.61 | 1.07 | 1.63
22 | TK22BR | 3.17 | 2.54 | 0.49 | 85.2 | 9.90
23 | TK22AC | 3.67 | 0.65 | 0.62 | 1.40 | 2.68
24 | TK22AR | 3.33 | 1.60 | 0.41 | 2.24 | 4.69
25 | TK22AV | 2.93 | 0.76 | 0.51 | 0.35 | 1.42
26 | TK36BC | 3.37 | 0.74 | 0.73 | 4.40 | 2.11
27 | TK37BC | 4.6 | 0.90 | 0.74 | 6.63 | 3.53
28 | TK37BR | 3.5 | 0.67 | 0.59 | 0.96 | 1.93
29 | TK37BV | 3.2 | 0.49 | 0.53 | 0.55 | 1.33

Further, to check for anisotropy in the meniscus, the ANOVA test is performed
on the body region, properties with circular, radial and vertical directions
as categories and properties are compared in different directions using
MATLAB’s default function ‘anova1’. For a few samples, properties were
abnormal compared to the remaining samples of the same category. It could be
due to the sample, experiment or fitting issue. These samples were detected as
outliers by the Matlab ’anova1’ function and were not considered for
comparison. The mean and the standard deviation of the groups of the samples
without considering the outliers are given in Table 3 and the results for the
ANOVA test are shown in Figure 6. In Figure 6(a), the aggregate modulus is
compared. The red line shows the sample’s median, and the black dashed line
limit gives the sample’s range. If the notched blue box of different samples
does not overlap, it is concluded that the true medians are different with 95%
confidence. If the probability against the null hypothesis ($p$-value) is less
than 0.05, it is concluded that the mean of the categories is different with
95% confidence. For aggregate modulus, the obtained $p$-value is 0.329,
indicating that the aggregate modulus is not different in different directions
and could possibly be treated as isotropic. The obtained $p$-value for
fractional order $\beta$ and permeability $\lambda_{\beta}$ are 0.001 and
0.0038, respectively. From Figures 6(b)-6(c), it is inferred that the
fractional order and permeability are relatively higher in circular directions
when compared to the other two directions, viz., vertical and radial. Further,
it can be considered similar in vertical and radial directions.

Table 3: Meniscus parameters variation across sample |  $\boldsymbol{M}\times$105 | $\boldsymbol{\beta}$ |  $\boldsymbol{\lambda_{\beta}}\times$10-12
---|---|---|---
| (Pa) |  |  (m2/Pa.s1-β)
Part | Mean | SD | Mean | SD | Mean | SD
Body Cir | 0.75 | 0.49 | 0.73 | 0.00 | 3.75 | 2.36
Body Rad | 0.64 | 0.14 | 0.64 | 0.06 | 1.25 | 0.40
Body Ver | 0.65 | 0.14 | 0.58 | 0.04 | 0.86 | 0.19
Anthorn Cir | 1.04 | 0.72 | 0.68 | 0.10 | 2.44 | 2.19
Anthorn Rad | 1.47 | 0.56 | 0.58 | 0.17 | 2.10 | 1.03
Anthorn Ver | 0.66 | 0.14 | 0.51 | 0.01 | 0.36 | 0.02

(a)

(b)

(c)

Figure 6: Material properties comparison in different orientations in the body
region(BC - Body circumferential, BR - Body Radial, BV - Body Vertical). (a)
Aggregate Modulus (M), (b) Fractional power $\beta$, (c) Permeability
$\lambda_{\beta}$

To study the anisotropic behaviour of the meniscus, the confined creep test is
performed as described in Section 2.1 with material properties taken from the
average values of the body region in circumferential, radial and vertical
directions. For the numerical study, the height of the sample is taken as 3mm
and pressure is increased to 0.07 MPa from zero and held constant. Figure 7
shows the pore pressure, displacement and the flux out of the sample from the
bottom (computed using Equation (9)) as a function of time for the three
different regions.

(a)

(b)

(c)

Figure 7: Flux, displacement and pore pressure for different directions
(Circumferential, Radial and Vertical)

###  3.2 Numerical modelling of fractional consolidation

Numerical modelling of confined creep experiments of the meniscus is done in
Abaqus using the parameters obtained from fitting the data given in Table 2.
Poisson’s ratio and Young’s modulus are the elastic parameters required for
Abaqus. As the confined compression tests depend only on the aggregate
modulus, Poisson’s ratio is arbitrarily assumed to be 0.3, and Young’s modulus
is calculated using the aggregate modulus and Poisson’s ratio. Owing to
symmetry, an axisymmetric model is considered for the numerical simulation.
Figure 8 shows an axisymmetric model employed for this study. The domain is
discretized with 4-noded bilinear quadrilateral elements (CAX4T) and a
structured mesh is used. Based on a systematic mesh convergence study, a mesh
size of 0.05mm was found to be adequate to model the behaviour. Axis-symmetric
boundary conditions are enforced on the left of the computational domain. The
displacements are restrained on the right and bottom faces. Zero pore pressure
condition is applied at the bottom for free fluid flow and pressure
$P_{A}=0.07MPa$ is applied at the top surface as a step load. For the
numerical simulation, a fixed time increment of $\Delta t=$ 0.1 s is used, and
the simulation is carried over for 400 s.

Figure 8: Axisymmetric model employed for the confined creep numerical
simulation

A comparison between the numerical, theoretical and experimental results is
shown in Figure 5 and it is inferred that a good agreement is seen between the
numerical, theoretical and experimental results. For the subsequent studies,
the numerical models are used for complex loading in further sections.

###  3.3 Fractional poroelastic model validation

To further validate the proposed poroelastic model, a confined compression
creep test with and without initial ramp loading and confined compression
stress relaxation is done numercally and compared with experiments. We also
compute the weight loss from the confined compression creep test.

#### Confined compression creep test with initial ramp loading

In the earlier study, the initial ramp present in the experiments was assumed
as a step load. However, to assess the impact of this assumption on the model
fitting parameters, we ran the FE model with the initial ramp (similar to that
in the experiments) for the creep test and the results from the numerical
simulation are compared with experiments. The creep test with ramp loading
within Abaqus is modelled by providing the amplitude of the load in the load
module with a ramp for a specific time period similar to experiments and then
kept constant for the rest of the simulation. Figure 9 shows the comparison of
numerical results with the experiments for a few samples, in particular for
samples TK11BC, TK11BR, TK11BV and TK16AR. The results of creep with ramp in
Figure 9 show slight deviation at the end of ramp loading caused due to this
assumption. But the deviation is small and can be neglected. Therefore we
conclude that the material properties extraction process as explained above is
correct.

(a) TK11 Body Circumferential

(b) TK11 Body Radial

(c) TK11 Body Vertical

(d) TK16 Anteriorhorn Radial

Figure 9: Validation by a creep with ramp test (a) Sample taken from the body
region of TK11 sample in the circumferential direction, (b) Sample taken from
the body region of TK11 sample in the radial direction, (c) Sample taken from
the body region of TK11 sample in the vertical direction, (d)Sample taken from
the anterior horn region of TK16 sample in the radial direction.

#### Weight loss from confined compression creep test

Weight loss for confined creep compression problem presented in Section 2.1
and compared with the experiments carried out by Bulle et al.[5]. Specific
weight of water at 37∘C is used($w_{s}=$ 997Kg/m3). Fluid flow is collected at
the bottom of the cylindrical sample with the area of the circular cross-
section $\boldsymbol{A}$ with diameter $d=3mm$. Material properties are taken
from the Table 2 which were obtained from fitting the fractional poroelastic
model with displacement data from the confined creep test. Weight loss with
time for some samples is shown in the Figure 10. Theoretical and experimental
results for a few samples match as shown in the Figure 10. Moreover, they do
not match well for some samples, as shown in Figure 1. This could be due to
the measurement error as it was carried out manually. The weight of the sample
is measured after every 75 seconds of creep by taking it out.

(a) TK16 Body Circumfrential

(b) TK16 Body Radial

(c) TK16 Body Vertical

(d) TK17 Anteriorhorn Radial

Figure 10: Comparison of theoretical Weight loss from the sample using
parameters obtained from fitting with the experiments (a) Sample taken from
the body region of TK16 sample in the circumferential direction, (b) Sample
taken from the body region of TK16 sample in the radial direction, (c) Sample
taken from the body region of TK16 sample in the vertical direction, (d)
Sample taken from the Anterior horn region of TK17 sample in the radial
direction.

#### Confined compression stress relaxation

The fractional poroelastic model is validated using the confined compression
stress relaxation tests performed by Bulle et al. [5]. Tests were performed on
the cylindrical sample of height $h$, confined on the sides, and the fluid can
flow freely from the bottom. Tests consist of preconditioning the sample with
10%$h$ compression at the ramp velocity 0.3%$h$/s. Then, five relaxation steps
with 2%$h$ compression at a ramp velocity of 0.3%$h$/s s were performed. In
all the cases, the force applied is measured with time. The finite element
modelling of this confined compression stress relaxation test is performed in
Abaqus similar to the creep test as mentioned in Section 3.2. Material
properties obtained from fitting the creep data were used. In this test, the
displacement is given as input, and the force is measured as output. Figure 11
shows a comparison of the FEM results with the experiments. The comparisons of
the model with the experiments are reasonably good.

(a) TK11 Body Circumferential

(b) TK11 Body Vertical

(c) TK11 Body Radial

Figure 11: Validation by comparing stresses from Abaqus with stress relaxation
test from experiments (a) Sample taken from the body region of TK11 sample in
the circumferential direction, (b) Sample taken from the body region of TK11
sample in the vertical direction, (c) Sample taken from the body region of
TK11 sample in the radial direction.

##  4 Discussion

We show that the fractional poroelastic model is more appropriate to model the
time-dependent behaviour of soft tissue presenting a hierarchically arranged
porous architecture. The model has been verified and validated. Our parameter
fittings show that the RMS error = 1.42$\times$10-5 for the fractional
poroelastic model is considerably lower than the one given by the classical
model that had an RMS error = 6.53$\times$10-5, see Figure 4. The theoretical
displacement curve starts from 0, whereas the experiment shows an initial
displacement. This is attributed to the assumption of undrained conditions at
the initial time and the assumption of incompressible material. We obtain a
value of fractional order $\beta$ in the range of $\beta=$ 0.51–0.73 as shown
in Table 3. From the numerical study, it can be inferred that the transport
phenomenon during the confined compression is ruled by a fractional version of
Darcy’s law. At the beginning of the confined compression test, the solid is
fully saturated, therefore the pore pressure reaches a maximum value (see
Figure 3(c)). This is also predicted by Terzaghi’s consolidation theory, in
line with the fact that the fluid pressure entirely carries the load. During
the consolidation process, the fluid flows out of the sample at a rate
depending on the anomalous permeability and on the order of the fractional
derivative (see Figure 3). The higher the value of $\beta$ the faster the
diffusion process. As the fluid flows out of the sample, the pore pressure
decreases and the solid starts deforming (see Figure 3(d)). Also, higher
values of $\beta$ imply a faster solid deformation process. We notice that the
aggregate modulus $M$ obtained from fitting mentioned in Table 2 is of the
same order as the literature [19, 20, 21, 22]. The anomalous permeability
$\lambda_{\beta}$ is in the order of $10^{-12}/10^{-13}$. From the statistical
analysis, it can be concluded that the meniscus in the body region is
elastically isotropic, and pore pressure diffusion is transversely isotropic
with symmetry in the vertical and radial directions. Higher hydraulic
permeability in the circular direction can be attributed to the fact that the
fibres are oriented in a circular direction in the body region, see Figure
1(b). To model the anisotropy in the anomalous transport phenomena, the
following anisotropic form of fractional Darcy’s law is proposed:

$\mathbf{J}_{p}=\underbrace{\begin{bmatrix}\lambda_{rad}&0&0\\\
0&\lambda_{cir}&0\\\
0&0&\lambda_{ver}\end{bmatrix}}_{\mathbf{\Lambda}}\underbrace{\begin{bmatrix}\prescript{}{0}{D}_{t}^{\beta_{rad}}&0&0\\\
0&\prescript{}{0}{D}_{t}^{\beta_{cir}}&0\\\
0&0&\prescript{}{0}{D}_{t}^{\beta_{ver}}\\\
\end{bmatrix}}_{\mathbf{\prescript{}{0}{D}_{t}^{\mathbf{\beta}}}}\begin{bmatrix}\nabla
p\end{bmatrix}$ (14)

where, $\mathbf{\Lambda}$ is the permeability tensor and
$\mathbf{\prescript{}{0}{D}_{t}^{\mathbf{\beta}}}$ is the fractional
derivative operator. We notice that for the body region, the anomalous
permeability tensor $\mathbf{\Lambda}$ is transversely isotropic with
$\lambda_{\rm cir}$ considerably higher than $\lambda_{\rm ver}$ and
$\lambda_{\rm rad}$. Similarly $\mathbf{\beta}$ tensor is transversely
isotropic with $\beta_{\rm cir}$ considerably higher than the $\beta_{\rm
rad}$ and $\beta_{\rm ver}$. Consequently, from the Figure 7(c), it is
observed that the fluid flow rate through the sample’s base during the
consolidation experiment in the circumferential direction is higher than the
other two directions. Pore pressure in the circumferential direction reaches a
steady state faster than the other two directions, as shown in the Figure 7(a)
because of the higher permeability in the circumferential direction. Based on
the mean values and the standard deviation in the anterior horn region, it is
opined that the aggregate modulus in the radial direction is higher than the
other two directions, viz., vertical and circular directions. It is inferred
that the permeability and the order of derivative are higher in the radial and
vertical directions than in the circumferential directions. The elastic tensor
results are transversely isotropic, with the modulus in the radial direction
being higher, and both $\mathbf{\Lambda}$ and $\mathbf{\beta}$ are
transversely isotropic with circumferential being lower. Different properties
in different directions show that the elastic, anomalous permeability and the
order of the derivative tensors are transversely isotropic. Furthermore,
properties vary with the anterior horn and body region, showing that it is not
homogeneous.

Figure 12: Directional and regional dependence of material properties in the
meniscus.

In summary (see Figure 12), we show that the body region which has a load-
bearing function exhibits a transversely isotropic behaviour due to the rate
of fluid flow being about three times higher (faster diffusion) in the
circumferential direction which is consistent with the preferential direction
of collagen channels. It explains the role of fluid pressure in sustaining the
load. Furthermore, we show that both the elastic and the permeability tensors
are transversely isotropic in the anterior horn. The aggregate modulus is
higher in the radial direction compared to the circumferential and vertical
directions. It explains the role of the anterior horn to be compliant in the
circumferential and radial directions to accommodate the kinematics of the
tissue.

##  5 Conclusions

The meniscus has a porous, hierarchical, multi-oriented fiber and bundle
structure filled with fluid that provides optimized load support and
lubrication properties. To understand its function, confined compression creep
tests were performed on different regions of meniscus tissue and in different
orientations. Biot’s theory with fractional Darcy’s law is used to find the
material properties. It is observed that, with the classical Darcy’s law,
fitting gives an RMS error of 6.53$\times$10-5, while the fractional Darcy’s
law gives an RMS error of 1.42$\times$10-5. It is shown that Biot’s theory
with fractional Darcy’s law is better suited for modelling the poroelastic
behaviour of the meniscus. Three material parameters, viz., aggregate
modulus($M$), fractional order($\beta$) and the
permeability($\lambda_{\beta}$) required for the fractional Biot’s theory were
found by fitting the theory with the confined compression creep experiment
results for all the samples. Material properties in different directions of
the body region are compared using ANOVA test methods. It is observed that the
aggregate modulus ($M_{\rm cir}=$ 75.4 $\pm$48.9 KPa, $M_{\rm rad}=$
63.7$\pm$14.1 KPa and $M_{\rm ver}=$ 65.3$\pm$13.7 KPa) is isotropic and the
permeability ($\lambda_{\beta_{\rm cir}}=$ 3.75$\times$10${}^{-12}\pm$
2.36$\times$10-12m2/Pa.s1-β, $\lambda_{\beta_{\rm rad}}=$
1.25$\times$10${}^{-12}\pm$ 4.04$\times$10-13m2/Pa.s1-β and
$\lambda_{\beta_{\rm ver}}=$ 8.61$\times$10${}^{-13}\pm$
1.94$\times$10-13m2/Pa.s1-β) and the fractional order ($\beta_{\rm cir}=$
0.73$\pm$0.00, $\beta_{\rm rad}=0.64\pm 0.06$ and $\beta_{\rm ver}=$ 0.58$\pm$
0.04) are transversely isotropic with properties in the circumferential
direction is greater than the vertical and radial directions.

Biot’s theory with fractional Darcy’s law is implemented numerically using the
FEM in Abaqus software using UMATHT and UMAT subroutines. This fractional
poroelastic theory is validated by using the material properties obtained from
fitting confined creep tests to model stress relaxation, creep with ramp, and
weight loss tests. The results are then compared with experimental results.
The comparison shows that Biot’s theory with fractional Darcy’s law is capable
of better characterising the poroelastic behaviour of the meniscus. Flux out
of the meniscus, which is difficult to measure experimentally, can easily be
computed efficiently using the current model. In future works, Biot’s theory
with fractional Darcy’s law will be extended to include anisotropy and use it
to understand the working of the meniscus in the knee joint.

##  6 Acknowledgements

O.B would like to acknowledge the European Union’s Horizon 2020 -EU.1.3.2. -
Nurturing excellence by means of cross-border and cross-sector mobility under
the Marie Skłodowska-Curie individual fellowship MSCA-IF-2017, MetaBioMec,
Grant agreement ID: 796405. The authors thank the Rizzoli Orthopaedic
Institute and, in particular, G. Marchiori and M. Berni for their invaluable
insights on the experiments.

## References

  * [1] G. Agustoni, F. P. Bonomo, S. P. A. Bordas, O. Barrera, High resolution micro-computed tomography reveals a network of collagen channels in the body region of the knee meniscus, Annals of Biomedical Engineering (2021).
  * [2] J. Maritz, G. Agustoni, K. Dragnevski, S. P. A. Bordas, O. Barrera, The functionally grading elastic and viscoelastic properties of the body region of the knee meniscus, Annals of Biomedical Engineering (2021).
  * [3] F. P. Bonomo, J. J. Gregory, O. Barrera, A procedure for slicing and characterizing soft heterogeneous and irregular-shaped tissue, Materials Today: Proceedings (2020).
  * [4] V. Vetri, K. Dragnevsk, M. Tkaczyk, M. Zingales, G. Marchiori, N. Lopomo, S. Zaffagnini, A. Bondi, J. Kennedy, D. Murray, O. Barrera, Advanced microscopy analysis of the micro-nanoscale architecture of human menisci., Scientific Reports (2019).
  * [5] R. Bulle, G. Alotta, G. Marchiori, M. Berni, N. F. Lopomo, S. Zaffagnini, S. P. A. Bordas, O. Barrera, The human meniscus behaves as a functionally graded fractional porous medium under confined compression conditions, Applied sciences 11 (20) (2021).
  * [6] G. A. Ateshian, J. A. Weiss, Anisotropic hydraulic permeability under finite deformation, Journal of Biomechanical Engineering 132 (11) (2010).
  * [7] H. Kurosawa, T. Fukubayashi, H. Nakajima, Load-bearing mode of the knee joint: physical behavior of the knee joint with or without menisci., Clinical Orthopaedics and Related Research 149 (1980) 283–290.
  * [8] N. Shrive, J. O’Connor, J. Goodfellow, Load-bearing in the knee joint., Clinical Orthopaedics and Related Research 131 (1978) 279–287.
  * [9] J. Waghorne, E. Elmukashfi, S. L. Giudice, V. Manuri, G. Pitarresi, O. Barrera, On the structure-function relationships of the meniscal tissue. a data driven approach to map soft tissues performance, To be submitted to Nature Material (2022).
  * [10] R. Readioff, R. Seil, C. Mouton, L. Marks, O. Barrera, An optimised patient-specific finite element model to study the influence of the intra-articular parameters on the contact mechanics in the knee, To be submitted to Journal of Biomechanical Engineering (2022).
  * [11] O. Barrera, et al., On the characteristics of natural hydraulic dampers: an image driven based approach to study the fluid flow behavior inside the human meniscal tissue., To be submitted to PNAS (2022).
  * [12] E. Elmukashfi, G. Marchiori, M. Berni, G. Cassiolas, N. F. Lopomo, H. Rappel, M. Girolami, O. Barrera, Model selection and sensitivity analysis in the biomechanics of soft tissues: A case study on the human knee meniscus, Advances in Applied Mechanics, Elsevier, 2022.
  * [13] A. van Kampen, The knee joint in sports medicine, International Orthopaedics 37 (2) (2013) 177–179.
  * [14] A. J. S. Fox, A. Bedi, S. A. Rodeo, The basic science of human knee menisci: structure, composition, and function, Sports Health 4 (4) (2012) 340–351.
  * [15] M. A. Soltz, G. A. Ateshian, Experimental verification and theoretical prediction of cartilage interstitial fluid pressurization at an impermeable contact interface in confined compression, Journal of Biomechanics 31 (1998) 927–934.
  * [16] M. A. Sweigart, C. F. Zhu, D. M. Burt, P. D. Deholl, C. M. Agrawal, T. O. Clanton, K.A.Athanasiou, Intraspecies and interspecies comparison of the compressive properties of the medial meniscus, Annals of Biomedical Engineering 32 (11) (2004) 1569–1579.
  * [17] O. Barrera, A unified modelling and simulation for coupled anomalous transport in porous media and its finite element implementation, Computational Mechanics 68 (2021) 1267–1282.
  * [18] A. M. Seitz, F. Galbusera, C. Krais, A. Ignatius, L. Dürselen, Stress-relaxation response of human menisci under confined compression conditions, Journal of the Mechanical Behavior of Biomedical Materials 26 (2013) 68 – 80.
  * [19] M. A. Sweigart, C. F. Zhu, D. M. Burt, P. D. Deholl, C. M. Agrawal, T. O. Clanton, K. A. Athanasiou, Intraspecies and interspecies comparison of the compressive properties of the medial meniscus, Annals of Biomedical Engineering 32 (11) (2004) 1569–1579.
  * [20] J. T. Moyer, R. Priest, T. Buman, A. C. Abraham, T. L. H. Donahue, Indentation properties and glycosaminoglycan content of human menisci in the deep zone, Acta Biomaterialia 9 (5) (2013) 6624–6629.
  * [21] A. M. Seitz, F. Galbusera, C. Krais, A. Ignatius, L. Durselen, Stress-relaxation response of human menisci under confined compression conditions, Journal of the Mechanical Behavior of Biomedical Materials 26 (2013) 68–80.
  * [22] H. N. Chia, M. L. Hull, Compressive moduli of the human medial meniscus in the axial and radial directions at equilibrium and at a physiological strain rate, Journal of Orthopaedic Research 26 (7) (2008) 951–956.
  * [23] I. Podlubny, Fractional Differential Equations, Academic Press, 1999.

## Appendix A

Using Equation (4a), the displacement field for the problem described in
Section 2.1 can be solved as follows:

$\epsilon_{zz}=\dfrac{\partial u_{z}}{\partial
z}=\dfrac{3}{3K+4G}\Big{(}-P_{A}+\alpha p\Big{)}$ (15)

Substituting Equation (6) in Equation (15) and integrating

$u(z,t)=\dfrac{3P_{A}}{3K+4G}\Bigg{[}-z+\gamma\alpha\sum\limits_{n=1,3}^{\infty}E_{1-\beta,1}\left(-{\dfrac{n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}}\right)b_{n}\sin\dfrac{n\pi
z}{2h}\Bigg{]}+d$ (16)

where, $b_{n}=\frac{8h}{(n\pi)^{2}}\left(-1\right)^{\frac{n-1}{2}}$. Applying
the boundary condition given by Section 2.2, we get:

$u(h,t)=\frac{3P_{A}}{3K+4G}\Bigg{[}-h+\gamma\alpha\sum\limits_{n=1,3}^{\infty}E_{1-\beta,1}\left(-{\frac{n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}}\right)\frac{8h}{(n\pi)^{2}}\left(-1\right)^{\frac{n-1}{2}}\sin\frac{n\pi}{2}\Bigg{]}+d=0$
(17)

Since,

$\Big{(}-1\Big{)}^{\frac{n-1}{2}}\sin\frac{n\pi}{2}=1\hskip
20.0pt\text{for}\qquad n=1,3,\cdots$ (18)

Therefore, the constant $d$ is

$d=-\dfrac{3P_{A}}{3K+4G}\Bigg{[}-h+\gamma\alpha\sum\limits_{n=1,3}^{\infty}E_{1-\beta,1}\left(-{\frac{n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}}\right)\dfrac{8h}{(n\pi)^{2}}\Bigg{]}$
(19)

Substituting Equation (19) in Equation (16), Displacement field equation is
obtained

$u(z,t)=\frac{3P_{A}}{3K+4G}\Bigg{[}(h-z)+\gamma\alpha\sum\limits_{n=1,3}^{\infty}E_{1-\beta,1}\left(-{\frac{n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}}\right)\dfrac{8h}{(n\pi)^{2}}\bigg{(}(-1)^{\frac{n-1}{2}}\sin\dfrac{n\pi
z}{2h}-1\bigg{)}\Bigg{]}$ (20)

## Appendix B

Weight loss from the sample over time W(t) for the problem described in
Section 2.1 can be solved using the fluid flux as:

$W(t)=W_{0}-\int_{0}^{t}\mathbf{J}_{p}\cdot w_{s}\mathbf{A}dt$ (21)

where $\mathbf{J}_{p}$ is the fluid flux given by Equation (2b), $w_{s}$ is
the specific weight, $\mathbf{A}$ is the cross-sectional area over which the
flux is calculated and $W_{0}$ being the initial weight of the sample.
Differentiating Equation (6) with respect to $z$, we get:

$\frac{\partial p}{\partial
z}=-P_{A}\gamma\sum_{n=1,3}^{\infty}E_{1-\beta,1}\bigg{[}-\frac{n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}\bigg{]}\dfrac{2}{h}(-1)^{\frac{n-1}{2}}\sin\bigg{(}\frac{n\pi
z}{2h}\bigg{)}$ (22)

Fluid can flow only at the bottom, as shown in Figure 2(a). Hence, the weight
loss is solved at $z=h$. Equation (22) at $z=h$ is:

$\dfrac{\partial p}{\partial
z}=-P_{A}\gamma\sum_{n=1,3}^{\infty}E_{1-\beta,1}\bigg{[}-\frac{n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}\bigg{]}\dfrac{2}{h}$
(23)

Using Equation (23) and Equation (2b) in weight loss Equation (21), we get:

$W(t)=W_{0}-\lambda_{\beta}w_{s}A\prescript{}{0}{D}_{t}^{\beta-1}\bigg{[}P_{A}\gamma\sum_{n=1,3}^{\infty}E_{1-\beta,1}\bigg{(}-\frac{n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}\bigg{)}\dfrac{2}{h}\bigg{]}$
(24)

Using the identity of the fractional derivative from [23], which is:

$\prescript{}{0}{D}_{t}^{\alpha}t^{\beta-1}E_{\mu,\beta}(\lambda
t^{\mu})=t^{\beta-1-\alpha}E_{\mu,\beta-\alpha}(\lambda t^{\mu})$ (25)

Equation (24) can be found as:

$W(t)=W_{0}-\lambda_{\beta}w_{s}AP_{A}\gamma
t^{1-\beta}\sum_{n=1,3}^{\infty}E_{1-\beta,2-\beta}\bigg{(}-\dfrac{n^{2}\pi^{2}\bar{\lambda}t^{1-\beta}}{4h^{2}}\bigg{)}\dfrac{2}{h}$
(26)

## Appendix C

S No. | Sample | M $\times 10^{5}(Pa)$ | $\lambda\times 10^{-13}(m^{2}/Pa.s)$ | RMS $\times 10^{-5}$
---|---|---|---|---
1 | TK11BC | 3.56 | 6.95 | 6.53
2 | TK11BR | 2.13 | 3.91 | 7.40
3 | TK11BV | 2.30 | 4.44 | 8.10
4 | TK16BC1 | 2.48 | 5.32 | 7.04
5 | TK16BR2 | 2.55 | 2.45 | 5.53
6 | TK16BV | 2.50 | 3.00 | 6.37
7 | TK16AC1 | 2.55 | 2.80 | 7.64
8 | TK16AR1 | 2.64 | 2.42 | 7.16
9 | TK16AV | 2.69 | 1.79 | 5.83
10 | TK16BC2 | 2.39 | 8.87 | 8.15
11 | TK16BR1 | 2.31 | 2.68 | 7.12
12 | TK16BV2 | 2.56 | 2.74 | 6.96
13 | TK16PV | 3.97 | 12.1 | 4.66
14 | TK17BC | 1.46 | 11.8 | 13.4
15 | TK17BR | 1.67 | 4.17 | 6.93
16 | TK17BV | 1.80 | 1.99 | 5.36
17 | TK17AC | 3.49 | 208 | 4.58
18 | TK17AR | 4.01 | 23.3 | 4.49
19 | TK18BC | 1.44 | 1.33e4 | 10.5
20 | TK18BR | 2.37 | 5.27 | 7.54
21 | TK18BV | 2.17 | 2.89 | 6.45
22 | TK22BR | 2.55 | 4.66e3 | 1.24
23 | TK22AC | 2.41 | 3.69 | 9.01
24 | TK22AR | 1.95 | 4.74 | 10.0
25 | TK22AV | 2.97 | 1.31 | 4.83
26 | TK36BC | 1.98 | 9.33 | 10.2
27 | TK37BC | 2.52 | 16.5 | 12.4
28 | TK37BR | 2.67 | 3.23 | 7.42
29 | TK37BV | 2.29 | 2.66 | 6.25
Table 1: Material properties obtained from fitting using classical Biot’s
theory

## Appendix D

(a) TK11 Body Circ

(b) TK11 Body Radial

(c) TK11 Body Vertical

(d) TK16 Anterior horn Radial

Figure 1: Comparison of theoretical weight loss from the sample using
parameters obtained from fitting with the experiments (a) Sample taken from
the body region of TK11 sample in the circumferential direction, (b) Sample
taken from the body region of TK11 sample in the radial direction, (c) Sample
taken from the body region of TK11 sample in the vertical direction, (d)
Sample taken from the Anthorn region of TK16 sample in the radial direction,