Helmholtz free energy function $\hat{\psi}_{g}(\rho_{s},T_{s})$ in (9).
Parameter values for $b$ can be determined from either isotropic or uniaxial
compression experiments (see Figure 3); however, calibration of this model
parameter requires experimental apparatuses that can reach and sustain
stresses greater than approximately $0.1$% of the elastic moduli (e.g.,
$E^{*}$, $G^{*}$, $K^{*}$) of the constituent solid that composes the
individual grains. Above this relative stress magnitude, the Hertzian contact
theory used to formulate the low-pressure elastic response,
$\hat{\psi}_{c}(\epsilon_{v}^{e},\epsilon_{s}^{e},B)$, begins to significantly
deviate from linear elastic behavior.

Under these loading conditions, the model can be approximated by the following
relationships between the measured volumetric strain $\epsilon_{v}$ and the
applied pressure or axial compression stress $p$:

$\displaystyle\epsilon_{0}^{e}$
$\displaystyle\approx\sqrt{\frac{4\rho_{0}p_{0}}{\bar{\rho}_{s}(1-\theta)p_{r}\bar{K}^{2}}},$
(48) $\displaystyle p$
$\displaystyle\approx\bigg{(}\frac{(1-\phi_{0})(1-\theta)p_{r}}{(1-\epsilon_{v})\rho_{0}}\bigg{)}\bar{K}^{2}(\epsilon_{0}^{e}+\epsilon_{v}-\epsilon_{0})^{2}+\bigg{(}\frac{1-\phi_{0}}{1-\epsilon_{v}}\bigg{)}^{2b+1}\rho_{0}C_{0}^{2}(\epsilon_{0}^{e}+\epsilon_{v}-\epsilon_{0}),$
$\displaystyle b$
$\displaystyle\approx\frac{\displaystyle\log\bigg{[}\frac{p}{\rho_{0}C_{0}^{2}(\epsilon_{0}^{e}+\epsilon_{v}-\epsilon_{0})}-\bigg{(}\frac{(1-\phi_{0})(1-\theta)p_{r}}{(1-\epsilon_{v})\rho_{0}}\bigg{)}\frac{\bar{K}^{2}}{\rho_{0}C_{0}^{2}}(\epsilon_{0}^{e}+\epsilon_{v}-\epsilon_{0})\bigg{]}}{\displaystyle
2\log\bigg{[}\frac{1-\phi_{0}}{1-\epsilon_{v}}\bigg{]}}-\frac{1}{2}.$

Here, $p_{0}$ and $\epsilon_{0}$ are the pressures and strains measured at the
onset of Stage III: shown in Sections 4.1 and 6.2. These measurements are
observer-dependent; however, they help isolate the high-pressure stress
response predicted by the model. Additionally, $\rho_{0}$ is the mass density
of the individual grains; $\theta$ is the _grading index_ , which has a
nominal value of 0.83 (Tengattini et al., 2016); $p_{r}$ is the nominal
reference pressure (frequently 1 kPa); $\bar{K}$ is the dimensionless,
reference bulk modulus; $\phi_{0}$ is the initial porosity of the sample
measured before compaction; and $C_{0}$ is bulk wave-speed of constituent
solid material.

Figure 10: (a) Data from one uniaxial compression test reported in Kuwik et
al. (2022) plotted in stress–strain space and compared with the approximation
of the model in (48). (b) Plot showing the estimated values for $b$ calculated
using (48) for each data point in Stage III of (a) with either a nominal value
for $\bar{K}$ (15340) and $\epsilon_{0}$ (20%) or samples from the range of
$\bar{K}$ values calculated in Figure 4b and initial strains between 20–25%.

Using (48), we calibrate the density–elasticity coupling parameter $b$ with
data from two common experimental geometries. For example, consider the model
parameters for Ottawa sand shown in Table 1. The value listed for $b$ was
determined using the pressure–strain measurements from one uniaxial
compression test reported in Kuwik et al. (2022): shown in Figure 10. The
parameter value of 10 was chosen to best match the observed behavior of Ottawa
sand.

##  Appendix B Numerical Integration Procedure

In this section, we briefly introduce the explicit time-integration procedure
used in Section 4.1 and Section 6 to update the effective granular stress
$\boldsymbol{\sigma}_{s}$ of a single material point, from current time
$t^{k}$ to updated time $t^{k+1}$. As inputs, this algorithm takes in the
material parameters listed in Table 1 and the relevant material point state at
time $t^{k}$: namely,

  1. (i)

the current relative breakage, $B^{k}$;

  2. (ii)

the current solid-phase, effective mass density, $\bar{\rho}_{s}^{k}$;

  3. (iii)

the current inelastic porosity, $\phi_{p}^{k}$;

  4. (iv)

the current left, elastic Cauchy–Green tensor, $\boldsymbol{B}^{ek}$;

  5. (v)

the current solid-phase temperature $T_{s}^{k}$;

  6. (vi)

the current solid-phase, specific internal energy, $\varepsilon_{s}^{k}$;

  7. (vii)

and the current effective granular stress, $\boldsymbol{\sigma}_{s}^{k}$.

Additionally, this algorithm requires an average deformation-rate (i.e., the
velocity gradient $\boldsymbol{L}_{s}^{k}=\nabla\boldsymbol{v}_{s}^{k}$)
experienced by the material during the time increment $\Delta
t^{k}=t^{k+1}-t^{k}$. As outputs, this algorithm updates the relevant material
point states and stresses: namely,

  1. (i)

the updated relative breakage, $B^{k+1}$;

  2. (ii)

the updated solid-phase, effective mass density, $\bar{\rho}_{s}^{k+1}$;

  3. (iii)

the updated inelastic porosity, $\phi_{p}^{k+1}$;

  4. (iv)

the updated left, elastic Cauchy–Green tensor, $\boldsymbol{B}^{ek+1}$;

  5. (v)

the updated solid-phase temperature $T_{s}^{k+1}$;

  6. (vi)

the updated solid-phase, specific internal energy, $\varepsilon_{s}^{k+1}$;

  7. (vii)

and the updated effective granular stress, $\boldsymbol{\sigma}_{s}^{k+1}$.

The stress update algorithm proceeds in three principle steps: (a) estimate
the material state and stresses at $t^{k+1}$, _assuming that no inelastic
deformation occurs during the time increment_ ; (b) check whether the material
is inelastically yielding at $t^{k+1}$ (i.e., $y_{1}\geq 0$, $y_{2}\geq 0$, or
$y_{3}\geq 0$); and (c) if the material is yielding, solve for $\lambda_{1}$,
$\lambda_{2}$, or $\lambda_{3}$ to ensure $y_{1}=0$, $y_{2}=0$, or $y_{3}=0$,
respectively. This method is sometimes referred to as a prediction–correction
method: (a) _predicts_ the material state at the end of the time increment —
sometimes called the _trial state_ — (b) checks whether the trial state is
valid; and (c) _corrects_ the material state in order to obey the model yield
conditions in (23)–(25). A complete outline of the stress update algorithm is
provided in Algorithms 1–3.

Using this algorithm, the entire state of a single, solid-phase material point
can be forward integrate from $t^{k}$ to $t^{k+1}$. In Section 4.1 and
Sections 6.1–6.3, this stress update algorithm is implemented in MATLAB to
integrate the material state over time in simple uniaxial or triaxial loading
conditions. In Section 6.4, this algorithm is implemented within the MPM
numerical simulation framework from Baumgarten and Kamrin (2019a).

Algorithm 1 Stress update algorithm: Part I.

1:

2:function
StressUpdate($B^{k},\bar{\rho}_{s}^{k},\phi_{p}^{k},\boldsymbol{B}^{ek},T_{s}^{k},\varepsilon_{s}^{k},\boldsymbol{L}_{s}^{k}$)

3: $\triangleright$ Update Mass Density:

4: $\bar{\rho}_{s}^{k+1}$ $:=$ $\bar{\rho}_{s}^{k}(1-\Delta t^{k}\
\text{tr}(\boldsymbol{L}_{s}^{k}))$ $\triangleright$ from Eq. (3)

5: $\triangleright$ Estimate Elastic Deformation w/o Plasticity:

6: $\boldsymbol{B}^{e*}$ $:=$ $\boldsymbol{B}^{ek}+\Delta
t^{k}(\boldsymbol{L}_{s}^{k}\boldsymbol{B}^{ek}+\boldsymbol{B}^{ek}\boldsymbol{L}_{s}^{k\top})$
$\triangleright$ Note:
$\boldsymbol{B}^{e}=\boldsymbol{F}^{e\top}\boldsymbol{F}^{e}$

7: $\triangleright$ Estimate Elastic Strain Invariants:

8: $\epsilon_{v}^{e}$ $:=$
$\tfrac{1}{2}\text{tr}(\boldsymbol{B}^{e*})-\tfrac{3}{2}$ $\triangleright$
from Eq. (8)

9: $\epsilon_{s}^{e}$ $:=$
$\sqrt{\tfrac{1}{6}\boldsymbol{B}^{e*}_{0}:\boldsymbol{B}^{e*}_{0}}$
$\triangleright$ $\boldsymbol{B}^{e*}_{0}$ denotes deviator of
$\boldsymbol{B}^{e*}$

10: $J^{e}$ $:=$ $\sqrt{\text{det}(\boldsymbol{B}^{e*})}$ $\triangleright$
from Eq. (8)

11: $\triangleright$ Check $y_{2}$ and $y_{3}$ Yield Functions:

12: if $J^{e}>1$ or $\bar{\rho}_{s}^{k+1}<\rho_{0}(1-\phi_{u}(1-B^{k})^{u})$
then

13: $\triangleright$ Material is Either Pulling Apart ($y_{2}=0$) or Unable to
Support Stresses ($y_{3}=0$)

14: $B^{k+1}$ $:=$ $B^{k}$

15: $\phi_{p}^{k+1}$ $:=$ $\phi_{p}^{k}$

16: $\boldsymbol{B}^{ek+1}$ $:=$ $\boldsymbol{1}$

17: $\varepsilon_{s}^{k+1}$ $:=$ $\varepsilon_{s}^{k}$

18: $T_{s}^{k+1}$ $:=$ $\varepsilon_{s}^{k+1}/c_{v}$

19: $\boldsymbol{\sigma}_{s}^{k+1}$ $:=$ $\boldsymbol{0}$

20: return
$\\{B^{k+1},\bar{\rho}_{s}^{k+1},\phi_{p}^{k+1},\boldsymbol{B}^{ek+1},T_{s}^{k+1},\varepsilon_{s}^{k+1},\boldsymbol{\sigma}_{s}^{k+1}\\}$

21: end if

22: $\triangleright$ Estimate Solid Volume Fraction and Solid Mass Density:

23: solve
$\phi_{s}=\bar{\rho}_{s}^{k+1}/\rho_{0}-{\phi_{s}}^{b+1}({J^{e}}^{-1}-1)$ for
$\phi_{s}$ $\triangleright$ from Eq. (16)

24: $\rho_{s}$ $:=$ $\bar{\rho}_{s}^{k+1}/\phi_{s}$ $\triangleright$ from Eq.
(1)

25: $\triangleright$ Estimate Scalar Stress Measures and Breakage Energy:

26: $\bar{p}$ $:=$ $(\bar{\rho}_{s}^{k+1}/\rho_{0})(1-\theta
B^{k})p_{r}(\bar{K}^{2}{\epsilon_{v}^{e}}^{2}+3\bar{G}\bar{K}{\epsilon_{s}^{e}}^{2})/4$
$\triangleright$ from Eq. (18)

27: $\bar{q}$ $:=$ $(\bar{\rho}_{s}^{k+1}/\rho_{0})(1-\theta
B^{k})3\bar{G}p_{r}(-\bar{K}\epsilon_{v}^{e}\epsilon_{s}^{e}/2+\sqrt{3\bar{G}\bar{K}}{\epsilon_{s}^{e}}^{2}/4)$
$\triangleright$ from Eq. (18)

28: $p^{*}$ $:=$
$\hat{p}_{H}(\rho_{s})[1-\Gamma_{0}(1-\rho_{0}/\rho_{s})/2]+\rho_{0}\Gamma_{0}(\hat{e}_{c}(\rho_{s})+c_{v}(T_{s}^{k}-T_{0}))$
$\triangleright$ from Eq. (19)

29: $E_{B}$ $:=$ $(\bar{\rho}_{s}^{k+1}/\rho_{0})\theta
p_{r}(-\bar{K}^{2}{\epsilon_{v}^{e}}^{3}/12-3\bar{G}\bar{K}\epsilon_{v}^{e}{\epsilon_{s}^{e}}^{2}/4+\bar{G}\sqrt{3\bar{G}\bar{K}}{\epsilon_{s}^{e}}^{3}/4)$
$\triangleright$ from Eq. (22)

30:

31:Continued in Algorithm 2…

Algorithm 2 Stress update algorithm: Part II.

32:

33:Continued from Algorithm 1…

34: $\triangleright$ Estimate Yield Stresses:

35: $\boldsymbol{\sigma}_{y}$ $:=$
$\bar{q}\boldsymbol{B}^{e*}_{0}\boldsymbol{B}^{e*}/(3\epsilon_{s}^{e})-\bar{p}\boldsymbol{B}^{e*}-\phi_{s}p^{*}\hat{C}(\phi_{s},J^{e})\boldsymbol{1}$
$\triangleright$ from Eq. (21)

