#  Production cross sections and angular distributions of neutron-rich rare
isotopes from 15 MeV/nucleon Kr-induced collisions: toward the r-process path

O. Fasoulaa, G.A. Souliotisa,∗, Y.K. Kwonb, K. Tshoob M. Veselskyc, S.J.
Yennellod and A. Bonaserad,e a Laboratory of Physical Chemistry, Department of
Chemistry, National and Kapodistrian University of Athens, Athens, Greece  b
The Rare Isotope Science Project (RISP), Institute for Basic Science, Daejeon,
Korea  c Institute of Experimental and Applied Physics, Czech Technical
University, Prague, Czech Republic  d Cyclotron Institute, Texas A&M
University, College Station, Texas, USA  e Laboratori Nazionali del Sud, INFN,
Catania, Italy  ∗ Corresponding author. Email<EMAIL_ADDRESS>
###### Abstract

We present our recent study of cross sections and angular distributions of
projectile fragments from heavy-ion reactions at beam energy of 15
MeV/nucleon. We studied the production cross sections and the angular
distributions of neutron-rich nuclides from collisions of a 86Kr (15
MeV/nucleon) beam with heavy targets (64Ni, 124Sn and 238U). Experimental data
from our previous work at Texas A & M were compared with model calculations.
Our calculations were based on a two-step approach: the dynamical stage of the
collision was described with, first, the phenomenological Deep-Inelastic
Transfer model (DIT) and, alternatively, with the microscopic Constrained
Molecular Dynamics model (CoMD). The de-excitation of the hot heavy projectile
fragments was performed with the Statistical Multifragmentation Model (SMM).
An overall good discription of the available data was obtained with the models
employed. Furthermore, we performed calculations with a radioactive beam of
92Kr (15 MeV/nucleon) interacting with a target of 238U. We observed that the
multinucleon transfer mechanism leads to extremely neutron-rich nuclides
toward and beyond the astrophysical r-process path.

##  1 Introduction

The exploration of the nuclear landscape toward the astrophysical r-process
path and the neutron drip-line are at the center of interest of the nuclear
physics community [1]. Essential to this development is the efficient
production of neutron-rich nuclides which constitutes a central issue in
current and upcoming rare isotope beam facilities worldwide (see, e.g., [2, 3,
4, 5, 6, 7, 8, 9, 10, 11]).

The traditional routes to produce neutron-rich nuclides are spallation,
fission and projectile fragmentation [12]. Spallation is an efficient
mechanism to produce rare isotopes for ISOL-type techniques [13]. Projectile
fission is appropriate in the region of asymmetric fission peaks of the light
and heavy fission fragments (see, e.g., [14]). Finally, projectile
fragmentation constitutes a general approach to produce exotic nuclei at beam
energies above 100 MeV/nucleon (see, e.g., [15, 16]). This approach is,
nevertheless, limited by the fact that optimum neutron excess in the fragments
is achieved by stripping the maximum possible number of protons (and a minimum
possible number of neutrons). Surpassing the limits of the traditional
approaches and reaching out to the neutron dripline is highly desirable at
present. Thus, the study of new synthesis routes constitutes a vigorous
endeavor of the nuclear science community.

Toward this end, to reach a high neutron-excess in the products, apart from
proton stripping, it is necessary to pickup neutrons from the target. Such a
possibility is offered by reactions of nucleon exchange at beam energies from
the Coulomb barrier [17, 18] to the Fermi energy (20–40 MeV/nucleon) [19, 20].
Detailed experimental data in this broad energy range are limited at present
[18, 21, 22]. In multinucleon transfer and deep-inelastic reactions near the
Coulomb barrier [18], the low velocities of the fragments and the wide angular
and ionic charge state distributions may limit the collection efficiency for
the most neutron-rich products. However, the reactions in the Fermi energy
regime (15–35 MeV/nucleon) combine the advantages of both low-energy (i.e.,
near and above the Coulomb barrier) and high-energy (i.e., above 100
MeV/nucleon) reactions. At this energy, the interaction of the projectile with
a neutron-rich target enhances the N/Z of the fragments, while the velocities
remain high enough to allow efficient in-flight collection and separation.

Our initial experimental studies of projectile fragments from 25 MeV/nucleon
reactions of 86Kr on 64Ni [19] and 124Sn [20] indicated substantial production
of neutron-rich fragments. Motivated by developments in several facilities
that will offer either very intense primary beams [4, 7, 9] at this energy
range or re-accelerated rare isotope beams [3, 4, 7, 8, 9], we continued our
experimental and theoretical studies at 15 MeV/nucleon [23, 29, 30, 31].

In this contribution, we present our study of the production cross sections
and the angular distributions of projectile fragments from collisions of a
86Kr (15 MeV/nucleon) beam with heavy targets. Data from our experimental work
at Texas A&M are compared with calculations based on either the
phenomenological deep-inelastic transfer (DIT) model, or the microscopic
constrained molecular dynamics model (CoMD). A overall good description of the
experimental data with DIT or CoMD was obtained, suggesting the possibility of
using the present theoretical models to predict the production of exotic
nuclei employing radioactive beams. As an example, we will show production
cross sections and rates of neutron-rich isotopes from a radioactive beam of
92Kr (15 MeV/nucleon) interacting with a 238U target.

##  2 Outline of Results and Comparisons

A detailed presentation of our previously obtained experimental data appear in
[23] in which the mass spectrometric measurements of production cross sections
of neutron-rich projectile fragments from the reactions of a 15 MeV/nucleon
86Kr beam with 64,58Ni and 124,112Sn targets were presented. We also mention
that the experimental data of 25 MeV/nucleon 86Kr-induced reactions are
described in detail in our articles [19, 20, 21, 22].

###  2.1 Study of Production Cross Sections

Figure 1:  Experimental mass distributions (open black circles) of elements
with Z = 37–30 from the reaction 86Kr(15 MeV/nucleon) + 124Sn [23] compared
with DIT/SMM calculations (yellow line) and CoMD/SMM calculations (red line).

In Fig. 1, we present the experimental mass distributions of projectile-like
fragments with Z = 37–30 from the reaction 86Kr (15 MeV/nucleon) + 124Sn [23]
and compare them with the calculations using the DIT code [24, 25] (solid
yellow line) and the CoMD code [27, 28] (solid red line). Both dynamical codes
were combined with the Statistical Multifragmentation Model (SMM) [26] used
for the de-excitation of the hot quasi-projectiles emerging after the
dynamical stage of the reaction. The results of the calculations are in
overall good agreement with the experimental data especially for isotopes
close to the projectile with Z = 36–32. We note that in some of the Kr
isotopes (A=86,89,90), the experimental cross sections are especially high
because of contamination from elastically scattered beam. The overestimation
of the cross sections for n-rich products with Z=31 and 30 from the CoMD
calculation may be related to issues of low excitation energy, as currently
calculated in CoMD [29].

Figure 2:  DIT/SMM calculated mass distributions of elements with Z = 37–30
from the reaction of a radioactive beam of 92Kr (15 MeV/nucleon) interacting
with a 238U target (dashed yellow line), compared to the DIT/SMM calculations
of the stable beam reaction 86Kr + 238U (solid yellow line). The experimental
mass distributions (open points) from 86Kr (15 MeV/nucleon) \+ 124Sn [23] are
also given for reference.

Furthermore, in Fig. 2 we show the DIT/SMM calculations of projectile-like
fragments from the reaction 86Kr(15 MeV/nucleon) + 238U (solid yellow line).
We have chosen the heaviest and most neutron-rich target available in order to
explore how far in N/Z we can go with this reaction. Also, we chose to study
the reaction with the 238U target, but with a radioactive beam of 92Kr(15
MeV/nucleon) (dashed yellow line). As there are no experimental data with the
238U target, we used the data from 86Kr (15 MeV/nucleon) + 124Sn for
reference. We observe that the distributions extend far to the neutron-rich
side especially when the radioactive beam is used. This tendency is pronounced
for isotopes rather close to the projectile (Z=34–36).

We point out that, for, e.g., bromine (Z=35), isotopes that have up to 15 more
neutrons (A = 96) than the corresponding stable isotope (A = 81) can be
obtained. Thus, by using neutron-rich radioactive beams, and via the mechanism
of peripheral multinucleon transfer, we have the possibility to produce
extremely neutron-rich nuclides. We note that similar observations were made
for the reaction with lighter projectiles (of Ar and Ca) on the 238U target in
our recent article [31]. We do not have experimental data for the 86Kr+238U
reaction, but we have plans to study it in the near future.

Figure 3:  Representation of DIT/SMM calculated production cross sections of
projectile fragments from the radioactive-beam reaction 92Kr (15 MeV/nucleon)
+ 238U on the Z–N plane. Cross section ranges are shown by open circles
according to the key. Closed squares show the stable isotopes. The purple line
shows the astrophysical r-process path and the green line shows the neutron
drip-line. The horizontal and vertical lines indicate, respectively, the Z and
N of the 92Kr projectile.

A comprehensive presentation of the DIT/SMM calculated production cross
sections of the projectile fragments from the above reaction on the Z vs N
plane is given in Fig. 3. In this figure, stable isotopes are represented by
closed squares, whereas fragments obtained by the 92Kr+238U reaction are given
by the open circles, with sizes corresponding to cross-section ranges
according to the figure key. The green line gives the location of the neutron
drip-line and the purple line indicates the expected path of the astrophysical
rapid neutron-capture process (r-process). Using this representation of the
previous calculation , we clearly observe that the neutron pickup products
from 92Kr reach and even exceed the path of the r-process near Z=30–36.

In Table I, we summarize the calculated cross sections and production rates of
several neutron-rich isotopes from the reaction of the radioactive beam of
92Kr (15 MeV/nucleon) with the 238U target. For the rate calculations, the
92Kr beam with intensity of 1.0 pnA (6.2$\times$109 particles/sec) is assumed
to interact with a 238U target of 20 mg/cm2 thickness. We see that it is
possible to produce extremely neutron-rich isotopes in these reactions with
the use of re-accelerated radioactive beams, such as 92Kr, that will become
available in upcoming rare-isotope facilities (e.g. [9, 10, 11]). Along these
lines, we wish to continue this work at the following facilities: a) the
Cyclotron Institute of Texas A&M University using the MARS separator, b) the
LNS/Catania using beams from the S800 Cyclotron and employing the MAGNEX
spectrometer, and c) the RISP facility with stable and radioactive beams from
the RAON accelerator complex using the KOBRA separator.

Table 1:  Cross sections and rate estimates (last column) of very neutron-rich isotopes from the reaction 92Kr (15 MeV/nucleon) + 238U. For the rates, a radioactive beam of 92Kr with intensity 1.0 particle nA (6.2$\times$109 particles/s) is assumed to interact with a 238U target of 20 mg/cm2 thickness.  Rare | Reaction | Cross | Rate (sec-1)
---|---|---|---
Isotope | Channel | Section (mb) |
93Kr | -0p+1n | 12 | 3.6$\times$103
94Kr | -0p+2n | 1.3 | 3.9$\times$102
95Kr | -0p+3n | 0.3 | 90
96Kr | -0p+4n | 0.05 | 15
92Br | -1p+1n | 0.8 | 2.4$\times$102
93Br | -1p+2n | 0.2 | 60
94Br | -1p+3n | 0.07 | 21
95Br | -1p+4n | 0.02 | 6
96Br | -1p+5n | 0.008 | 2
90Se | -2p+0n | 0.25 | 75
91Se | -2p+1n | 0.14 | 42
92Se | -2p+2n | 0.05 | 15
93Se | -2p+3n | 0.02 | 6

###  2.2 Study of Angular Distributions

After the description of the production cross sections, which have been the
main focus of our work so far, we initiated a detailed study of the angular
distributions of the projectile-like fragments. In regards to experimental
data, we refer again to our previous work [23] with the MARS spectrometer. As
described in [23], the data were obtained at two angle settings, by sending
the beam on the target at the appropriate angle in the MARS target chamber.
The angle settings were: the 4.0°setting and the 7.5°setting, covering the
polar angular ranges of $\Delta\theta$=2.2°–5.8°, and
$\Delta\theta$=5.6°–9.2°, respectively. We remind that approptiate integration
of these data provided the total production cross sections that are reported
in [23] and presented in Figs. 1,2.

We start our presentation with the angular distributions of selected
projectile fragments from the reaction of 86Kr+124Sn at 15 MeV/nucleon
(Figures 4,5 and 6). For this system, the grazing angle is $\theta_{gr}\simeq$
9.0°. (A verical dashed line indicates this angle in the figures.) We, thus,
expect that the 7.5°angle setting is the appropriate one to efficiently
collect the quasi-elastic products, as has been discussed in [23, 29].

The experimental data of the differential cross sections for each nuclide
consist of two points: one at $\theta$=4.0°and the other at $\theta$=7.5°,
each plotted with horizontal error bars indicating the polar angular
acceptance of $\Delta\theta$ = 3.6°). (The vertical error bars indicate, as
usual, the statistical uncertainties). The data are compared with calculations
employing the DIT model (yellow symbols) and the CoMD model, both followed by
the SMM model, as also shown in the previous comparisons for the cross
sections. The calculations were binned and presented in angular steps of two
degrees.

