# Morphometric description of strength and degradation in porous media

A. Guével<EMAIL_ADDRESS>H. Rattez  E. Veveakis  Civil and
Environmental Engineering Department, Duke University, Durham, NC 27708-0287,
USA  Institute of Mechanics, Materials and Civil Engineering (IMMC),
Université Catholique de Louvain, B-1348 Louvain-la-Neuve, Belgium

###### Abstract

The influence of the microstructural geometry on the behavior of porous media
is widely recognized, particularly in geomaterials, but also in biomaterials
and engineered materials. Recent advances in imaging techniques, such as X-ray
microcomputed tomography, and in modeling make it possible to capture the
exact morphometry of the microstructure with high precision. However, most
existing continuum theories only partially account for the morphometry. We
propose here a unifying approach to link the strength of porous materials with
the necessary and sufficient microstructural information, using Minkowski
functionals, as per Hadwiger’s theorem. A morphometric strength law is
inferred from synthetic microstructures with a wide range of porosities and
heterogeneities, through qualitative 2D phase-field simulations. Namely, the
damage is modeled at the microstructural level by tracking the solid-pore
interfaces under mechanical loading. The strength is found to be best
described by an exponential function of the morphometers, thus generalizing
early works on metals and ceramics. We then show that the predictiveness of
this relationship extends to real porous media, including rocks and bones.

###### keywords:

porous media, morphometry, strength, phase-field modeling, degradation,
damage, Minkowski functionals

††journal: Journal of the Mechanics and Physics of Solids

##  1 Introduction

Porous media represent a wide range of materials but also a tremendous
challenge to be fully understood and harnessed. Among them, geomaterials,
stemming from millions of years of transformations under harsh conditions,
represent a particularly complex subclass, inasmuch as processes in these
media are multiphysics and multiscales. Recent advances in geosciences found,
however, that this complexity may boil down to the great heterogeneity and
stochasticity of geomaterials’ microstructures. For instance, pressure
solution (Niemeijer et al. (2009); Croizé et al. (2013); van den Ende et al.
(2019); Guével et al. (2020)), strain localization (Vardoulakis and Sulem
(1995); Kawamoto et al. (2018)), frictional instabilities (Rattez et al.
(2018a, b)), fault reactivation (Veveakis et al. (2014); Lesueur et al.
(2020)), and granular flow (Buscarnera and Einav (2021)) largely depend on the
microstructural geometry, or morphometry. The same conclusion holds for
engineered porous materials as well, such as ceramics (Salvini et al. (2018))
and energetic materials (Chun et al. (2020)). Biomaterials, forming a third
subclass of porous media (Huyghe et al. (2002)), are also increasingly studied
in the light of their microstructures, in particular in bones mechanics (Augat
and Schorlemmer (2006); Wachter et al. (2002)).

Therefore, modeling the microstructural dynamics is a crucial step towards
better understanding the macroscopic behavior of porous media. In metallurgy,
this has been achieved with phase-field modeling (Allen and Cahn (1979);
Provatas and Elder (2010); Bhattacharyya et al. (2019)), which tracks the
grains’ interfaces in finite elements, and in geosciences, with Discrete
Element Modeling (O’Sullivan (2011); Kawamoto et al. (2018)), which tracks the
grains as discrete elements, and recently also with phase-field modeling
(Guével et al. (2020)). That said, due to computational limitations, relying
solely on microscopic simulations for large scales is unrealistic. On the
other hand, the existing constitutive macroscopic laws rely heavily on
destructive experimental calibration, which can be unfeasible because of the
limited availability of materials or the impossibility to reproduce the
environmental conditions. For example, obtaining rocks from high depths can be
prohibitively expensive or even impossible, and geological time scales are
hardly reproducible in laboratory. As for biomaterials, such as bones, they
are best studied in vivo (see Fragogeorgi et al. (2019) e.g.).

The correlation between the strength of materials and their microstructural
geometry was first theorized, through the foundation of poromechanics, by Biot
(1941), and for soils specifically, by Terzaghi (1943). This effort was
intensified following the catastrophic failure of welded Allies ships through
brittle fracture during the Second World War, with focus on metals and
ceramics. The strength of brittle polycrystalline metals was quantitatively
addressed by Orowan (1949), proposing what would ultimately become the
Hall–Petch effect (Hall (1951); Petch (1953)). The Hall–Petch relationship
asserts that the yield stress, or similarly the flow stress, scales as the
inverse square root of the mean grain size. However, Li et al. (2016) recently
showed, upon critically reviewing a large amount of experimental data, that
other types of laws, including inverse exponential, turn out as satisfying.
Furthermore, Ryshkewitch (1953) experimentally found that the strength of
porous brittle polycrystalline materials, such as ceramics, is an inverse
exponential function of the porosity. This was confirmed more recently with
porosity-controlled experiments and modeling by Liu et al. (2017). Building
upon the previous results from non-porous and porous brittle materials,
Knudsen (1959) inferred a multiplicative dependence of the strength on grain
size and porosity (Eq. (7) therein), which was checked experimentally for
ceramics at different temperatures. These efforts culminated to the advent of
modern poromechanics (see Coussy (2004); Dormieux et al. (2006)), reconciling
continuum mechanics with Biot’s theory, and relevant to a wide range of
disciplines, including geo, bio and material sciences. However, the strength
of porous media was mostly studied with respect to porosity and grain size.
According to morphometry theory and in particular Hadwiger’s theorem (see
Armstrong et al. (2019) and Section 2.1), two more descriptors of the
microstructure, which we refer to as morphometers, are needed to fully account
for the morphometric effects. On these bases, we propose to investigate a
morphometric strength law of the form:

$\sigma_{f}=\sigma_{f}^{*}f(M_{0},M_{1},M_{2},M_{3})$ (1)

where we choose to measure the strength via a flow stress $\sigma_{f}$ (see
Section 3.1), the $M_{i}$ (i=0,…,3) denote the four required morphometers, the
star quantities are reference values depending on the material and
environmental conditions, and $f$ is a multiplicative function of the $M_{i}$,
following Knudsen (1959). In particular, $M_{0}$ denotes the porosity and
$M_{2}$ the grain size. Our objective is to investigate the possible
application of this law to porous materials in general, including geomaterials
and biomaterials. Thereby, we subsume the microscopic mechanisms responsible
for failure under the generic umbrella of damage. For geomaterials, this
encompasses debonding, dissolution, cracks and breakage. Recent investigations
in geosciences on the influence of the morphometry notably include the work of
Zhang et al. (2016), Wetzel et al. (2021) and Buscarnera and Einav (2021). For
biomaterials, in particular for bones, damage comprises dissolution, as in
osteoporosis, and cracks, particularly in cortical bones (Augat and
Schorlemmer (2006)). While the focus in biomechanics is mostly on the
relationship between elastic moduli and morphometers, some results corroborate
the conclusions inferred for the previous subclasses of porous materials, such
as the exponential relationship between strength and porosity (see Wachter et
al. (2002) and A). However, a mathematically-consistent law such as Eq. (1)
has not been proposed yet, accounting for the necessary and sufficient
microstructural information.

In order to numerically investigate the dependency of the strength of porous
materials on the morphometers and specify Eq. (1), we control the variability
of the morphometry by using synthetic microstructures (SMs). This also allows
to run a sufficiently large ensemble of simulations to be statistically
meaningful and account for the stochastic nature of real microstructures (Chun
et al. (2020)). Interest in synthesizing microstructures is intensifying in
material sciences, owing to recent computational advances. In particular,
machine learning enables to create highly realistic, tunable microstructures
(Mosser et al. (2017); Chun et al. (2020)). While the latter represents the
state of the art, we have used, as a first step, a more straightforward and
much less computationally-demanding approach (see Section 2.2).

In all, we aim here at bridging the gap between microscopic and macroscopic
modeling, by calibrating macroscopic laws with morphometric parameters
upscaled from simulations on synthetic microstructures, hinging, as much as
possible, on non-destructive methods. We should then check that the
predictiveness of this relationship extends to real porous media, by using CT
scans. While complex multiphysics processes like pressure solution may still
require explicit modeling of the microstructure (Guével et al. (2020)), we
suggest that as far as the influence of the microstructural geometry is
concerned, it is possible to upscale only the essential morphometric
information.

##  2 Methods

###  2.1 Upscaling of the microstructural information

The morphometers $M_{i}$ in Eq. (1) are appropriately described by the
Minkowski functionals of the domain formed by the grains $\Omega$ (Armstrong
et al. (2019)). Hadwiger’s theorem guarantees that a microstructure is fully
described by $d+1$ Minkowski functionals, where $d$ is the microstructure’s
dimension, in the sense that any other descriptor that is additive, motion-
invariant and conditionally continuous would be a linear combination of those
functionals (see Hadwiger (1951) and also Klain (1995) for a short proof). In
3D, the 4 needed Minkowski functionals are (Armstrong et al. (2019))

$M_{0}(\Omega)=\int_{\Omega}\text{d}V,$ (2)

the total volume of the grains,

$M_{1}(\Omega)=\int_{\partial\Omega}\text{d}S,$ (3)

their total surface area,

$M_{2}(\Omega)=\int_{\partial\Omega}\left(1/r_{1}+1/r_{2}\right)\text{d}S,$
(4)

their total mean curvature, where $r_{1}$ and $r_{2}$ denote the principal
radii of curvature of the surface element $\text{d}S$, and

$M_{3}(\Omega)=\int_{\partial\Omega}\left(1/r_{1}r_{2}\right)\text{d}S=4\pi\chi(\Omega),$
(5)

their total Gaussian curvature, directly related to the Euler characteristic
$\chi$ by the Gauss-Bonnet theorem (Armstrong et al. (2019)). In practice, we
will use the porosity $n$ as a measure of $M_{0}$ and the mean grain size for
$M_{2}$ , the latter approximating the inverse of the average mean curvature.
In this contribution, we will restrict our results to $d=2$ as a first step,
and therefore use only 3 morphometers to constrain the strength law Eq. (1).

###  2.2 Synthetic microstructures

To find the function $f(M_{i})$ in Eq. (1), we run a large amount of
simulations on SMs, which are generated using the Python open-source package
PoreSpy 111http://porespy.org (Gostick et al. (2019)). The starting point is a
random noise (see Fig. 1a), that is, a 2D array $R[i,j]$
($i=0,\dots,N-1$,$j=0,\dots,N-1$) of dimension $N$ of random values between 0
and 1. It is then convoluted with a Gaussian filter (see Fig. 1b)

$G[i,j]=e^{-\frac{i^{2}+j^{2}}{2s^{2}}},$ (6)

which standard deviation $s=\frac{N}{40h}$, or blurriness, is inversely
proportional to a heterogeneity parameter $h$, and where 40 is a scaling
factor; the heterogeneity of the output is then independent from the image
size. The filtered output $F[i,j]$ is thus calculated as follows (see Gonzalez
and Woods (2008) e.g.):

$F[i,j]=R[i,j]*G[i,j]=\sum_{u=0}^{N-1}\sum_{v=0}^{N-1}R[u,v]G[i-u,j-v],$ (7)

where $*$ denotes the convolution operator. Finally, after uniformizing the
blurred noise $F[i,j]$ between 0 and 1, the array is binarized into $B[i,j]$
(see Fig. 1c) by using the porosity $n$ as the threshold, yielding the final
SM. The two input parameters controling the generation of the SMs are thus the
heterogeneity $h$ and the porosity $n$. We will see that $h$ controls the
narrowness and skewness of the morphometers distributions.

\begin{overpic}[width=130.08731pt]{random_noise.png}
\put(-8.0,87.0){\color[rgb]{0,0,0} $a)$} \end{overpic}

\begin{overpic}[width=130.08731pt]{gaussian_blur_sigma167_blob3.png}
\put(-8.0,87.0){\color[rgb]{0,0,0} $b)$} \end{overpic}

\begin{overpic}[width=130.08731pt]{poro03_blob3.png}
\put(-8.0,87.0){\color[rgb]{0,0,0} $c)$} \end{overpic}