36: $p_{y}$ $:=$ $-\tfrac{1}{3}\text{tr}(\boldsymbol{\sigma}_{y})$

37: $q_{y}$ $:=$
$\sqrt{\tfrac{3}{2}\boldsymbol{\sigma}_{y0}:\boldsymbol{\sigma}_{y0}}.$
$\triangleright$ $\boldsymbol{\sigma}_{y0}$ denotes deviator of
$\boldsymbol{\sigma}_{y}$

38: $\triangleright$ Check $y_{1}$ Yield Function:

39: $y_{1}$ $:=$ $E_{B}(1-B^{k})^{2}/E_{c}+q_{y}^{2}/(Mp_{y})^{2}-1$
$\triangleright$ from Eqs. (23) and (28)

40: if $y_{1}<0$ then

41: $\triangleright$ Material is Not Yielding

42: $B^{k+1}$ $:=$ $B^{k}$

43: $\phi_{p}^{k+1}$ $:=$ $\phi_{p}^{k}$

44: $\boldsymbol{B}^{ek+1}$ $:=$ $\boldsymbol{B}^{e*}$

45: $\boldsymbol{\sigma}_{s}^{k+1}$ $:=$
$\bar{q}\boldsymbol{B}^{e*}_{0}\boldsymbol{B}^{e*}/(3\epsilon_{s}^{e})-\bar{p}\boldsymbol{B}^{e*}-\phi_{s}p^{*}\hat{A}(\phi_{s},J^{e})\boldsymbol{1}$
$\triangleright$ from Eq. (17)

46: $\varepsilon_{s}^{k+1}$ $:=$ $\varepsilon_{s}^{k}+\tfrac{1}{2}(\Delta
t/\bar{\rho}_{s}^{k+1})(\boldsymbol{\sigma}_{s}^{k}+\boldsymbol{\sigma}_{s}^{k+1}):\boldsymbol{L}_{s}^{k}$
$\triangleright$ from Eq. (5)

47: $T_{s}^{k+1}$ $:=$
$(\varepsilon_{s}^{k+1}-\hat{\psi_{c}}(\epsilon_{v}^{e},\epsilon_{s}^{e},B^{k+1})-\hat{e}_{c}(\rho_{s}))/c_{v}$
$\triangleright$ from Eqs. (18) and (19)

48: return
$\\{B^{k+1},\bar{\rho}_{s}^{k+1},\phi_{p}^{k+1},\boldsymbol{B}^{ek+1},T_{s}^{k+1},\varepsilon_{s}^{k+1},\boldsymbol{\sigma}_{s}^{k+1}\\}$

49: else

50: $\triangleright$ Material is Yielding ($y_{1}=0$): Calculate $\lambda_{1}$
such that $y_{1}=0$.

51: while $y_{1}>0$ do

52: $\triangleright$ Estimate $\lambda_{1}$ (e.g., with Bisection or
Newton–Raphson Method)

53: $\lambda_{1}$ $\leftarrow$ best estimate

54: $\triangleright$ Estimate Inelastic Deformation Rates

55: $\dot{B}^{k}$ $:=$
$\lambda_{1}\frac{(1-B^{k})^{2}}{E_{c}}\cos^{2}(\omega)$ $\triangleright$ from
Eqs. (26) and (28)

56: $\xi_{v}^{pk}$ $:=$
$\lambda_{1}\frac{E_{B}(1-B^{k})^{2}}{E_{c}}\frac{-p_{y}}{(p_{y}^{2}+q_{y}^{2})}\text{sin}^{2}(\omega)+\lambda_{1}M_{d}\frac{q_{y}}{(Mp_{y})^{2}}$
$\triangleright$ from Eqs. (26) and (28)

57: $\xi_{s}^{pk}$ $:=$
$\lambda_{1}\frac{E_{B}(1-B^{k})^{2}}{E_{c}}\frac{q_{y}}{(p_{y}^{2}+q_{y}^{2})}\text{sin}^{2}(\omega)+\lambda_{1}\frac{q_{y}}{(Mp_{y})^{2}}$
$\triangleright$ from Eqs. (26) and (28)

58: $\boldsymbol{\tilde{D}}^{pk}$ $:=$
$\tfrac{3}{2}\xi_{s}^{p}\boldsymbol{\sigma}_{y0}/q_{y}+\tfrac{1}{3}\xi_{v}^{p}\boldsymbol{1}$
$\triangleright$ from Eq. (20)

59:

60:Continued in Algorithm 3…

Algorithm 3 Stress update algorithm: Part III.

61:

62:Continued from Algorithm 2…

63: $\triangleright$ Update Estimated Material State

64: $B^{k+1}$ $:=$ $B^{k}+\Delta t^{k}\dot{B}^{k}$

65: $\phi_{p}^{k+1}$ $:=$ $\phi_{p}^{k}+\Delta
t^{k}\xi_{v}^{pk}(1-\phi_{p}^{k})$ $\triangleright$ from Eq. (12)

66: $\boldsymbol{B}^{ek+1}$ $:=$ $\boldsymbol{B}^{e*}-\Delta
t^{k}(\boldsymbol{\tilde{D}}^{pk}\boldsymbol{B}^{ek}+\boldsymbol{B}^{ek}\boldsymbol{\tilde{D}}^{pk})$

67: $\epsilon_{v}^{e}$ $:=$
$\tfrac{1}{2}\text{tr}(\boldsymbol{B}^{ek+1})-\tfrac{3}{2}$ $\triangleright$
from Eq. (8)

68: $\epsilon_{s}^{e}$ $:=$
$\sqrt{\tfrac{1}{6}\boldsymbol{B}^{ek+1}_{0}:\boldsymbol{B}^{ek+1}_{0}}$
$\triangleright$ $\boldsymbol{B}^{ek+1}_{0}$ denotes deviator of
$\boldsymbol{B}^{ek+1}$

69: $J^{e}$ $:=$ $\sqrt{\text{det}(\boldsymbol{B}^{ek+1})}$ $\triangleright$
from Eq. (8)

70: $\triangleright$ Update Volume Fraction and Mass Density from Eqs. (1) and
(16)

71: solve
$\phi_{s}=\bar{\rho}_{s}^{k+1}/\rho_{0}-{\phi_{s}}^{b+1}({J^{e}}^{-1}-1)$ for
$\phi_{s}$

72: $\rho_{s}$ $:=$ $\bar{\rho}_{s}^{k+1}/\phi_{s}$

73: $\triangleright$ Update Scalar Stress Measures from Eqs. (18), (19) and
(22)

74: $\bar{p}$ $:=$ $(\bar{\rho}_{s}^{k+1}/\rho_{0})(1-\theta
B^{k})p_{r}(\bar{K}^{2}{\epsilon_{v}^{e}}^{2}+3\bar{G}\bar{K}{\epsilon_{s}^{e}}^{2})/4$

75: $\bar{q}$ $:=$ $(\bar{\rho}_{s}^{k+1}/\rho_{0})(1-\theta
B^{k})3\bar{G}p_{r}(-\bar{K}\epsilon_{v}^{e}\epsilon_{s}^{e}/2+\sqrt{3\bar{G}\bar{K}}{\epsilon_{s}^{e}}^{2}/4)$

76: $p^{*}$ $:=$
$\hat{p}_{H}(\rho_{s})[1-\Gamma_{0}(1-\rho_{0}/\rho_{s})/2]+\rho_{0}\Gamma_{0}(\hat{e}_{c}(\rho_{s})+c_{v}(T_{s}^{k}-T_{0}))$

77: $E_{B}$ $:=$ $(\bar{\rho}_{s}^{k+1}/\rho_{0})\theta
p_{r}(-\bar{K}^{2}{\epsilon_{v}^{e}}^{3}/12-3\bar{G}\bar{K}\epsilon_{v}^{e}{\epsilon_{s}^{e}}^{2}/4+\bar{G}\sqrt{3\bar{G}\bar{K}}{\epsilon_{s}^{e}}^{3}/4)$

78: $\triangleright$ Update Yield Stresses and $y_{1}$ Yield Function

79: $\boldsymbol{\sigma}_{y}$ $:=$
$\bar{q}\boldsymbol{B}^{ek+1}_{0}\boldsymbol{B}^{ek+1}/(3\epsilon_{s}^{e})-\bar{p}\boldsymbol{B}^{ek+1}-\phi_{s}p^{*}\hat{C}(\phi_{s},J^{e})\boldsymbol{1}$
$\triangleright$ from Eq. (21)

80: $p_{y}$ $:=$ $-\tfrac{1}{3}\text{tr}(\boldsymbol{\sigma}_{y})$

81: $q_{y}$ $:=$
$\sqrt{\tfrac{3}{2}\boldsymbol{\sigma}_{y0}:\boldsymbol{\sigma}_{y0}}.$
$\triangleright$ $\boldsymbol{\sigma}_{y0}$ denotes deviator of
$\boldsymbol{\sigma}_{y}$

82: $y_{1}$ $:=$ $E_{B}(1-B^{k})^{2}/E_{c}+q_{y}^{2}/(Mp_{y})^{2}-1$
$\triangleright$ from Eqs. (23) and (28)

83: end while

84: $\triangleright$ Update Final Stress and Energy State

85: $\boldsymbol{\sigma}_{s}^{k+1}$ $:=$
$\bar{q}\boldsymbol{B}^{ek+1}_{0}\boldsymbol{B}^{ek+1}/(3\epsilon_{s}^{e})-\bar{p}\boldsymbol{B}^{ek+1}-\phi_{s}p^{*}\hat{A}(\phi_{s},J^{e})\boldsymbol{1}$
$\triangleright$ from Eq. (17)

86: $\varepsilon_{s}^{k+1}$ $:=$ $\varepsilon_{s}^{k}+\tfrac{1}{2}(\Delta
t/\bar{\rho}_{s}^{k+1})(\boldsymbol{\sigma}_{s}^{k}+\boldsymbol{\sigma}_{s}^{k+1}):\boldsymbol{L}_{s}^{k}$
$\triangleright$ from Eq. (5)

87: $T_{s}^{k+1}$ $:=$
$(\varepsilon_{s}^{k+1}-\hat{\psi_{c}}(\epsilon_{v}^{e},\epsilon_{s}^{e},B^{k+1})-\hat{e}_{c}(\rho_{s}))/c_{v}$
$\triangleright$ from Eqs. (18) and (19)

88: return
$\\{B^{k+1},\bar{\rho}_{s}^{k+1},\phi_{p}^{k+1},\boldsymbol{B}^{ek+1},T_{s}^{k+1},\varepsilon_{s}^{k+1},\boldsymbol{\sigma}_{s}^{k+1}\\}$

89: end if

90:end function

##  Appendix C Notes on Model Theory

In this section, we provide additional details about the thermodynamic and
mechanical theories used in formulating the model. The model developed in this
work is constructed within the framework of continuum mechanics, using the
fluid–sediment mixture theories of Bedford and Drumheller (1983) and Jackson
(2000); the kinematic theory of Coleman and Noll (1963); and the breakage
mechanics theory of Einav (2007). Here, we provide brief overviews of these
theories, including derivation of equations presented in the main body of this
work.

###  C.1 Derivation of Governing Equations

We begin this summary by considering the system of governing equations
presented in (3)–(6). These are derived from basic conservation laws — namely,
(i) mass conservation, (ii) momentum conservation, (iii) energy conservation,
and (iv) the imbalance of entropy — as applied to mixtures of fluids and
granular sediments. In this brief derivation we primarily consider the
representative volume of mixed material, $\Omega$, shown in Figure 1a. We
assume that the grains within this volume are composed of a solid material
with _true_ density $\rho_{s}$, and we assume that the pore space between the
grains is occupied by a fluid material with _true_ density $\rho_{f}$. This
representative volume of material, $\Omega$, can therefore be decomposed into
a solid volume, $\Omega_{s}$, and a fluid volume, $\Omega_{f}$, such that
$\Omega=\Omega_{s}\cup\Omega_{f}$ (see Figure 1b). The fraction of the total
volume of $\Omega$ that is occupied by solid material (i.e., $\Omega_{s}$) is
defined as the _solid volume fraction_ , $\phi_{s}$, and the fraction occupied
by fluid material (i.e., $\Omega_{f}$) is defined as the _fluid volume
fraction_ or _porosity_ , $\phi_{f}$.

Figure 1: Illustration of a representative volume element (RVE) of fluid-
saturated, granular material. The representative volume $\Omega$ shown in (a)
is separable into a solid components $\Omega_{s}$ and a fluid component
$\Omega_{f}$ shown in (b). These two components of the RVE, their external
boundaries $\partial\Omega_{s}$ and $\partial\Omega_{f}$, and their interior
surface $\partial\Omega^{*}$ are used to construct _effective_ material
fields, which are continuous in the volume $\Omega$ shown in (c).

Further, we assume that the all of the material within the volume $\Omega$ has
a well-defined _true_ velocity, $\boldsymbol{v}^{*}$; specific internal
energy, $\varepsilon^{*}$; specific entropy $s^{*}$; and temperature $T^{*}$.
If we allow the true material fields to have both spatial and temporal
variations, we may define:

$\boldsymbol{v}^{*}=\boldsymbol{v}^{*}(\boldsymbol{x},t),\quad\varepsilon^{*}=\varepsilon^{*}(\boldsymbol{x},t),\quad
s^{*}=s^{*}(\boldsymbol{x},t),\quad\text{and}\quad
T^{*}=T^{*}(\boldsymbol{x},t),$ (49)