Three different of calculations were performed with the CoMD model. a) The
”standard” calculation (red symbols) using standard parameters of the CoMD
code, as in our previous works [29, 31]. b) A calculation with pairing (blue
symbols). For this calculation, we introduced a phenomenological description
of n-n and p-p pairing, according to the recent article [32]. In this
description, for two neutrons (or protons ) that have anti-parallel spins and
energies near the Fermi energy, an additional attractive interaction was
introduced with strength adjusted to reproduce the empirical pairing energy of
the Bethe-Weizsacker equation. c) a calculation in which the compressibility
of cold symmetric nuclear matter was increased to K=272 MeV (purple symbols).
We remind that in the ”standard” CoMD calculation, the parameters of the
effective interaction correspond to a compressibility of K=200 MeV.

Figure 4:  Angular distributions of neutron rich isotopes of 36Kr and 38Sr
from the reaction of 86Kr (15 MeV/nucleon) with 124Sn. Black squares (with
horizontal errorbars): experimental data [23]. Symbols connected with dotted
lines: calculations as follows. Yellow symbols: DIT. Red symbols: standard
CoMD. Blue symbols: CoMD with pairing. Purple symbols: CoMD with
compressibility K=272 MeV.

In the upper half of Fig. 4, we show the angular distributions for Krypton
(Z=36). For A=86, we see that the data point at $\theta$=7.5°is higher than
that at $\theta$=4.0°, as expected. This step seems to be described by the DIT
calculation, despite the larger value of the latter at $\theta$=8°and
$\theta$=10°. The CoMD calculations (in all three variations) appear to
describe the shape of the angular distribution rather well; the CoMD
calculation with K=272 is however lower that the other calculations. For A=85
and A=84, there is some agreement of the calculations with the data. However,
for A=88 (+2n channel) the calculations cannot describe the step in the data.
This may be related to enhanced neutron-pair transfer, possibly beyond our
tentative pairing calculation.

In the lower half of Fig. 4, we show the angular distributions for Strontium
(Z=38, +2p channels). For A=88, we note that the +2p channel does not peak at
the grazing angle as the +2n channel (Fig. 8d). Moreover, for A=86, the
removal of neutrons leads to a flat distribution in the data. The addition of
neutrons, however, leads to angular distributions peaked at $\theta_{gr}$.
This feature cannot be described by DIT or CoMD, and requires further
investigation.

Figure 5:  Angular distributions of neutron rich isotopes of 35Br and 34Se
from the reaction of 86Kr (15 MeV/nucleon) with 124Sn. Black squares (with
horizontal errorbars): experimental data [23]. Symbols connected with dotted
lines: calculations as follows. Yellow symbols: DIT. Red symbols: standard
CoMD. Blue symbols: CoMD with pairing. Purple symbols: CoMD with
compressibility K=272 MeV.

In the upper half of Fig. 5, we show the angular distributions of several
isotopes of bromine (Z=35). For A=85 (-1p channel), the calculations describe
rather well the $\theta$=7.5°data point, but are higher than the
$\theta$=4.0°data point. Interestingly, the CoMD calculation with K=272 MeV
appears to describe both data points well. Again, the removal of neutrons
(A=83, -1p-2n) results in a nearlyflat angular distribution that seems to be
described by the calculations. Concerning the neutron pickup nuclides (A=86,
-1p+1n) and (A=87, -1p+2n), the experimental data are peaked at $\theta_{gr}$
and cannot be described by the DIT and CoMD calculations. We think that for
A=87 (-1p+2n), the experimental data may be contaminated from background of
elastically scattered beam. Moreover, especially for A=86, we speculate that
part of the discrepancy may be due to a contribution of a single charge
exchange (SCE) process that, of course, cannot be described by DIT or CoMD.

In the lower half of Fig. 5, we show the angular distributions of several
isotopes of Selenium (Z=34) (-2p channels). Observations similar to those of
the previous figure pertain here. For A=82 (-2p-2n), the experimental angular
distribution is rather flat and is well described by the calculations.
However, for A=84 (-2p), the experimental angular distribution is again peaked
at $\theta_{gr}$ and cannot be adequately described by the calculations. The
situation for the angular distributions of A=85 (-2p+1n) and A=86 (-2p+2n) is
similar to the previous ones. We may attribute the discrepancy for A=85 in
part to a contributing SCE process along with a +1n process. For A=86
(-2p+2n), we speculate that, a DCE (double charge exchange) process might also
contribute, along with two successive SCE processes. The subject of potential
contributions of SCE and DCE processes in these heavy-ion reactions below the
Fermi energy is currently under investigation by our group via systematic
studies of the momentum distributions of the relevant nuclides, and may shed
light to interesting aspects of the reaction mechanism(s).

In Fig. 6, we show the angular distributions for several isotopes of Arsenic
(Z=33, -3p channels) in the upper half, and Germanium (Z=32, -4p channels) in
the lower half. The angular distributions of these nuclides reveal features
analogous to those of the nuclides shown in the previous figures.
Specifically, the ”cold” proton removal channels (-3p and -4p, respectively)
are peaked at $\theta_{gr}$ (only one point in the data), whereas the DIT and
CoMD calculations show, as before, a monotonic decrease with increasing angle.
Furthermore, we notice a transition from a $\theta_{gr}$ peaked angular
distribution to a flat distribution between A=82 to A=80 for As, and, between
A=80 to A=78 for Ge.

Figure 6:  Angular distributions of neutron rich isotopes of 33As and 32Ge
from the reaction of 86Kr (15 MeV/nucleon) with 124Sn. Black squares (with
horizontal errorbars): experimental data [23]. Symbols connected with dotted
lines: calculations as follows. Yellow symbols: DIT. Red symbols: standard
CoMD. Blue symbols: CoMD with pairing. Purple symbols: CoMD with
compressibility K=272 MeV.

We now present the angular distributions of selected projectile fragments from
the reaction of 86Kr (15 MeV/nucleon) + 64Ni in Figs. 7 and 8. For this
system, the grazing angle is $\theta_{gr}\simeq$ 6.0°. (A vertical dashed line
indicates this angle in Figs. 7 and 8). We expect that the 4.0°angular setting
is the optimum for the efficient collection of quasi-elastic products [23].
The general feature of the experimental angular distributions is that the
4.0°data points are higher than the corresponding 7.5°ones. In the upper half
of Fig. 7, we show the angular distributions of several isotopes of Krypton
(Z=36). For A=86, the data point at $\theta$=4.0°is, as expected, higher than
that at $\theta$=7.5°. This step seems to be described by the calculations,
which are however lower than the data. For A=85 and A=84, there is a better
agreement of the calculations with the data. However, again in the two-neutron
pickup channel, A=88, the calculations are substantially lower than the
experimental data. In the lower part of Fig. 7, we present the angular
distributions of several isotopes of Strontium (Z=38, +2p channels). For these
isotopes, the calculations appear to describe the decreasing trend of the data
with increasing angle.

Figure 7:  Angular distributions of neutron rich isotopes of 36Kr and 38Sr
from the reaction of 86Kr (15 MeV/nucleon) with 64Ni. Black squares (with
horizontal errorbars): experimental data [23]. Symbols connected with dotted
lines: calculations as follows. Yellow symbols: DIT. Red symbols: standard
CoMD. Blue symbols: CoMD with pairing. Purple symbols: CoMD with
compressibility K=272 MeV.

In the upper part of Fig. 8, we show the angular distributions of several
isotopes of Bromine (Z=35, -1p channels). For A=85 (-1p), there is a
satisfactory description of the data by the calculations. This is also true
for the neutron removal channel, A=83 (-1p-2n ). As we move to the neutron
pickup channels, A=86 (-1p+1n), and A=87 (-1p+2n), the calculations appear to
be substantially lower than the experimental data. In the lower part of Fig.
8, we show the angular distributions of several nuclides of Selenium (Z=34,
-2p channels). Observations similar to those of the upper half of figure
pertain here. Specifically, for A=84 (-2p) the calculations appear to follow
the data to some extent. Moving to the neutron removal channel, for A=82
(-2p-2n) the calculations describe the data fairly well. Finally, for the
neutron pickup channels, A=85 (-2p+1n), and A=86 (-2p+2n) the calculations
again follow the trend of the data, but are lower than them.

As a general remark of the angular distribution study, we point out that a
rather satisfactory overall agreement of the calculations with the
experimental data was achieved. Several discrepancies were noticed and are the
subject of further investigation. As expected, the angular distributions of
the Kr+Sn system are more extended than the Kr+Ni system, which are more
forward focused due to inverse kinematics.

Referring to the Kr+Sn system, the DIT calculations are rather similar in most
cases to the CoMD calculations, especially to the standard CoMD calculation
and the CoMD calculation with the paring. Specifically, the CoMD calculation
with the pairing seems to be slightly higher than the standard CoMD
calculation for large angles (where, of course, there are no experimental
data). Finally, the effect of the compressibility in the CoMD code was
investigated by performing calculations with K=272 MeV, along with the
standard value of K=200 MeV in the other CoMD calculations. We observed that
with K=272, the CoMD calculated angular distributions are peaked more toward
the $\theta_{gr}$, as compared to those with the ”standard” value of K=200
MeV, and are thus in better agreement with the data. However, they have the
tendency to be lower that the other CoMD calculations.

Figure 8:  Angular distributions of neutron rich isotopes of 35Br and 34Se)
from the reaction of 86Kr (15 MeV/nucleon) with 64Ni. Black squares (with
horizontal errorbars): experimental data [23]. Symbols connected with dotted
lines: calculations as follows. Yellow symbols: DIT. Red symbols: standard
CoMD. Blue symbols: CoMD with pairing. Purple symbols: CoMD with
compressibility K=272 MeV.

##  3 Summary and Conclusions

In this contribution, we presented our recent studies of production cross
sections and angular distributions of projectile fragments from 86Kr-induced
reactions at 15 MeV/nucleon. We studied the cross sections and the angular
distributions of neutron-rich nuclides from collisions of a 86Kr (15
MeV/nucleon) beam with several heavy targets (64Ni, 124Sn and 238U).
Experimental data from our previous work at Texas A&M were compared with model
calculations. The calculations were based on a two-step approach: the
dynamical stage of the collision was described with either the Deep-Inelastic
Transfer model (DIT), or with the microscopic Constrained Molecular Dynamics
model (CoMD). The de-excitation of the hot heavy projectile fragments was
performed with the Statistical Multifragmentation Model (SMM). We first
studied the production cross sections of neutron-rich nuclides from collisions
of a 86Kr (15 MeV/nucleon) beam with 124Sn and 238U. We also proceeded with
calculations with a radioactive beam of 92Kr (15 MeV/nucleon) with 238U and
observed that the multinucleon transfer mechanism leads to very neutron-rich
nuclides in the mass region A=80–100 toward and beyond the r-process path.

Subsequently, we initiated a study of the angular distributions of projectile
fragments from the reactions of 86Kr (15 MeV/nucleon) with targets of 64Ni and
124Sn. We compared our experimental data at the two angle settings of the MARS
spectrometer, namely at 4.0°and at 7.5°, with detailed calculations using the
DIT and CoMD models. Three variants of the CoMD calculations were performed:
the standard calculation (as in all our previous works), b) a calculation
including a pairing interaction, and c) a calculation with compressibility
K=272 for symmetric nuclear matter. A rather satisfactory agreement of the
calculations with the experimental data was obtained. Several discrepancies
were noticed and are the subject of further investigation.

Our current studies indicate that heavy-ion reactions below the Fermi energy
can be exploited as an effective route to access extremely neutron-rich
isotopes toward the r-process path and the neutron drip-line. In this vein,
future experiments in several accelerator facilities may be planned that will
enable a variety of nuclear studies in unexplored regions of the nuclear
chart.

## References

  * [1] J. Erler et al, Nature 486, 509 (2011).
  * [2] D. F. Geesaman, C. K. Gelbke, R. V. F. Janssens, B. M. Sherrill, Ann. Rev. Nucl. Part. Sci. 56, 53 (2006)
  * [3] FRIB main page: www.frib.msu.edu
  * [4] GANIL main page: www.ganil.fr
  * [5] GSI main page: www.gsi.de
  * [6] RIBF main page: www.rarf.riken.go.jp/Eng/facilities/RIBF.html
  * [7] ATLAS main page: www.phy.anl.gov/atlas/facility/index.html
  * [8] EURISOL main page: www.eurisol.org
  * [9] RISP main page: www.risp.re.kr/eng/pMainPage.do
  * [10] K. Tshoo, Y. K. Kim, Y. K. Kwon et al, Nucl. Instrum. Methods B 317, 242 (2013).
  * [11] K. Tshoo, H. Chae, J. Park, J.Y. Moon, Y.K. Kwon, G.A. Souliotis et al, Nucl. Instrum. Methods B 376, 188 (2016).
  * [12] Y. Blumenfeld, T. Nilsson and P. Van Duppen, Phys. Scr. T152 014023 (2013).
  * [13] A. Kelić, M. V. Ricciardi, K. -H. Schmidt, BgNS Transactions, 13, 98 (2009).
  * [14] H. Alvarez-Pol et al., Phys. Rev. C 82, 041602 (2010).
  * [15] O. B. Tarasov et al., Phys. Rev. C 80, 034609 (2009).
  * [16] S. Lukyanov et al., Phys. Rev. C 80, 014609 (2009).
  * [17] V. V. Volkov, Phys. Rep. 44, 93 (1978).
  * [18] L. Corradi, G. Pollarolo, S. Szilner, J. Phys. G 36, 113101 (2009).
  * [19] G. A. Souliotis et al., Phys. Lett. B 543, 163 (2002).
  * [20] G. A. Souliotis et al., Phys. Rev. Lett. 91, 022701 (2003).
  * [21] G. A. Souliotis et al., Nucl. Instrum. Methods B 204 166 (2003).
  * [22] G. A. Souliotis et al., Nucl. Instum. Methods B 266, 4692 (2008).
  * [23] G. A. Souliotis et al., Phys. Rev. C 84, 064607 (2011).
  * [24] L. Tassan-Got and C. Stephan, Nucl. Phys. A 524, 121 (1991).
  * [25] M. Veselsky and G.A. Souliotis, Nucl. Phys. A 765, 252 (2006).
  * [26] J. Bondorf et al., Phys. Rep. 257, 133 (1995).
  * [27] M. Papa et al., Phys. Rev. C 64, 024612 (2001).
  * [28] M. Papa et al, J. Comp. Phys. 208, 403 (2005).
  * [29] P.N. Fountas, G.A. Souliotis, M. Veselsky and A. Bonasera, Phys. Rev. C 90, 064613 (2014).
  * [30] N. Vonta, G.A. Souliotis et al., Phys. Rev. C 92, 064611 (2016).
  * [31] A. Papageorgiou, G.A. Souliotis et al., Journal of Physics G 45, 095105 (2018).
  * [32] C. Agodi, G. Giuliani, F. Cappuzzello, A. Bonasera, et al., Physical Review C 97, 034616 (2018).