\begin{overpic}[width=130.08731pt]{watershed.png}
\put(-8.0,87.0){\color[rgb]{0,0,0} $d)$} \end{overpic}

Figure 1: a) Random noise $R[x,y]$, i.e. array of random values between 0 and
1. b) Gaussian filter $F[x,y]$ applied onto the random noise, with a standard
deviation of 1.76, or equivalently, $h=3$. c) Resulting binarized SM $B[x,y]$
upon applying the threshold $n=0.3$. d) SM segmented into 110 grains via
watershed algorithm.

To restrict the study to realistic microstructures, $n$ is varied from $0.15$
to $0.4$ and $h$ from $2$ to $5$, creating an ensemble of 42 SMs (see B). This
ensemble forms the training batch, from which the dependence of the flow
stress on the morphometry will be inferred. Throughout this work, the
dimensions of the synthetic SMs are restricted to $200\times 200$ pixels, as a
compromise between sufficiently realistic microstructures and reasonable mesh
sizes, proportional to this resolution. It is then possible to obtain the
distributions of the morphometers $M_{0}$, $M_{1}$, $M_{2}$, i.e. of the pore
size, perimeter and grain size, respectively, starting with the marker-based
watershed segmentation algorithm introduced by Gostick (2017). The values of
the filter size and of the Gaussian blur therein must be chosen carefully to
maximize the number of segmented regions (grains or pores). In particular,
since the optimal values for those parameters vary with the grain/pore size,
they are expected to vary with $h$. Namely, we found that the filter size
varies from 2 to 9 pixels for our range of SMs, whereas the blur varies from
0.2 to 0.8. Upon segmenting the pores, their sizes are obtained by determining
the equivalent disks for every pores. Upon segmenting the grains, their
perimeter is obtained via an edge-finding algorithm optimized for fast
computing introduced by Benkrid et al. (2000). The grain sizes are obtained
similarly to the pore sizes. The last morphometer $M_{3}$, the Euler
characteristic, is not concerned by distributions since it is a topological
quantity. While the distributions of $M_{0}$ and $M_{1}$ are useful when
comparing the SMs with real microstructures (see Section 3.4), only the grain
size distribution is calculated for every SMs, in order to determine the mean
grain size $M_{2}$. The latter is taken as the mean of the lognormal fit of
the grain size distribution, as explained in Section 3.4. $M_{0}$ and $M_{1}$
are directly obtained from calculating the porosity, i.e. the ratio of the
black pixels over the total number of pixels, and the total perimeter,
respectively. Finally, the Euler characteristic $M_{3}$ is calculated from an
integral geometry formula in discretized space. In the current 2D restriction,
$M_{3}$ is simply the number of objects minus the number of holes. We refer to
the documentation of PoreSpy and skimage available online for further details
(see C).

###  2.3 Damage phase-field modeling of the microstructure

The evolution of the microstructure is modeled via a phase field $\varphi$
differentiating the pores from the grains, as introduced in Guével et al.
(2020), and describing the damage of the grains. As such, the phase field
measures the changes of morphometry, and thus controls the evolution of the
morphometers (see Section 3.6). For all simulations, we use the open-source
multiphysics finite-element platform MOOSE (Permann et al. (2020)). Following
Guével et al. (2020) and assuming that the grains cannot heal, we consider the
following Allen–Cahn equation, coupled with elastic mechanics:

$\begin{cases}\tau\dot{\varphi}=\kappa\triangle\varphi-
Gg^{\prime}(\varphi)-\lambda(\bm{\epsilon})h^{\prime}(\varphi),\\\
\nabla\cdot\bm{C}(\varphi)\bm{\epsilon}=0,\end{cases}$ (8)

where $\tau$ is the relaxation time, $\kappa$ the interfacial energy
coefficient, $G$ the double-well’s height, $\bm{\epsilon}$ the strain tensor,
$\triangle$ and $\nabla\cdot$ denote the Laplacian and divergence operators,
the dot denotes the time derivative and the prime the derivative with respect
to $\varphi$. The polynomials $g(\varphi)=\varphi^{2}(1-\varphi)^{2}$ and
$h(\varphi)=\varphi^{2}(3-2\varphi)$ correspond to a double-well potential and
an interpolation function, respectively. Furthermore,
$\lambda(\bm{\epsilon})=1/2(\bm{C}_{g}-\bm{C}_{p})\bm{\epsilon}\cdot\bm{\epsilon}$
denotes the difference of elastic energy between the grains and pores phases.
Therein, $\bm{C}_{g}$ and $\bm{C}_{p}$ denote the elastic tensor of the grains
and the pores, respectively. Finally,
$\bm{C}(\varphi)=(1-h(\varphi))\bm{C}_{p}+h(\varphi)\bm{C}_{g}$ is the elastic
tensor of the mixture. In doing so, we apply the Voigt-Taylor homogenization
scheme (see Ammar et al. (2009) e.g.), which consists in interpolating the
partial stresses of the grains and pores phases while assuming homogeneous
strains, thus recovering a familiar poromechanics formulation. The damage of
the grains occurs when the elastic energy is large enough to tilt the double-
well potential beyond the saddle-node bifurcation point (Guével et al.
(2020)).

This model stems from the coupling of the microscale mechanics (at the grains
scale) and the macroscale mechanics (mixture of grains and pores). Namely, Eq.
(8)-1 stems from the micro-momentum balance (Fried and Gurtin (1993); Gurtin
(1996))

$\nabla\cdot\bm{\xi}+\pi=0,$ (9)

where $\pi$ is the microforce (scalar), energy-conjugate of $\varphi$,
$\bm{\xi}$ is the microstress (vector), energy-conjugate of $\nabla\varphi$,
whereas Eq. (8)-2 corresponds to the usual macro-momentum balance. They are
coupled via the generalized relaxation equation (Guével et al. (2020)),
equivalent to the second law of thermodynamics as described by the following
dissipation inequality

$-\dot{\psi}-\pi\dot{\varphi}+\bm{\xi}\cdot\nabla\dot{\varphi}-\bm{\sigma}:\dot{\bm{\epsilon}}\leq
0,$ (10)

where $\psi$ denotes the free energy, described in its usual Landau form
(Landau (1937))

$\psi(\varphi,\nabla\varphi,\bm{\epsilon})=Gg(\varphi)+(1-h(\phi))\frac{1}{2}\bm{\epsilon}\cdot\bm{C}_{p}\bm{\epsilon}+h(\varphi)\frac{1}{2}\bm{\epsilon}\cdot\bm{C}_{g}\bm{\epsilon}+\frac{\kappa}{2}\lvert\nabla\phi\rvert^{2},$
(11)

and $\bm{\sigma}=\bm{C}(\varphi)\bm{\epsilon}$ is the Cauchy stress tensor,
resulting from the interpolation of the partial stresses. Therein, the state
variables are $\varphi$, describing the normal variations of the grains
interface, $\nabla\varphi$, describing the tangential variations of the grains
interface, and the macroscopic (elastic) strain tensor $\bm{\epsilon}$. For
simplicity, we consider $\nabla\varphi$ and $\bm{\epsilon}$ to be non-
dissipative, thus neglecting viscous effects; the viscous effect associated
with the phase field would delay phase changes as shown in Guével et al.
(2020).

Here, damage is not specified but includes dissolution, debonding and
microcracks. In particular, the final system of equations (8) is similar to
the one used for modeling fractures at the continuum scale (see Kuhn and
Müller (2010) e.g.). The fundamental difference is that our model is applied
directly at the grains scale, whereas phase-field damage models have so far
been applied to the continuum scale. Hence, instead of the usual
differentiation between the intact and the damaged phases, our model
differentiates between the grains and the pores phases, where the growth of
the pores phase at the expense of the grains phase embodies the
microstructural damage. We distinguish the latter from its upscaled
manifestation at the continuum scale, which we call the degradation (see
Section 3.6). In particular, while the damaged phase in the macroscopic theory
is not energetic, here, the pores phase is associated with an elastic energy
(albeit not allowing shearing), yielding a mixture of elastic energies, as
described above. Therefore, the only differences in the form of our phase-
field equation with the continuum damage phase-field modeling are the presence
of an elastic energy for the phase $\varphi=0$ and the degrees of the
polynomials representing the potential $g(\varphi)$ and the interpolation
function $h(\varphi)$. Our model can also be compared to the one used for
grain growth in (non-porous) metals under mechanical loading (see Tonks and
Millett (2011) e.g.), where each grain of a polycrystal is represented by an
order parameter. Unlike for metals, the pores phase in porous media such as
geomaterials plays a crucial role, so that it is given an independent role in
our model. The main difference between modeling porous media and metals is
thus that we discriminate the pores from the skeleton, without
differentiating, as a first step, the different grains, inasmuch as they are
not necessarily distinguishable within the skeleton, depending on the
materials; for instance, as discussed in the following, unlike a sandstone, a
cortical bone can hardly be decomposed into grains. In practice, this
translates into dropping the grain-to-grain interfacial energy term present in
the polycrystals phase-field model. In all, modeling the damage of the
microstructure of porous media is performed here as a combination of the
concepts used in continuum damage and polycrystals phase-field modeling.

Upon choosing a characteristic length $l_{0}$, a characteristic time
$t_{0}=\tau/G$ and a characteristic specific energy $G=1GPa$, we write (8) in
dimensionless form as follows:

$\begin{cases}\dot{\varphi}=\hat{\kappa}\triangle\varphi-g^{\prime}(\varphi)-\hat{\lambda}(\bm{\epsilon})h^{\prime}(\varphi),\\\
\nabla\cdot\hat{\bm{C}}(\varphi)\bm{\epsilon}=0,\end{cases}$ (12)

where $\hat{\kappa}=\kappa/Gl_{0}^{2}$ is the (dimensionless) interfacial
group, $\hat{\lambda}=\lambda/G$ is the activation energy group and
$\hat{\bm{C}}=\bm{C}/G$. The dimensionless derivatives are noted similarly and
the hat notation for dimensionless quantities is dropped in the following. The
evolution of the damage via $\varphi$ is fully determined by $\lambda$, upon
fixing $\kappa$ to an appropriately small value (see Section 3.3). We also
restrict our attention for simplicity to 2D problems, so that the elastic
energy $\lambda$ is fully determined by the two Lamé parameters. As explained
in Guével et al. (2020), the grains are considered as solids of very low
porosity, and the pores as shear-free solids that are much more deformable.
Namely, for geomaterials, we choose the first and second Lamé parameters of
the grains both to be $30GPa$, and for the pores $1GPa$ and 0, respectively.
Further major assumptions include small strains and that the assembly of
grains considered is representative of the material. To focus on the effect of
the microstructure on the strength, we will first perform simulations in
tension with fixed lateral boundaries. Then, to determine the effect of the
confinement, we will perform axisymmetric biaxial loadings. Unless mentioned
otherwise, the stress used in output throughout this work is the vertical
stress averaged on the top boundary of the domain. The associated strain is
the vertical displacement of the top boundary divided by the initial length of
the square domain.

##  3 Model’s performance

###  3.1 Measure of the strength

In solid mechanics, the word ”strength” has a broad meaning. Most often, it
corresponds to the maximum value that the stress can attain (peak stress).
Alternative measures of strength include the yield stress, or even the value
of stress at a given value of strain when the material is undergoing prolonged
hardening (Sari et al. (2020)). Since the latter is the case in the present
study, the strength is defined throughout this work as a value of the vertical
stress for a certain value of a post-yield deformation attributed to a given
material. The latter is determined for the SMs a posteriori, after obtaining
the stress-strain output for all the training batch, as the minimum yield
strain. Namely, as shown in Fig. 2 for 16 of the 42 SMs, we find that
$\epsilon_{f}=11.5\%$ is the minimum strain required for all the SMs to enter
yielding, that is, so that the stress-strain response is non-linear. The flow
stress $\sigma_{f}$ is the stress corresponding to this value of
$\epsilon_{f}$. It will be determined similarly for the real microstructures.