for all spatial points $\boldsymbol{x}$ in the volume $\Omega$ and all times
$t$. From these _true_ material fields, we can define _effective_ material
fields — namely, a solid and fluid effective mass density, $\bar{\rho}_{s}$
and $\bar{\rho}_{f}$; a homogenized solid and fluid velocity,
$\boldsymbol{v}_{s}$ and $\boldsymbol{v}_{f}$; homogenized specific internal
energies, $\varepsilon_{s}$ and $\varepsilon_{f}$; homogenized specific
entropies, $s_{s}$ and $s_{f}$; and average temperatures, $T_{s}$ and $T_{f}$
— which also have spatial and temporal variation:

$\bar{\rho}_{s}=\bar{\rho}_{s}(\boldsymbol{x},t),\quad\boldsymbol{v}_{s}=\boldsymbol{v}_{s}(\boldsymbol{x},t),\quad\varepsilon_{s}=\varepsilon_{s}(\boldsymbol{x},t),\quad
s_{s}=s_{s}(\boldsymbol{x},t),\quad\text{and}\quad
T_{s}=T_{s}(\boldsymbol{x},t),$ (50)

and,

$\bar{\rho}_{f}=\bar{\rho}_{f}(\boldsymbol{x},t),\quad\boldsymbol{v}_{f}=\boldsymbol{v}_{f}(\boldsymbol{x},t),\quad\varepsilon_{f}=\varepsilon_{f}(\boldsymbol{x},t),\quad
s_{f}=s_{f}(\boldsymbol{x},t),\quad\text{and}\quad
T_{f}=T_{f}(\boldsymbol{x},t).$ (51)

For the problems considered in this work, it is generally infeasible to model
the _true_ material behaviors directly. Therefore, these _effective_ material
fields are constructed to represent volume and mass averages of conserved
material quantities over an ensemble of RVEs (e.g., see Figure 1c):

$\displaystyle\int_{\Omega}\bar{\rho}_{s}\ dv\equiv\int_{\Omega_{s}}\rho_{s}\
dv,\quad$ $\displaystyle\int_{\Omega}\bar{\rho}_{f}\
dv\equiv\int_{\Omega_{f}}\rho_{f}\ dv,$ (52)
$\displaystyle\int_{\Omega}\bar{\rho}_{s}\boldsymbol{v}_{s}\
dv\equiv\int_{\Omega_{s}}\rho_{s}\boldsymbol{v}^{*}\ dv,\quad$
$\displaystyle\int_{\Omega}\bar{\rho}_{f}\boldsymbol{v}_{f}\
dv\equiv\int_{\Omega_{f}}\rho_{f}\boldsymbol{v}^{*}\ dv,$
$\displaystyle\int_{\Omega}\bar{\rho}_{s}s_{s}\
dv\equiv\int_{\Omega_{s}}\rho_{s}s^{*}\ dv,\quad$
$\displaystyle\int_{\Omega}\bar{\rho}_{f}s_{f}\
dv\equiv\int_{\Omega_{f}}\rho_{f}s^{*}\ dv,$
$\displaystyle\int_{\Omega}\bar{\rho}_{s}T_{s}\
dv\equiv\int_{\Omega_{s}}\rho_{s}T^{*}\ dv,\quad$
$\displaystyle\int_{\Omega}\bar{\rho}_{f}T_{f}\
dv\equiv\int_{\Omega_{f}}\rho_{f}T^{*}\ dv,$

and finally,

$\displaystyle\int_{\Omega}\bar{\rho}_{s}(\varepsilon_{s}+\tfrac{1}{2}\boldsymbol{v}_{s}\cdot\boldsymbol{v}_{s})\
dv\equiv\int_{\Omega_{s}}\rho_{s}(\varepsilon^{*}+\tfrac{1}{2}\boldsymbol{v}^{*}\cdot\boldsymbol{v}^{*})\
dv,\quad\text{and}$ (52, cont.)
$\displaystyle\int_{\Omega}\bar{\rho}_{f}(\varepsilon_{f}+\tfrac{1}{2}\boldsymbol{v}_{f}\cdot\boldsymbol{v}_{f})\
dv\equiv\int_{\Omega_{f}}\rho_{f}(\varepsilon^{*}+\tfrac{1}{2}\boldsymbol{v}^{*}\cdot\boldsymbol{v}^{*})\
dv.$

If the integrating volumes are chosen to be _large_ relative to the size of
the grains but _small_ relative to the scale of any details of interest (e.g.,
impactor geometries), these _effective_ material fields should be smoothly
varying, continuous, and retain important details about the properties of the
deforming constituent materials.

This construction of the effective material fields ensures that conservation
of mass, momentum, and energy in the _effective_ continua will correspond to
conservation of mass, momentum, and energy in the _true_ solid and fluid
constituents of the mixture. By analogy, we also assume that this feature
extends to the imbalance of entropy across scales; however, a more thorough
analysis of how changing integrating volumes affects the evaluation of entropy
is needed. What remains, then, is to answer the question: _how do fluxes at
the boundary of the volume, $\partial\Omega$, govern the change of mass,
momentum, and energy in the volume, $\Omega$?_

To answer this question, we need to define the boundaries of the spatial
integrating volumes above — namely, $\partial\Omega$, $\partial\Omega_{s}$,
$\partial\Omega_{f}$, and $\partial\Omega^{*}$ in Figure 1. The exterior
boundary of the spatial volume $\Omega$ is $\partial\Omega$, and it can be
decomposed into a solid boundary $\partial\Omega_{s}$ and a fluid boundary
$\partial\Omega_{f}$, with
$\partial\Omega=\partial\Omega_{s}\cup\partial\Omega_{f}$. (Note that
$\partial\Omega$, $\partial\Omega_{s}$, and $\partial\Omega_{f}$ are assumed
to remain fixed in space.) The _interior boundary_ , $\partial\Omega^{*}$,
denotes the boundary between the solid and fluid materials, and is assumed to
be _co-moving with the particles_.

If we assume that the _true_ constituent materials are non-reactive and non-
volatile (i.e, there is no combustion or vaporization at the interior grain
surfaces associated with $\partial\Omega^{*}$), then the conservation of mass
for each material within a representative volume $\Omega$ can be expressed in
the following form:

$\displaystyle\frac{\partial}{\partial t}\int_{\Omega_{s}}{\rho_{s}\
dv}+\int_{\partial\Omega_{s}}{\rho_{s}\boldsymbol{v}^{*}\cdot\hat{\boldsymbol{n}}\
da}=0,$ (53) $\displaystyle\frac{\partial}{\partial
t}\int_{\Omega_{f}}{\rho_{f}\
dv}+\int_{\partial\Omega_{f}}{\rho_{f}\boldsymbol{v}^{*}\cdot\hat{\boldsymbol{n}}\
da}=0,$

with $\hat{\boldsymbol{n}}$ the outward facing surface normal vector. This
statement of mass conservation (i.e., the change of mass within the volume
must be balanced by the flux of mass across its boundary) can be re-expressed
in terms of the _effective_ material fields defined above, _if we assume a
similarity between the volume integral identities in ( 52) and the boundary
integrals of sufficiently large RVEs_:

$\displaystyle\frac{\partial}{\partial
t}\int_{\Omega}{\bar{\rho}_{s}dv}+\int_{\partial\Omega}{\bar{\rho}_{s}\boldsymbol{v}_{s}\cdot\hat{\boldsymbol{n}}da}=0,$
(54) $\displaystyle\frac{\partial}{\partial
t}\int_{\Omega}{\bar{\rho}_{f}dv}+\int_{\partial\Omega}{\bar{\rho}_{f}\boldsymbol{v}_{f}\cdot\hat{\boldsymbol{n}}da}=0.$

Since these expressions should be true for any representative volume,
$\Omega$, we can apply the divergence theorem and isolate the strong-form mass
conservation rules:

$\displaystyle\frac{\partial\bar{\rho}_{s}}{\partial
t}+\text{div}(\bar{\rho}_{s}\boldsymbol{v}_{s})=0,$ (55)
$\displaystyle\frac{\partial\bar{\rho}_{f}}{\partial
t}+\text{div}(\bar{\rho}_{f}\boldsymbol{v}_{f})=0.$

Together with (2), these final expressions can be shown to be equivalent to
(3).

Similarly, the rate of change of the momentum for each material within a
representative volume $\Omega$ can be expressed in the following form:

$\displaystyle\frac{\partial}{\partial
t}\int_{\Omega_{s}}{\rho_{s}\boldsymbol{v}^{*}\
dv}+\int_{\partial\Omega_{s}}{\rho_{s}\boldsymbol{v}^{*}(\boldsymbol{v}^{*}\cdot\hat{\boldsymbol{n}})\
da}=\int_{\Omega_{s}}{\rho_{s}\boldsymbol{g}dv}-\int_{\partial\Omega^{*}}\boldsymbol{t}(\hat{\boldsymbol{n}})\
da+\int_{\partial\Omega_{s}}\boldsymbol{t}(\hat{\boldsymbol{n}})\ da,$ (56)
$\displaystyle\frac{\partial}{\partial
t}\int_{\Omega_{f}}{\rho_{f}\boldsymbol{v}^{*}\
dv}+\int_{\partial\Omega_{f}}{\rho_{f}\boldsymbol{v}^{*}(\boldsymbol{v}^{*}\cdot\hat{\boldsymbol{n}})\
da}=\int_{\Omega_{f}}{\rho_{f}\boldsymbol{g}\
dv}+\int_{\partial\Omega^{*}}\boldsymbol{t}(\hat{\boldsymbol{n}})\
da+\int_{\partial\Omega_{f}}\boldsymbol{t}(\hat{\boldsymbol{n}})\ da,$

with $\boldsymbol{t}(\hat{\boldsymbol{n}})$ the surface traction vector, a
function of the outward pointing surface normal vector $\hat{\boldsymbol{n}}$;
and $\boldsymbol{g}$ the gravitational acceleration vector. This statement of
momentum conservation (i.e., the change of momentum within the volume must be
balanced by flux of momentum across its boundary and the external forces
acting on it) can be re-expressed in terms of the effective material fields in
(52), if we introduce _effective_ internal stresses and interaction forces to
account for the tractions along $\partial\Omega$ and $\partial\Omega^{*}$,
respectively:

$\displaystyle\int_{\partial\Omega}{\boldsymbol{\Sigma}_{s}\hat{\boldsymbol{n}}\
da}\equiv\int_{\partial\Omega_{s}}{\boldsymbol{t}(\hat{\boldsymbol{n}})-\rho_{s}\boldsymbol{\delta
v}_{s}(\boldsymbol{\delta v}_{s}\cdot\hat{\boldsymbol{n}})\ da},$ (57)
$\displaystyle\int_{\partial\Omega}{\boldsymbol{\Sigma}_{f}\hat{\boldsymbol{n}}\
da}\equiv\int_{\partial\Omega_{f}}{\boldsymbol{t}(\hat{\boldsymbol{n}})-\rho_{f}\boldsymbol{\delta
v}_{f}(\boldsymbol{\delta v}_{f}\cdot\hat{\boldsymbol{n}})\ da},$
$\displaystyle\int_{\Omega}(\boldsymbol{f}_{b}+\boldsymbol{f}_{d})dv\equiv\int_{\partial\Omega^{*}}\boldsymbol{t}(\hat{\boldsymbol{n}})da.$

Here $\boldsymbol{\Sigma}_{s}$ is the solid Cauchy stress;
$\boldsymbol{\Sigma}_{f}$ is the fluid Cauchy stress; $\boldsymbol{f}_{b}$ is
the _buoyant_ interaction force; $\boldsymbol{f}_{d}$ is the _drag_
interaction force; and $\boldsymbol{\delta v}_{s}$ and $\boldsymbol{\delta
v}_{f}$ are the _fluctuational_ components of the _true_ material velocities.
These velocity fluctuations are defined as $\boldsymbol{\delta
v}_{s}\equiv\boldsymbol{v}^{*}-\boldsymbol{v}_{s}$ and $\boldsymbol{\delta
v}_{f}\equiv\boldsymbol{v}^{*}-\boldsymbol{v}_{s}$. Importantly, these
fluctuations give rise to the _tortuous stresses_ in the fluid and _thermal
stress_ in the granular solid (e.g., see Batchelor, 1970; Haff, 1983; Coussy,
2004). Substituting these expressions — along with the effective material
identities in (52) — into the momentum balance integrals in (56), we may apply
the divergence theorem to find:

$\displaystyle\frac{\partial}{\partial
t}\int_{\Omega}\bar{\rho}_{s}\boldsymbol{v}_{s}\
dv+\int_{\Omega}\text{div}\big{(}\bar{\rho}_{s}\boldsymbol{v}_{s}\otimes\boldsymbol{v}_{s}\big{)}\
dv=\int_{\Omega}\big{(}\bar{\rho}_{s}\boldsymbol{g}-\boldsymbol{f_{b}}-\boldsymbol{f_{d}}+\text{div}(\boldsymbol{\Sigma}_{s})\big{)}\
dv,$ (58) $\displaystyle\frac{\partial}{\partial
t}\int_{\Omega}\bar{\rho}_{f}\boldsymbol{v}_{f}\
dv+\int_{\Omega}\text{div}\big{(}\bar{\rho}_{f}\boldsymbol{v}_{f}\otimes\boldsymbol{v}_{f}\big{)}\
dv=\int_{\Omega}\big{(}\bar{\rho}_{f}\boldsymbol{g}+\boldsymbol{f_{b}}+\boldsymbol{f_{d}}+\text{div}(\boldsymbol{\Sigma}_{f})\big{)}\
dv,$