# Energetic Formulation of Large-Deformation Poroelasticity

Mina Karimi<EMAIL_ADDRESS>Department of Civil and Environmental
Engineering, Carnegie Mellon University  Mehrdad Massoudi  National Energy
Technology Laboratory, Pittsburgh PA 15236-0940  Noel Walkington  Center for
Nonlinear Analysis, Department of Mathematical Sciences, Carnegie Mellon
University  Matteo Pozzi  Department of Civil and Environmental Engineering,
Carnegie Mellon University  Kaushik Dayal<EMAIL_ADDRESS>Department of
Civil and Environmental Engineering, Carnegie Mellon University  Center for
Nonlinear Analysis, Department of Mathematical Sciences, Carnegie Mellon
University  Department of Materials Science and Engineering, Carnegie Mellon
University

###### Abstract

The modeling of coupled fluid transport and deformation in a porous medium is
essential to predict the various geomechanical process such as CO2
sequestration, hydraulic fracturing, and so on. Current applications of
interest, for instance, that include fracturing or damage of the solid phase,
require a nonlinear description of the large deformations that can occur. This
paper presents a variational energy-based continuum mechanics framework to
model large-deformation poroelasticity. The approach begins from the total
free energy density that is additively composed of the free energy of the
components. A variational procedure then provides the balance of momentum,
fluid transport balance, and pressure relations. A numerical approach based on
finite elements is applied to analyze the behavior of saturated and
unsaturated porous media using a nonlinear constitutive model for the solid
skeleton. Examples studied include the Terzaghi and Mandel problems; a gas-
liquid phase-changing fluid; multiple immiscible gases; and unsaturated
systems where we model injection of fluid into soil. The proposed variational
approach can potentially have advantages for numerical methods as well as for
combining with data-driven models in a Bayesian framework.

††preprint: To appear in International Journal for Numerical and Analytical
Methods in Geomechanics††preprint: DOI: https://doi.org/10.1002/nag.3326

##  1 Introduction

The modeling of fluid transport in porous deformable media, and analyzing the
coupled deformation-transport behavior, is of importance for a variety of
engineering and scientific applications, ranging from the modeling of natural
and engineered biological materials to geomechanics, e.g. [25, 54, 33, 50, 52,
21, 6, 55, 57, 39, 26, 24, 59, 43, 53, 58, 28, 34, 52, 31, 13] and numerous
other works. It has therefore been the focus of numerous formulations,
starting from the seminal work of Biot [10, 11]. In this paper, we formulate a
general finite-deformation model for poroelasticity based on a variational
approach.

The classical linear theory proposed by Biot has been applied to study a wide
range problems, e.g. [19, 60, 47, 20]. However, the many assumptions of the
linear theory, including linear elastic response for the solid phase,
infinitesimal strain, and constant permeability are too restrictive to model
the complex nonlinear response of many porous systems, such as polymers,
damaged geomaterials with fracture, and soft tissues.

Therefore, there have been many efforts to generalize the theory to the
nonlinear setting, with early works starting from Biot [12]. More recently,
examples of nonlinear large-deformation poroelastic models include models for
elastomeric and biological materials [41]; multiphase continuum formulations
presented in [27] and [32] for saturated fields reinforced by elastic fibers
and polymeric gels, respectively; and for computational geomechanics, such as
multiphase finite element formulations for fully and partially saturated
porous geomaterials, by Borja and coworkers [38, 14, 49]. Another older
approach is based on mixture theory, e.g. [17], which leads to a formidable
set of equations.

The key contribution of this paper is a variational derivation that starts
from a free energy, and uses variational principles to obtain the governing
equations, in contrast to other works in the literature. For instance, a
closely-related body of work is by Borja and coworkers [16, 14], in which they
use a thermodynamic approach based on an energy and also the balance laws of
continuum mechanics; they formulate a three-phase continuum mixture theory to
analyze a partially saturated porous media. They use a free energy function of
the medium is defined based on the solid strain, displacement vector, and
pressure of each phase. In a related approach, [36] formulated an energy-based
method for unsaturated porous medium where they used the energy equation which
is developed by [8].

Our approach is instead to directly use the energy and associated variational
principles to derive the governing equations. Similar approaches have been
used in geomechanics to model solid materials with a fluid phase, i.e. two-
phase porous medium models, e.g. [7, 56, 5]. Variational approaches for
poromechanics are relatively recent, and appear to begin with [29, 18]. While
we follow the overall structure of these works in our formulation, these have
been restricted to simpler settings, e.g. saturated systems, two-phase
systems, and so on. In this work, as in those works, the governing equations
are derived from energy minimization. However, [18] assume a multiplicative
decomposition of the deformation gradient into an elastic part and a
poromechanical part – corresponding to the deformation caused by fluid
transport – whereas we do not use this assumption.

This work presents a variational energy-based model for poromechanics, which
can be applied to both saturated and unsaturated porous mediums, containing
compressible constituents. The proposed model is amenable to the use of
arbitrary energy density functions for both the fluid and solid phases. In
particular, it can model unsaturated porous media with compressible
constituents, and multiple immiscible fluid phases with compressible and
incompressible constituents. Our strategy is to begin with the free energy
density function, and use a variational approach to obtain the conservation of
mass and momentum, and the pressure equality equations. An important advantage
of the energetic formulation and the minimization structure is the advantage
for robust numerical schemes using finite elements and other methods [51]. In
addition, it can enable efficient numerical approaches to combining physics-
based models with data-driven strategies [35].

#### Organization.

Section 2 provides the general formulation, starting with the kinematics and
then using a variational approach to develop the general governing equations.
Sections 3 discusses the constitutive models for saturated and unsaturated
systems, which includes a discussion on porous media with single and multiple
fluids. Section 4 examines the relation between the referential and current
volume fractions for different assumptions on the fluid and solid responses.
Section 5 provides numerical solutions of benchmark problems proposed by
Terzaghi and Mandel. Sections 6, 7, and 8 provides numerical examples to
investigate the settings of a phase-changing (gas to liquid) fluid, multiple
compressible fluids, and unsaturated systems.

##  2 Model Formulation

We assume that our system consists of $\displaystyle N$ phases or components
indexed by $\displaystyle i$, with $\displaystyle i=1$ corresponding to the
solid skeleton phase, and $\displaystyle i=2,\ldots,N$ corresponding to the
fluid phases. We follow the motion of the solid skeleton with respect to the
reference configuration, i.e., we use a Lagrangian description. Therefore, we
use material time derivatives. Once we obtain the governing equations, they
can be pushed forward to the current configuration.

The gradient and divergence operators in the reference and current
configurations will be denoted
$\displaystyle\nabla_{0},\nabla,\operatorname{div}_{0},\operatorname{div}$
respectively. In general, quantities with the subscript $\displaystyle 0$ will
refer to referential quantities. We use $\displaystyle\frac{\ \mathrm{d}\ }{\
\mathrm{d}t}$ to denote the material time derivative following the motion of
the solid skeleton.

###  2.A Kinematics

Figure 1: Representative volume element of the porous medium, and relation of
reference and current configurations.

We consider a porous medium $\displaystyle\Omega$, containing solid and fluid
components. We consider a two-scale picture motivated by homogenization.
Namely, there is a macroscopic deformation $\displaystyle{\mathbold
x}({\mathbold x}_{0},t)$, where $\displaystyle{\mathbold x}$ is the location
at time $\displaystyle t$ of a material particle located at
$\displaystyle{\mathbold x}_{0}$ in the reference configuration. To each
macroscale location, i.e., for a given $\displaystyle{\mathbold x}_{0}$, there
corresponds a microscale representative volume element (RVE), which accounts
for the geometry of the porous skeleton and other microscale features (Figure
1). Therefore, at the macroscopic homogenized, an RVE is equivalent to a
material point.

We define:

  1. 1.

the macroscale deformation gradient tensor $\displaystyle{\mathbold
F}:=\frac{\partial{\mathbold x}}{\partial{\mathbold x}_{0}}$, the Jacobian
$\displaystyle J:=\det{\mathbold F}$, and the displacement
$\displaystyle{\mathbold u}({\mathbold x}_{0},t):={\mathbold x}({\mathbold
x}_{0},t)-{\mathbold x}_{0}$;

  2. 2.

the densities $\displaystyle P_{i}$ and $\displaystyle P_{0i}$ as the mass of
component $\displaystyle i$ per unit (current and referential respectively)
volume of the mixture;

  3. 3.

the true densities $\displaystyle\rho_{i}$ and $\displaystyle\rho_{0i}$ as the
mass of component $\displaystyle i$ per unit (current and referential
respectively) volume occupied by component $\displaystyle i$;

  4. 4.

the volume fractions $\displaystyle\phi_{i}$ and $\displaystyle\phi_{0i}$ of
phase $\displaystyle i$ in the current and reference respectively.

We highlight that we will use a Lagrangian description for all of these
quantities. That is, while $\displaystyle P_{i},\rho_{i},\phi_{i}$ all
describe quantities in the current configuration, they will be written as
functions of $\displaystyle{\mathbold x}_{0}$.

It follows from these definitions that:

$\displaystyle\displaystyle P_{0i}$
$\displaystyle\displaystyle=\phi_{0i}\rho_{0i}$ (1a)
$\displaystyle\displaystyle P_{i}$
$\displaystyle\displaystyle=\phi_{i}\rho_{i}$ (1b)

A central assumption is that the entire RVE is mapped by
$\displaystyle{\mathbold F}$ as an affine transformation. Further, the mass of
each phase contained in an RVE is the same in the reference and current
configurations111We note, however, that the mass of the fluid phases in both
the reference and current configuration can change in time due to mass
transport.. Together, these imply the following relations between quantities
in the reference and current configurations:

$\displaystyle\displaystyle\text{Differential volume: }\
\mathrm{d}\Omega_{0}\mapsto\ \mathrm{d}\Omega=J\ \mathrm{d}\Omega_{0}$ (2a)
$\displaystyle\displaystyle\text{Density: }P_{0i}\mapsto P_{i}=J^{-1}P_{0i}$
(2b) $\displaystyle\displaystyle\text{True density:
}\rho_{0i}\mapsto\rho_{i}=\frac{\phi_{0i}}{\phi_{i}}J^{-1}\rho_{0i}$ (2c)

The relation between the referential and current versions of other quantities
such as the volume fraction will depend on assumptions such as
incompressibility, and will be discussed in later sections.

We notice that $\displaystyle\phi_{0i}$ is a fixed referential quantity that
depends only on the geometry of the reference RVE, while
$\displaystyle\phi_{i}$ depends on the deformation. Therefore, the latter is a
suitable quantity to appear as an argument of the total free energy below.
Further, $\displaystyle P_{0i}$ is not a fixed referential quantity, because
the mass of the fluid phases in the RVE can change due to transport.
Therefore, $\displaystyle P_{0i}$ is a suitable quantity to appear as an
argument of the total free energy.

###  2.B Energetics

The total free energy of the porous medium consisting of a solid skeleton and
$\displaystyle N-1$ fluid phases is given by:

$E\left[{\mathbold
x},\\{P_{0i}\\},\\{\phi_{i}\\}\right]=\int_{\Omega_{0}}\left(W_{0}\left(\nabla_{0}{\mathbold
x},\\{P_{0i}\\},\\{\phi_{i}\\}\right)-\psi_{0}({\mathbold
x},\\{P_{0i}\\})\right)\ \mathrm{d}\Omega_{0}-\int_{\partial A}{\mathbold
t}_{0}^{*}\cdot{\mathbold x}\ \mathrm{d}S-\sum_{i=2}^{N}\int_{\partial
B_{i}}\eta^{*}_{0i}P_{0i}\ \mathrm{d}S$ (3)