Figure 2: Stress-strain responses for 16 of the 42 SMs
($n\in\\{0.15,0.20,0.30,0.40\\}$, $h\in\\{2,3,4,5\\}$) until
$\epsilon_{f}=11.5\%$, showing that all the SMs reach post-yield for this
value. Changes from dotted to full lines correspond to an increase in porosity
$n$, and increases in grey intensity to an increase in heterogeneity $h$. The
rest of the SMs corroborates this assertion but are not shown here for
visibility purposes.

###  3.2 Mesh convergence

Before obtaining any meaningful results, it is important to study the mesh
convergence. The initial conditions of the simulations are indeed determined
by digitizing the input image, with respect to the chosen mesh dimensions. In
particular, smaller microstructural dimensions require a finer mesh.
Therefore, we expect to use finer meshes for SMs of smaller grain sizes. We
refer to B for the summary of the grain sizes of our training batch. To assess
the mesh convergence, we compare the values of the flow stress for the
different SMs, to find that, indeed, the mesh dimensions increase with the
grain size (see Fig. 3). We then pick the minimum mesh size required for mesh
convergence within a 2% error. In practice, using Fig. 3, we use the mesh
sizes shown in Tab. 1. Note that the grain sizes are measured in pixels where
$1px=l_{0}/N$, $l_{0}$ is the characteristic length of the problem taken as
the image size and $N$ the image dimension. When $l_{0}$ is determined in the
case of real microstructures, the length measurements will be given in $mm$
(or $\mu m$) in the following.

Figure 3: Evolution of the minimum mesh size (number of elements per side of the domain) required for mesh convergence with respect to the grain size. Mean grain size (px) | Mesh size (elts)
---|---
4-6 | $360\times 360$
6-7 | $350\times 360$
7-8 | $300\times 300$
8-13 | $280\times 280$
Table 1: Number of elements used for the meshes in our simulations for
different mean grain sizes.

###  3.3 Influence of the interfacial coefficient $\kappa$

While the activation energy $\lambda(\bm{\epsilon})$, coupled with the
mechanical energy, is the main drive for the phase change, the interfacial
coefficient $\kappa$ quantifies its diffusive character. When
$\kappa\gg\lambda$, that is, when the activation energy is negligible compared
to the interfacial energy, the process is purely diffusive, and therefore
nonphysical in the current context. Conversely, when $\kappa\rightarrow 0$,
the process falls back on the corresponding sharp interface problem, thereby
losing the regularization provided by the diffuse interfaces, characteristic
of phase-field modeling. We thus choose the value of $\kappa$ not too large to
avoid spurious interfacial diffusion, but not too low to keep diffuse
interfaces. The former would typically induce the closure of the smallest
pores, independently from the mechanical loading. To do so, we perform
simulations for different values of $\kappa$, while the microstructure and the
other parameters are fixed. To be conservative, we choose the SM with the
smallest mean grain size in our training batch, namely for $(n=0.4,h=5)$.
Indeed, spurious diffusion will occur first for the smallest microstructural
length scales. Note that more elaborated techniques exist to circumvent
spurious kinetic effects (see Tourret et al. (2021) and references therein).
We choose $\kappa=10^{-5}$ since for $\kappa\geq 10^{-4}$, spurious diffusion
occurs, whereas for $\kappa\leq 10^{-6}$, the results are similar to the case
where $\kappa=0$ (see Fig. 4). Specifically, it can be seen that for
$\kappa\geq 10^{-4}$ the smallest pores close under diffusion of the grains
phase, as opposed to the case where $\kappa\leq 10^{-5}$. The limit case
$\kappa\gg\lambda$ of pure diffusion is also shown in Fig. 4, where, after
only one timestep, the microstructure mixes into a homogeneous phase where
$\varphi\approx 0.6$. Thereupon, the mechanical response is that of a
homogeneous elastic solid, which homogenized elastic modulus can be estimated
around 6.

Figure 4: Mechanical response of the SM $(n=0.4,h=5)$ for different values of
$\kappa$, including the limit cases $\kappa=0$ (sharp interface) and
$\kappa\gg\lambda$ (pure diffusion). Note the spurious diffusion in the case
$\kappa=10^{-4}$, leading to the closure of the smallest pores, as opposed to
the case $\kappa=10^{-5}$.

###  3.4 Resemblance of the SMs with real microstructures

Before drawing realistic results from the SMs, let us also check how
representative they are of real microstructures. We find that the SMs resemble
real microstructures through the distributions of the three first
morphometers. We find that the grain size distributions of the SMs are skewed
towards smaller grains (i.e. right-tailed), as expected from thresholding the
initial symmetric Gaussian noise. The heterogeneity $h$ controls their
skewness (i.e. deviation from a normal distribution) and narrowness (see Fig.
5). The distributions of SMs with larger heterogeneity $h$ are less skewed but
narrower.

Figure 5: Evolution of the grain size distribution with the heterogeneity $h$
for a fixed porosity $n=0.2$ for a $200\times 200$ pixels resolution, showing
that higher values of $h$ lead to less skewed but narrower distributions.
Indeed, the skewness values, calculated as Pearson’s moment coefficient of
skewness, are found here to be 2.75 for $h=2$, 1.96 for $h=3$, 1.07 for $h=4$,
0.95 for $h=5$.

More precisely, the SMs exhibit lognormal distributions (see Fig. 6),
similarly to a large variety of porous media. This type of distribution can be
found in geomaterials (see Hwang and Powers (2003), Marks and Einav (2015) and
references therein), as well as engineered porous materials (see Liu et al.
(2017) e.g.). For direct comparison, we consider a Mt Simon sandstone,
presented in Kohanpur et al. (2020), and digitalized on
digitalrocksportal.org, of average porosity $27.1\%$. Each sides of the CT
scans have a resolution of 1200 pixels measuring each $2.8\mu m$. To imitate a
given slice of porosity $27.4\%$ (Fig. 6a), we generate a SM (Fig. 6e) of same
porosity and resolution, and of heterogeneity $h$ determined to match the mean
grain size. The distributions of the three first morphometers for the real and
synthetic microstructure are compared in Fig. 6 and Tab. 2, and checked to be
both well fitted by lognormal distributions.

\begin{overpic}[width=91.0631pt]{mtsimon890.png}
\put(-8.0,85.0){\color[rgb]{0,0,0} $a)$} \end{overpic}

\begin{overpic}[width=91.0631pt]{PSD_mtsimon890.png}
\put(1.0,65.0){\color[rgb]{0,0,0} $b)$} \end{overpic}

\begin{overpic}[width=91.0631pt]{perim_mtsimon890.png}
\put(1.0,65.0){\color[rgb]{0,0,0} $c)$} \end{overpic}

\begin{overpic}[width=91.0631pt]{GSD_mtsimon890.png}
\put(1.0,65.0){\color[rgb]{0,0,0} $d)$} \end{overpic}

\begin{overpic}[width=91.0631pt]{SM_poro027_blob39.png}
\put(-8.0,85.0){\color[rgb]{0,0,0} $e)$} \end{overpic}

\begin{overpic}[width=91.0631pt]{PSD_SM.png} \put(1.0,65.0){\color[rgb]{0,0,0}
$f)$} \end{overpic}

\begin{overpic}[width=91.0631pt]{perim_SM.png}
\put(1.0,65.0){\color[rgb]{0,0,0} $g)$} \end{overpic}

\begin{overpic}[width=91.0631pt]{GSD_SM.png} \put(1.0,65.0){\color[rgb]{0,0,0}
$h)$} \end{overpic}

Figure 6: a) CT scan of Mt Simon sandstone segmented into 733 grains, of porosity $M_{0}=27.4\%$ (mean pore size $41.51\mu m$), total perimeter $M_{1}=190.11mm$ (mean perimeter $1.02mm$) , mean grain size $M_{3}=79.80\mu m$ and Euler number $M_{4}=-351$, along with the associated lognormally-fitted distributions in b), c), d) respectively. e) SM generated for $n=27.4\%$ and $h=3.9$ with resolution $1200\times 1200$ pixels, matching the sandstone microstructure, segmented into 717 grains, of porosity $M_{0}=27.4\%$ (mean pore size $47.49\mu m$), total perimeter $M_{1}=167.17mm$ (mean perimeter $1.06mm$) , mean grain size $M_{3}=79.26\mu m$ and Euler number $M_{4}=-202$, along with the associated distributions in f), g), h) respectively. | Sandstone | Synthetic | Error (%)
---|---|---|---
Porosity $M_{0}$ (%) | 27.4 | 27.4 | 0
Perimeter $M_{1}$ (mm) | 190 | 167 | 12
Mean grain size $M_{2}$ ($\mu$m) | 79.8 | 79.3 | 0.68
Mean pore size ($\mu$m) | 41.5 | 47.5 | 14
Mean perimeter (mm) | 1.02 | 1.06 | 3.9
Table 2: Comparison of the morphometers of the Mt Simon sandstone with the
synthetic analogue, along with the mean values of their lognormal
distributions.

###  3.5 Determinism of the SMs

Since the generation of the SMs starts from a random noise, it is not a priori
deterministic. Indeed, varying the random seeds for a given set of input
parameters $n$ and $h$ yields different microstructures. To quantify the
influence of this discrepancy on the mechanical response, we compare the
numerical results for $n=0.3$, $h=3$ as an example, with 10 different seeds
(see D). Namely, we calculate the set of flow stresses $\sigma_{f}$ (for
$\epsilon_{f}=11.5\%$) for the different seeds. We find that the standard
deviation of this set divided by its mean is around $4\%$, which we consider
as an approximate measure of the potential error in the following predictive
fitting.

###  3.6 Morphometers evolution

Numerical results are shown in Fig. 7 for the SM $n=0.2$ and $h=5$, from the
initial SM input, which digitization provides the initial conditions for the
phase field $\varphi$, to the (post-yield) flow stress state. The associated
distributions of the volumetric strain and of the phase-field rate illustrate
the causality between mechanical loading and damage. Upon tensile loading, the
pores open up, leading to damage events responsible for failure, namely the
coalescence of neighboring pores, well known possible microstructural
manifestation of strain hardening (Pardoen and Hutchinson (2000)). This can
represent debonding or microcracks. In this particular example, the main
damage event (maximum of $\lvert\dot{\varphi}\rvert$) responsible for the flow
stress state occurs in the center (see black oval in Fig. 7).

Figure 7: a) SM generated for $n=0.2$ and $h=5$. b) Digitization of the SM
used as input for $\varphi$. c) Output for $\varphi$ after tensile loading
post-yield at a certain flow stress for a vertical strain of $11.5\%$. d)
Volumetric strain distribution for the same flow stress, concentrated in a
central zone (black oval). e) Distribution of $\dot{\varphi}$ for the same
flow stress, negative since the pore phase $\varphi=0$ is produced, with
maximum absolute value also in the central zone, coinciding with the
concentration of mechanical energy.

As a macroscopic descriptor of the microscopic damage, we define the relative
degradation $D_{r}$, similarly to breakage mechanics (Einav (2007)), so that
$D_{r}=0$ corresponds to an undamaged material and $D_{r}=100\%$ to a final
state of damage. The latter is arbitrarily defined as the flow state (see
Section 3.1). Thus, $D_{r}$ is a relative measure of degradation, and not a
state variable, as it depends on the initial and final state considered, and
is calculated as the following upscaling function:

$D_{r}=\frac{\int_{\bar{\Omega}}\varphi(\epsilon)\text{d}\epsilon-\int_{\bar{\Omega}}\varphi(\epsilon_{0})\text{d}\epsilon}{\int_{\bar{\Omega}}\varphi(\epsilon_{f})\text{d}\epsilon-\int_{\bar{\Omega}}\varphi(\epsilon_{0})\text{d}\epsilon},$
(13)