with $\otimes$ denoting the tensor product operator.

Following the analyses in Drumheller (2000), Baumgarten et al. (2021), and
Baumgarten (2021), we may assume that the tortuous stresses and thermal
stresses mentioned above are small and that there is a smoothly varying fluid
pressure, $p_{f}$, within the volume $\Omega$. Using these two assumptions,
the Cauchy stresses and buoyant interaction force can be re-expressed as
follows:

$\boldsymbol{\Sigma}_{s}=\boldsymbol{\sigma}_{s}-\phi_{s}p_{f}\boldsymbol{1},\quad\boldsymbol{\Sigma}_{f}=\boldsymbol{\tau}_{f}-\phi_{f}p_{f}\boldsymbol{1},\quad\text{and}\quad\boldsymbol{f}_{b}=p_{f}\nabla\phi_{f},$
(59)

where $\boldsymbol{\sigma}_{s}$, $\boldsymbol{\tau}_{f}$, and $p_{f}$ denote
the effective granular stress tensor, the effective fluid shear stress tensor,
and the fluid pore pressure described in the main text. Substituting these
into the integral expressions in (58), and recognizing that these expressions
should be true for any representative volume $\Omega$, we can isolate the
strong-form momentum conservation rules:

$\displaystyle\frac{\partial\bar{\rho}_{s}\boldsymbol{v}_{s}}{\partial
t}+\text{div}\big{(}\bar{\rho}_{s}\boldsymbol{v}_{s}\otimes\boldsymbol{v}_{s}\big{)}=\bar{\rho}_{s}\boldsymbol{g}-\boldsymbol{f_{d}}+\text{div}\boldsymbol{\sigma}_{s}-\phi_{s}\nabla
p_{f},$ (60)
$\displaystyle\frac{\partial\bar{\rho}_{f}\boldsymbol{v}_{f}}{\partial
t}+\text{div}\big{(}\bar{\rho}_{f}\boldsymbol{v}_{f}\otimes\boldsymbol{v}_{f}\big{)}=\bar{\rho}_{f}\boldsymbol{g}+\boldsymbol{f_{d}}+\text{div}\boldsymbol{\tau}_{f}-\phi_{f}\nabla
p_{f}.$

Together with (2) and (3), these expressions can be shown to be equivalent to
(4).

The final conservation expressions that we derive in this section are for the
conservation of total energy. If we assume that the _true_ constituent
materials have a well defined internal heat flow rate, $\boldsymbol{q}^{*}$,
and a well-defined rate of heat generation, $q^{*}$, associated with the
absolute temperature of the materials, $T^{*}$, then we can express
conservation of energy within a representative volume as follows:

$\displaystyle\frac{\partial}{\partial
t}\int_{\Omega_{s}}\rho_{s}\big{(}\varepsilon^{*}+\tfrac{1}{2}\boldsymbol{v}^{*}\cdot\boldsymbol{v}^{*}\big{)}\
dv=$
$\displaystyle-\int_{\partial\Omega_{s}}{\rho_{s}\big{(}\varepsilon^{*}+\tfrac{1}{2}\boldsymbol{v}^{*}\cdot\boldsymbol{v}^{*}\big{)}(\boldsymbol{v}^{*}\cdot\hat{\boldsymbol{n}})\
da}$ (61)
$\displaystyle-\int_{\partial\Omega_{s}}(\boldsymbol{q}^{*}\cdot\hat{\boldsymbol{n}})\
da+\int_{\partial\Omega^{*}}(\boldsymbol{q}^{*}\cdot\hat{\boldsymbol{n}})\ da$
$\displaystyle+\int_{\Omega_{s}}q^{*}\ dv$
$\displaystyle+\int_{\partial\Omega_{s}}(\boldsymbol{t}(\hat{\boldsymbol{n}})\cdot\boldsymbol{v}^{*})\
da-\int_{\partial\Omega^{*}}(\boldsymbol{t}(\hat{\boldsymbol{n}})\cdot\boldsymbol{v}^{*})\
da$
$\displaystyle+\int_{\Omega_{s}}(\rho_{s}\boldsymbol{g}\cdot\boldsymbol{v}^{*})\
dv,$

and,

$\displaystyle\frac{\partial}{\partial
t}\int_{\Omega_{f}}\rho_{f}\big{(}\varepsilon^{*}+\tfrac{1}{2}\boldsymbol{v}^{*}\cdot\boldsymbol{v}^{*}\big{)}\
dv=$
$\displaystyle-\int_{\partial\Omega_{f}}{\rho_{f}\big{(}\varepsilon^{*}+\tfrac{1}{2}\boldsymbol{v}^{*}\cdot\boldsymbol{v}^{*}\big{)}(\boldsymbol{v}^{*}\cdot\hat{\boldsymbol{n}})\
da}$ (62)
$\displaystyle-\int_{\partial\Omega_{f}}(\boldsymbol{q}^{*}\cdot\hat{\boldsymbol{n}})\
da-\int_{\partial\Omega^{*}}(\boldsymbol{q}^{*}\cdot\hat{\boldsymbol{n}})\ da$
$\displaystyle+\int_{\Omega_{f}}q^{*}\ dv$
$\displaystyle+\int_{\partial\Omega_{f}}(\boldsymbol{t}(\hat{\boldsymbol{n}})\cdot\boldsymbol{v}^{*})\
da+\int_{\partial\Omega^{*}}(\boldsymbol{t}(\hat{\boldsymbol{n}})\cdot\boldsymbol{v}^{*})\
da$
$\displaystyle+\int_{\Omega_{f}}(\rho_{f}\boldsymbol{g}\cdot\boldsymbol{v}^{*})\
dv.$

These statements of energy conservation — i.e., the change in _total energy_
within the volume must be balanced by flux of _total energy_ across its
boundary, the external work acted on it, the flow of heat across its boundary,
and the generation of heat within it — can be re-expressed in terms of
_effective_ material fields in (52), if we introduce _effective_ rates of heat
flow and heat generation:

$\displaystyle\int_{\partial\Omega}(\boldsymbol{q}_{s}\cdot\hat{\boldsymbol{n}})\
da\equiv\int_{\partial\Omega_{s}}\big{(}(\boldsymbol{q}^{*}\cdot\hat{\boldsymbol{n}})+\rho_{s}\varepsilon^{*}(\boldsymbol{\delta
v}_{s}\cdot\hat{\boldsymbol{n}})\big{)}\ da,$ (63)
$\displaystyle\int_{\partial\Omega}(\boldsymbol{q}_{f}\cdot\hat{\boldsymbol{n}})\
da\equiv\int_{\partial\Omega_{f}}\big{(}(\boldsymbol{q}^{*}\cdot\hat{\boldsymbol{n}})+\rho_{f}\varepsilon^{*}(\boldsymbol{\delta
v}_{f}\cdot\hat{\boldsymbol{n}})\big{)}\ da,$

and,

$\displaystyle\int_{\Omega}q_{s}\ dv\equiv\int_{\Omega_{s}}q^{*}\ dv,$ (63,
cont.) $\displaystyle\int_{\Omega}q_{f}\ dv\equiv\int_{\Omega_{f}}q^{*}\ dv,$
$\displaystyle\int_{\Omega}q_{i}\
dv\equiv-\int_{\partial\Omega^{*}}(\boldsymbol{q}^{*}\cdot\hat{\boldsymbol{n}})\
da.$

Additionally, if we assume that the magnitude of the local velocity
fluctuations (i.e., $\boldsymbol{\delta v}_{s}$ and $\boldsymbol{\delta
v}_{f}$) are small and that the interior grain boundaries (associated with
$\partial\Omega^{*}$) generally move with the velocity of the solid particles,
we can use the following approximation for the internal work acting on the
interior boundary $\partial\Omega^{*}$:

$\int_{\partial\Omega^{*}}(\boldsymbol{t}(\hat{\boldsymbol{n}})\cdot\boldsymbol{v}^{*})\
da\approx\int_{\Omega}\bigg{(}(\boldsymbol{f}_{d}+\boldsymbol{f}_{d})\cdot\boldsymbol{v}_{s}-\phi_{s}p_{f}\text{div}(\boldsymbol{v}_{s})+\frac{\phi_{s}p_{f}}{\rho_{s}}\bigg{(}\frac{d^{s}\rho_{s}}{dt}\bigg{)}\bigg{)}\
dv.$ (64)

(Note, the first term in this approximation derives from the definition of
$\boldsymbol{f}_{d}$ and $\boldsymbol{f}_{b}$ above; the second term derives
from the pressure-volume work acting on particles intersecting the boundary
$\partial\Omega$; and the third term derives from the pressure-volume work
acting on interior particles with changing _true_ densities.)

Substituting the identities in (63) and the approximation in (64) into the
integral expressions in (62), we may apply the divergence theorem to find
_effective_ energy conservation expressions:

$\displaystyle\frac{\partial}{\partial
t}\int_{\Omega}\bar{\rho}_{s}\big{(}\varepsilon_{s}+\tfrac{1}{2}\boldsymbol{v}_{s}\cdot\boldsymbol{v}_{s}\big{)}dv=$
$\displaystyle-\int_{\Omega}{\text{div}\big{(}\bar{\rho}_{s}\big{(}\varepsilon_{s}+\tfrac{1}{2}\boldsymbol{v}_{s}\cdot\boldsymbol{v}_{s}\big{)}\boldsymbol{v}_{s}\big{)}\
dv}$ (65)
$\displaystyle+\int_{\Omega}\text{div}(\boldsymbol{\sigma}_{s}\boldsymbol{v}_{s}-\boldsymbol{q}_{s})\
dv$
$\displaystyle+\int_{\Omega}\bigg{(}q_{s}-q_{i}+(\bar{\rho}_{s}\boldsymbol{g}-\boldsymbol{f_{d}}-\phi_{s}\nabla
p_{f})\cdot\boldsymbol{v}_{s}-\frac{\phi_{s}p_{f}}{\rho_{s}}\frac{d^{s}\rho_{s}}{dt}\bigg{)}\
dv,$

and,

$\displaystyle\frac{\partial}{\partial
t}\int_{\Omega}\bar{\rho}_{f}\big{(}\varepsilon_{f}+\tfrac{1}{2}\boldsymbol{v}_{f}\cdot\boldsymbol{v}_{f}\big{)}dv=$
$\displaystyle-\int_{\Omega}{\text{div}\big{(}\bar{\rho}_{f}\big{(}\varepsilon_{f}+\tfrac{1}{2}\boldsymbol{v}_{f}\cdot\boldsymbol{v}_{f}\big{)}\boldsymbol{v}_{f}\big{)}\
dv}$ (66)
$\displaystyle+\int_{\Omega}\text{div}(\boldsymbol{\tau}_{f}\boldsymbol{v}_{f}-\boldsymbol{q}_{f})\
dv$
$\displaystyle+\int_{\Omega}(q_{f}+q_{i}+(\bar{\rho}_{f}\boldsymbol{g}-\phi_{f}\nabla
p_{f})\cdot\boldsymbol{v}_{f}+\boldsymbol{f_{d}}\cdot\boldsymbol{v}_{s})\ dv$
$\displaystyle-\int_{\Omega}\bigg{(}\frac{\phi_{s}p_{f}}{\rho_{s}}\frac{d^{s}\rho_{s}}{dt}+p_{f}\text{div}(\phi_{f}\boldsymbol{v}_{f}+\phi_{s}\boldsymbol{v}_{s})\bigg{)}\
dv.$

Through some careful manipulation, the final integrand above can be simplified
using (1)–(3) to $-\phi_{f}p_{f}/\rho_{f}(d^{f}\rho_{f}/dt)$. Since these
expressions should be valid for any representative volume $\Omega$, we can
isolate the strong-form energy conservation rules:

$\displaystyle\frac{\partial\bar{\rho}_{s}E_{s}}{\partial
t}+\text{div}(\bar{\rho}_{s}E_{s}\boldsymbol{v}_{s})=$ $\displaystyle\
\text{div}(\boldsymbol{\sigma}_{s}\boldsymbol{v}_{s}-\boldsymbol{q}_{s})+\bar{\rho}_{s}\boldsymbol{g}\cdot\boldsymbol{v}_{s}-\boldsymbol{f}_{d}\cdot\boldsymbol{v}_{s}-\phi_{s}\nabla
p_{f}\cdot\boldsymbol{v}_{s}-\frac{\phi_{s}p_{f}}{\rho_{s}}\frac{d^{s}\rho_{s}}{dt}$
(67) $\displaystyle+q_{s}-q_{i},$
$\displaystyle\frac{\partial\bar{\rho}_{f}E_{f}}{\partial
t}+\text{div}(\bar{\rho}_{f}E_{f}\boldsymbol{v}_{f})=$ $\displaystyle\
\text{div}(\boldsymbol{\tau}_{f}\boldsymbol{v}_{f}-\boldsymbol{q}_{f})+\bar{\rho}_{f}\boldsymbol{g}\cdot\boldsymbol{v}_{f}+\boldsymbol{f}_{d}\cdot\boldsymbol{v}_{s}-\phi_{f}\nabla
p_{f}\cdot\boldsymbol{v}_{f}-\frac{\phi_{f}p_{f}}{\rho_{f}}\frac{d^{f}\rho_{f}}{dt}$
$\displaystyle+q_{f}+q_{i},$