where $\displaystyle W_{0}$ is the (Helmholtz) energy per unit referential
volume, and $\displaystyle\psi_{0}$ is the force potential per unit
referential volume. The boundary conditions correspond to specified mechanical
traction $\displaystyle{\mathbold t}^{*}$ on the solid skeleton over some part
of the boundary $\displaystyle\partial A$, and specified referential chemical
potential $\displaystyle\eta^{*}_{0i}$ over some part of the boundary
$\displaystyle\partial B_{i}$ for each fluid phase.

We will assume that the energy density $\displaystyle W_{0}$ is additively
composed of the energy densities of the components:

$W_{0}=\phi_{0s}W_{0s}({\mathbold
F})+\sum_{i=2}^{N}\phi_{0i}W_{0i}(J,P_{0i},\phi_{i})$ (4)

where $\displaystyle W_{0i}$ is the energy density of phase $\displaystyle i$
per unit referential volume.

For the solid phase, it is standard to assume that $\displaystyle W_{0_{1}}$
is only a function of $\displaystyle{\mathbold F}$, and that $\displaystyle
P_{0_{1}}$ and $\displaystyle\phi_{1}$ will not appear.

For the fluid phases, we are typically given a free energy density (per unit
current volume) that is a function only of the true density of the
corresponding phase. Denoting this given free energy density as $\displaystyle
W_{i}(\rho_{i})$, we write the energy density $\displaystyle W_{0i}$ in terms
of the arguments of $\displaystyle E$:

$W_{0i}=\frac{J\phi_{i}}{\phi_{0i}}W_{i}\left(J^{-1}P_{0i}/\phi_{i}\right)$
(5)

where we have used that $\displaystyle\rho_{i}=J^{-1}P_{0i}/\phi_{i}$ from (1)
and (2).

###  2.C Balance of Mass

Following [29, 19], we define the (referential) chemical potential
$\displaystyle\eta_{0}$ as the variation of $\displaystyle E$ with respect to
the referential density for the fluid phases in the bulk:

$\eta_{0i}=-\delta_{P_{0i}}E=-\frac{\partial W_{0}}{\partial
P_{0i}}+\frac{\partial\psi_{0}}{\partial P_{0i}}$ (6)

where we use the notation $\displaystyle\delta_{P_{0i}}E$ to denote the
variational (functional) derivative of $\displaystyle E$ with respect to
$\displaystyle P_{0i}$. The boundary condition that appears from the
variational principle corresponds to $\displaystyle\eta_{0i}=\eta^{*}_{0i}$
specified on $\displaystyle\partial B_{i}$; this is related, but not directly
equivalent, to specifying the pressure in the general setting.

The velocity in the reference configuration – relative to the solid skeleton –
of the $\displaystyle i$-th fluid phase is assumed to be proportional to the
gradient of the chemical potential:

${\mathbold v}_{0i}={\mathbold K}_{i}\nabla_{0}\eta_{0i}={\mathbold
K}_{i}\left(-\nabla_{0}\frac{\partial W_{0}}{\partial
P_{0i}}+\nabla_{0}\frac{\partial\psi_{0}}{\partial P_{0i}}\right)$ (7)

where $\displaystyle{\mathbold K}_{i}$ is the positive-definite second-order
permeability tensor of the porous medium in the reference configuration.
Following [7], we can push forward (7) as follows:

${\mathbold v}_{i}=J^{-1}{\mathbold F}{\mathbold v}_{0i}={\mathbold
k}_{i}\left(-\nabla\frac{\partial W_{0}}{\partial
P_{0i}}+\nabla\frac{\partial\psi_{0}}{\partial P_{0i}}\right)$ (8)

where $\displaystyle{\mathbold k}_{i}$ is the permeability tensor of the
porous medium in the current configuration, $\displaystyle{\mathbold
K}_{i}=J{\mathbold F}^{-1}{\mathbold k}_{i}{\mathbold F}^{-T}$ [38], using
$\displaystyle\nabla_{0}={\mathbold F}^{T}\nabla$. In an isotropic medium, the
permeability tensor can be written:

${\mathbold k}_{i}=\frac{\tilde{k}_{i}}{g}{\mathbold I}$ (9)

where $\displaystyle g$ is the gravitational acceleration, and
$\displaystyle\tilde{k}_{i}$ is the hydraulic conductivity, which is defined
as $\displaystyle\tilde{k}_{i}=\frac{\kappa}{\gamma_{i}}\rho_{i}g$.
$\displaystyle\kappa$ is the true permeability of the solid;
$\displaystyle\gamma_{i}$ is the dynamic viscosity of the fluid;
$\displaystyle\rho_{i}$ is the fluid density; and $\displaystyle{\mathbold I}$
is the second-order identity tensor.

Therefore, the fluid flux vector for the $\displaystyle i$-th phase is:

${\mathbold q}_{i}=P_{i}{\mathbold v}_{i}=-{\mathbold
k}_{i}\left(P_{i}\nabla\frac{\partial W_{0}}{\partial
P_{0i}}-P_{i}\nabla\frac{\partial\psi_{0}}{\partial P_{0i}}\right)$ (10)

The conservation of mass for each fluid phase in the current configuration can
be written as follows:

$-\int_{\partial\Omega}{\mathbold q}_{i}\cdot{\mathbold n}\
\mathrm{d}S=\frac{\ \mathrm{d}\ }{\ \mathrm{d}t}\left(\int_{\Omega}P_{i}\
\mathrm{d}V\right)\Rightarrow-J\operatorname{div}{\mathbold q}_{i}=\frac{\
\mathrm{d}\ }{\ \mathrm{d}t}P_{i}J$ (11)

where the differential form follows from the integral form by arbitrariness of
$\displaystyle\Omega$.

Specializing to a uniform gravitational field, the referential potential can
be written as $\displaystyle\psi_{0}={\mathbold b}_{0}\cdot{\mathbold x}$,
where $\displaystyle{\mathbold b}_{0}$ is the body force per unit referential
volume. This can be written in terms of the densities and volume fractions of
the individual phases:

${\mathbold b}_{0}:=P_{0}{\mathbold
g}=\left(\sum_{i=1}^{N}P_{0i}\right){\mathbold g}$ (12)

where $\displaystyle P_{0}$ is the total mass per unit referential volume, and
$\displaystyle{\mathbold g}$ is the force per unit mass due to gravity. Then,
we can simplify the potential term in (10) to:

$\frac{\partial\psi_{0}}{\partial P_{0i}}={\mathbold g}\cdot{\mathbold
x}\Rightarrow\nabla\frac{\partial\psi_{0}}{\partial P_{0i}}={\mathbold g}$
(13)

Consequently, the fluid flux vector in (10)) can be rewritten as follows:

${\mathbold q}_{i}=-{\mathbold k}_{i}\left(P_{i}\nabla\frac{\partial
W_{0}}{\partial P_{0i}}-P_{i}{\mathbold g}\right)$ (14)

###  2.D Balance of Momentum

Setting to zero the variation of $\displaystyle E$ with respect to
$\displaystyle{\mathbold x}$ yields the force equilibrium equation in the
referential and current forms:

$\operatorname{div}_{0}{\mathbold
T}+\frac{\partial\psi_{0}}{\partial{\mathbold x}}=\mathbold{$\displaystyle
0$}\Leftrightarrow\operatorname{div}\mathbold{\sigma}+J^{-1}\frac{\partial\psi_{0}}{\partial{\mathbold
x}}=\bf 0$ (15)

where $\displaystyle{\mathbold T}:=\frac{\partial W_{0}}{\partial{\mathbold
F}}$ is the first Piola stress tensor;
$\displaystyle\mathbold{\sigma}=\frac{1}{J}{\mathbold T}{\mathbold F}^{T}$ is
the Cauchy stress tensor; the body force is defined through the derivative of
$\displaystyle\frac{\partial\psi_{0}}{\partial{\mathbold x}}$; and the
corresponding boundary conditions from the variation are
$\displaystyle{\mathbold T}\hat{\mathbold n}_{0}={\mathbold t}_{0}^{*}$.

If we specialize to a uniform gravitational field, the body force in the
reference $\displaystyle{\mathbold b}_{0}$ is given by the expression in (12).
In the current configuration, we have

${\mathbold b}:=P{\mathbold g}=\left(\sum^{N}_{i=1}P_{i}\right){\mathbold g}$
(16)

where $\displaystyle P$ is the total mass per unit current volume. We note
that $\displaystyle{\mathbold b}_{0}$ and $\displaystyle{\mathbold b}$ are
related by the expression $\displaystyle{\mathbold b}=J^{-1}{\mathbold
b}_{0}$.

###  2.E Balance of Fluid Pressure

We first notice that the variables $\displaystyle\phi_{i}$ are not independent
when we take variations; they must satisfy the constraint
$\displaystyle\sum^{N}_{i=1}\phi_{i}=1$. We can alternatively write this as:

$\sum^{N}_{i=1}\phi_{i}=1\Rightarrow\phi_{2}=1-\phi_{s}-\phi_{3}-\cdots-\phi_{i}-\cdots-\phi_{N}$
(17)

We will eliminate $\displaystyle\phi_{2}$ in terms of the other
$\displaystyle\phi_{i}$. Further, we notice that $\displaystyle\phi_{s}$ will
be prescribed; similarly, some of the $\displaystyle\phi_{i}$ could be
prescribed, e.g., due to incompressibility.

Setting to zero the variation of $\displaystyle E$ with respect to
$\displaystyle\phi_{i}$, and using the form from (4), gives:

$\frac{\partial W_{0}}{\partial\phi_{i}}=\phi_{02}\frac{\partial
W_{02}}{\partial\phi_{2}}\frac{\partial\phi_{2}}{\partial\phi_{i}}+\phi_{0i}\frac{\partial
W_{0i}}{\partial\phi_{i}}=0\Rightarrow-\phi_{02}\frac{\partial
W_{02}}{\partial\phi_{2}}+\phi_{0i}\frac{\partial W_{0i}}{\partial\phi_{i}}=0$
(18)

Since this is true for every $\displaystyle i$, it follows that:

$\phi_{02}\frac{\partial W_{02}}{\partial\phi_{2}}=\phi_{03}\frac{\partial
W_{03}}{\partial\phi_{3}}=\cdots=\phi_{0N}\frac{\partial
W_{0N}}{\partial\phi_{N}}$ (19)

We next recall that the derivative of the Helmholtz free energy with respect
to volume, keeping temperature and mass fixed, is the pressure $\displaystyle
p$. Rewriting this in terms of the density $\displaystyle\rho$ – which is
inversely proportional to the volume when the mass is fixed – and in terms of
the Helmholtz free energy density per unit volume, we have the following
relation between the fluid pressure $\displaystyle p$ and the Helmholtz free
energy density:

$p=W(\rho)-\rho\frac{\ \mathrm{d}W}{\ \mathrm{d}\rho}$ (20)

This relation holds for each fluid phase individually.

Applying the equilibrium relations above to the form assumed in (5), we get
that:

$\phi_{0i}\frac{\partial
W_{0i}}{\partial\phi_{i}}=J\left(W_{i}\left(J^{-1}P_{0i}/\phi_{i}\right)-\frac{J^{-1}P_{0i}}{\phi_{i}}\frac{\
\mathrm{d}W_{i}}{\ \mathrm{d}\rho_{i}}\right)=Jp_{i}$ (21)

Since $\displaystyle J$ depends only on $\displaystyle{\mathbold x}_{0}$, it
follows that the pressures in each phase $\displaystyle p_{i}$ at a given
$\displaystyle{\mathbold x}_{0}$ must be equal.

We highlight that the balance of pressure provides a local condition that must
be satisfied pointwise and does not require the solution of PDE, in contrast
to the balances of mass and momentum.

###  2.F Compatibility with the Second Law of Thermodynamics

To ensure that the proposed model is consistent with the second law of
thermodynamics, we first notice that in an isothermal setting it reduces to
the requirement that the dissipation is non-negative for any process [30].

Following [46, 4, 40, 22, 3, 1, 2], we begin by taking the time derivative of
the free energy from (3) to get:

$\frac{\ \mathrm{d}E}{\
\mathrm{d}t}=\int_{\Omega_{0}}\sum_{i}\left(\frac{\partial W_{0}}{\partial
P_{0i}}-\frac{\partial\psi_{0}}{\partial P_{0i}}\right)\frac{\
\mathrm{d}P_{0i}}{\ \mathrm{d}t}\ \mathrm{d}\Omega_{0}$ (22)

where we have used that $\displaystyle\frac{\ \mathrm{d}\phi_{i}}{\
\mathrm{d}t}=0$ and $\displaystyle\frac{\ \mathrm{d}{\mathbold x}}{\
\mathrm{d}t}=\bf 0$ per our quasistatic model. We have assumed that the
boundary terms vanish for simplicity, but they can easily be accounted for
following, e.g., [46, 1].

From Section 2.C, we can push back the balance of mass to the reference
configuration to obtain:

$-J\operatorname{div}{\mathbold q}_{i}=\frac{\ \mathrm{d}}{\
\mathrm{d}t}P_{i}J\Rightarrow-\operatorname{div}_{0}{\mathbold q}_{0i}=\frac{\
\mathrm{d}}{\ \mathrm{d}t}P_{0i}$ (23)

and the corresponding flux for each fluid phase is given by:

${\mathbold q}_{0i}=-P_{i}{\mathbold K}_{i}\nabla_{0}\left(\frac{\partial
W_{0}}{\partial P_{0i}}-\frac{\partial\psi_{0}}{\partial P_{0i}}\right)$ (24)

Using these expressions, and imposing the condition that
$\displaystyle{\mathbold K}$ is a positive-definite tensor, we can rewrite
$\displaystyle\frac{\ \mathrm{d}E}{\ \mathrm{d}t}$ as follows and use
integration-by-parts:

$\begin{split}\frac{\ \mathrm{d}E}{\
\mathrm{d}t}&=\int_{\Omega_{0}}\sum_{i}\left(\frac{\partial W_{0}}{\partial
P_{0i}}-\frac{\partial\psi_{0}}{\partial
P_{0i}}\right)\operatorname{div}_{0}\left[P_{i}{\mathbold
K}_{i}\nabla_{0}\left(\frac{\partial W_{0}}{\partial
P_{0i}}-\frac{\partial\psi_{0}}{\partial P_{0i}}\right)\right]\
\mathrm{d}\Omega_{0}\\\
&=-\int_{\Omega_{0}}\sum_{i}P_{i}\nabla_{0}\left(\frac{\partial
W_{0}}{\partial P_{0i}}-\frac{\partial\psi_{0}}{\partial
P_{0i}}\right)\cdot{\mathbold K}_{i}\nabla_{0}\left(\frac{\partial
W_{0}}{\partial P_{0i}}-\frac{\partial\psi_{0}}{\partial P_{0i}}\right)\
\mathrm{d}\Omega_{0}\\\ &\leq 0\end{split}$ (25)

which shows that $\displaystyle E$ is non-increasing in the proposed model, as
required by the second law of thermodynamics.

##  3 Models for Saturated and Unsaturated Systems

###### Single Compressible Fluid in a Saturated Medium.

The energy density is:

$W_{0}({\mathbold F},P_{0f},\phi_{f})=(1-\phi_{0f})W_{0s}({\mathbold
F})+\phi_{0f}W_{0f}(J,P_{0f},\phi_{f})$ (26)

where we have used $\displaystyle N=2$ (solid $\displaystyle s$ and fluid
$\displaystyle f$). ∎

###### Multiple Immiscible Compressible Fluids.

The energy density is:

$W_{0}({\mathbold
F},\\{P_{0i}\\}^{N}_{i=2},\\{\phi_{i}\\}^{N}_{i=2})=\phi_{0s}W_{0s}({\mathbold
F})+\sum_{i=2}^{N}\phi_{0i}W_{0i}(J,P_{0i},\phi_{i})$ (27)

where the subscript $\displaystyle s$ denotes quantities corresponding to the
solid. ∎

###### Unsaturated Systems.

Consider an unsaturated medium with a solid phase, a incompressible fluid
phase (e.g. water), and an compressible fluid phase (e.g. air). The state
variables for the energy are $\displaystyle{\mathbold
F},P_{0I},P_{0c},\phi_{I},\phi_{c}$, where the subscripts $\displaystyle c$
and $\displaystyle I$ denote quantities corresponding to the compressible
fluid and incompressible fluid respectively.

For the incompressible fluid phase, we require that the true density in the
current configuration has a fixed value $\displaystyle\tilde{\rho}_{I}$:

$\frac{P_{0I}}{J\phi_{I}}=\tilde{\rho}_{I}$ (28)

We notice that the energy is constant for an incompressible fluid: all states
that satisfy incompressibility have the same energy, and all other states are
inadmissible (i.e., have infinite energy). Since this energy is a constant, it
will not appear once we take variations and consequently we can set it to
$\displaystyle 0$ for simplicity.

Therefore, we can write the total energy density as:

$W_{0}({\mathbold
F},P_{0I},P_{0c},\phi_{I},\phi_{c})=\phi_{0s}W_{0s}({\mathbold
F})+\phi_{0c}W_{0c}(J,\phi_{c},P_{0c})$ (29)

We highlight some issues. First, if we further assume that the solid skeleton
is incompressible and affinely deformed by $\displaystyle{\mathbold F}$, this
imposes the additional constraint $\displaystyle J=1$. Second, we notice that
the incompressibility constraint is generally not satisfied in the reference
configuration, i.e. $\displaystyle\rho_{0I}\neq\tilde{\rho}_{I}$ in general.
Third, it is sometimes reasonable to neglect the energy of the compressible
fluid (e.g., air at low pressure in rock); in that case, we need only consider
the energy of the solid. However, we emphasize that the problem still contains
the influence of the fluid through the mass transport equations (11). ∎

##  4 Models for Referential and Current Volume Fractions

###### Incompressible Solid with Compressible Pore Space.

Figure 2: Volume change for an incompressible solid phase with a compressible
fluid.

Consider a porous medium with an incompressible solid skeleton and
compressible fluid in the pore space (Figure 2). From incompressibility of the
solid skeleton – meaning that the true current density $\displaystyle\rho_{s}$
has a given fixed value – and the fact that the mass of the solid skeleton is
conserved because there is no mass transport of the solid, we have that the
volume occupied by the solid phase in the reference and current is constant.
Therefore, the current and referential volume fractions of the solid phase are
related by $\displaystyle\phi_{s}=J^{-1}\phi_{0s}$. Using that the fluid
volume fractions in the reference and current are given by
$\displaystyle\phi_{0f}=1-\phi_{0s}$ and $\displaystyle\phi_{f}=1-\phi_{s}$,
we have that $\displaystyle\phi_{f}=1-J^{-1}(1-\phi_{0f})$. ∎

###### Affinely Deformed Solid with Compressible Fluid.

Figure 3: Volume change for a compressible solid phase with a compressible
fluid.

Consider a porous medium with a compressible solid phase and a compressible
fluid phase (Figure 3). Let $\displaystyle\phi_{0f}$ and
$\displaystyle\phi_{0s}$ be the referential volume fractions of fluid and
solid phases, respectively. The deformation of the solid phase is assumed to
be affine with the macroscopic deformation gradient $\displaystyle{\mathbold
F}$. This implies that $\displaystyle\phi_{0s}=\phi_{s}$ as follows. We have
the following relation between current and reference volume elements for the
entire RVE as well as for the solid skeleton:

$\displaystyle\displaystyle\ \mathrm{d}\Omega=J\ \mathrm{d}\Omega_{0}$ (30)
$\displaystyle\displaystyle\ \mathrm{d}\Omega_{s}=J\ \mathrm{d}\Omega_{0s}$

Further, from the definition of the volume fraction of the solid skeleton, we
have:

$\displaystyle\displaystyle\ \mathrm{d}\Omega_{s}=\phi_{s}\
\mathrm{d}\Omega=\phi_{s}J\ \mathrm{d}\Omega_{0}$ (31)
$\displaystyle\displaystyle\ \mathrm{d}\Omega_{0s}=\phi_{0s}\
\mathrm{d}\Omega_{0}$

therefore, by substituting last equation in second equation, we can derive:

$\ \mathrm{d}\Omega_{0}=J\phi_{0s}\ \mathrm{d}\Omega_{0}$ (32)

We can therefore conclude that $\displaystyle\phi_{0s}=\phi_{s}$.

Further, since $\displaystyle\phi_{f}=1-\phi_{s}$ and
$\displaystyle\phi_{0f}=1-\phi_{0s}$, it follows that
$\displaystyle\phi_{0f}=\phi_{f}$. ∎

###### Compressible Solid with Compressible and Incompressible Fluids.

Figure 4: Volume change for a compressible solid phase with incompressible and
compressible fluid phases.

Consider a porous medium with compressible solid skeleton and air and
incompressible water (Figure 4). Let $\displaystyle\phi_{0s}$,
$\displaystyle\phi_{0I}$, and $\displaystyle\phi_{0c}$ be the reference volume
fractions of compressible solid, incompressible fluid, and compressible fluid
phases, respectively. Assuming that the solid phase deforms affinely, we have
that $\displaystyle\phi_{s}=\phi_{0s}$. For the incompressible fluid phase, we
require that the true density in the current configuration have a fixed value
$\displaystyle\tilde{\rho}_{I}$, giving
$\displaystyle\phi_{I}=\frac{P_{0I}}{J\tilde{\rho}_{I}}$. Using that the
volume fractions must sum to $\displaystyle 1$, we have
$\displaystyle\phi_{c}=1-\phi_{0s}-\frac{P_{0I}}{J\tilde{\rho}_{I}}$. ∎

##  5 Model Verification with Terzaghi’s and Mandel’s Problems

In this section, we verify the model by investigating its response in the
model consolidation problems proposed by Terzaghi and Mandel. The general form
of the strain energy density for a saturated porous medium with an
incompressible solid phase and compressible fluid phase is:

$W_{0}=(1-\phi_{0f})W_{0s}({\mathbold F})+\phi_{0f}W_{0f}(J,P_{0f},\phi_{f})$
(33)

Following section 2.D, the first Piola stress tensor can be derived as:

${\mathbold T}=(1-\phi_{0f})\frac{\partial W_{0s}}{\partial{\mathbold
F}}+\phi_{0f}\frac{\partial W_{0f}}{\partial J}J{\mathbold F}^{-T}$ (34)

Using (5) and (21), we can write:

$\frac{\partial W_{0f}}{\partial J}=-p\frac{\phi_{f}}{\phi_{0f}}$ (35)

Thus, the Cauchy stress tensor can be written as follows:

$\mathbold{\sigma}=(1-\phi_{0f})J^{-1}{\mathbold F}^{T}\frac{\partial
W_{0s}}{\partial{\mathbold F}}-\phi_{f}p{\mathbold I}$ (36)

Equation (36) shows that the total Cauchy stress is obtained by a summation of
the solid and fluid partial stresses
$\displaystyle\mathbold{\sigma}=\mathbold{\sigma}_{s}+\mathbold{\sigma}_{f}$,
where $\displaystyle\mathbold{\sigma}_{f}=-\phi_{f}p{\mathbold I}$ is an
isotropic stress tensor. In addition, the general expression of the Cauchy
effective stress tensor is
$\displaystyle\mathbold{\sigma}^{\prime}=\mathbold{\sigma}+p{\mathbold I}$,
which can be written as:

$\mathbold{\sigma}^{\prime}=(1-\phi_{f})\left({\mathbold F}^{T}\frac{\partial
W_{0s}}{\partial{\mathbold F}}+p{\mathbold I}\right)$ (37)

Assuming that the deformation of the solid phase is small and that the
response is isotropic, we linearize it to get:

$\mathbold{\sigma}=(1-\phi_{f})\left(2\mu\mathbold{\epsilon}+\lambda\operatorname{tr}(\mathbold{\epsilon}){\mathbold
I}\right)-\phi_{f}p{\mathbold I}$ (38)

where $\displaystyle\operatorname{tr}$ denotes the trace operator. We then
define the renormalized Lame constants
$\displaystyle\tilde{\mu}=(1-\phi_{0f})\mu$ and
$\displaystyle\tilde{\lambda}=(1-\phi_{0f})\lambda$. We define
$\displaystyle\tilde{p}=\phi_{f}p$ and $\displaystyle p$ as the total pressure
and the intrinsic fluid pressure, respectively.

Following section 2.C, the conservation of mass for the fluid phase can be
written as:

$\frac{\ \mathrm{d}}{\ \mathrm{d}t}P_{f}J=-J\operatorname{div}{\mathbold q}$
(39)

Following [15], the bulk modulus of the fluid phase can be written:

$K_{f}=\rho_{f}\frac{\ \mathrm{d}{p}}{\ \mathrm{d}\rho_{f}}$ (40)

We recall further that:

$P_{f}=\phi_{f}\rho_{f}$ (41)

Using (40) and (41), we can rewrite the conservation of mass for the fluid and
solid phases respectively as:

$\displaystyle\displaystyle\frac{\ \mathrm{d}J}{\
\mathrm{d}t}\frac{\phi_{f}}{J}+\frac{\ \mathrm{d}\phi_{f}}{\
\mathrm{d}t}+\frac{\ \mathrm{d}p}{\
\mathrm{d}t}\frac{\phi_{f}}{K_{f}}=-\frac{1}{\rho_{f}}\operatorname{div}{\mathbold
q}$ (42) $\displaystyle\displaystyle\frac{\ \mathrm{d}J}{\
\mathrm{d}t}\frac{\phi_{s}}{J}+\frac{\ \mathrm{d}\phi_{s}}{\ \mathrm{d}t}=0$
(43)

The second equation, for the solid phase, is derived in the same way as the
equation (42) for the fluid phase, but also using that there is no transport
of the solid phase and that the solid phase is incompressible (the bulk
modulus $\displaystyle K_{s}\to\infty$).

Adding (42) and (43) yields the balance of mass of the porous system as
follows:

$\frac{\ \mathrm{d}J}{\ \mathrm{d}t}\frac{1}{J}+\frac{\ \mathrm{d}p}{\
\mathrm{d}t}\frac{\phi_{f}}{K_{f}}=-\frac{1}{\rho_{f}}\operatorname{div}{\mathbold
q}$ (44)

Ignoring the effect of gravity and using (2) and (6), it follows that:

$\displaystyle\displaystyle\eta_{0}=\frac{\partial W_{f}(\rho)}{\partial\rho}$
(45) $\displaystyle\displaystyle\frac{\partial p}{\partial\eta_{0}}=-P_{f}$
(46)

Using these equations and considering small deformations, (44) can be written
as:

$\frac{\ \mathrm{d}\operatorname{tr}(\mathbold{\epsilon})}{\
\mathrm{d}t}+\frac{\ \mathrm{d}p}{\
\mathrm{d}t}\frac{\phi_{f}}{K_{f}}=\frac{\kappa}{\gamma}\nabla^{2}{p}$ (47)

###  5.A Terzaghi’s one-dimensional consolidation

In this section, we investigate the one-dimensional Terzaghi’s problem. The
geometry of the sample is shown in Fig. 5, and the numerical values of the
model parameters are listed in Table 1. We consider a fully saturated sample
with porosity $\displaystyle\phi_{f}=$\displaystyle 0.375$$ that is loaded
instantaneously at $\displaystyle t=$\displaystyle 0$$ by a vertical load
$\displaystyle w=$\displaystyle 10\text{\,}\mathrm{kPa}$$ at the top boundary,
and then remains constant in time. This sets up a mechanical stress in the
sample.

The fluid flux is allowed to drain out of the top surface and is zero on the
bottom surface. The pressure boundary conditions of are:

$\text{ Top BC: }\tilde{p}(z=0)=0;\quad\text{ Bottom BC:
}\left.\frac{\partial\tilde{p}}{\partial z}\right|_{z=h}=0$ (48)

The closed-form solution to this problem can be written as [19]:

$\displaystyle\displaystyle\bar{p}(\bar{z},\bar{t})=\sum_{n=0}^{\infty}\frac{4}{\pi(2n+1)}\sin\left(\frac{(2n+1)\pi}{2}\bar{z}\right)\exp\left(-\frac{(2n+1)^{2}\pi^{2}}{4}\bar{t}\right)$
(49)
$\displaystyle\displaystyle\tilde{p}=\frac{wK_{f}}{K_{f}+\phi_{f}(\tilde{\lambda}+2\tilde{\mu})}\bar{p};\quad
z=h\bar{z};\quad
t=\frac{h^{2}\kappa}{\gamma}\left(\frac{K_{f}(\tilde{\lambda}+2\tilde{\mu})}{K_{f}+\phi_{f}(\tilde{\lambda}+2\tilde{\mu})}\right)\bar{t}$
(50)

Figure 5: Geometry of Terzaghi’s problem. Property | Value
---|---
Solid phase Lame constant, $\displaystyle\tilde{\lambda}$ | $\displaystyle 40\text{\,}\mathrm{MPa}$
Solid phase Lame constant, $\displaystyle\tilde{\mu}$ | $\displaystyle 40\text{\,}\mathrm{MPa}$
Hydraulic conductivity, $\displaystyle\tilde{k}$ | $\displaystyle 10\text{\times}{10}^{-8}\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$
Fluid bulk modulus, $\displaystyle K_{f}$ | $\displaystyle 2270\text{\,}\mathrm{MPa}$
Table 1: Properties of solid and fluid phases.

Fig. 6 shows the comparison of the closed-form and numerical results at
different times. Due to the compressibility of the pore space, the fluid
drains gradually from the top boundary of specimen; therefore the pore
pressure value decreases and the effective stress increases as time increases.
As time increases, the pore pressure tends to $\displaystyle 0$ and the solid
partial stress tends to the external load.

Figure 6: Variation of pressure with height for Terzaghi’s problem.

###  5.B Mandel’s two-dimensional consolidation

In this section, we discuss a two-dimensional consolidation of a porous
sample. The geometry of the specimen is shown in Fig. 7, and the material
properties and soil porosity are assumed the same as the previous example. We
consider a specimen with infinite length in the $\displaystyle z$ direction
which is between two horizontal plates (the top and bottom boundaries).

A constant load $\displaystyle w=$\displaystyle 10\text{\,}\mathrm{kPa}$$ is
applied in the vertical direction at $\displaystyle t=$\displaystyle 0$$, and
remains constant thereafter. The left and right boundaries are traction-free.

For the fluid flow, we assume that there is no flux at the top and bottom
boundaries, and the pressure is $\displaystyle 0$ on the left and right
boundaries (located at $\displaystyle x=\pm a/2$).

The closed-form solution for the pressure distribution is [19]:

$\displaystyle\displaystyle\bar{p}(\bar{x},\bar{t})=2\sum_{n=1}^{\infty}\frac{\cos(\alpha_{n}\bar{x})-\cos\alpha_{n}}{\alpha_{n}-\sin\alpha_{n}\cos\alpha_{n}}\sin\alpha_{n}\exp(-\alpha_{n}^{2}\bar{t})$
(51)
$\displaystyle\displaystyle\tilde{p}=\frac{1}{3}B(1+\nu_{u})w\bar{p};\quad
x=a\bar{x};\quad
t=\frac{a^{2}\kappa}{\gamma}\left(\frac{K_{f}(\tilde{\lambda}+2\tilde{\mu})}{K_{f}+\phi_{f}(\tilde{\lambda}+2\tilde{\mu})}\right)\bar{t}$
(52)

where $\displaystyle\alpha_{n}$ is the solution of the following equation:

$\frac{\tan\alpha_{n}}{\alpha_{n}}=\frac{1-\nu}{\nu_{u}-\nu};\quad\text{ where
}\nu=\frac{\tilde{\lambda}}{2(\tilde{\lambda}+\tilde{\mu})};\quad\nu_{u}=\frac{3\nu+B(1-2\nu)}{3-B(1-2\nu)};\quad
B=\frac{K_{f}}{\phi_{f}K_{u}};\quad
K_{u}=\lambda+\frac{2\tilde{\mu}}{3}+\frac{K_{f}}{\phi_{f}}$ (53)

Figure 7: Geometry and the boundary conditions for Mandel’s problem.

Once the load is applied, the fluid pressure decreases at the vicinity of left
and right boundaries due to drainage out of the sides. Due to relatively slow
drainage, the pressure decrease cannot be observed immediately in the whole
domain, while the total load remains constant. Therefore, the pore pressure
increases beyond the initial value at the center of specimen, which is the
Mandel effect. The closed-form solution and our numerical results are shown in
Fig. 8 for various times; we notice a signature of the Mandel effect at
$\displaystyle t=$\displaystyle 0.5\text{\,}\mathrm{s}$$ wherein we have the
maximum value of pore pressure in the central region of the specimen.
Similarly, Fig. 9 shows the pore pressure distribution at $\displaystyle
t=$\displaystyle 1.0\text{\,}\mathrm{s}$$. We see clearly the increase of pore
pressure beyond the initial pressure at the center, and it goes to
$\displaystyle 0$ on the right and left boundaries.

Figure 8: Variation of pore pressure in the horizontal direction at different
times. Figure 9: Pressure distribution for Mandel’s problem at $\displaystyle
t=$\displaystyle 1\text{\,}\mathrm{s}$$.

The variation of non-dimensional pore pressure ($\displaystyle\bar{p}$) at
different locations (from Fig. 7) vs. logarithmic time is plotted in Fig. 10.
The log scale makes more prominent the pore pressure increase above its
initial value shortly after applying the vertical load due to the Mandel
effect. Also, the comparison of these three graphs at points $\displaystyle
A$, $\displaystyle B$, and $\displaystyle C$, clearly shows that the specimen
experiences the highest amount of early times pressure increase at the central
points (point $\displaystyle A$), and the Mandel effect vanishes at the
vicinity of left and right boundaries.

Figure 10: Pore pressure vs. time at different distances from the center of
specimen.

##  6 Van der Waals compressible fluid

In this section, we examine a porous medium containing a single fluid phase
that is modeled as a Van der Waals (VdW) gas, motivated by experimental
approaches such as [37]. A key feature of the VdW model is that it allows for
a change from gas to liquid, and provides a simple model for gases such as
carbon dioxide (CO2). An important challenge in simulating gas injection in
porous media is to predict the gas-liquid phase transformation of the injected
gas under high pressure in the vicinity of the injection point. The VdW model
for the fluid phase enables us to investigate this complex behavior.

The Helmholtz free energy density per unit current volume of a VdW gas is
given by the expression:

$W_{f}(\rho)=c\rho RT\left(1-\log(cRT)\right)-\rho
RT\log\left(\rho-b\right)-a\rho^{2}$ (54)

where $\displaystyle c$, $\displaystyle a$, $\displaystyle b$, $\displaystyle
R$ and $\displaystyle S_{0}$ are constants, and $\displaystyle T$ is the
temperature [42].

Following (5), we rewrite this as:

$W_{0f}(J,P_{0},\phi)=\frac{J\phi}{\phi_{0}}W\left(J^{-1}P_{0}/\phi\right)=c\frac{P_{0}}{\phi}RT\left(1-\log(cRT)\right)-\frac{P_{0}}{\phi}RT\log\left(\frac{\phi
J}{P_{0}}-b\right)-a\frac{P_{0}^{2}}{\phi^{2}J}$ (55)

where we have used that $\displaystyle\phi=\phi_{0}$ from Section 4 (affinely
deformed solid with compressible fluid).

We use a Neo-Hookean model for the solid phase:

$W_{0s}({\mathbold F})=\frac{\mu}{2}(\operatorname{tr}({\mathbold
F}^{T}{\mathbold F})-2)-\mu\log J+\frac{\lambda}{2}(\log J)^{2}$ (56)

Using (6) and (10), we find the chemical potential and fluid flux:

$\displaystyle\displaystyle\eta_{0}$
$\displaystyle\displaystyle=\left(cRT\left(1-\log(cRT)\right)-RT\log\left(\frac{1}{\rho}-b\right)+\frac{RT}{1-b\rho}-2a\rho\right)$
(57) $\displaystyle\displaystyle{\mathbold q}$
$\displaystyle\displaystyle=-\phi{\mathbold
k}\left(\frac{RT\nabla\rho}{1-b\rho}+\frac{bRT\rho\nabla\rho}{(1-b\rho)^{2}}-2a\rho\nabla\rho\right)$
(58)

where $\displaystyle\rho=\frac{P}{\phi}$ is the true current density of the
fluid.

The first Piola stress can be obtained from the energy density:

${\mathbold T}=(1-\phi)\left(\mu{\mathbold F}-\mu{\mathbold
F}^{-T}+\lambda\log{J}{\mathbold F}^{-T}\right)-P_{0}\left(\frac{RT{\mathbold
F}^{-T}}{1-b\rho}-a\rho{\mathbold F}^{-T}\right)$ (59)

The geometry and boundary conditions are shown in Fig. 11, and the material
properties of soil and CO2 gas are listed in table 2. This provides a
simplified model of CO2 injection through a well. The numerical calculation
uses soil porosity of $\displaystyle\phi_{0f}=$\displaystyle 0.2$$ and a
constant flux $\displaystyle{\mathbold q}=$\displaystyle
4.4\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{2}\,\mathrm{s}$$
$\displaystyle\left($\displaystyle
2.225\text{\times}{10}^{3}\text{\,}\mathrm{L}\mathrm{/}\mathrm{m}^{2}\,\mathrm{s}$\right)$
is assumed at the injection point.

Figure 11: Geometry and boundary conditions of the specimen.

Property | Value
---|---
Solid phase Lame constant, $\displaystyle\lambda$ | $\displaystyle 57.7\text{\,}\mathrm{MPa}$
Solid phase Lame constant, $\displaystyle\mu$ | $\displaystyle 38.46\text{\,}\mathrm{MPa}$
Hydraulic conductivity, $\displaystyle\tilde{k}$ | $\displaystyle 0.1\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$
Gas constant, $\displaystyle R$ |  $\displaystyle 8.32$ $\displaystyle{\mathrm{{}^{3}}\mathrm{m}.\mathrm{Pa}}\mathrm{/}{\mathrm{K}.\mathrm{mol}}$
CO2 constant, $\displaystyle a$ |  $\displaystyle 0.364$ $\displaystyle\mathrm{Pa}\,\mathrm{m}^{6}\mathrm{/}\mathrm{mol}^{2}$
CO2 constant, $\displaystyle b$ |  $\displaystyle 42.67\text{\times}{10}^{-6}$ $\displaystyle\mathrm{m}^{3}\mathrm{/}\mathrm{mol}$
Critical temperature, $\displaystyle T_{c}$ | $\displaystyle 303.4\text{\,}\mathrm{K}$

Table 2: Properties of the solid phase and VdW gas (corresponding to CO2, from
[9])

.

Fig. 12 illustrates the variation of pore pressure vs. inverse fluid density
($\displaystyle\rho^{-1}$), along the vertical direction below the center of
injection (shown by the dashed line in Fig 11). As shown in Fig. 12, the
simulation is repeated at different temperatures: below, above, and close to
the critical temperature of CO2. As expected, near the center of injection,
shortly after start of fluid injection, we see a gas-liquid phase
transformation due to the high pressures when the temperatures are below the
critical temperatures ($\displaystyle T<T_{c},~{}T=$\displaystyle
270\text{\,}\mathrm{K}$$ and $\displaystyle T=$\displaystyle
285\text{\,}\mathrm{K}$$).

Figure 12: Pore pressure variation vs. inverse fluid density at different
temperatures at $\displaystyle t=$\displaystyle 20\text{\,}\mathrm{s}$$.