where ${\bar{\Omega}}$ denotes the total spatial domain, $\epsilon_{0}=0$,
$\epsilon_{f}$ is the strain when $\sigma_{f}$ is reached, so that
$\int_{\bar{\Omega}}\varphi(\epsilon)\text{d}\epsilon$ computes the domain’s
portion occupied by the grains, since $\varphi=1$ characterizes the grains and
$\varphi=0$ the pores. The evolution of the relative degradation of the SM of
Fig. 7 with respect to the strain is shown in Fig. 8, along with the
corresponding stress-strain curve and the evolution of the 4 morphometers;
although only 3 morphometers are required in 2D to fully constrain the
strength law, we find it informative to include all of them here. As expected,
the tensile loading opens up the pores, which increases the porosity
$n=M_{0}$. Moreover, as the neighboring pores merge with each other towards
failure (see Fig. 7), the number of pores decreases (by around 9%) and
therefore, the Euler characteristic $\chi=M_{3}$ decreases in absolute value.
As a consequence, the total perimeter $M_{1}$ decreases and the mean grain
size $M_{2}$ increases, although both by less than 5%, as opposed to $M_{0}$
and $M_{3}$ varying by around 20%. The observed morphometers evolution is
consistent with the evolution of strength found later in Section 4.1 with
respect to the initial morphometers values. Indeed, the material weakens upon
loading, and therefore can be seen as a strong material initially and as a
weaker one at the end of the loading. Consistently, we find that the strength
decreases with the porosity $M_{0}$ and the mean grain size $M_{2}$, but
increases with the total perimeter $M_{1}$ and the absolute value of the Euler
characteristic $M_{3}$.

\begin{overpic}[width=173.44534pt]{stress_degradation_strain.pdf}
\put(-2.0,58.0){\color[rgb]{0,0,0} $a)$}
\end{overpic}\begin{overpic}[width=173.44534pt]{morphometers_evolution.pdf}
\put(-2.0,58.0){\color[rgb]{0,0,0} $b)$} \end{overpic}

Figure 8: a) Stress-strain curve (in blue) corresponding to the evolution of
the SM presented in Fig. 7, along with the degradation-strain curve (in
orange). b) Evolution of the morphometers with the degradation, whose values
are normalized with their final values. Note that the values of the Euler
characteristic $M_{3}$ is negative before normalization.

##  4 Results

###  4.1 Scaling of flow stress with morphometers in SMs

We now gather the results of all the simulations on the 42 SMs (see B), to
find the best-fit function $f(M_{i})$ relating the flow stress to the
morphometers, via the least squares method. In the present 2D case, we only
need 3 morphometers, chosen as $M_{0}$, $M_{1}$ and $M_{2}$ for instance.
Following our initial assumption of a multiplicative dependence of the
strength on the morphometers, as a generalization of Knudsen’s theory (Knudsen
(1959)), we find that the best fit is exponential, with an adjusted $r^{2}$
coefficient of 0.9995 and an average error of 2.49%, as compared with other
potential candidates (see Tab. 3). In the present 2D case, the strength law
Eq. (1) thus reads

$\sigma_{f}=\sigma_{f}^{*}e^{-\frac{M_{0}-M_{0}^{*}}{M_{0}^{*}}+\frac{M_{1}-M_{1}^{*}}{M_{1}^{*}}-\frac{M_{2}-M_{2}^{*}}{M_{2}^{*}}}.$
(14)

Notably, the fitting coefficient for the porosity $1/M_{0}^{*}=6.6$ is found
to be in the range of values found experimentally by Knudsen (1959). The
prediction is found to be as accurate when replacing $M_{2}$ by $M_{3}$. We
check that a Hall-Petch type of fitting with respect to the mean grain size
$M_{2}$ is also satisfying, while keeping the exponential dependence for
$M_{0}$ and $M_{1}$. We find that the flow stress $\sigma_{f}$ increases with
the total perimeter $M_{1}$ and with the absolute value of the Euler
characteristic $M_{3}$. We also check that it decreases with the porosity
$M_{0}$ and with the grain size $M_{2}$, as already well known. For a given
random seed, the prediction is thus found to be erroneous by less than 3% on
average. Then, since choosing a seed amounts to an error of approximatively
4.33% (see Section 3.5), we infer that the overall predictive error is
$\sqrt{2.49^{2}+4.33^{2}}\approx 5\%$. To use the latter formula, we have
assumed that the error propagation is multiplicative, since the choice of the
seed and of the morphometry is concomitant. This error estimation relies,
however, on assuming that the stochastic error of 4.33% is valid for all the
SMs, which is not here verified.

Best fit |  | Exponential
---
($M_{0}$, $M_{1}$, $M_{2}$)
| Exponential
---
($M_{0}$, $M_{1}$, $M_{3}$)
Power | Hall-Petch | Linear
Adjusted $r^{2}$ | 0.9995 | 0.9995 | 0.9974 | 0.9990 | 0.9780
Average error (%) | 2.49 | 2.47 | 5.96 | 3.70 | 15.03
Fitting function | $\alpha_{1}e^{\alpha_{2}M_{0}+\alpha_{3}M_{1}+\alpha_{4}M_{2}}$ | $\alpha_{1}e^{\alpha_{2}M_{0}+\alpha_{3}M_{1}+\alpha_{4}M_{3}}$ | $\alpha_{1}M_{0}^{\alpha_{2}}M_{1}^{\alpha_{3}}M_{2}^{\alpha_{4}}$ | $\alpha_{1}e^{\alpha_{2}M_{0}+\alpha_{3}M_{1}}M_{2}^{-0.5}$ | $\alpha_{1}M_{0}+\alpha_{2}M_{1}+\alpha_{3}M_{2}$
Table 3: Different fitting models with associated adjusted $r^{2}$ coefficient
and average error, showing that the best fitting function is exponential.

The dependence of the strength on the porosity and grain size was already
established to be exponential, as recovered here, or similarly, a power law
(see Knudsen (1959) and references therein). However, there are few results,
if any, relating the strength to the two other morphometers, the surface area
$M_{2}$ and the Euler number $M_{4}$, which could have corroborated our
results.

###  4.2 Dependence of the strength on the confinement

While the focus of this work concerns the morphometry, the strength can depend
on other factors such as the confining pressure $P$, which we shall briefly
discuss here. Such dependence can be included in the reference stress
$\sigma_{f}^{*}$. We determine it by performing simulations of axisymmetric
biaxial compressions for different values of $P$, for a given microstructure
with values of porosity and heterogeneity in the middle range of our training
batch ($n=0.3$, $h=3$), in plane strain conditions. To constrain the range of
confining pressures, we estimate the preconsolidation pressure $P_{c}$ by
performing an isotropic compression (see E). By reading the yield stress value
in Fig. 16, we find an approximate value of $P_{c}=2$. We then plot the
obtained flow stresses with respect to the confining pressures (see Fig. 9),
where the best fit function is a parabola. This recovers the parabolic failure
criterion, recently suggested and experimentally verified (see Yuan et al.
(2020), Wang et al. (2019), Singh and Rao (2005)). This also captures the
transition from a linear behavior at low confining pressures to a pressure-
insensitive behavior for high confinement (Yuan et al. (2020)).

Figure 9: The points represent numerical results for the flow stress versus
the confining pressure in axisymmetric biaxial compression. The curve
corresponds to the fit $-0.0513p^{2}+0.181p+3.14$, for which the adjusted
$r^{2}$ is exactly 1.

###  4.3 Strength prediction from real microstructures

To assess the validity of the morphometric strength law Eq. (14), we predict
the strength of real materials from CT scans, and compare it with phase-field
simulations as previously. For that, we employ the following procedure. 1) We
choose a set of $m$ CT scans $C^{1},\dots,C^{m}$ of a given material and pick
one of them as the reference $C^{*}$. 2) We binarize the $C^{i}$ and calculate
the morphometers for each $C^{i}$ as done previously for the SMs. In
particular, this determines the reference values for the morphometers
$M_{i}^{*}$. 3) We run one simulation for the reference CT scan $C^{*}$ to
determine the remaining reference value $\sigma_{f}^{*}$. 4) We deduce the
$\sigma_{f}^{i}$ for the remaining CT scans using the strength law. 5) We run
the simulations for the remaining CT scans for comparison.

####  4.3.1 Geomaterials

Starting with geomaterials, we use CT scans222digitalrocksportal.org of Mt
Simon and Berea sandstones, two benchmark rocks. The Mt Simon sandstone was
introduced in Section 3.4. The Berea sandstone, of average porosity $17.9\%$,
is presented in Lucas-Oliveira et al. (2020), with a digitization resolution
of $2.25\mu m$ per pixel. The resolution of the CT scans of both sandstones is
reduced to $300\times 300$ pixels to reduce the mesh size and have reasonable
computation times, in particular yielding erroneously high porosity values in
Tab. 4. This does not hinder the present qualitative study since the different
parts of the Berea sandstone remain comparable with each others. Using the
same model Eq. (12) and parameters values, we follow the procedure described
above. The results are gathered in Fig. 12 and detailed in Fig. 10 and Tab. 4
for the example of a reference CT scan and a predicted flow stress for a
weaker part with higher porosity. The predictions show a good agreement with
an average error of 2.04% for the Mt Simon sandstone and 3.88% for the Berea
sandstone. We must note that the choice of the reference CT scan $C^{*}$ was
made judiciously to minimize this error. For instance, for the Mt Simon
sandstone, while most choices of $C^{*}$ yielded an average error around 2 or
3%, one in particular yielded an average error of 25%. This seems to pertain
to the usual problem in modeling geomaterials, that is, finding a
representative elementary domain. Indeed, it is sensible that in order to
deduce the properties of a material’s domain, one must use a representative
domain as a reference.

\begin{overpic}[width=433.62pt]{sandstones.pdf}
\put(0.0,52.0){\color[rgb]{0,0,0} $a)$} \put(31.0,52.0){\color[rgb]{0,0,0}
$b)$} \put(0.0,24.0){\color[rgb]{0,0,0} $c)$}
\put(31.0,24.0){\color[rgb]{0,0,0} $d)$} \put(60.0,40.0){\color[rgb]{0,0,0}
$e)$} \end{overpic} Figure 10: a) CT scans of a Mt Simon sandstone (layer 108)
and b) another layer of the same sandstone (890) with 8% increase of porosity
with similar other morphometers, hence weaker. c) a) CT scans of a Berea
sandstone (layer 000) and b) another layer of the same sandstone (940) with 7%
increase of porosity with similar other morphometers, hence weaker. e)
Mechanical response of the 4 sandstones microstructures.

####  4.3.2 Bones

To showcase the general predictive power of the morphometric strength law Eq.
(14), we perform the same predictive procedure for bones, which are porous
materials exhibiting a wide range of porosity and microstructural architecture
different from that of geomaterials. In particular, the example of bones
enables us to address anisotropic microstructures. Typical values of
morphometers for geomaterials and bones are shown in Tab. 4 for comparison. As
far as the microstructure is concerned, there are two main kinds of bones, the
cortical bone and the trabecular bone (see Cowin (2002) for a general review
of bone mechanics). The former occupies the perimeter of the bone and has a
low porosity ranging typically from 5 to 25%. The latter occupies the inside
of the bone and has a scaffold microstructure, with a much higher porosity
ranging from 50 to 90%. The main morphometer varying from a bone to another is
the porosity. Indeed, the creation of porosity, because of age or disease, is
the main sign of bone degradation. Thus, we will use our strength law to
predict the strength of bones with varying porosities, given a reference bone,
typically a healthy bone with a normal porosity.

For the trabecular bone, we used the CT scan in Fig. 1 of Neumann et al.
(2018), providing a transverse section. We choose the healthy microstructure
as the reference, to predict the strength of the arthritic one (see Fig. 11).
As for the cortical bone, failing to find satisfying input images, we model it
via SMs with anisotropy to resemble real cortical microstructures in their
longitudinal section (see Fig. 9 in Cooper et al. (2004) and Fig. 1 in Granke
et al. (2011)). For that, the heterogeneity $h$ is chosen as an array $(1,4)$,
meaning that the pores are elongated in the vertical direction four times more
than in the horizontal one (see Fig. 11). We then generate different SMs with
varying porosities (3, 5, 7, 10, 15 and 20%) and fixed anisotropy.