with
$E_{s}\equiv\varepsilon_{s}+\tfrac{1}{2}\boldsymbol{v}_{s}\cdot\boldsymbol{v}_{s}$
and
$E_{f}\equiv\varepsilon_{f}+\tfrac{1}{2}\boldsymbol{v}_{f}\cdot\boldsymbol{v}_{f}$,
the specific total energies for each effective material phase. Together with
(2), (3), and (4), this expression can be shown to be identical to the model
equations in (5).

This general approach to deriving the governing equations is also performed
for the imbalance of entropy associated with the model expression in (6);
however, we do not report this derivation here. Further (and more detailed)
discussion and derivations of this system of governing equations can be found
in Drumheller (2000), Jackson (2000), Appendix A of Baumgarten et al. (2021),
and Chapter 2 of Baumgarten (2021).

###  C.2 Derivation of Granular Kinematics

Here, we continue the discussion of the proposed model theory by considering
the physical connection between material motions and deformations within
granular materials, following the basic kinematic principles described in
Gurtin et al. (2010). In this analysis, we consider both the true granular
solid (shown in Figure 2a) alongside the effective continuum material that we
use to represent it (shown in Figure 2b). The purpose of this discussion is to
clarify the elastic-plastic decomposition of material deformations used in the
formulation of the model.

True Granular Solid

(a)

Effective Continuum Material

(b)

Figure 2: Reference and deformed continuum bodies associated with (a) the true
granular solid and (b) the effective continuum material that we use to
represent it. The motion functions
$\boldsymbol{\chi}^{*}(\boldsymbol{X}^{*},t)$ and
$\boldsymbol{\chi}_{s}(\boldsymbol{X}_{s},t)$ map points in the true reference
body ($\boldsymbol{X}^{*}$) and effective reference body
($\boldsymbol{X}_{s}$) to spatial points ($x$) in the true and effective
deformed bodies, respectively. In the effective continuum material, it is
useful to describe the deformation gradient, $\boldsymbol{F}$, which maps an
infinitesimal vector $\boldsymbol{dX}_{s}$ in the effective reference body
$\mathcal{B}$ to a vector $\boldsymbol{dx}_{s}$ in the effective deformed body
$\mathcal{B}_{t}$.

First, consider the _true_ granular material shown in Figure 2a. We identify
the region of space occupied by a _reference configuration_ of this material
as the _true reference body_ , $\mathcal{B}^{*}$. This reference configuration
is commonly taken to be a hypothetical _undeformed_ configuration of the
material that is later moved into a _deformed_ configuration. A position
within this region of space is referred to as a _true material point_ , which
is denoted by $\boldsymbol{X}^{*}$ and represents the location of a _material
particle_.

Similarly, we identify the region of space that is actually occupied by the
_true_ granular material — at some time $t$ — as the _true deformed body_ ,
$\mathcal{B}^{*}_{t}$. Each material point $\boldsymbol{X}^{*}$ in
$\mathcal{B}^{*}$ can be mapped to a unique spatial point $\boldsymbol{x}$ in
$\mathcal{B}^{*}_{t}$ using the _motion function_
$\boldsymbol{\chi}^{*}(\boldsymbol{X}^{*},t)$. The time rate of change of the
motion function — i.e., the rate of change of the position $\boldsymbol{x}$ of
a material point $\boldsymbol{X}^{*}$ — is the _true_ velocity of the granular
solid, $\boldsymbol{v}^{*}$:

$\boldsymbol{v}^{*}=\frac{\partial\boldsymbol{\chi}^{*}(\boldsymbol{X}^{*},t)}{\partial
t}.$ (68)

We also consider the effective continuum material, shown in Figure 2b, that we
use to represent the true granular solid. Using the true reference body
$\mathcal{B}^{*}$, we construct an _effective reference body_ $\mathcal{B}$
through the homogenization procedure described in the previous section. We
refer to a position within this region as an _effective_ material point,
$\boldsymbol{X}_{s}$, which denotes the location of an _effective_ material
particle.

Ideally, the motion and deformation of the effective continuum body
$\mathcal{B}$ should represent the aggregate motion of the true continuum body
$\mathcal{B}^{*}$. This can be accomplished by introducing a second motion
function: the _effective_ motion function
$\boldsymbol{\chi}_{s}(\boldsymbol{X}_{s},t)$, which maps an effective
material point $\boldsymbol{X}_{s}$ in $\mathcal{B}$ to a unique spatial point
$\boldsymbol{x}$ in the _effective deformed body_ $\mathcal{B}_{t}$. Since we
are interested in the aggregate motion of the material in the _neighborhood_
of $\boldsymbol{X}_{s}$, we identify the time rate of change of this motion
function — i.e., the rate of change of the position $\boldsymbol{x}$ mapped to
by a particular effective material point $\boldsymbol{X}_{s}$ — as the
homogenized solid velocity field, $\boldsymbol{v}_{s}$:

$\boldsymbol{v}_{s}=\frac{\partial\boldsymbol{\chi}_{s}(\boldsymbol{X}_{s},t)}{\partial
t}.$ (69)

Importantly, the _effective_ material point $\boldsymbol{X}_{s}$ _does not_
identify a true material point (denoted by $\boldsymbol{X}^{*}$). Rather,
$\boldsymbol{X}_{s}$ denotes a reference position in $\mathcal{B}$ that we use
to follow the aggregate motion and deformation of the material in the
_neighborhood_ of $\boldsymbol{X}_{s}$. This distinction is important. Since
the velocity of the true granular material $\boldsymbol{v}^{*}$ and the
homogenized velocity field $\boldsymbol{v}_{s}$ are _not_ identical, there can
be significant differences between the local deformations of individual solid
grains (in $\mathcal{B}^{*}_{t}$) and the effective deformation of _material
neighborhoods_ (in $\mathcal{B}_{t}$). Additionally, it is possible that _all_
of the _true_ material that is initially in the neighborhood of the effective
material point $\boldsymbol{X}_{s}$ may be exchanged for an equivalent mass of
material from outside of this neighborhood, through regular granular
rearrangement (e.g., see particle size-segregation phenomena in Liu et al.,
2023).

_How, then, can we predict the evolution of the effective granular stress,
$\boldsymbol{\sigma}_{s}$, without direct knowledge about the microscopic
deformations and motions of the individual grains that compose the granular
solid?_ In solid mechanics, it is common to consider the fundamental
connection between material strains — a measure of the local material
deformation — and material stresses. This connection arises due to the
physical link between the motion of atoms and molecules in the material and
the position-dependent nature of inter-atomic and inter-molecular forces.
Since we have chosen to represent the true granular solid with an effective
continuum material, we cannot use specific information about the true velocity
field $\boldsymbol{v}^{*}(\boldsymbol{x},t)$ or the motion function
$\boldsymbol{\chi}^{*}(\boldsymbol{X}^{*},t)$ to determine these true,
microscopic deformations. Here, we are limited to describing material
deformations in aggregate using the homogenized velocity field
$\boldsymbol{v}_{s}(\boldsymbol{x},t)$, the effective motion function
$\boldsymbol{\chi}_{s}(\boldsymbol{X}_{s},t)$, and their respective gradients.

To characterize this aggregate deformation, we consider the motion of the
effective continuum material in the _neighborhood_ of an effective material
point $\boldsymbol{X}_{s}$. To do this, we introduce an infinitesimal vector
$\boldsymbol{dX}_{s}$ in the effective reference body $\mathcal{B}$, shown in
Figure 2b, which points from the effective material point $\boldsymbol{X}_{s}$
to an arbitrary material point within its neighborhood:
$\boldsymbol{X}_{s}+\boldsymbol{dX}_{s}$. In principle, this second material
point can be mapped to a spatial position $\boldsymbol{x}$ using the motion
function $\boldsymbol{\chi}_{s}(\boldsymbol{X}_{s}+\boldsymbol{dX}_{s},t)$.
The deformation of this material neighborhood can therefore be characterized
by relating a _relative position vector_ $\boldsymbol{dx}_{s}$ (in
$\mathcal{B}_{t}$) to the infinitesimal vector $\boldsymbol{dX}_{s}$ (in
$\mathcal{B}$), where:

$\boldsymbol{dx}_{s}=\boldsymbol{\chi}_{s}(\boldsymbol{X}_{s}+\boldsymbol{dX}_{s},t)+\boldsymbol{\chi}_{s}(\boldsymbol{X}_{s},t).$
(70)

As in standard solid mechanics, we relate these vector quantities using the
_effective deformation gradient_ $\boldsymbol{F}$:

$\boldsymbol{F}=\frac{\partial\boldsymbol{\chi}_{s}(\boldsymbol{X}_{s},t)}{\partial\boldsymbol{X}_{s}},\quad\text{and}\quad\boldsymbol{dx}_{s}=\boldsymbol{F}\
\boldsymbol{dX}_{s}.$ (71)

Following the derivation in Gurtin et al. (2010), the time rate of change of
the deformation gradient can be expressed as follows:

$\frac{d^{s}\boldsymbol{F}}{dt}=\boldsymbol{L}_{s}\boldsymbol{F},\quad\text{with}\quad\boldsymbol{L}_{s}\equiv\nabla\boldsymbol{v}_{s}.$
(72)

Conceptually, the deformation gradient tensor $\boldsymbol{F}$ maps the vector
between neighboring effective material points ($\boldsymbol{dX}_{s}$) in the
reference body $\mathcal{B}$ to a relative position vectos
($\boldsymbol{dx}_{s}$) in the deformed body $\mathcal{B}_{t}$. If
$\boldsymbol{F}$ is the identity tensor $\boldsymbol{1}$, then the effective
material in this neighborhood would be _effectively unstretched_ and
_undeformed_.

Recall, however, that $\boldsymbol{X}_{s}$ and its associated mapping function
$\boldsymbol{\chi}_{s}(\boldsymbol{X}_{s},t)$ _do not_ identify the motion of
the _true_ material that composes the individual grains. Similarly,
$\boldsymbol{F}$ describes how the material in the neighborhood of
$\boldsymbol{X}_{s}$ deforms when subjected to the aggregate velocity field
$\boldsymbol{v}_{s}(\boldsymbol{x},t)$. Since the _true_ velocity
$\boldsymbol{v}^{*}$ of the granular solid and the _effective_ velocity
$\boldsymbol{v}_{s}$ of the aggregate material are _not_ identical, there can
be significant differences between the local deformations of individual solid
grains and the effective deformations described by $\boldsymbol{F}$. In
particular, since the grains in the material are not fixed to their neigbors,
it is possible that the material around $\boldsymbol{X}_{s}$ might exhibit
significant distortion and strain (i.e., $\boldsymbol{F}\neq\boldsymbol{1}$)
_without any corresponding deformation within the individual grains that
compose this neighborhood_. In other words, _except for particular
arrangements and loading histories_ , the tensor $\boldsymbol{F}$ is generally
_insufficient_ for characterizing the deformation of the true granular solid —
which, in turn, gives rise to the effective granular stress,
$\boldsymbol{\sigma}_{s}$.

To overcome this limitation, we consider the Kroner–Lee decomposition of the
deformation gradient into an _inelastic_ or _plastic_ component,
$\boldsymbol{F}^{p}$, and an _elastic_ component $\boldsymbol{F}^{e}$.

$\boldsymbol{F}=\boldsymbol{F}^{e}\boldsymbol{F}^{p}.$ (73)

In plasticity theory, this decomposition of material deformation represents a
conceptual mapping of the material in the neighborhood of an effective
material point $\boldsymbol{X}_{s}$ first into an _intermediate space_ with
the _plastic distortion_ tensor $\boldsymbol{F}^{p}$ (Figure 3c,d) and then
into the final _deformed space_ with the _elastic distortion_ tensor
$\boldsymbol{F}^{e}$ (Figure 3d,e). Here we use $\boldsymbol{F}^{p}$ to
represent _any_ aggregate granular motion that does not result in _elastic_
deformations within the individual grains (e.g., granular rearrangement, grain
fragments sliding along fracture planes, sub-grain plasticity, etc.).
$\boldsymbol{F}^{e}$, on the other hand, describes _any_ aggregate granular
motion associated with microscopic, elastic deformations of the grains. Since
the effective deformation gradients $\boldsymbol{F}$ is already disconnected
from the _true_ physical motion of the individual grains (through the
homogenization procedure described in the previous section), we must assume
that this type of decomposition of $\boldsymbol{F}$ is admissible.

Additionally, we assume that $\boldsymbol{F}^{e}$ and $\boldsymbol{F}^{p}$
follow a similar evolution rule to $\boldsymbol{F}$:

$\frac{d^{s}\boldsymbol{F}^{e}}{dt}=\boldsymbol{L}^{e}\boldsymbol{F}^{e},\quad\text{and}\quad\frac{d^{s}\boldsymbol{F}^{p}}{dt}=\boldsymbol{L}^{p}\boldsymbol{F}^{p},\quad\text{with}\quad\boldsymbol{L}_{s}=\boldsymbol{L}^{e}+\boldsymbol{F}^{e}\boldsymbol{L}^{p}\boldsymbol{F}^{e-1}=\boldsymbol{L}^{e}+\boldsymbol{\tilde{D}}^{p}.$
(74)