Fig. 13 plots the variation of pore pressure with distance from the injection
point along the vertical direction. Near the center of injection, we have the
highest value of pore pressure which decreases gradually as we move away from
the injection point. We can expect that if the temperature is just below
$\displaystyle T_{c}$, the injected fluid will undergo a phase transformation
in the high pressure region close to the center of injection. Therefore, it is
possible to find the range of parameters such as injection flux and injection
temperature to avoid the phase transformation. We highlight that we do not
consider heat transfer effects in any of these calculations.

Figure 13: Pore pressure variation vs. distance from the injection point at
$\displaystyle t=20s$ for different temperatures.

To further examine the effect of the pore pressure distribution around the
injection point, the predicted pore pressure value at different times is
plotted in Fig 14. This simulation uses a flux $\displaystyle
2.2\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{2}\,\mathrm{s}$
$\displaystyle\left($\displaystyle
1.112\text{\times}{10}^{3}\text{\,}\mathrm{L}\mathrm{/}\mathrm{m}^{2}\,\mathrm{s}$\right)$
and temperature $\displaystyle T=$\displaystyle 320\text{\,}\mathrm{K}$$ at
the injection point. The plot clearly shows the increase of pore pressure with
time.

Figure 14: Pore pressure vs. distance at different times.

##  7 Multicomponent Fluid Phase

In this section, we apply the model to a porous medium with two immiscible
ideal gases, motivated by experimental approaches such as [48].

The Helmholtz free energy density per unit current volume of an ideal gas is:

$W_{i}(\rho_{i})=-\rho_{i}RT\left(1+\frac{3}{2}\log{\left(\frac{3}{2}RT\right)}-\log{(\rho_{i}\xi_{i})}\right)$
(60)

where $\displaystyle R$ and $\displaystyle\xi_{i}$ are constants and
$\displaystyle T$ is the temperature [42]. Following (5), we rewrite the
energy density of the gases as:

$W_{0i}(J,P_{0i},\phi_{i})=\frac{J\phi}{\phi_{0i}}W(J^{-1}P_{0i}/\phi_{i})=-\frac{P_{0i}}{\phi_{i}}RT\left(1+\frac{3}{2}\log\left(\frac{3}{2}RT\right)-\log\left(\frac{P_{0i}\xi_{i}}{J\phi_{i}}\right)\right)$
(61)

The strain energy density of the solid phase is modeled as neo-Hookean as in
(56).

By applying (6) and (10), we find the chemical potential and flux of the fluid
phases:

$\displaystyle\displaystyle\eta_{0i}$
$\displaystyle\displaystyle=-RT\left(\frac{3}{2}\log\left(\frac{3}{2}RT\right)-\log\left(\rho_{i}\xi_{i}\right)\right)$
(62) $\displaystyle\displaystyle{\mathbold q}_{i}$
$\displaystyle\displaystyle=-\phi_{i}{\mathbold k}_{i}RT\nabla\rho_{i}$ (63)

The volume fraction of the solid phase is fixed. Consequently, the unknown
variables of system are $\displaystyle{\mathbold F}$, $\displaystyle P_{i0}$
and $\displaystyle\phi_{i}$. Therefore, the balance of mass for the fluid
phases and the relation between fluid densities can be found using (19):

$\displaystyle\displaystyle\frac{\ \mathrm{d}}{\ \mathrm{d}t}P_{i}J$
$\displaystyle\displaystyle=J\operatorname{div}\left({\mathbold
k}_{i}RT\phi_{i}\nabla\rho_{i}\right)$ (64)
$\displaystyle\displaystyle\frac{P_{1}}{\phi_{1}}$
$\displaystyle\displaystyle=\frac{P_{2}}{\phi_{2}}$ (65)

We assume a porous medium with porosity $\displaystyle\phi_{0f}=$\displaystyle
0.2$$, and relation between the viscosity of two fluids is
$\displaystyle\gamma_{1}=2\gamma_{2}$, and the geometry and boundary
conditions are shown in Fig. 15. The properties of fluids and soil are given
in Table 3.

At the beginning of the injection ($\displaystyle t=0$), the domain is
saturated with gas $\displaystyle 1$; we then begin injecting gas
$\displaystyle 2$ with the flux $\displaystyle q_{2}=$\displaystyle
0.44\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{2}\,\mathrm{s}$$
$\displaystyle\left($\displaystyle
2.225\text{\times}{10}^{2}\text{\,}\mathrm{L}\mathrm{/}\mathrm{m}^{2}\,\mathrm{s}$\right)$
at the injection point. Fig. 16 shows the pore pressure distribution of both
gases around the injection point at $\displaystyle t=$\displaystyle
550\text{\,}\mathrm{s}$$. Fig. 16 (left) shows that the injection of gas
$\displaystyle 2$ into the porous medium causes the density of gas
$\displaystyle 1$ ($\displaystyle P_{1}$) to decrease around the injection
point, and therefore the pore pressure of gas $\displaystyle 1$ decreases in
this region. Further, Fig. 16 (right) shows the increase of the pore pressure
due to gas $\displaystyle 2$ in regions close to the injection point.

Figure 15: Geometry and boundary conditions. Property | Value
---|---
Solid phase Lame modulus, $\displaystyle\lambda$ | $\displaystyle 57.7\text{\,}\mathrm{MPa}$
Solid phase Lame modulus, $\displaystyle\mu$ | $\displaystyle 38.6\text{\,}\mathrm{MPa}$
Hydraulic conductivity, $\displaystyle\tilde{k}_{1}$ | $\displaystyle 0.1\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$
Hydraulic conductivity, $\displaystyle\tilde{k}_{2}$ | $\displaystyle 0.2\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$
Gas constant, $\displaystyle R$ | $\displaystyle 8.32\text{\,}{\mathrm{{}^{3}}\mathrm{m}.\mathrm{Pa}}\mathrm{/}{\mathrm{K}.\mathrm{mol}}$
Table 3: Properties of solid and fluid phases. Figure 16: Pore pressure
distribution of gas $\displaystyle 1$ (left) and gas $\displaystyle 2$ (right)
at $\displaystyle t=$\displaystyle 550\text{\,}\mathrm{s}$$.

Due to the compressibility of fluids, the volume fraction of each gas changes
with time and distance from the injection point. Fig. 17 shows the variation
of volume fraction of each gas with vertical location below the injection
point, at time ($\displaystyle t=$\displaystyle 550\text{\,}\mathrm{s}$$).
After $\displaystyle 550\text{\,}\mathrm{s}$ from the start of the injection,
the volume fraction of gas $\displaystyle 1$ close to the injection point
decreases to $\displaystyle 0.001$ and the volume fraction of gas
$\displaystyle 2$ increases to $\displaystyle 0.199$. Since the medium
initially is saturated with gas $\displaystyle 1$, we find that the volume
fraction of gas $\displaystyle 1$ tends to $\displaystyle 0.2$ far from the
injection point.

Figure 17: Volume fraction of each fluid component vs. distance from the
injection point at $\displaystyle t=$\displaystyle 550\text{\,}\mathrm{s}$$.

##  8 Unsaturated Medium with an Incompressible Fluid

In this section, we discuss an unsaturated porous medium consisting of a
compressible solid phase, incompressible fluid (liquid) phase, and a
compressible fluid (gas) phase, motivated by experimental approaches such as
[23].

Using the ideal gas model for the compressible fluid phase, the total energy
density can be written as:

$W_{0}({\mathbold F},P_{0g},P_{0l})=\phi_{s}W_{0s}({\mathbold
F})-P_{0g}RT\left(1+\frac{3}{2}\log\left(\frac{3}{2}RT\right)-\log\left(\frac{P_{0g}\xi}{J(1-\phi_{s})-\frac{P_{0l}}{\rho_{l}}}\right)\right)$
(66)

Using (6), the chemical potential of gas and liquid can be derived as:

$\displaystyle\displaystyle\eta_{0g}$
$\displaystyle\displaystyle=-RT\left(\frac{3}{2}\log\left(\frac{3}{2}RT\right)-\log\left(\frac{\rho_{0g}\xi}{J(1-\phi_{s})-\frac{P_{0l}}{\rho_{l}}}\right)\right)$
(67) $\displaystyle\displaystyle\eta_{0l}$
$\displaystyle\displaystyle=\frac{\rho_{g}RT}{\rho_{l}}=\frac{P_{g}RT}{(1-\phi_{s}-\phi_{l})\rho_{l}}$
(68)

Further, the flux of gas and liquid can be written as:

$\displaystyle\displaystyle{\mathbold q}_{g}$
$\displaystyle\displaystyle=-\phi_{g}{\mathbold k}_{g}RT\nabla\rho_{g}$ (69)
$\displaystyle\displaystyle{\mathbold q}_{l}$
$\displaystyle\displaystyle=-\phi_{l}{\mathbold k}_{l}RT\nabla\rho_{l}$ (70)

The geometry is the same as Section 7, and the properties of the solid and gas
phases are listed in Table 4. We consider a porous medium with initial solid
volume fraction $\displaystyle\phi_{0s}=0.9$, saturated with the ideal gas
($\displaystyle\phi_{0g}=0.1$ and $\displaystyle p_{g}=$\displaystyle
3300\text{\,}\mathrm{Pa}$$ at $\displaystyle t=0$). We then inject the liquid
with the flux $\displaystyle q_{l}=$\displaystyle
200\text{\,}\mathrm{L}\mathrm{/}\mathrm{m}^{2}\,\mathrm{s}$$.

Property | Value
---|---
Solid phase Lame constant, $\displaystyle\lambda$ | $\displaystyle 57.7\text{\,}\mathrm{MPa}$
Solid phase Lame constant, $\displaystyle\mu$ | $\displaystyle 38.6\text{\,}\mathrm{MPa}$
True permeability, $\displaystyle\kappa$ |  $\displaystyle 1.8\text{\times}{10}^{-7}$ $\displaystyle\mathrm{m}^{2}$
Gas constant, $\displaystyle R$ | $\displaystyle 8.32\text{\,}{\mathrm{{}^{3}}\mathrm{m}.\mathrm{Pa}}\mathrm{/}{\mathrm{K}.\mathrm{mol}}$
Gas viscosity, $\displaystyle\gamma_{g}$ |  $\displaystyle 1.8\text{\times}{10}^{-5}$ $\displaystyle\mathrm{Pa}\text{\,}\mathrm{s}\text{/}$
Liquid density, $\displaystyle\rho_{l}$ |  $\displaystyle 1000$ $\displaystyle\mathrm{kg}\text{/}{\mathrm{m}}^{3}$
Liquid viscosity, $\displaystyle\gamma_{l}$ | $\displaystyle 0.001\text{\,}\mathrm{Pa}\text{\,}\mathrm{s}$
Table 4: Properties of the solid and fluid phases.

Fig. 18 shows the volume fraction distributions for the gas and liquid phases.
As expected, injecting the liquid into the porous medium with gas causes the
liquid to displace the gas in the pore spaces, and the volume fraction of gas
decreases in the vicinity of the injection point. Also, the pore pressure
distribution due to the liquid injection is shown in Fig. 19; we see that the
highest value of pore pressure is in the vicinity of the injection point, and
tends gradually to the far-field pressure of the gas.

Figure 18: Volume fraction of liquid (left) and gas (right) around the
injection point at $\displaystyle t=$\displaystyle 200\text{\,}\mathrm{s}$$.
Figure 19: The pore pressure distribution around the injection point at
$\displaystyle t=$\displaystyle 200\text{\,}\mathrm{s}$$.

##  9 Concluding Remarks

In this study, we have proposed a variational energy-based continuum mechanics
framework in the large deformation regime to model coupled fluid transport and
deformation in a porous deformable medium. The free energy density function of
the porous medium is defined as additively composed of the free-energies of
the solid and fluid phases of the medium. The variational structure provides
for the standard use of the Finite Element Method to solve the proposed model;
in general, variational approaches can provide important advantages for
numerical solution [45].

The proposed formulation was tested on two benchmark problems, namely
Terzaghi’s 1-d and Mandel’s 2-d consolidation problems. Further, the model
enabled the modeling of porous media containing a single (gas-liquid phase
transforming) and multiple immiscible fluid phases. Furthermore, the method
was used to model the behavior of an unsaturated porous medium. As a specific
application of the proposed model, we investigated fluid injection into soil.

An important advantage of the variational formulation, that we are currently
investigating and will be reported in the future, is in enabling an overall
variational structure to the problem of combining a physics-based approach
with a data-driven approach through a Bayesian framework. In addition, plastic
effects are important in many realistic situations; it would therefore be
useful to extend this model to incorporate powerful variational models for
plasticity that have been proposed, e.g. [61, 44].

## Author Contributions

All authors contributed to the formulation of the model. Karimi, Walkington,
and Dayal contributed to the numerical implementation. All authors contributed
to the writing of the paper.

###### Acknowledgements.

We thank the National Science Foundation for support through XSEDE resources
provided by Pittsburgh Supercomputing Center. Mina Karimi acknowledges
financial support from the Scott Institute. Noel Walkington acknowledges
financial support from NSF (DMS 1729478 and DMS 2012259). Kaushik Dayal
acknowledges financial support from NSF (CMMI MOMS 1635407, DMREF 1628994),
ARO (MURI W911NF-19-1-0245), ONR (N00014-18-1-2528), BSF (2018183), and an
appointment to the National Energy Technology Laboratory sponsored by the U.S.
Department of Energy. This work was funded (in part) by the Dowd Fellowship
from the College of Engineering at Carnegie Mellon University; we thank Philip
and Marsha Dowd for their financial support and encouragement. We thank Basant
Sharma for alerting to us to typos in a previous version.