We use the same model as previously, but with different elastic moduli for the
skeleton phase. The same problem of determining the moduli of the skeleton
phase holds as for geomaterials. Indeed, elastic moduli are usually measured
for the mixture of the skeleton and pores phases. As explained in Section
(2.3), we have assimilated the skeleton phase of geomaterials to a rock with
very low porosity. Similarly, we assume the solid (mineral) phase of bones to
be analogous to a bone of very low porosity, namely a cortical bone, for which
the Lamé parameters can be found to be around $10GPa$ and $7GPa$ respectively
(see Lai et al. (2015)). We choose the flow stress for the bones as
$\sigma_{f}=25\%$. For the trabecular bone, since we use transverse sections,
the loading is compressive. As for the cortical bone, modeled longitudinally,
we choose a tensile loading. In agreement with Morgan et al. (2018), we find
that for a given post-yield vertical strain, the trabecular microstructure, in
compression, has a softening response, whereas the cortical SM, in tension, is
limited to a hardening response (see Fig. 11). For the synthetic cortical
bones, we show for instance the results for $n=5\%$ as a healthy bone and
$n=10\%$ as the arthritic bone to be predicted, corresponding to real porosity
values (see Fig. 6 in Pauly et al. (2015) for the lateral tibia). Finally, we
show a relatively accurate prediction of the strength from our morphometric
law (see Fig. 12), with an error of 4.2% for the trabecular bone (detailed in
Tab. 4), and an average error of 1.3% for the synthetic trabecular bones. We
note that for the latter, the accurate prediction required using two different
reference microstructures, one for $n<10\%$ and one for $n\geq 10\%$,
indicating that the material may have two distinct mechanical behaviors for
those two ranges of porosities. We note also that we have used $M_{3}$ instead
of $M_{2}$ for convenience, as determining the grain size distributions turns
out less accurate for microstructures with notably low or high porosities and
may not bear any physical meaning.

\begin{overpic}[width=433.62pt]{bones.pdf} \put(-1.0,52.0){\color[rgb]{0,0,0} $a)$} \put(30.0,52.0){\color[rgb]{0,0,0} $b)$} \put(-1.0,24.0){\color[rgb]{0,0,0} $c)$} \put(30.0,24.0){\color[rgb]{0,0,0} $d)$} \put(59.0,40.0){\color[rgb]{0,0,0} $e)$} \end{overpic} Figure 11: a) CT scans of a healthy and b) arthritic trabecular bone’s microstructure (from Neumann et al. (2018)). c) Synthetic microstructures of a healthy and d) arthritic cortical bone. e) Mechanical response of the 4 bones. | $M_{0}(\%)$ |  $M_{1}$ (mm) |  $M_{2}$ ($\mu$m) | $M_{3}$ |  $\sigma_{f}$ (numerics) |  $\sigma_{f}$ (predicted) | Error (%)
---|---|---|---|---|---|---|---
Mt Simon sandstone | 20.6 | 142 | 118 | / | 3.71 | * | /
Weaker MS sandstone | 22.2 | 145 | 121 | / | 3.39 | 3.59 | 6.09
Berea sandstone | 37.7~ | 253 | 7.96 | / | 1.33 | * | /
Weaker B. sandstone | 40.2~ | 253 | 8.35 | / | 1.23 | 1.26 | 2.91
Healthy trabecular bone | 29.5 | 451 | / | -423 | 1.51 | * | /
Arthritic trabecular bone | 41.1 | 552 | / | -265 | 0.986 | 0.944 | 4.2
Healthy cortical bone | 5 | 37.4 | / | -67 | 5.23 | * | /
Arthritic cortical bone | 10 | 73.1 | / | -91 | 4.41 | 4.36 | 1.08
Table 4: Details of the prediction of the strength of the Mt Simon and Berea
sandstones and of the trabecular and cortical bones, for one reference and one
prediction. *reference CT scans. ~Porosities larger than real values from
lowering initial resolution.

In all, the average errors between predicted and simulated flow stresses are
all below 5%, the expected maximum prediction error calculated in Section 4.1.

Figure 12: Predicted flow stress with the morphometric strength law Eq. (14)
versus the flow stress obtained from phase-field simulations, for SMs,
geomaterials (Mt Simon and Berea sandstones) and biomaterials (trabecular
bones and synthetic cortical bones).

##  5 Conclusion

We have suggested from qualitative 2D numerical simulations on synthetic and
real microstructures the possibility to fully account for the morphometry of
porous media in their macroscopic properties, such as the strength. To do so,
we have used a damage phase-field modeling of the microstructure, capturing
the exact microstructural geometry and encapsulating causes of damage such as
debonding, dissolution and microcracks. The necessary and sufficient
microstructural information is then upscaled in the form of morphometers
defined from Minkowski functionals, as per Hadwiger’s theorem. In the context
of porous media, the four morphometers are chosen as the porosity, the total
surface area of the skeleton, the mean grain size and the Euler
characteristic; only three of them are required in 2D. We have inferred from a
wide range of synthetic microstructures, with various porosities and
heterogeneities, that the strength is best described by an exponential
function of the morphometers. Those generated microstructures were checked to
resemble real microstructures via the lognormal distribution of their
morphometers. The exponential strength law was also checked to be accurate for
two subclasses of porous media, geomaterials (sandstones) and biomaterials
(bones), from digitalized CT scans. The strength prediction only requires the
initial values of the morphometers for a given material’s sample, along with
the reference strength and morphometers values from another microstructure
made of the same material.

The overarching goal is to predict the behavior of porous media with the
minimal amount of data possible and in particular, minimizing the use of
destructive tests in favor of remote sensing. This is of interest in
geosciences where samples may be inaccessible and in biosciences where
materials are best studied in vivo. The following steps in this endeavor will
include 3D modeling and most importantly, experimental corroboration of the
morphometric strength law. Then, the influence of environmental factors, like
temperature and acidity, may be incorporated in a similar fashion. Building on
aforementioned works, machine learning may also be used to ensure the realism
of the synthetic microstructures. Ultimately, classical constitutive laws,
such as viscoplasticity, may be re-explored to explicitly include the
morphometric dependence. A major impediment to such program, albeit common to
any modeling of heterogeneous microstructures, is that the reference sample
must be representative of the material. This holistic description opens new
avenues for modeling porous media in general, including geomaterials,
biomaterials and engineered materials, at least as far as the microstructural
geometry is concerned. It is indeed of paramount interest to draw bridges
between the progress in different disciplines concerned with the same
problematics.

##  6 Acknowledgement

The authors Alexandre Guével and Manolis Veveakis gratefully acknowledge the
support of US DOE grant DE-NE0008746.

##  Appendix A Exponential relationship between strength and porosity for
bones

Figure 13: Exponential fit of the yield stress evolution with porosity for
cortical bones, adapted from the experimental results of Wachter et al.
(2002).

##  Appendix B Summary table

Figure 14: Summary table gathering the morphometers of the SMs, with
corresponding flow stresses obtained numerically with phase-field simulations,
versus the flow stresses predicted with the morphometric strength law, along
with the associated prediction error. The predicted flow stress is found to be
best predicted with the exponential function
$\sigma_{f}=12.37342^{-6.64471e+0.938539p-0.951469g}$, obtained from
Mathematica’s function NonlinearModelFit

##  Appendix C Calculations of the morphometers

The morphometers are calculated via the Python libraries PoreSpy et skimage,
as described in the following table. The segmentation of the SMs is performed
via the SNOW algorithm introduced in Gostick (2017), using the function
filters.snow_partitioning in PoreSpy.

Morphometer | Library | Function
---|---|---
Porosity $M_{0}$ | PoreSpy | metrics.porosity
Perimeter $M_{1}$ | skimage | measure.perimeter
Grain size $M_{2}$ | PoreSpy |  pore.diameter*
Euler number $M_{3}$ | skimage | measure.euler_number
Table 5: * The function pore.diameter was initially defined (and named) to
calculate the sizes of the pores obtained from segmentation. However, it can
be applied to calculate the grain size when applied on the segmentation of the
grains.

##  Appendix D Determinism of the SMs

Figure 15: Mechanical response for different random seeds for the SM $n=0.3$,
$h=3$.

##  Appendix E Preconsolidation

Figure 16: Stress-strain curve of the isotropic compression simulation
determining the preconsolidation stress, taken as the yield stress of this
curve and found to be around 2.

## References