Here we decompose the effective velocity gradient $\boldsymbol{L}_{s}$ into an
_elastic distortion rate_ $\boldsymbol{L}^{e}$ and a _inelastic distortion
rate_ $\boldsymbol{L}^{p}$. However, since the decomposition of
$\boldsymbol{F}$ is multiplicative, we may assume that $\boldsymbol{L}^{p}$ is
symmetric and transform it into the same space as $\boldsymbol{L}_{s}$, which
we refer to as the _effective inelastic deformation rate_
$\boldsymbol{\tilde{D}}^{p}$:

$\boldsymbol{\tilde{D}}^{p}=\boldsymbol{F}^{e}\boldsymbol{L}^{p}\boldsymbol{F}^{e-1}.$
(75)

In this way, the decomposition of $\boldsymbol{L}_{s}$ in (74) becomes
_additive_.

Together, $\boldsymbol{F}$, $\boldsymbol{F}^{e}$, and $\boldsymbol{F}^{p}$
characterize the aggregate deformation of the effective continuum material
shown in Figure 2b. $\boldsymbol{F}^{e}$ represents the elastic deformations
within the _true_ granular material that give rise to the effective granular
stress $\boldsymbol{\sigma}_{s}$, and $\boldsymbol{F}^{p}$ represents all
other aggregate deformations (e.g., granular rearrangement, fragmentation,
sub-grain plasticity, etc.). However, this picture of material deformation is
incomplete without some model for $\boldsymbol{L}^{e}$ or
$\boldsymbol{\tilde{D}}^{p}$, which govern the evolution of these two
deformation measures. Although the material deformation rate
$\boldsymbol{L}_{s}$ is determined purely by aggregate material motion,
$\boldsymbol{L}^{e}$ and $\boldsymbol{\tilde{D}}^{p}$ must be predicted in
some other way.

###  C.3 Breakage Mechanics Theory

In addition to the elastic-plastic deformation theory above, the constitutive
model proposed in this work also uses the theory of breakage mechanics (Einav,
2007). This modeling approach incorporates the micro-mechanics of particle
fragmentation through the auxiliary variable $B$, called the _relative
breakage_. This internal variable provides additional information about the
_true_ granular solid _that is not otherwise captured_ by the deformation
gradient $\boldsymbol{F}$ or its components $\boldsymbol{F}^{e}$ and
$\boldsymbol{F}^{p}$.

Originally introduced in Hardin (1985), $B$ measures pulverization in granular
sediments through changes in their aggregate _particle size distribution_ —
i.e., crushing generally causes larger grains to fracture into smaller
particles. This particle size distribution is frequently reported using the
_cumulative size distribution function_ , $F(d)$, which describes the mass
fraction of material that is finer than a given sieve diameter $d$. An
illustration of several cumulative size distributions are shown in Figure 3.
Since $F(d)$ represents the mass fraction that would pass through a particular
sieve size, it is easily measured experimentally.

Figure 3: Illustration of cumulative size distribution function, $F(d)$, for
an idealized, mono-disperse granular material (right); a partially-pulverized
granular material (center); and a hypothetical, fully-pulverized granular
material (right). The hypothetical distribution for a fully pulverized
sediment is taken to be a fractal distribution:
$F(d)=(d/d_{\text{max}})^{0.4}$, (see Einav, 2007).

In order to characterize a distribution of particle sizes using the single
internal variable $B$, Einav (2007) proposed the following reduction:

$B\equiv\frac{\int_{d_{\text{min}}}^{d_{\text{max}}}\big{(}F(x)-F_{0}(x)\big{)}x^{-1}dx}{\int_{d_{\text{min}}}^{d_{\text{max}}}\big{(}F_{u}(x)-F_{0}(x)\big{)}x^{-1}dx},\quad\text{with}\quad
B\in[0,1].$ (76)

Here $d_{\text{min}}$ and $d_{\text{max}}$ denote the minimum and maximum
grain diameters in a given sample, respectively; $F_{0}(d)$ denotes the
_initial_ cumulative size distribution of a sample before any particles have
fragmented; and $F_{u}(d)$ denotes the _hypothetical_ fractal distribution
that the pulverized sediment approaches as more and more particles fragment.
In this work, we take the _initial_ cumulative size distribution to be a
_hypothetical_ , mono-disperse distribution at $d=d_{\text{max}}$:

$F_{0}(d)=H(d-d_{\text{max}}),$ (77)

with $H(x)$ the Heaviside function. Additionally, following Einav (2007), we
let the _ultimate_ cumulative size distribution be the _hypothetical_ fractal
distribution beginning at $d=d_{\text{max}}$:

$F_{u}(d)=\bigg{(}\frac{d}{d_{\text{max}}}\bigg{)}^{(3-\alpha)},\quad\text{with}\quad\alpha=2.6.$
(78)

Using this approach, information about the distribution of grain sizes can be
incorporated into the model through the Helmholtz free energy function,
$\psi_{s}$, in (9); the compression yield function, $y_{1}$, in (23); the
inelastic deformation rates in (26); and the limiting porosities,
$\phi_{\text{max}}$ and $\phi_{\text{min}}$ from Rubin and Einav (2011).
Further (and more detailed) discussion of this modeling theory can be found in
Einav (2007), Nguyen and Einav (2009), and Rubin and Einav (2011).

###  C.4 Coleman–Noll Procedure and Expressions for $\hat{A}$ and $\hat{C}$

We now turn to thermodynamic restrictions that are specific to the model
proposed in this work. In particular, we consider the model restrictions
imposed by the selection of specific functional forms for the Helmholtz free
energies, $\psi_{s}$ and $\psi_{f}$ in (9). This process is frequently
referred to as the Coleman–Noll procedure (see Coleman and Noll, 1963) and it
determines the thermodynamic requirements of the model shown in (11), (17),
(21), and (31). In this section, we provide a brief derivation of these model
terms, including explicit forms of the functions $\hat{A}(\phi_{s},J^{e})$ and
$\hat{C}(\phi_{s},J^{e})$ used in the model — shown in (86).

We begin this derivation by recalling the definition of the _effective_
Helmholtz free energies:

$\psi_{s}\equiv\varepsilon_{s}-T_{s}s_{s},\quad\text{and}\quad\psi_{f}\equiv\varepsilon_{f}-T_{f}s_{f},$
(79)

with $\varepsilon_{s}$, $\varepsilon_{f}$, $T_{s}$, $T_{f}$, $s_{s}$, and
$s_{f}$ defined in (52). If we take the governing equations in (5) and (6) as
given, then the evolution rule for the free energies in (79) can expressed in
the following form:

$\displaystyle\bar{\rho}_{s}\frac{d^{s}\psi_{s}}{dt}\leq\boldsymbol{\sigma}_{s}:\nabla\boldsymbol{v}_{s}+\frac{\phi_{s}p_{f}}{\rho_{s}}\bigg{(}\frac{d^{s}\rho_{s}}{dt}\bigg{)}-\bar{\rho}_{s}s_{s}\bigg{(}\frac{d^{s}T_{s}}{dt}\bigg{)}+\frac{(T_{s}-T_{f})q_{i}}{T_{f}}-\frac{\boldsymbol{q}_{s}\cdot\nabla
T_{s}}{T_{s}}$ , (80)
$\displaystyle\bar{\rho}_{f}\frac{d^{f}\psi_{f}}{dt}\leq\boldsymbol{\tau}_{f}:\nabla\boldsymbol{v}_{f}+\frac{\phi_{f}p_{f}}{\rho_{f}}\bigg{(}\frac{d^{f}\rho_{f}}{dt}\bigg{)}-\bar{\rho}_{f}s_{f}\bigg{(}\frac{d^{f}T_{f}}{dt}\bigg{)}+\frac{(T_{s}-T_{f})q_{i}}{T_{s}}-\frac{\boldsymbol{q}_{f}\cdot\nabla
T_{f}}{T_{f}}$
$\displaystyle+\boldsymbol{f}_{d}\cdot(\boldsymbol{v}_{s}-\boldsymbol{v}_{f}).$

In the first step of the Coleman–Noll procedure, we assert that since (5) and
(6) must be true for mixtures with _any arbitrary_ temperature fields, $T_{s}$
and $T_{f}$, the inequalities in (80) can only be satisfied if the following
inequalities are also satisfied _independently_ :

$\boldsymbol{q}_{s}\cdot\nabla T_{s}\leq 0,\quad\boldsymbol{q}_{f}\cdot\nabla
T_{f}\leq 0,\quad\text{and}\quad q_{i}(T_{s}-T_{f})\geq 0.$

Additionally, a similar analysis of (80) with respect to _any arbitrary_
difference in material velocities, $\boldsymbol{v}_{s}-\boldsymbol{v}_{f}$,
requires that the final part of (11) also be satisfied _independently_ :

$\boldsymbol{f}_{d}\cdot(\boldsymbol{v}_{s}-\boldsymbol{v}_{f})\geq 0.$

Substituting (11) into (80) yields the model equations presented in (10):

$\displaystyle\bar{\rho}_{s}\frac{d^{s}\psi_{s}}{dt}$
$\displaystyle=\boldsymbol{\sigma}_{s}:\nabla\boldsymbol{v}_{s}+\frac{\phi_{s}p_{f}}{\rho_{s}}\bigg{(}\frac{d^{s}\rho_{s}}{dt}\bigg{)}-\bar{\rho}_{s}s_{s}\bigg{(}\frac{d^{s}T_{s}}{dt}\bigg{)}-D_{s},$
$\displaystyle\bar{\rho}_{f}\frac{d^{f}\psi_{f}}{dt}$
$\displaystyle=\boldsymbol{\tau}_{f}:\nabla\boldsymbol{v}_{f}+\frac{\phi_{f}p_{f}}{\rho_{f}}\bigg{(}\frac{d^{f}\rho_{f}}{dt}\bigg{)}-\bar{\rho}_{f}s_{f}\bigg{(}\frac{d^{f}T_{f}}{dt}\bigg{)}-D_{f}.$

with $D_{s}\geq 0$ and $D_{f}\geq 0$ the positive rates of _mechanical
dissipation_.

We continue this analysis by recalling the functional forms for $\psi_{s}$ and
$\psi_{f}$ in (9) — namely,
$\psi_{s}=\hat{\psi}_{c}(\epsilon_{v}^{e},\epsilon_{s}^{e},B)+\hat{\psi}_{g}(\rho_{s},T_{s})$
and $\psi_{f}=\hat{\psi}_{f}(\rho_{f},T_{s})$. If we substitute these
expressions into (10), apply the chain rule, and group terms on the right-hand
side we find:

$\displaystyle 0=$
$\displaystyle-\bar{\rho}_{s}\bigg{(}\frac{\partial\hat{\psi}_{c}}{\partial\epsilon_{v}^{e}}\frac{d^{s}\epsilon_{v}^{e}}{dt}+\frac{\partial\hat{\psi}_{c}}{\partial\epsilon_{s}^{e}}\frac{d^{s}\epsilon_{s}^{e}}{dt}+\frac{\partial\hat{\psi}_{c}}{\partial
B}\frac{d^{s}B}{dt}+\frac{\partial\hat{\psi}_{g}}{\partial\rho_{s}}\frac{d^{s}\rho_{s}}{dt}+\frac{\partial\hat{\psi}_{g}}{\partial
T_{s}}\frac{d^{s}T_{s}}{dt}\bigg{)}+\boldsymbol{\sigma}_{s}:\nabla\boldsymbol{v}_{s}$
(81)
$\displaystyle+\frac{\phi_{s}p_{f}}{\rho_{s}}\bigg{(}\frac{d^{s}\rho_{s}}{dt}\bigg{)}-\bar{\rho}_{s}s_{s}\bigg{(}\frac{d^{s}T_{s}}{dt}\bigg{)}-D_{s},$
$\displaystyle 0=$
$\displaystyle-\bar{\rho}_{f}\bigg{(}\frac{\partial\psi_{f}}{\partial\rho_{f}}\frac{d^{f}\rho_{f}}{dt}+\frac{\partial\psi_{f}}{\partial
T_{f}}\frac{d^{f}T_{f}}{dt}\bigg{)}+\boldsymbol{\tau}_{f}:\nabla\boldsymbol{v}_{f}+\frac{\phi_{f}p_{f}}{\rho_{f}}\bigg{(}\frac{d^{f}\rho_{f}}{dt}\bigg{)}-\bar{\rho}_{f}s_{f}\bigg{(}\frac{d^{f}T_{f}}{dt}\bigg{)}-D_{f}.$

Additionally, recalling the elastic strain definitions in (8) — namely,
$\boldsymbol{E}^{e}=\tfrac{1}{2}(\boldsymbol{F}^{e\top}\boldsymbol{F}^{e}-\boldsymbol{1})$,
$\epsilon_{v}^{e}=\text{tr}(\boldsymbol{E}^{e})$, and
$\epsilon_{s}^{e}=\sqrt{\tfrac{2}{3}\boldsymbol{E}_{0}:\boldsymbol{E}_{0}}$ —
along with the definition of the left Cuachy–Green tensor
($\boldsymbol{B}^{e}=\boldsymbol{F}^{e}\boldsymbol{F}^{e\top}$) we can use
(74) to re-express (81) as follows:

$\displaystyle 0=$
$\displaystyle-\bar{\rho}_{s}\frac{\partial\hat{\psi}_{c}}{\partial\epsilon_{v}^{e}}(\boldsymbol{B}^{e}:\boldsymbol{L}^{e})-\frac{\bar{\rho}_{s}}{3\epsilon_{s}^{e}}\frac{\partial\hat{\psi}_{c}}{\partial\epsilon_{s}^{e}}(\boldsymbol{B}^{e}_{0}\boldsymbol{B}^{e}:\boldsymbol{L}^{e})-\bar{\rho}_{s}\frac{\partial\hat{\psi}_{c}}{\partial
B}\frac{d^{s}B}{dt}$ (82)
$\displaystyle-\bar{\rho}_{s}\frac{\partial\hat{\psi}_{g}}{\partial\rho_{s}}\frac{d^{s}\rho_{s}}{dt}-\bar{\rho}_{s}\frac{d^{s}T_{s}}{dt}\bigg{(}\frac{\partial\hat{\psi}_{g}}{\partial
T_{s}}+s_{s}\bigg{)}+\boldsymbol{\sigma}_{s}:\nabla\boldsymbol{v}_{s}+\frac{\phi_{s}p_{f}}{\rho_{s}}\bigg{(}\frac{d^{s}\rho_{s}}{dt}\bigg{)}-D_{s},$
$\displaystyle 0=$
$\displaystyle-\frac{\phi_{f}}{\rho_{f}}\frac{d^{f}\rho_{f}}{dt}\bigg{(}\rho_{f}^{2}\frac{\partial\psi_{f}}{\partial\rho_{f}}-p_{f}\bigg{)}-\bar{\rho}_{f}\frac{d^{f}T_{f}}{dt}\bigg{(}\frac{\partial\psi_{f}}{\partial
T_{f}}+s_{f}\bigg{)}+\boldsymbol{\tau}_{f}:\nabla\boldsymbol{v}_{f}-D_{f}.$

At this stage, we make a _significant simplifying assumption_ : when the
_true_ solid phase density is changing (i.e., $d^{s}\rho_{s}/dt\neq 0$) the
contribution of the pressure-volume work exerted by the fluid on the solid is
_negligible_ (i.e., $\rho_{s}^{2}\ \partial\hat{\psi}_{g}/\partial\rho_{s}\gg
p_{f}$). Using this assumption, we therefore isolate two _independent_
expressions for the imbalance of free energy:

$\displaystyle 0=$
$\displaystyle-\bar{\rho}_{s}\frac{\partial\hat{\psi}_{c}}{\partial\epsilon_{v}^{e}}(\boldsymbol{B}^{e}:\boldsymbol{L}^{e})-\frac{\bar{\rho}_{s}}{3\epsilon_{s}^{e}}\frac{\partial\hat{\psi}_{c}}{\partial\epsilon_{s}^{e}}(\boldsymbol{B}^{e}_{0}\boldsymbol{B}^{e}:\boldsymbol{L}^{e})-\bar{\rho}_{s}\frac{\partial\hat{\psi}_{c}}{\partial
B}\frac{d^{s}B}{dt}-\bar{\rho}_{s}\frac{\partial\hat{\psi}_{g}}{\partial\rho_{s}}\frac{d^{s}\rho_{s}}{dt}$
(83)
$\displaystyle-\bar{\rho}_{s}\frac{d^{s}T_{s}}{dt}\bigg{(}\frac{\partial\hat{\psi}_{g}}{\partial
T_{s}}+s_{s}\bigg{)}+\boldsymbol{\sigma}_{s}:\nabla\boldsymbol{v}_{s}-D_{s},$

and,

$0=-\frac{\phi_{f}}{\rho_{f}}\frac{d^{f}\rho_{f}}{dt}\bigg{(}\rho_{f}^{2}\frac{\partial\psi_{f}}{\partial\rho_{f}}-p_{f}\bigg{)}-\bar{\rho}_{f}\frac{d^{f}T_{f}}{dt}\bigg{(}\frac{\partial\psi_{f}}{\partial
T_{f}}+s_{f}\bigg{)}+\boldsymbol{\tau}_{f}:\nabla\boldsymbol{v}_{f}-D_{f}.$
(84)

Let us continue this analysis by first focusing on (83). In particular, we may
substitute $d^{s}\rho_{s}/dt$ using (8) and (74) by recalling the _mechanical
model_ for the solid density in (14):
$\rho_{s}=\rho_{0}(1+\hat{\alpha}(\phi_{s})(J^{e-1}-1))$. Through careful
manipulation of this equation, the time-rate of change of the _true_ solid
mass density can be expressed as,

$\frac{d^{s}\rho_{s}}{dt}=-\rho_{s}\bigg{(}\frac{\hat{\alpha}(\phi_{s})\boldsymbol{L}^{e}:\boldsymbol{1}+\phi_{s}\frac{\partial\hat{\alpha}}{\partial\phi_{s}}(1-J^{e})\boldsymbol{L}_{s}:\boldsymbol{1}}{J^{e}+\hat{\alpha}(\phi_{s})(1-J^{e})+\phi_{s}\frac{\partial\hat{\alpha}}{\phi_{s}}(1-J^{e})}\bigg{)},\quad\text{with}\quad\hat{\alpha}(\phi_{s})=\phi_{s}^{b}.$
(85)

Further, if we introduce $\hat{A}(\phi_{s},J^{e})$ and
$\hat{C}(\phi_{s},J^{e})$ as follows,

$\hat{A}(\phi_{s},J^{e})\equiv\frac{\phi_{s}^{b}+b\phi_{s}^{b}(1-J^{e})}{J^{e}+(b+1)\phi_{s}^{b}(1-J^{e})},\quad\text{and}\quad\hat{C}(\phi_{s},J^{e})\equiv\frac{\phi_{s}^{b}}{J^{e}+(b+1)\phi_{s}^{b}(1-J^{e})},$
(86)

then (83) can be re-expressed as,

$\displaystyle 0=$
$\displaystyle\bigg{(}\boldsymbol{\sigma}_{s}-\bar{\rho}_{s}\frac{\partial\hat{\psi}_{c}}{\partial\epsilon_{v}^{e}}\boldsymbol{B}^{e}-\frac{\bar{\rho}_{s}}{3\epsilon_{s}^{e}}\frac{\partial\hat{\psi}_{c}}{\partial\epsilon_{s}^{e}}\boldsymbol{B}^{e}_{0}\boldsymbol{B}^{e}+\phi_{s}\rho_{s}^{2}\frac{\partial\hat{\psi}_{g}}{\partial\rho_{s}}\hat{A}(\phi_{s},J^{e})\boldsymbol{1}\bigg{)}:\boldsymbol{L}^{e}$
(87)
$\displaystyle+\bigg{(}\boldsymbol{\sigma}_{s}+\phi_{s}\rho_{s}^{2}\frac{\partial\hat{\psi}_{g}}{\partial\rho_{s}}\big{(}\hat{A}(\phi_{s},J^{e})-\hat{C}(\phi_{s},J^{e})\big{)}\boldsymbol{1}\bigg{)}:\boldsymbol{\tilde{D}}^{p}-\bar{\rho}_{s}\frac{\partial\hat{\psi}_{c}}{\partial
B}\frac{d^{s}B}{dt}-D_{s}$
$\displaystyle-\bar{\rho}_{s}\frac{d^{s}T_{s}}{dt}\bigg{(}\frac{\partial\hat{\psi}_{g}}{\partial
T_{s}}+s_{s}\bigg{)}.$

Here, a significant number of algebraic steps have been omitted; however, for
brevity, we believe these omissions are necessary.

Now, we may conclude the Coleman–Noll procedure by asserting that (84) and
(87) must remain true for _any arbitrary_ densities, temperatures, and rates
of deformation. This necessitates that _each_ of the following expressions
must be true _independently_ :

$\displaystyle\boldsymbol{\sigma}_{s}=\frac{\bar{q}}{3\epsilon_{s}^{e}}\boldsymbol{B}^{e}_{0}\boldsymbol{B}^{e}-\bar{p}\boldsymbol{B}^{e}-\phi_{s}p^{*}\hat{A}(\phi_{s},J^{e})\boldsymbol{1},$
(88)
$\displaystyle\quad\text{with}\quad\bar{p}=-\bar{\rho}_{s}\frac{\partial\hat{\psi}_{c}}{\partial\epsilon_{v}^{e}},\quad\bar{q}=\bar{\rho}_{s}\frac{\partial\hat{\psi}_{c}}{\partial\epsilon_{s}^{e}},\quad
p^{*}=\rho_{s}^{2}\frac{\partial\hat{\psi}_{g}}{\partial\rho_{s}};$
$\displaystyle
D_{s}=\boldsymbol{\sigma}_{y}:\boldsymbol{\tilde{D}}^{p}+E_{B}\frac{d^{s}B}{dt}\geq
0,$ $\displaystyle\quad\text{with}\quad
E_{B}=-\rho_{s}\frac{\partial\hat{\psi}_{c}}{\partial
B},\quad\boldsymbol{\sigma}_{y}=\frac{\bar{q}}{3\epsilon_{s}^{e}}\boldsymbol{B}^{e}_{0}\boldsymbol{B}^{e}-\bar{p}\boldsymbol{B}^{e}-\phi_{s}p^{*}\hat{C}(\phi_{s},J^{e})\boldsymbol{1};$
$\displaystyle p_{f}=\rho_{f}^{2}\frac{\partial\psi_{f}}{\partial\rho_{f}};$
$\displaystyle D_{f}=\boldsymbol{\tau}_{f}:\nabla\boldsymbol{v}_{f}\geq 0;$
$\displaystyle s_{s}=-\frac{\partial\hat{\psi}_{g}}{\partial
T_{s}};\quad\text{and}\quad s_{f}=-\frac{\partial\psi_{f}}{\partial T_{f}}.$

In (88), we see that the definition of the effective stresses,
$\boldsymbol{\sigma}_{s}$, $\boldsymbol{\tau}_{f}$, and $p_{f}$ _are not_
arbitrary, they are a direct consequence of the specific choice of Helmholtz
free energy functions in (9) and enforcement of the first and second laws of
thermodynamics.

###  C.5 Proof of Model Dissipation

The final component of the modeling theory that we discuss in this section
concerns proving that the model obeys the thermodynamic requirements listed in
(88). In particular, we are interested in showing that,

$D_{s}=\boldsymbol{\sigma}_{y}:\boldsymbol{\tilde{D}}^{p}+E_{B}\frac{d^{s}B}{dt}\geq
0.$

The other components of (88) are designed to be true by construction. However,
because the model is _non-associative_ — i.e., the inelastic yield functions
in (23)—(25) _do not_ uniquely determine the inelastic deformation rate,
$\boldsymbol{\tilde{D}}^{p}$ — it is important to validate that the model is
thermodynamically sound by substituting (20) into (88) as follows:

$D_{s}=\boldsymbol{\sigma}_{y}:\bigg{(}\frac{3\xi_{s}^{p}}{2q_{y}}\boldsymbol{\sigma}_{y0}+\tfrac{1}{3}\big{(}\xi_{v}^{p}+\xi_{2}^{p}+\xi_{3}^{p}\big{)}\boldsymbol{1}\bigg{)}+E_{B}\frac{d^{s}B}{dt}.$
(89)

Here
$q_{y}=\sqrt{\tfrac{3}{2}\boldsymbol{\sigma}_{y0}:\boldsymbol{\sigma}_{y0}}$,
and $\xi_{s}^{p}$, $\xi_{v}^{p}$, $\xi_{2}^{p}$, and $\xi_{3}^{p}$ are defined
in (26)–(30). We can re-express (89) by applying the associative property of
the tensor contraction operator (:), as follows:

$D_{s}=\frac{3\xi_{s}^{p}}{2q_{y}}\boldsymbol{\sigma}_{y0}:\boldsymbol{\sigma}_{y0}+\tfrac{1}{3}\big{(}\xi_{v}^{p}+\xi_{2}^{p}+\xi_{3}^{p}\big{)}\text{tr}(\boldsymbol{\sigma}_{y})+E_{B}\frac{d^{s}B}{dt}.$
(90)

Recalling that $p_{y}=-\tfrac{1}{3}\text{tr}(\boldsymbol{\sigma}_{y})$, we
arrive at the simple expression for mechanical dissipation from (22):

$D_{s}=q_{y}\xi_{s}^{p}-p_{y}(\xi_{v}^{p}+\xi_{2}^{p}+\xi_{3}^{p})+E_{B}\frac{d^{s}B}{dt}.$
(91)

To prove that this expression for $D_{s}$ satisfies the condition that
$D_{s}\geq 0$, we consider four independent conditions:

  1. (i)

no yielding (i.e., $y_{1}<0$, $y_{2}<0$, and $y_{3}<0$);

  2. (ii)

dense compression yielding (i.e., $y_{1}=0$, $y_{2}<0$, $y_{3}<0$);

  3. (iii)

dense tensile yielding (i.e., $y_{1}<0$, $y_{2}=0$, $y_{3}<0$);

  4. (iv)

and dilute compression yielding (i.e., $y_{1}<0$, $y_{2}<0$, $y_{3}=0$).

The first case above is trivial: in the absence of yielding (23)–(25) restrict
$\xi_{s}^{p}$, $\xi_{v}^{p}$, $\xi_{2}^{p}$, and $\xi_{3}^{p}$ to zero, along
with $d^{s}B/dt=0$. Together, this ensures that $D_{s}=0$. For the other three
conditions, however, we must substitute in the yield equations in (23)–(25)
along with the flow rules in (26)–(30) into (91) directly.