## Data Availability Statement

A version of the code developed for this work is available at
github.com/minakari/Poromechanics.

## References

  * [1] Vaibhav Agrawal and Kaushik Dayal.  A dynamic phase-field model for structural transformations and twinning: Regularized interfaces with transparent prescription of complex kinetics and nucleation. part i: Formulation and one-dimensional characterization.  Journal of the Mechanics and Physics of Solids, 85:270–290, 2015\.
  * [2] Vaibhav Agrawal and Kaushik Dayal.  A dynamic phase-field model for structural transformations and twinning: Regularized interfaces with transparent prescription of complex kinetics and nucleation. part ii: Two-dimensional characterization and boundary kinetics.  Journal of the Mechanics and Physics of Solids, 85:291–307, 2015\.
  * AD [17] Vaibhav Agrawal and Kaushik Dayal.  Dependence of equilibrium griffith surface energy on crack speed in phase-field models for fracture coupled to elastodynamics.  International Journal of Fracture, 207(2):243–249, 2017.
  * AK [06] Rohan Abeyaratne and James K Knowles.  Evolution of phase transitions: a continuum theory.  Cambridge University Press, 2006.
  * Ana [15] Lallit Anand.  A derivation of the theory of linear poroelasticity from chemoelasticity.  Journal of Applied Mechanics, 82(11), 2015.
  * AZ [15] Gianluca Alaimo and Massimiliano Zingales.  Laminar flow through fractal porous materials: the fractional-order transport equation.  Communications in Nonlinear Science and Numerical Simulation, 22(1-3):889–902, 2015.
  * BA [95] Ronaldo I Borja and Enrique Alarcón.  A mathematical framework for finite strain elastoplastic consolidation part 1: Balance laws, variational formulation, and linearization.  Computer Methods in Applied Mechanics and Engineering, 122(1-2):145–171, 1995.
  * BC [96] Lynn Schreyer Bennethum and John H Cushman.  Multiscale, hybrid mixture theory for swelling systems—i: balance laws.  International Journal of Engineering Science, 34(2):125–145, 1996\.
  * Bie [07] Andreas Bielinski.  Numerical simulation of CO2 sequestration in geological formations.  2007\.
  * Bio [35] Maurice A Biot.  Le problem de la consolidation des matieres argileuses sous une charge.  Annaies de la Societe Scientifique de Bruxelles, pages 110–113, 1935.
  * Bio [41] Maurice A Biot.  General theory of three-dimensional consolidation.  Journal of applied physics, 12(2):155–164, 1941.
  * Bio [73] Maurice A Biot.  Nonlinear and semilinear rheology of porous solids.  Journal of Geophysical Research, 78(23):4924–4937, 1973.
  * BOKCW [20] Quan M Bui, Daniel Osei-Kuffuor, Nicola Castelletto, and Joshua A White.  A scalable multigrid reduction framework for multiphase poromechanics of heterogeneous media.  SIAM Journal on Scientific Computing, 42(2):B379–B396, 2020.
  * Bor [04] Ronaldo I Borja.  Cam-clay plasticity. part v: A mathematical framework for three-phase deformation and strain localization analyses of partially saturated porous media.  Computer methods in applied mechanics and engineering, 193(48-51):5301–5338, 2004.
  * Bor [05] Ronaldo I Borja.  Conservation laws for three-phase partially saturated granular media.  In Unsaturated Soils: Numerical and Theoretical Approaches, pages 3–14. Springer, 2005.
  * Bor [06] Ronaldo I Borja.  On the mechanical energy and effective stress in saturated and unsaturated porous continua.  International Journal of Solids and Structures, 43(6):1764–1786, 2006.
  * Bow [82] Ray M Bowen.  Compressible porous media models by use of the theory of mixtures.  International Journal of Engineering Science, 20(6):697–735, 1982\.
  * CA [10] Shawn A Chester and Lallit Anand.  A coupled theory of fluid permeation and large deformations for elastomeric materials.  Journal of the Mechanics and Physics of Solids, 58(11):1879–1906, 2010.
  * Cou [04] Olivier Coussy.  Poromechanics.  John Wiley & Sons, 2004.
  * Cow [99] Stephen C Cowin.  Bone poroelasticity.  Journal of biomechanics, 32(3):217–238, 1999.
  * DJB [20] Ravindra Duddu, Stephen Jiménez, and Jeremy Bassis.  A non-local continuum poro-damage mechanics model for hydrofracturing of surface crevasses in grounded glaciers.  Journal of Glaciology, 66(257):415–429, 2020.
  * dMPBD [18] Robert Buarque de Macedo, Hossein Pourmatin, Timothy Breitzman, and Kaushik Dayal.  Disclinations without gradients: A nonlocal model for topological defects in liquid crystals.  Extreme Mechanics Letters, 23:29–40, 2018.
  * DPS [93] DE Dria, GA Pope, and Kamy Sepehrnoori.  Three-phase gas/oil/brine relative permeabilities measured under co2 flooding conditions.  SPE reservoir engineering, 8(02):143–150, 1993.
  * DSF [10] Fernando P Duda, Angela C Souza, and Eliot Fried.  A theory for species migration in a finitely strained solid with application to polymer network swelling.  Journal of the Mechanics and Physics of Solids, 58(4):515–529, 2010\.
  * Ehl [09] Wolfgang Ehlers.  Challenges of porous media models in geo-and biomechanical engineering including electro-chemically active polymers and gels.  International Journal of Advances in Engineering Sciences and Applied Mathematics, 1(1):1–24, 2009.
  * FEN [19] Zhiqiang Fan, Peter Eichhubl, and Pania Newell.  Basement fault reactivation by fluid injection into sedimentary reservoirs: Poroelastic effects.  Journal of Geophysical Research: Solid Earth, 124(7):7354–7369, 2019.
  * FG [12] Salvatore Federico and Alfio Grillo.  Elasticity and permeability of porous fibre-reinforced materials under large deformations.  Mechanics of Materials, 44:58–71, 2012.
  * FJK+ [21] Marvin Fritz, Prashant K. Jha, Tobias Köppl, J. Tinsley Oden, Andreas Wagner, and Barbara Wohlmuth.  Modeling and simulation of vascular tumors embedded in evolving capillary networks.  Computer Methods in Applied Mechanics and Engineering, 384:113975, 2021.
  * Gaj [10] Alessandro Gajo.  A general approach to isothermal hyperelastic modelling of saturated porous media at finite strains with compressible solid constituents.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 466(2122):3061–3087, 2010.
  * GFA [10] Morton E Gurtin, Eliot Fried, and Lallit Anand.  The mechanics and thermodynamics of continua.  Cambridge University Press, 2010.
  * HG [79] Majid Hassanizadeh and William G Gray.  General conservation equations for multi-phase systems: 1. averaging procedure.  Advances in water resources, 2:131–144, 1979.
  * HZZS [08] Wei Hong, Xuanhe Zhao, Jinxiong Zhou, and Zhigang Suo.  A theory of coupled diffusion and large deformation in polymeric gels.  Journal of the Mechanics and Physics of Solids, 56(5):1779–1793, 2008.
  * JA [20] Wencheng Jin and Chloé Arson.  Fluid-driven transition from damage to fracture in anisotropic porous media: a multi-scale xfem approach.  Acta Geotechnica, 15(1):113–144, 2020.
  * JJ [14] Birendra Jha and Ruben Juanes.  Coupled multiphase flow and poromechanics: A computational model of pore pressure effects on fault slip and earthquake triggering.  Water Resources Research, 50(5):3776–3808, 2014.
  * KMDP [22] Mina Karimi, Mehrdad Massoudi, Kaushik Dayal, and Matteo Pozzi.  High-dimensional nonlinear bayesian inference of poroelastic fields from pressure data.  submitted, 2022.
  * KPC [07] Natalie Kleinfelter, Moongyu Park, and John H Cushman.  Mixture theory and unsaturated flow in swelling soils.  Transport in porous media, 68(1):69–89, 2007.
  * KPZB [12] Samuel CM Krevor, Ronny Pini, Lin Zuo, and Sally M Benson.  Relative permeability and trapping of co2 and water in sandstone rocks at reservoir conditions.  Water resources research, 48(2), 2012.
  * LBR [04] Chao Li, Ronaldo I Borja, and Richard A Regueiro.  Dynamics of porous media at finite strain.  Computer methods in applied mechanics and engineering, 193(36-38):3837–3870, 2004.
  * LHA [18] Xilin Lü, Maosong Huang, and José E Andrade.  Modeling the static liquefaction of unsaturated sand containing gas bubbles.  Soils and foundations, 58(1):122–133, 2018.
  * MD [14] Jason Marshall and Kaushik Dayal.  Atomistic-to-continuum multiscale modeling with long-range electrostatic interactions in ionic solids.  Journal of the Mechanics and Physics of Solids, 62:137–162, 2014\.
  * MDW [16] Christopher W MacMinn, Eric R Dufresne, and John S Wettlaufer.  Large deformations of a soft porous material.  Physical Review Applied, 5(4):044020, 2016.
  * MM [09] Ingo Müller and Wolfgang H Müller.  Fundamentals of thermodynamics and applications: with historical annotations and many citations from Avogadro to Zermelo.  Springer Science & Business Media, 2009.
  * MN [12] Maruti Kumar Mudunuru and KB Nakshatrala.  A framework for coupled deformation–diffusion analysis with application to degradation/healing.  International journal for numerical methods in engineering, 89(9):1144–1170, 2012.
  * OP [04] M Ortiz and A Pandolfi.  A variational cam-clay theory of plasticity.  Computer Methods in Applied Mechanics and Engineering, 193(27-29):2645–2666, 2004.
  * OR [12] John Tinsley Oden and Junuthula Narasimha Reddy.  Variational methods in theoretical mechanics.  Springer Science & Business Media, 2012.
  * PF [90] Oliver Penrose and Paul C Fife.  Thermodynamically consistent models of phase-field type for the kinetic of phase transitions.  Physica D: Nonlinear Phenomena, 43(1):44–62, 1990.
  * Rud [01] JW Rudnicki.  Coupled deformation-diffusion effects in the mechanics of faulting and failure of geomaterials.  Appl. Mech. Rev., 54(6):483–502, 2001.
  * Sar [66] AM Sarem.  Three-phase relative permeability measurements by unsteady-state method.  Society of Petroleum Engineers Journal, 6(03):199–205, 1966.
  * SB [14] Xiaoyu Song and Ronaldo I Borja.  Mathematical framework for unsaturated flow in the finite deformation range.  International Journal for Numerical Methods in Engineering, 97(9):658–682, 2014.
  * Sel [13] Antony PS Selvadurai, editor.  Mechanics of poroelastic media, volume 35.  Springer Science & Business Media, 2013.
  * SF [73] Gilbert Strang and George Fix.  An Analysis of the Finite Element Method.  Prentice-Hall, 1973.
  * Sim [92] Bruce R Simon.  Multiphase poroelastic finite element models for soft tissue structures.  Appl. Mech. Rev., 45, 1992.
  * SKH+ [19] Shriram Srinivasan, Satish Karra, Jeffrey Hyman, Hari Viswanathan, and Gowri Srinivasan.  Model reduction for fractured porous media: a machine learning approach for identifying main flow pathways.  Computational Geosciences, 23(3):617–629, 2019.
  * SLL+ [20] Zhichen Song, Xiang Li, José J Lizárraga, Lianheng Zhao, and Giuseppe Buscarnera.  Spatially distributed landslide triggering analyses accounting for coupled infiltration and volume change.  Landslides, pages 1–14, 2020.
  * SOS [13] WaiChing Sun, Jakob T Ostien, and Andrew G Salinger.  A stabilized assumed deformation gradient finite element formulation for strongly coupled poromechanical simulations at finite strain.  International Journal for Numerical and Analytical Methods in Geomechanics, 37(16):2755–2788, 2013.
  * SWB [16] Shabnam J Semnani, Joshua A White, and Ronaldo I Borja.  Thermoplasticity and strain localization in transversely isotropic materials based on anisotropic critical state plasticity.  International Journal for Numerical and Analytical Methods in Geomechanics, 40(18):2423–2449, 2016.
  * TA [20] Arash Dahi Taleghani and Milad Ahmadi.  Thermoporoelastic analysis of artificially fractured geothermal reservoirs: A multiphysics problem.  Journal of Energy Resources Technology, 142(8), 2020.
  * TW [18] Konstantina Trivisa and Franziska Weber.  Analysis and simulations on a model for the evolution of tumors under the influence of nutrient and drug application.  SIAM Journal on Numerical Analysis, 56(1):542–569, 2018.
  * vSAAGC [20] Virginia von Streng, Rami Abi-Akl, Bianca Giovanardi, and Tal Cohen.  Morphogenesis and proportionate growth: A finite element investigation of surface growth with coupled diffusion.  arXiv preprint arXiv:2005.10747, 2020.
  * Wan [00] Herbert F Wang.  Theory of linear poroelasticity with applications to geomechanics and hydrogeology, volume 2.  Princeton University Press, 2000.
  * WMO [06] Kerstin Weinberg, Alejandro Mota, and Michael Ortiz.  A variational constitutive model for porous metal plasticity.  Computational Mechanics, 37(2):142–152, 2006.