## References

  * Allen and Cahn (1979) Allen, S.M., Cahn, J.W., 1979\.  A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening.  Acta Metallurgica 27, 1085–1095.  doi:10.1016/0001-6160(79)90196-2.
  * Ammar et al. (2009) Ammar, K., Appolaire, B., Cailletaud, G., Forest, S., 2009\.  Combining phase field approach and homogenization methods for modelling phase transformation in elastoplastic media.  Revue européenne de mécanique numérique 18, 485–523.  doi:10.3166/ejcm.18.485-523.
  * Armstrong et al. (2019) Armstrong, R.T., McClure, J.E., Robins, V., Liu, Z., Arns, C.H., Schlüter, S., Berg, S., 2019.  Porous Media Characterization Using Minkowski Functionals: Theories, Applications and Future Directions.  Transport in Porous Media 130, 305–335.  doi:10.1007/s11242-018-1201-4.
  * Augat and Schorlemmer (2006) Augat, P., Schorlemmer, S., 2006\.  The role of cortical bone and its microstructure in bone strength.  Age and Ageing 35, 27–31.  doi:10.1093/ageing/afl081.
  * Benkrid et al. (2000) Benkrid, K., Crookes, D., Benkrid, A., 2000.  Design and fpga implementation of a perimeter estimator.  Proceedings of the Irish Machine Vision and Image Processing Conference , 51–57.
  * Bhattacharyya et al. (2019) Bhattacharyya, S., Sahara, R., Ohno, K., 2019.  A first-principles phase field method for quantitatively predicting multi-composition phase separation without thermodynamic empirical parameter.  Nature Communications 10, 1–2.  doi:10.1038/s41467-019-11248-z.
  * Biot (1941) Biot, M.A., 1941.  General theory of three-dimensional consolidation.  Journal of Applied Physics 12, 155–164.
  * Buscarnera and Einav (2021) Buscarnera, G., Einav, I., 2021\.  The mechanics of brittle granular materials with coevolving grain size and shape.  Proceedings the Royal Society A 477\.  doi:20201005.
  * Chun et al. (2020) Chun, S., Roy, S., Nguyen, Y.T., Choi, J.B., Udaykumar, H.S., Baek, S.S., 2020.  Deep learning for synthetic microstructure generation in a materials-by-design framework for heterogeneous energetic materials.  Scientific Reports 10, 1–15.  doi:10.1038/s41598-020-70149-0.
  * Cooper et al. (2004) Cooper, D.M., Matyas, J.R., Katzenberg, M.A., Hallgrimsson, B., 2004\.  Comparison of Microcomputed Tomographic and Microradiographic Measurements of Cortical Bone Porosity.  Calcified Tissue International 74, 437–447.  doi:10.1007/s00223-003-0071-z.
  * Coussy (2004) Coussy, O., 2004.  Poromechanics.  Wiley.  doi:10.1002/0470092718.
  * Cowin (2002) Cowin, S.C. (Ed.), 2002.  Bone Mechanics Handbook.  CRC Press.  doi:10.1115/1.1579463.
  * Croizé et al. (2013) Croizé, D., Renard, F., Gratier, J.P., 2013.  Compaction and Porosity Reduction in Carbonates: A Review of Observations, Theory, and Experiments, in: Advances in Geophysics. Elsevier Inc.. volume 54, pp. 181–238.  doi:10.1016/B978-0-12-380940-7.00003-2.
  * Dormieux et al. (2006) Dormieux, L., Kondo, D., Ulm, F.J., 2006.  Microporomechanics.  John Wiley and Sons Inc.
  * Einav (2007) Einav, I., 2007.  Breakage mechanics-Part I: Theory.  Journal of the Mechanics and Physics of Solids 55, 1274–1297.  doi:10.1016/j.jmps.2006.11.003.
  * van den Ende et al. (2019) van den Ende, M.P., Niemeijer, A.R., Spiers, C.J., 2019.  Influence of Grain Boundary Structural Evolution on Pressure Solution Creep Rates.  Journal of Geophysical Research: Solid Earth , 1–30doi:10.1029/2019JB017500.
  * Fragogeorgi et al. (2019) Fragogeorgi, E.A., Rouchota, M., Georgiou, M., Velez, M., Bouziotis, P., Loudos, G., 2019\.  In vivo imaging techniques for bone tissue engineering.  Journal of Tissue Engineering 10\.  doi:10.1177/2041731419854586.
  * Fried and Gurtin (1993) Fried, E., Gurtin, M.E., 1993\.  Continuum theory of thermally induced phase transitions based on an order parameter.  Physica D: Nonlinear Phenomena 68, 326–343.  doi:10.1016/0167-2789(93)90128-N.
  * Gonzalez and Woods (2008) Gonzalez, R.C., Woods, R.E., 2008\.  Digital image processing.  Prentice Hall.
  * Gostick et al. (2019) Gostick, J., Khan, Z., Tranter, T., Kok, M., Agnaou, M., Sadeghi, M., Jervis, R., 2019.  PoreSpy: A Python Toolkit for Quantitative Analysis of Porous Media Images.  Journal of Open Source Software 4, 1296.  doi:10.21105/joss.01296.
  * Gostick (2017) Gostick, J.T., 2017.  Versatile and efficient pore network extraction method using marker-based watershed segmentation.  Physical Review E 96, 1–15.  doi:10.1103/PhysRevE.96.023307.
  * Granke et al. (2011) Granke, M., Grimal, Q., Saïed, A., Nauleau, P., Peyrin, F., Laugier, P., 2011\.  Change in porosity is the major determinant of the variation of cortical bone elasticity at the millimeter scale in aged women.  Bone 49, 1020–1026.  doi:10.1016/j.bone.2011.08.002.
  * Guével et al. (2020) Guével, A., Rattez, H., Veveakis, E., 2020.  Viscous phase-field modeling for chemo-mechanical microstructural evolution: application to geomaterials and pressure solution.  International Journal of Solids and Structures 207, 230–249.  doi:10.1016/j.ijsolstr.2020.09.026.
  * Gurtin (1996) Gurtin, M.E., 1996.  Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance.  Physica D: Nonlinear Phenomena 92, 178–192.  doi:10.1155/S1110757X03204083.
  * Hadwiger (1951) Hadwiger, H., 1951.  The deformation and ageing of mild steel: III Discussion of results.  Proceedings of the Physical Society. Section B 64, 747–753.
  * Hall (1951) Hall, E.O., 1951.  The deformation and ageing of mild steel: III Discussion of results.  Proceedings of the Physical Society. Section B 64, 747–753.  doi:10.1088/0370-1301/64/9/303.
  * Huyghe et al. (2002) Huyghe, J.M., Loon, R.V., Kemenade, P.M.V., Smit, T.H., Perfusion, B., 2002.  We all are porous media.  Poromechanics II , 17–28.
  * Hwang and Powers (2003) Hwang, S.I., Powers, S.E., 2003\.  Lognormal Distribution Model for Estimating Soil Water Retention Curves for Sandy Soils.  Soil Science 168, 156–166.  doi:10.1097/01.ss.0000058888.60072.e3.
  * Kawamoto et al. (2018) Kawamoto, R., Andò, E., Viggiani, G., Andrade, J.E., 2018\.  All you need is shape: Predicting shear banding in sand with LS-DEM.  Journal of the Mechanics and Physics of Solids 111, 375–392.  doi:10.1016/j.jmps.2017.10.003.
  * Klain (1995) Klain, D.A., 1995.  A Short Proof Of Hadwiger’S Characterization Theorem.  Mathematika 42, 329–339.  doi:10.1112/S0025579300014625.
  * Knudsen (1959) Knudsen, F.P., 1959.  Dependence of Mechanical Strength of Brittle Polycrystalline Specimens on Porosity and Grain Size.  Journal of the American Ceramic Society 42, 376–387.  doi:10.1111/j.1151-2916.1959.tb13596.x.
  * Kohanpur et al. (2020) Kohanpur, A.H., Rahromostaqim, M., Valocchi, A.J., Sahimi, M., 2020\.  Two-phase flow of CO2-brine in a heterogeneous sandstone: Characterization of the rock and comparison of the lattice-Boltzmann, pore-network, and direct numerical simulation methods.  Advances in Water Resources 135.  doi:10.1016/j.advwatres.2019.103469.
  * Kuhn and Müller (2010) Kuhn, C., Müller, R., 2010\.  A continuum phase field model for fracture.  Engineering Fracture Mechanics 77, 3625–3634.  doi:10.1016/j.engfracmech.2010.08.009.
  * Lai et al. (2015) Lai, Y.S., Chen, W.C., Huang, C.H., Cheng, C.K., Chan, K.K., Chang, T.K., 2015\.  The effect of graft strength on knee laxity and graft in-situ forces after posterior cruciate ligament reconstruction.  PLoS ONE 10.  doi:10.1371/journal.pone.0127293.
  * Landau (1937) Landau, L.D., 1937.  On the theory of phase transitions.  Zh. Eks. Teor. Fiz. 7, 19–32.  doi:10.1038/138840a0.
  * Lesueur et al. (2020) Lesueur, M., Poulet, T., Veveakis, M., 2020.  Permeability Hysteresis From Microchannels Opening During Dissolution/Reprecipitation Cycle.  Geophysical Research Letters 47, 1–8.  doi:10.1029/2020GL088674.
  * Li et al. (2016) Li, Y., Bushby, A.J., Dunstan, D.J., 2016.  The Hall-Petch effect as a manifestation of the general size effect.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472.  doi:10.1098/rspa.2015.0890.
  * Liu et al. (2017) Liu, D., Šavija, B., Smith, G.E., Flewitt, P.E., Lowe, T., Schlangen, E., 2017\.  Towards understanding the influence of porosity on mechanical and fracture behaviour of quasi-brittle materials: experiments and modelling.  International Journal of Fracture 205, 57–72.  doi:10.1007/s10704-017-0181-7.
  * Lucas-Oliveira et al. (2020) Lucas-Oliveira, E., Araujo-Ferreira, A.G., Trevizan, W.A., Coutinho dos Santos, B.C., Bonagamba, T.J., 2020.  Sandstone surface relaxivity determined by NMR T2 distribution and digital rock simulation for permeability evaluation.  Journal of Petroleum Science and Engineering 193, 107400.  doi:10.1016/j.petrol.2020.107400.
  * Marks and Einav (2015) Marks, B., Einav, I., 2015.  A mixture of crushing and segregation: The complexity of grainsize in natural granular flows.  Geophysical Research Letters 42, 274–281.  doi:10.1002/2014GL062470.
  * Morgan et al. (2018) Morgan, E.F., Unnikrisnan, G.U., Hussein, A.I., 2018.  Bone Mechanical Properties in Healthy and Diseased States.  Annual Review of Biomedical Engineering 20, 119–143.  doi:10.1146/annurev-bioeng-062117-121139.
  * Mosser et al. (2017) Mosser, L., Dubrule, O., Blunt, M.J., 2017.  Reconstruction of three-dimensional porous media using generative adversarial neural networks.  Physical Review E 96.  doi:10.1103/PhysRevE.96.043309, arXiv:1704.03225.
  * Neumann et al. (2018) Neumann, A., Haschka, J., Kleyer, A., Schuster, L., Englbrecht, M., Berlin, A., Figueiredo, C.P., Simon, D., Muschitz, C., Kocijan, R., Resch, H., Rech, J., Schett, G., 2018.  Cortical bone loss is an early feature of nonradiographic axial spondyloarthritis.  Arthritis Research and Therapy 20, 1–11.  doi:10.1186/s13075-018-1620-1.
  * Niemeijer et al. (2009) Niemeijer, A., Elsworth, D., Marone, C., 2009.  Significant effect of grain size distribution on compaction rates in granular aggregates.  Earth and Planetary Science Letters 284, 386–391.  doi:10.1016/j.epsl.2009.04.041.
  * Orowan (1949) Orowan, E., 1949.  Fracture and strength of solids.  Reports on Progress in Physics 12, 185–232.
  * O’Sullivan (2011) O’Sullivan, C., 2011.  Particulate Discrete Element Modelling: a Geomechanics Perspective.  Spon Press.
  * Pardoen and Hutchinson (2000) Pardoen, T., Hutchinson, J.W., 2000\.  Extended model for void growth and coalescence.  Journal of the Mechanics and Physics of Solids 48, 2467–2512.  doi:10.1016/S0022-5096(00)00019-3.
  * Pauly et al. (2015) Pauly, H.M., Larson, B.E., Coatney, G.A., Button, K.D., DeCamp, C.E., Fajardo, R.S., Haut, R.C., Haut Donahue, T.L., 2015\.  Assessment of cortical and trabecular bone changes in two models of post-traumatic osteoarthritis.  Journal of orthopaedic research 33, 1835–1845.  doi:10.1002/jor.22975.
  * Permann et al. (2020) Permann, C.J., Gaston, D.R., Andrš, D., Carlsen, R.W., Kong, F., Lindsay, A.D., Miller, J.M., Peterson, J.W., Slaughter, A.E., Stogner, R.H., Martineau, R.C., 2020.  MOOSE: Enabling massively parallel multiphysics simulation.  SoftwareX 11, 100430\.  doi:10.1016/j.softx.2020.100430.
  * Petch (1953) Petch, N.J., 1953.  The Cleavage Strength of Polycrystals.  Journal of the Iron and Steel Institute 174, 25–28.
  * Provatas and Elder (2010) Provatas, N., Elder, K., 2010\.  Phase-Field Methods in Material Science and Engineering.  Wiley-VCH , 312doi:10.1002/9783527631520.
  * Rattez et al. (2018a) Rattez, H., Stefanou, I., Sulem, J., 2018a.  The importance of Thermo-Hydro-Mechanical couplings and microstructure to strain localization in 3D continua with application to seismic faults. Part I: Theory and linear stability analysis.  Journal of the Mechanics and Physics of Solids 115, 54–76.  doi:10.1016/j.jmps.2018.03.004.
  * Rattez et al. (2018b) Rattez, H., Stefanou, I., Sulem, J., Veveakis, M., Poulet, T., 2018b.  The importance of Thermo-Hydro-Mechanical couplings and microstructure to strain localization in 3D continua with application to seismic faults. Part II: Numerical implementation and post-bifurcation analysis.  Journal of the Mechanics and Physics of Solids 115, 1–29.  doi:10.1016/j.jmps.2018.03.003.
  * Ryshkewitch (1953) Ryshkewitch, E., 1953.  Compression Strength of Porous Sintered Alumina and Zirconia.  Journal of the American Ceramic Society 36, 65–68.  doi:10.1016/b978-0-408-15960-9.50076-0.
  * Salvini et al. (2018) Salvini, V.R., Pandolfelli, V.C., Spinelli, D., 2018.  Mechanical Properties of Porous Ceramics.  Recent Advances in Porous Ceramics doi:10.5772/intechopen.71612.
  * Sari et al. (2020) Sari, M., Alevizos, S., Poulet, T., Lin, J., Veveakis, M., 2020.  A visco-plastic framework for interface processes in sedimentary reservoir rocks at HPHT conditions.  Geomechanics for Energy and the Environment 22, 100165.  doi:10.1016/j.gete.2019.100165.
  * Singh and Rao (2005) Singh, M., Rao, K.S., 2005.  Bearing Capacity of Shallow Foundations in Anisotropic Non-Hoek–Brown Rock Masses.  Journal of Geotechnical and Geoenvironmental Engineering 131, 1014–1023.  doi:10.1061/(asce)1090-0241(2005)131:8(1014).
  * Terzaghi (1943) Terzaghi, K., 1943.  Theoretical Soil Mechanics.  John Wiley and Sons Inc.  doi:10.1680/geot.1964.14.1.1.
  * Tonks and Millett (2011) Tonks, M., Millett, P., 2011\.  Phase field simulations of elastic deformation-driven grain growth in 2D copper polycrystals.  Materials Science and Engineering A 528, 4086–4091.  doi:10.1016/j.msea.2011.02.007.
  * Tourret et al. (2021) Tourret, D., Liu, H., LLorca, J., 2021.  Phase-field modeling of microstructure evolution: Recent applications, perspectives and challenges.  Progress in Materials Science , 100810doi:10.1016/j.pmatsci.2021.100810, arXiv:2104.09915.
  * Vardoulakis and Sulem (1995) Vardoulakis, I., Sulem, J., 1995\.  Bifurcation analysis in geomechanics.  doi:10.1201/9781482269383.
  * Veveakis et al. (2014) Veveakis, E., Poulet, T., Alevizos, S., 2014.  Thermo-poro-mechanics of chemically active creeping faults. 1: Theory and steady state considerations.  Journal of Geophysical Research: Solid Earth 119, 4583–4605.  doi:10.1002/2013JB010071.
  * Wachter et al. (2002) Wachter, N.J., Krischak, G.D., Mentzel, M., Sarkar, M.R., Ebinger, T., Kinzl, L., Claes, L., Augat, P., 2002\.  Correlation of bone mineral density with strength and microstructural parameters of cortical bone in vitro.  Bone 31, 90–95.  doi:10.1016/S8756-3282(02)00779-2.
  * Wang et al. (2019) Wang, D.J., Tang, H., Shen, P., Cai, Y., 2019.  A Parabolic Failure Criterion for Transversely Isotropic Rock: Modification and Verification.  Mathematical Problems in Engineering 2019\.  doi:10.1155/2019/8052560.
  * Wetzel et al. (2021) Wetzel, M., Kempka, T., Kühn, M., 2021.  Diagenetic trends of synthetic reservoir sandstone properties assessed by digital rock physics.  Minerals 11, 1–21.  doi:10.3390/min11020151.
  * Yuan et al. (2020) Yuan, C., Guo, Y., Wang, W., Cao, L., Fan, L., Huang, C., 2020.  Study on “Triaxial Loading-Unloading-Uniaxial Loading” and Microscopic Damage Test of Sandstone.  Frontiers in Earth Science 8, 1–11.  doi:10.3389/feart.2020.00078.
  * Zhang et al. (2016) Zhang, Y.D., Buscarnera, G., Einav, I., 2016.  Grain size dependence of yielding in granular soils interpreted using fracture mechanics, breakage mechanics and Weibull statistics.  Geotechnique 66, 149–160.  doi:10.1680/jgeot.15.P.119.