We begin with case (ii): $y_{1}=0$. If $y_{1}=0$, then it can be shown that
$y_{2}<0$ and $y_{3}<0$, and accordingly $\xi_{2}^{p}=\xi_{3}^{p}=0$.
Following (23), $y_{1}=0$ also implies that,

$\frac{E_{B}(1-B)^{2}}{E_{c}}+\frac{q_{y}^{2}}{(Mp_{y})^{2}}=1.$

Returning to (91), we can substitute the flow equations in (26), along with
$\xi_{s}^{p}=\xi_{3}^{p}=0$, to find:

$\displaystyle D_{s}=\lambda_{1}\bigg{(}$
$\displaystyle\frac{E_{B}(1-B)^{2}}{E_{c}}\frac{q_{y}^{2}}{(p_{y}^{2}+q_{y}^{2})}\sin^{2}(\omega)+\frac{q_{y}^{2}}{(Mp_{y}^{2})}$
$\displaystyle+\frac{E_{B}(1-B)^{2}}{E_{c}}\frac{p_{y}^{2}}{(p_{y}^{2}+q_{y}^{2})}\sin^{2}(\omega)-M_{d}\frac{q_{y}}{M^{2}p_{y}}$
$\displaystyle+\frac{E_{B}(1-B)^{2}}{E_{c}}\cos^{2}(\omega)\bigg{)}.$

Together with (23), this expression simplifies to:

$D_{s}=\lambda_{1}\bigg{(}1-\frac{M_{d}q_{y}}{M^{2}p_{y}}\bigg{)}\quad\text{if}\quad
y_{1}=0.$ (92)

Since, in general, the model requires that $q_{y}<Mp_{y}$ and $M_{d}<M$ — see
(28) — it is simple to show that (92) satisfies the requirement that
$D_{s}\geq 0$.

Next, we consider case (iii) above: $y_{2}=0$. If $y_{2}=0$, then it can be
shown that $y_{1}<0$ and $y_{3}<0$, and accordingly
$\xi_{s}^{p}=\xi_{v}^{p}=\xi_{3}^{p}=0$ and $d^{s}B/dt=0$. Following (24),
$y_{2}=0$ also implies that $p_{y}=0$. Substitution of these equalities into
(91) yields the trivial result: $D_{s}=0$ when $y_{2}=0$.

Finally, we consider case (iv) above: $y_{3}=0$. If $y_{3}=0$, then it can be
shown that $y_{1}<0$ and $y_{2}<0$, and accordingly
$\xi_{s}^{p}=\xi_{v}^{p}=\xi_{2}^{p}=0$ and $d^{s}B/dt=0$. Following (25),
$y_{3}=0$ also implies that $\lambda_{3}\leq 0$. If we substitute the flow
equation in (30) into (91), we find:

$D_{s}=\lambda_{3}p_{y},\quad\text{if}\quad y_{3}=0.$ (93)

Since the model requires that $\lambda_{3}\geq 0$ and $p_{y}\geq 0$ — by (24)
— it is simple to show that (93) satisfies the requirement that $D_{s}\geq 0$.

##  Appendix D Explicit Expressions for Model EOS Functions

In this section, we provide additional details about the model equations used
for the solid pressure term, $p^{*}$, in (17) and (21) and the fluid pore
pressure, $p_{f}$. The purpose of this section is to provide brief
descriptions of these model equations alongside the abbreviated explicit
equations of state. The equations selected for these two components of the
model derive from existing models in the literature. For this reason, we
direct the reader to other sources for further reading (e.g., Mie, 1903;
Grüneisen, 1912; Tillotson, 1962; Brundage, 2013; Tonge and Ramesh, 2016).

###  D.1 Mie–Grüneisen EOS

In the model proposed in this work, the solid pressure term, $p^{*}$, is used
to represent the response of a granular sediment in compression _when
approaching zero porosity_ — i.e., $\phi_{f}=0$ and $\phi_{s}=1$. From (17),
we have defined $p^{*}$ as follows:

$p^{*}=\rho_{s}^{2}\frac{\partial\hat{\psi}_{g}}{\partial\rho_{s}}.$

At this stage, we may define $p^{*}$ using _any thermodynamically consistent_
model — i.e., a model with either an explicit or implicit equation for
$\hat{\psi}_{g}(\rho_{s},T_{s})$ in (9). In this work, we have chosen the
Mie–Grüneisen EOS to define this relationship.

The Mie–Grüneisen EOS is a common model for solids in the shock-compression
literature (e.g., see Heuzé, 2012; Tonge and Ramesh, 2016), which utilizes a
simple pressure-energy relationship:

$p^{*}-p_{0}=\rho_{s}\Gamma(\varepsilon^{*}-\varepsilon_{0}).$ (94)

Here, $\varepsilon^{*}$ denotes the _true_ specific internal energy of the
material and generally has two contributions: a thermal component
($c_{v}T^{*}$) and a cold component ($\hat{e}_{c}(\rho_{s})$). For the model
proposed in this work, we adopt a similar definition for the granular internal
energy associated with $p^{*}$:

$\varepsilon_{g}=\hat{e}_{c}(\rho_{s})+c_{v}T_{s}.$ (95)

Following the work of Drumheller (1998), we let the specific entropy
associated with $p^{*}$ take the following form:

$s_{g}=c_{v}\ln(T_{s}/T_{0})-c_{v}\Gamma_{0}(1-\rho_{0}/\rho_{s})+s_{0},$ (96)

where $T_{0}$ and $\rho_{0}$ are reference temperatures and densities at zero
pressure and $s_{0}$ is the non-zero reference entropy. Together, these define
the specific Helmholtz free energy function, $\hat{\psi}_{g}(\rho_{s},T_{s})$:

$\hat{\psi}_{g}(\rho_{s},T_{s})=\hat{e}_{c}(\rho_{s})+c_{v}T_{s}-T_{s}\big{(}c_{v}\ln(T_{s}/T_{0})-c_{v}\Gamma_{0}(1-\rho_{0}/\rho_{s})+s_{0}\big{)}.$
(97)

The final component of the Mie–Grüneisen EOS is the specific functional form
of $\hat{e}_{c}(\rho_{s})$, which we adopt from Tonge and Ramesh (2016):

$\hat{e}_{c}(\rho_{s})=e^{\Gamma_{0}(1-J)}\big{(}c_{v}\hat{T}_{H}(J)-c_{v}T_{0}\big{)}+c_{v}T_{0},\quad\text{with}\quad
J=\rho_{0}/\rho_{s}.$ (98)

Here, $\hat{T}_{H}(J)$ is the temperature Hugoniot, as a function of the
specific volume ratio, $\rho_{0}/\rho_{s}$. If the corresponding pressure
Hugoniot function, $\hat{p}_{H}(J)$, is known, then $\hat{T}_{H}(J)$ can be
expressed as,

$c_{v}\hat{T}_{H}(J)=-\int_{1}^{J}\frac{\hat{p}_{H}(J^{\prime})}{\rho_{s}}\bigg{[}1-\frac{\Gamma_{0}}{2}(1-J^{\prime})\bigg{]}e^{-\Gamma_{0}(1-J^{\prime})}dJ^{\prime}.$
(99)

For the first-order Mie–Grüneisen EOS, the pressure Hugoniot,
$\hat{p}_{H}(J)$, is given by the following piece-wise function:

$\hat{p}_{H}(J)=\begin{cases}\displaystyle\frac{\rho_{0}C_{0}^{2}(1-J)}{(1-S_{0}(1-J))^{2}}&\quad\text{if}\quad
J<1,\\\\[10.00002pt] \rho_{0}C_{0}^{2}(1-J)&\quad\text{if}\quad J\geq
1,\end{cases}\quad\text{with}\quad J=\rho_{0}/\rho_{s}.$ (100)

Substituting these explicit functions into (97), we arrive at the expression
for $p^{*}$ from (19)

$p^{*}=\rho_{s}^{2}\frac{\partial\hat{\psi}_{g}}{\partial\rho_{s}}=\hat{p}_{H}\bigg{(}\frac{\rho_{0}}{\rho_{s}}\bigg{)}\bigg{[}1-\frac{\Gamma_{0}}{2}\bigg{(}1-\frac{\rho_{0}}{\rho_{s}}\bigg{)}\bigg{]}+\rho_{0}\Gamma_{0}\big{(}\hat{e}_{c}(\rho_{s})+c_{v}(T_{s}-T_{0})\big{)}$

###  D.2 Tillotson EOS

In the model proposed in this work, the fluid pore pressure, $p_{f}$,
describes the response of the pore fluid to compression and tension. Not only
does $p_{f}$ partially determine the motion of the pore fluid, it also drives
the motion of the solid material, through normal tractions along the boundary
$\partial\Omega^{*}$ in Figure 1. From (31), we have defined $p_{f}$ as
follows:

$p_{f}=\rho_{f}^{2}\frac{\partial\psi_{f}}{\partial\rho_{f}}.$

At this stage, we may define $p_{f}$ using _any thermodynamically consistent_
model, and here we choose the Tillotson EOS.

The Tillotson EOS (Tillotson, 1962) is another common model for solids and
liquids in the shock-compression literature (e.g., see Brundage, 2013), which
is slightly more sophisticated than the Mie–Grüneisen EOS. In addition to the
compressive and thermal stresses, the Tillotson EOS also incorporates
vaporization, which may occur in the pore fluid at sufficiently high impact
velocities. As with the Mie–Grüneisen EOS, the specific internal energy of the
fluid $\varepsilon_{f}$ generally has two contributions: a thermal component
($c_{vf}(T_{f}-T_{0})$) and a cold component ($\hat{e}_{cf}(\rho_{f})$). For
the model proposed in this work, we follow the implementation of Brundage
(2013):

$\varepsilon_{f}=\hat{e}_{cf}(\rho_{f})+c_{vf}(T_{f}-T_{0}).$ (101)

However, unlike the Mie–Grüneisen EOS, the Tillotson EOS _does not_ have a
simple closed form expression for the entropy of the pore fluid. Rather, the
Tillotson EOS defines the pressure $p_{f}$ as a function of $\rho_{f}$ and
$\varepsilon_{f}$, with the entropy implicitly determined by these equations.

Additional details about the Tillotson EOS may be found in Tillotson (1962)
and Brundage (2013). Here we provide the full form of the model equations used
in this work:

$\hat{p}_{f1}(\rho_{f},\varepsilon_{f})=\bigg{[}a_{f}+\frac{b_{f}}{\varepsilon_{f}/(E_{0}\eta^{2})+1}\bigg{]}\rho_{f}\varepsilon_{f}+A_{f}\mu+B_{f}\mu^{2},$
(102)
$\hat{p}_{f2}(\rho_{f},\varepsilon_{f})=a_{f}\rho_{f}\varepsilon_{f}+\bigg{[}\frac{b_{f}\rho_{f}\varepsilon_{f}}{\varepsilon_{f}/(E_{0}\eta^{2})+1}+A_{f}\mu
e^{-\beta_{f}(\rho_{0}/\rho_{f}-1)}\bigg{]}e^{-\alpha_{f}(\rho_{0}/\rho_{f}-1)^{2}},$
(103)

and,

$\hat{p}_{f3}(\rho_{f},\varepsilon_{f})=\bigg{[}a_{f}+\frac{b_{f}}{\varepsilon_{f}/(E_{0}\eta^{2})+1}\bigg{]}\rho_{f}\varepsilon_{f}+A_{f}\mu,$
(104)

with $\eta=\rho_{f}/\rho_{0}$ and $\mu=\eta-1$. From $\hat{p}_{f1}$,
$\hat{p}_{f2}$, and $\hat{p}_{f3}$, we construct the fluid pore pressure
function $\hat{p}_{f}(\rho_{f},\varepsilon_{f})$ as follows:

$\hat{p}_{f}(\rho_{f},\varepsilon_{f})=\begin{cases}\hat{p}_{f1}(\rho_{f},\varepsilon_{f}),&\text{if}\quad\rho_{f}\geq\rho_{f0},\\\\[10.00002pt]
\hat{p}_{f1}(\rho_{f},\varepsilon_{f}),&\text{if}\quad\rho_{f0}>\rho_{f}\geq\rho_{\text{IV}},\
\varepsilon_{f}\leq E_{\text{IV}},\\\\[10.00002pt]
\displaystyle\frac{(\varepsilon_{f}-E_{\text{IV}})\hat{p}_{f1}+(E_{\text{CV}}-\varepsilon_{f})\hat{p}_{f2}}{E_{\text{CV}}-E_{\text{IV}}},&\text{if}\quad\rho_{f0}>\rho_{f}\geq\rho_{\text{IV}},\
E_{\text{CV}}>\varepsilon_{f}\geq E_{\text{IV}},\\\\[10.00002pt]
\hat{p}_{f2}(\rho_{f},\varepsilon_{f}),&\text{if}\quad\rho_{0}>\rho_{f},\
\varepsilon_{f}\geq E_{\text{CV}},\\\\[10.00002pt]
\hat{p}_{f3}(\rho_{f},\varepsilon_{f}),&\text{if}\quad\rho_{\text{IV}}>\rho_{f},\
E_{\text{CV}}>\varepsilon_{f}.\end{cases}$ (105)