# The Effect of Epigenetic Blocking on Dynamic Multi-Objective Optimisation
Problems

Sizhe Yuen  1 University of SouthamptonMaritime EngineeringSouthamptonUnited
Kingdom<EMAIL_ADDRESS>,  Thomas H.G. Ezard  1 National Oceanography
CentreOcean and Earth ScienceSouthamptonUnited Kingdom<EMAIL_ADDRESS>and
Adam J. Sobey  1 The Alan Turing InstituteMarine and Maritime Group, Data-
centre EngineeringLondonUnited Kingdom<EMAIL_ADDRESS>
###### Abstract.

Hundreds of Evolutionary Computation approaches have been reported. From an
evolutionary perspective they focus on two fundamental mechanisms: cultural
inheritance in Swarm Intelligence and genetic inheritance in Evolutionary
Algorithms. Contemporary evolutionary biology looks beyond genetic
inheritance, proposing a so-called “Extended Evolutionary Synthesis”. Many
concepts from the Extended Evolutionary Synthesis have been left out of
Evolutionary Computation as interest has moved towards specific
implementations of the same general mechanisms. One such concept is epigenetic
inheritance, which is increasingly considered central to evolutionary
thinking. Epigenetic mechanisms allow quick non- or partially-genetic
adaptations to environmental changes. Dynamic multi-objective optimisation
problems represent similar circumstances to the natural world where fitness
can be determined by multiple objectives (traits), and the environment is
constantly changing.

This paper asks if the advantages that epigenetic inheritance provide in the
natural world are replicated in dynamic multi-objective optimisation problems.
Specifically, an epigenetic blocking mechanism is applied to a state-of-the-
art multi-objective genetic algorithm, MOEA/D-DE, and its performance is
compared on three sets of dynamic test functions, FDA, JY, and UDF. The
mechanism shows improved performance on 12 of the 16 test problems, providing
initial evidence that more algorithms should explore the wealth of epigenetic
mechanisms seen in the natural world.

Genetic algorithms, multi-objective optimisation, dynamic optimisation,
epigenetics

††copyright: rightsretained††conference: The Genetic and Evolutionary
Computation Conference; July 9–13; Boston, MA, USA††ccs: Theory of computation
Evolutionary algorithms††ccs: Computing methodologies Genetic algorithms

##  1\. Borrowing from the Extended Evolutionary Synthesis

The Modern Synthesis (Huxley, 1943), a combination of Darwin and Wallace’s
ideas of natural selection (Darwin and Wallace, 1858; Darwin, 1859), and
Mendel’s principles of inheritance (Bateson et al., 1902), has been an
important inspiration for the concepts used in Evolutionary Algorithms.
Despite the number of approaches available, the core inspiration is often just
genetic inheritance. However, modern evolutionary theory has since continued
to explore the mechanisms of evolution, extending the Modern Synthesis to
include concepts of non-genetic inheritance such as epigenetics, parental
effects, multilevel selection, and cultural inheritance in a portfolio
proposed as the Extended Evolutionary Synthesis (Pigliucci, 2007; Pigliucci
and Muller, 2010).

Epigenetic mechanisms in evolutionary theory alter DNA expression, leading to
a change in phenotype without a change in the underlying genotype (Dolinoy,
2008). This “phenotypic plasticity” leads to a faster rate of change to
quickly adapt to changes in the environment, and the ability to revert changes
if the environmental conditions do not activate the epigenetic mechanism. This
is to allow adaptation to a natural world that is changing, optimising
organismal fitness without altering its underlying genotype. Dynamic problems
reflect a similar set of challenges by changing the optimal Pareto set or
Pareto front over time (Farina et al., 2004). This paper therefore explores
if, and, if so, how epigenetics might be used within Evolutionary Algorithms
to improve performance. Due to the faster-than-generational adaptation
capabilities of epigenetics, these are first explored in a dynamic multi-
objective problem context.

There is limited literature exploring epigenetic mechanisms in genetic
algorithms; those that exist are focused on static single-objective problems.
Results are nevertheless consistent with expectations based on evolutionary
theory. The Knapsack problem was used to determine the performance of the
cytosine methylation epigenetic process when added to a traditional genetic
algorithm (Chrominski et al., 2020), resulting in a 25-30% reduction in the
number of generations needed to reach an optimal solution. A constant
probability of methylation is used, so the rate at which genes are blocked do
not react to dynamic changes based on external factors or progression of the
search. In epigenetic models, constant rates of variation are non-specific to
the environment and act similar to mutation, while varying rates are directed
by environmental cues (Jablonka et al., 1992).

##  2\. Epigenetic blocking mechanism

There are three forms of epigenetic transfer possible: mitotic, germline, and
experience-dependent (Vaughn et al., 2007; Champagne, 2008). Germline and
experience-dependent transfer passes down epigenetic marks that direct the
epigenetic process for future generations, while mitotic transfer only
propagates changes in the same generation. For an initial exploration of
epigenetic mechanisms, a simple form of genetic blocking is chosen for
inspiration based on mitotic transfer. A probabilistic blocking mechanism is
used to block some variables in each individual from being changed during
crossover. The parameters of the epigenetic process are not inherited across
multiple generations (Figure 1). A simple mechanism allows the properties to
be analysed and to more clearly understand the effect of an epigenetic process
on the performance of an algorithm. The simplicity also allows the mechanism
to be adapted to any Evolutionary Algorithm by altering the crossover method
without adding complicated features.

The mechanism has a probability to trigger during the reproduction stage for
every parent without bias towards fitness as blocking both the fitter and
less-fit parents have merits. Blocking fitter parents reduces stagnation of
the population, i.e. maintains diversity, should the variables or objectives
of the dynamic problem changes, while blocking less fit parents increases the
convergence of the population through more rapid selection at the variable
level.

Figure 1. The blocking mechanism where some variables are blocked from
carrying over to the next generation.

By varying the probability and the number of genes that are blocked, and
controlling the duration of each dynamic cycle in a problem, the impact of
epigenetic blocking can be analysed and compared to a baseline algorithm, and
to a constant probability of triggering the mechanism.

##  3\. Experimental setup

The MOEA/D-DE algorithm (Li and Zhang, 2009) with the re-initialisation
strategy outlined in (Biswas et al., 2014) was chosen as the base algorithm
for benchmarking due to its strong performance on dynamic optimisation
problems. The fast converging nature of MOEA/D lets it adapt quickly in
changing environments (Jiang and Yang, 2017). A population size of 500 is
used.

Three variants of the epigenetic blocking mechanism are investigated by
changing two parameters: the probability for the mechanism to trigger, and the
number of variables blocked in the process (block size). A summary of the
three variants is as follows:

  * •

E \- with a constant probability of 0.1 and a constant block size of 6.

  * •

EIB \- with a constant probability of 0.1 and a varying block size from 1 up
to the number of variables in the problem.

  * •

EIP \- with a varying probability from 0 up to a maximum of 0.8 with a
constant block size.

The gradual increase in either probability or block size is intended to
increase the convergence of the population, as blocking more prevents diverse
changes. The maximum probability is limited to 0.8 to prevent stagnation where
the blocking occurs too often. Increases to probability are rounded to 0.01,
increasing every 2 generations to the maximum probability. The block size
increases are rounded to the nearest whole number and depends on the number of
variables. For example, a problem with 30 variables and a population of 500
will increase the block size every 4 generations.

Table 1. Hyperparameters table Algorithm | Hyperparameters
---|---
MOEA/D-DE |  $P(m)=1/D$, $P(x)=0.9,CR=1,F=0.5$
E | $P(b)=0.1,s=6$
EIB | $P(b)=0.1,s=\frac{evals}{max\\_evals}*num\\_vars$
EIP | $P(b)=\frac{evals}{max\\_evals}*0.8,s=6$

Details of these hyperparameters are shown in in Table 1. $P(b)$ is the
probability for the mechanism to trigger and $s$ is the block size. The
Inverted Generational Distance (IGD) (Sato et al., 2004) metric is used to
show the performance of the epigenetic variants to the MOEA/D-DE algorithm.

The FDA (Farina et al., 2004), JY (Jiang and Yang, 2017), and UDF (Biswas et
al., 2014) benchmark functions are chosen to test the performance of the
epigenetic blocking mechanism. The FDA and UDF functions are considered to be
simpler to solve as they are based on the multi-objective ZDT (Deb et al.,
2002) problems while the JY problems are more complex (Jiang et al., 2020),
including elements such as linkage between variables and multiple knee points.
For the benchmark problem properties, $\tau$ is set at 5 and $nT$ at 10. This
gives 5 generations before the problem changes and 10 distinct steps. In total
this gives 100 generations (50,000 iterations) to complete a full cycle back
to the original variables and objectives of the problem. 2 full cycles are
benchmarked with 20 independent runs for analysis.

##  4\. Performance on dynamic multi-objective optimisation problems

To determine how the changes in the block rate and the probability of change
affect the performance, the three variants are compared by taking the average
IGD every generation to match the rate at which the probability varies. The
difference between each variant and the baseline are then compared to
demonstrate the improvement in performance.

###  4.1. Performance of epigenetic approaches against each other

Table 2. The total % difference for a two full dynamic cycles. The p-values
from a Wilcoxon signed-rank test is shown in brackets, bold indicates the best
performing variant. The lighter blue boxes indicate an improvement in
performance; darker red boxes indicate a decrease.

Problem | E | EIB | EIP
---|---|---|---
FDA1 | -102 (0.000) | -14 (0.027) | -188 (0.015)
FDA2 | 399 (0.000) | 400 (0.000) | 316 (0.000)
FDA3 | 163 (0.000) | 189 (0.016) | 810 (0.000)
JY1 | 657 (0.000) | 504 (0.000) | 917 (0.000)
JY2 | 193 (0.012) | 41 (0.783) | 470 (0.000)
JY3 | 669 (0.000) | 1164 (0.000) | 422 (0.000)
JY5 | -231 (0.000) | -28 (0.518) | 226 (0.002)
JY6 | -19 (0.209) | -15 (0.164) | -30 (0.020)
JY7 | 96 (0.871) | 318 (0.029) | 391 (0.001)
JY8 | 16 (0.108) | 29 (0.028) | 37 (0.001)
UDF1 | -181 (0.031) | 170 (0.013) | -35 (0.687)
UDF2 | -52 (0.323) | 15 (0.728) | -40 (0.024)
UDF3 | 61 (0.000) | -17 (0.000) | 35 (0.586)
UDF4 | 124 (0.003) | 413 (0.000) | -71 (0.027)
UDF5 | -2 (0.217) | 122 (0.262) | -120 (0.287)
UDF6 | 206 (0.000) | 180 (0.000) | 128 (0.000)

The total percentage difference between the baseline and each variant are
summarised in Table 2. Positive values represent the epigenetic mechanism
performing better than the baseline algorithm and negative values a decrease.
Values highlighted have a p-value lower than 0.05 from a Wilcoxon signed-rank
test and are statistically significant at the 95% level.

The performance summary of the three variants is as follows:

  * •

E \- Positive performance on 8 out of 16 problems, best on 2 problems.

  * •

EIB \- Positive performance on 10 out of 16 problems, best on 4 problems.

  * •

EIP \- Positive performance on 9 out of 16 problems, best on 6 problems.

The presence of epigenetic blocking generally improves the performance
compared to the baseline. It is only the FDA1 case where all of the epigenetic
enhanced algorithms struggle. Both EIB and EIP generally outperform the E
variant, where the blocking rate and probability of blocking remain static.
This is expected as using a constant probability and block size does not allow
the mechanism to react dynamically to changes in the problem.

EIP performs best on the most problems, 6, with 5 of these 6 problems in the
JY problem set. The performance suggests that EIP works best on more complex
problems with its aggressive use of the epigenetic mechanic compared to EIB
and E. However, EIP doesn’t perform well on the UDF problem set, considered to
be the simplest problems that are usually dominated by convergence, with 3
non-detectable results and 2 negative results where it performs worse. The
aggressive behaviour of this variant leads to larger positive and negative
results making it less robust.

In comparison, the EIB variant shows more consistency across the problem sets
with only 1 significant negative result and is able to achieve a best result
on each problem set. Although EIP can find stronger positive results in some
problems, its inconsistency and poor performance on the UDF problem set shows
that it is more suitable for complex problems, whereas EIB finds positive
improvement across a larger range of problems. Both of these approaches
outperform the E mechanism, indicating that increasing blocking over the
lifecycle was beneficial.

###  4.2. Performance against the baseline

The total difference in performance between the baseline and EIB over the 2
cycles is shown in Figure 2. The difference on each two generation interval is
summed with dark red denoting when the baseline performs better and light blue
denoting when EIB performs better. In problems such as JY5 and JY6, both
algorithms have parts of the search where they perform well, leading to no
detectable difference as can be seen in the p-value in Table 2. In problems
with statistically significant negative performance such as FDA1, it is not
overwhelmingly negative, with a positive improvement 40% of the time. The
results show strong performance overall against the baseline, indicating the
strength of the epigenetic mechanism in solving dynamic problems.

(a) Category I problems

(b) Category II problems

(c) Category III problems

Figure 2. Comparison of EIB against the baseline MOEA/D-DE algorithm

##  5\. Discussion

The addition of an epigenetic blocking mechanism is able to significantly
improve the performance of the MOEA/D-DE algorithm on a range of dynamic
problems.

Both the EIB and EIP variants performed better on more problems compared to E,
suggesting that varying the rates at which the mechanism triggers and the
variables that are blocked is more effective than keeping the parameters
static. From evolutionary theory, this behaviour is expected in a dynamic
environment with epigenetic variations that are guided rather than random
(Jablonka et al., 1992). The varying rates give more control over the
convergence and diversity of the populations. The increased consistency from
EIB compared to EIP shows the effectiveness of blocking more variables, even
at a slower rate. Convergence is increased with more variables blocked, hence
it performs well in general and on convergence based problems like UDF. EIP
blocks fewer variables but more often, finding success when the correct
variables are blocked, but worse performance on the wrong variables. More
diversity can be retained with EIP as only a small number of variables are
blocked compared to blocking a majority of variables in EIB. The increased
diversity would explain EIP’s better performance on the complex JY problems,
and it would be expected for EIP to perform well on complex real world
problems where diversity plays a larger role than convergence.

The epigenetic mechanism shows more sensitivity to Pareto set changes because
the variables are directly blocked. Better performance is observed on Category
II and III problems where there are Pareto front changes, and worse
performance on the Category I problem where there are only Pareto set changes.
A change in the Pareto front is analogous to a change in objective (fitness
optimum) in the natural world, which epigenetic processes can more quickly
adapt to. A change in the Pareto set is more complicated and can be related to
the concept of plasticity (Price et al., 2003), where the environment
influences the developmental stage of phenotypes.

An important advantage of this epigenetic mechanism is the ability to include
the mechanism into any crossover method of a genetic algorithm. The mechanism
sits on top of the crossover method, altering it slightly by blocking the
variables. In the natural world, epigenetics also sits on top of existing
genetic mechanisms to allow for quick adaptation.

##  6\. Conclusion

This paper demonstrates the impact of including a simple epigenetic blocking
mechanism into an existing evolutionary algorithm, and the advantages
epigenetics can provide to algorithms solving dynamic problems. The epigenetic
mechanism improves upon the baseline MOEA/D-DE algorithm on 12 of the 16
dynamic multi-objective test problems, with a conclusively negative result on
only 1 problem. Increasing the number of variables blocked helps increase
convergence, giving consistent improvement on all categories of test problems.
Increasing the rate of which the mechanism is activated is more effective on
the complex JY problems, suggesting a retention of diversity in the
population.

There are further epigenetic mechanisms and features to look into. In nature,
epigenetics enable fast adaptations triggered by environmental cues, and the
epigenetic method itself is inherited by future generations. The next step in
the proposed blocking mechanism is to include environmental cues to direct how
the mechanism should be triggered, and a form of inheritance to direct which
variables should be blocked.

## References

  * (1)
  * Bateson et al. (1902) William Bateson, Gregor Mendel, and William Morton Wheeler. 1902\.  _Mendel’s principles of heredity; a defence_.  Cambridge,University press. 236 pages.
  * Biswas et al. (2014) Subhodip Biswas, Swagatam Das, Ponnuthurai Suganthan, and Carlos Coello. 2014.  Evolutionary multiobjective optimization in dynamic environments: A set of novel benchmark functions.  _Proceedings of the 2014 IEEE Congress on Evolutionary Computation, CEC 2014_ , 3192–3199.  https://doi.org/10.1109/CEC.2014.6900487
  * Champagne (2008) Frances A Champagne. 2008\.  Epigenetic mechanisms and the transgenerational effects of maternal care.  _Front Neuroendocrinol_ 29 (2008), 386,397.  Issue 3.
  * Chrominski et al. (2020) Kornel Chrominski, Magdalena Tkacz, and Mariusz Boryczka. 2020\.  Epigenetic Modification of Genetic Algorithm. In _Computational Science – ICCS 2020_. Springer International Publishing, Cham, 267–278.
  * Darwin (1859) Charles Darwin. 1859\.  _On the Origin of Species by Means of Natural Selection_.  Murray, London.
  * Darwin and Wallace (1858) Charles Darwin and Alfred Wallace. 1858.  On the Tendency of Species to form Varieties; and on the Perpetuation of Varieties and Species by Natural Means of Selection.  _Journal of the Proceedings of the Linnean Society of London. Zoology_ 3, 9 (1858), 45–62.
  * Deb et al. (2002) Kalyanmoy Deb, L. Thiele, Marco Laumanns, and Eckart Zitzler. 2002\.  Scalable multi-objective optimization test problems.  _Proc. of 2002 IEEE Congress on Evolutionary Computation_ 1 (06 2002), 825–830.
  * Dolinoy (2008) Dana C Dolinoy. 2008\.  The agouti mouse model: an epigenetic biosensor for nutritional and environmental alterations on the fetal epigenome.  _Nutrition reviews_ 66 Suppl 1 (2008), S7 – S11.
  * Farina et al. (2004) M. Farina, K. Deb, and P. Amato. 2004.  Dynamic multiobjective optimization problems: test cases, approximations, and applications.  _IEEE Transactions on Evolutionary Computation_ 8, 5 (2004), 425–442.  https://doi.org/10.1109/TEVC.2004.831456
  * Grudniewski and Sobey (2021) Przemyslaw A. Grudniewski and Adam J. Sobey. 2021.  cMLSGA: A Co-Evolutionary Multi-Level Selection Genetic Algorithm for Multi-Objective Optimization.  _CoRR_ abs/2104.11072 (2021).  arXiv:2104.11072
  * Huxley (1943) J. Huxley. 1943\.  _Evolution, the modern synthesis_.  New York ; London: Harper & brothers.
  * Jablonka et al. (1992) Eva Jablonka, Michael Lachmann, and Marion J. Lamb. 1992\.  Evidence, mechanisms and models for the inheritance of acquired characters.  _Journal of Theoretical Biology_ 158, 2 (1992), 245–268.  https://doi.org/10.1016/S0022-5193(05)80722-2
  * Jiang et al. (2020) Shouyong Jiang, Marcus Kaiser, Shengxiang Yang, Stefanos Kollias, and Natalio Krasnogor. 2020\.  A Scalable Test Suite for Continuous Dynamic Multiobjective Optimization.  _IEEE Transactions on Cybernetics_ 50, 6 (2020), 2814–2826.  https://doi.org/10.1109/TCYB.2019.2896021
  * Jiang and Yang (2017) Shouyong Jiang and Shengxiang Yang. 2017.  Evolutionary Dynamic Multiobjective Optimization: Benchmarks and Algorithm Comparisons.  _IEEE Transactions on Cybernetics_ 47, 1 (2017), 198–211.  https://doi.org/10.1109/TCYB.2015.2510698
  * Li and Zhang (2009) H. Li and Q. Zhang. 2009.  Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II.  _IEEE Transactions on Evolutionary Computation_ 13, 2 (2009), 284–302.
  * Pigliucci (2007) Massimo Pigliucci. 2007\.  Do We Need an Extended Evolutionary Synthesis?  _Evolution_ 61, 12 (2007), 2743–2749.
  * Pigliucci and Muller (2010) Massimo Pigliucci and Gerd Muller. 2010.  _Evolution – the Extended Synthesis_.  MIT Press.
  * Price et al. (2003) Trevor Price, Anna Qvarnström, and Darren Irwin. 2003\.  The role of phenotypic plasticity in driving genetic evolution.  _Proceedings. Biological sciences / The Royal Society_ 270 (08 2003), 1433–40.
  * Sato et al. (2004) Hiroyuki Sato, Hernán E. Aguirre, and Kiyoshi Tanaka. 2004.  Local dominance using polar coordinates to enhance multiobjective evolutionary algorithms.  _Proceedings of the 2004 Congress on Evolutionary Computation (IEEE Cat. No.04TH8753)_ 1 (2004), 188–195 Vol.1.
  * Stolfi and Alba (2018) Daniel H. Stolfi and Enrique Alba. 2018.  Epigenetic algorithms: A New way of building GAs based on epigenetics.  _Information Sciences_ 424 (2018), 250 – 272.
  * Vaughn et al. (2007) Matthew W Vaughn, Miloš Tanurdžić, Zachary Lippman, Hongmei Jiang, Robert Carrasquillo, Pablo D Rabinowicz, Neilay Dedhia, W. Richard McCombie, Nicolas Agier, Agnès Bulski, Vincent Colot, R.W Doerge, and Robert A Martienssen. 2007.  Epigenetic Natural Variation in Arabidopsis thaliana.  _PLOS Biology_ 5, 7 (06 2007), 1–13.
  * Yuen et al. (2021) Sizhe Yuen, Thomas H. G. Ezard, and Adam J. Sobey. 2021\.  Epigenetic opportunities for Evolutionary Computation.  arXiv:2108.04546 [cs.NE]