# Free energy calculations along entropic pathways III. Nucleation of
capillary bridges and bubbles.

Caroline Desgranges and Jerome Delhommelle  Department of Chemistry,
University of North Dakota, Grand Forks ND 58202

###### Abstract

Using molecular simulation, we analyze the capillary condensation and
evaporation processes for Argon confined in a cylindrical nanopore. For this
purpose, we define the entropy of the adsorbed fluid as a reaction coordinate
and determine the free energy associated with both processes along entropic
pathways. For capillary condensation, we identify a complex free energy
profile resulting from the multi-stage nature of this phenomenon. We find
capillary condensation to proceed through the nucleation of a liquid bridge
across the nanopore, followed by its expansion throughout the pore to give
rise to the stable phase of high density. In the case of capillary
evaporation, the free energy profile along the entropy pathway also exhibits
different regimes, corresponding to the initial destabilization of the layered
structure of the fluid followed by the formation, and subsequent expansion, of
a bubble across the nanopore.

##  I Introduction

The stability of liquid bridges between solid surfaces, as well as the
mechanism by which they form, is a central phenomenon in interface science and
adhesion Derjaguin (1992); Coasne and Pellenq (2004); Saugey et al. (2005);
Major et al. (2006); Casanova et al. (2007); Puibasset (2008); Edison and
Monson (2013). The formation of these bridges, as well as the nucleation of
bubbles, actually occurs in a wide range of systems that span many length
scales, from colloidal systems to nanosized capillaries Szoszkiewicz and Riedo
(2005); Bruschi et al. (2010); Greiner et al. (2010). This has been shown e.g.
by atomic force microscopy and surface force apparatus experiments, which have
revealed how liquid contacts can form and snap off, paving the way for
applications in nanotribology and nanolithography Lin et al. (1994); Willett
et al. (2000); Gogotsi et al. (2001); He et al. (2001); Heuberger et al.
(2001); Patel et al. (2002); Maeda and Israelachvili (2002); Jang et al.
(2003); Weeks et al. (2005); Berim and Ruckenstein (2008). In recent years,
the use of molecular modeling and simulation to provide a molecular level
understanding of this phenomenon has drawn considerable interest Thommes et
al. (2006); Fortini and Dijkstra (2006); Horikawa et al. (2011); Monson
(2012); Gommes (2012); Mszr et al. (2013); Zeng et al. (2014); Hiratsuka et
al. (2016). Density functional theory calculations Restagno et al. (2000);
Talanquer and Oxtoby (2001); Ustinov and Do (2005); Men et al. (2009), as well
as molecular dynamics Zhang and Chakrabarti (1994); Gelb and Gubbins (1997);
Landman (1998); Leung and Luzar (2000); Yasuoka et al. (2000) and Monte Carlo
simulations Liu and Grest (1991); Gac et al. (1994); Gelb et al. (1999); Bock
and Schoen (1999); Bolhuis and Chandler (2000); Stroud et al. (2001); Liu and
Monson (2006); Mota and Esteves (2007); Winkler et al. (2010); Nguyen et al.
(2011); Gor et al. (2012); Siderius and Shen (2013); van Meel et al. (2015)
have started to shed light on the formation of solid-like and liquid-like
junctions in nanopores and on the dynamics of cavitation between solid
surfaces. In particular, Neimark and coworkers Kornev et al. (2002); Neimark
et al. (2003); Vishnyakov and Neimark (2003a, b); Neimark and Vishnyakov
(2005a, b, c) have examined the capillary condensation process in cylindrical
nanopores and have shown that the scenario proposed by Everett and Haynes
Everett and Haynes (1972) for a macroscopic capillary could also extend to
nanosized capillaries. This suggested a very complex and intriguing pathway
for the capillary condensation process, involving a series of structural
changes in the adsorbed fluid and the formation of a liquid bridge across the
pore section as an intermediate towards the high density adsorbed phase.

Condensation and evaporation processes are challenging processes to model and
simulate. This is due to the large free energy barrier of nucleation that the
system has to overcome to complete the phase transition Oxtoby and Evans
(1988); ten Wolde and Frenkel (1998); Ten Wolde et al. (1999); Shen and
Debenedetti (1999); Desgranges and Delhommelle (2011); Patel et al. (2011);
Desgranges and Delhommelle (2014); Lauricella et al. (2015). Simulating the
pathways underlying condensation and evaporation in nanoscale capillaries
requires suitable simulation techniques that allow for the sampling of these
rare events by driving the formation of the new phase. In this work, we focus
on elucidating both the capillary condensation and evaporation pathways and on
shedding light on the relation between the structure of the nanoconfined fluid
and its entropy. For this purpose, we extend the recently developed $\mu VT-S$
simulation method Desgranges and Delhommelle (2016a, b) to study these two
phenomena. The $\mu VT-S$ approach is implemented within the grand-canonical
ensemble, in which the chemical potential $\mu$, the temperature $T$ and the
volume $V$ of the system are held fixed. As shown in the case of the
nucleation of liquid droplets for single component systems Desgranges and
Delhommelle (2016a) and binary mixtures Desgranges and Delhommelle (2016b),
this method provides a direct way to calculate the entropy of the system
throughout the nucleation process and, therefore, to use it as a reaction
coordinate for such activated processes. In this work, we carry out $\mu VT-S$
simulations with the aim of (i) simulating the adsorption and desorption of
the fluid in the grand-canonical statistical ensemble, that is especially
suited to study adsorption phenomena, and (ii) using the entropy of the
confined fluid as the reaction coordinate to elucidate the two pathways
underlying nanoscale capillary condensation and evaporation. To achieve this,
we impose a target value for the entropy of the system, through an umbrella
sampling bias potential Torrie and Valleau (1977). Then, by varying the value
of the entropy of the system, we are able to sample the entire pathway
connecting the two nanoconfined phases, starting from the metastable confined
vapor and ending with the nanoconfined liquid in the case of capillary
condensation. Doing so, we are able to identify the different intermediate
structures, specifically bridges and bubbles, that form in the nanopores
during the capillary condensation and evaporation processes.

The paper is organized as follows. We start by discussing how we extend the
$\mu VT-S$ approach to study the condensation and evaporation processes of
Argon in MCM-41 silica mesoporous molecular sieves. We also present the
potential models used for the adsorbate and its interaction with the
cylindrical nanopore, and provide the technical details for the simulations.
We then examine the capillary condensation process for Argon in the nanopore
and determine the free energy barrier associated with the process. We carry
out a detailed structural analysis to identify the formation of the liquid
bridge along this pathway. We also discuss the capillary evaporation process
that occurs during the desorption of the confined fluid and apply the same
analyses as for the reverse pathway. We finally draw our main conclusions in
the last section.

##  II Models

We study the capillary condensation and evaporation of Argon in a cylindrical
pore of 10 atomic diameters. The geometry is defined as follows. The
cylindrical pore is aligned with the $z$-axis, and we define the length, or
lateral dimension, as $L_{z}$. Simulations are carried out at the boiling
point for Argon ($T=87.3$ K) within a nanopore found in MCM-41 silica
mesoporous molecular sieves Vishnyakov and Neimark (2003b). We model Argon
with a Lennard-Jones potential and take the following parameters for the
exclusion diameter $\sigma=3.4$ Å and for the potential well-depth
$\epsilon/k_{B}=119.8$ K. The interactions between the adsorbed Argon atoms
and the nanopore are modeled with the following functional form Neimark et al.
(2003); Vishnyakov and Neimark (2003a, b); Neimark and Vishnyakov (2005a, b,
c)

$\begin{array}[]{ll}U_{sf}(r,R)=&\\\
\pi^{2}\rho_{s}\epsilon_{sf}\sigma_{sf}^{2}\left\\{{63\over
32}\left[{{R-r}\over\sigma_{sf}}\left(1+{r\over
R}\right)\right]^{-10}F\left[{-9\over 2},{-9\over 2};1;\left({r\over
R}\right)^{2}\right]-3\left[{{R-r}\over\sigma_{sf}}\left(1+{r\over
R}\right)\right]^{-4}F\left[{-3\over 2},{-3\over 2};1;\left({r\over
R}\right)^{2}\right]\right\\}&\\\ \end{array}$ (1)

in which $r$ is the radial coordinate of the Ar atom adsorbed in the pore, $R$
is the pore radius (here $5\sigma$), $\rho_{s}$ is the surface density of
adsorption centers and $F(\alpha,\beta;\gamma;\delta)$ is the hypergeometric
series. The parameters for the solid-fluid interactions were taken as
$\rho_{s}\epsilon_{sf}=2253$ K/nm2 and $\sigma_{sf}=3.17$ Å. This functional
form accurately models the interaction between adsorbate and the structureless
cylindrical layer of adsorption centers on the pore wall as shown by
Ravikovitch et al. Ravikovitch et al. (1997). Following Vishnyakov and Neimark
Vishnyakov and Neimark (2003b), we carry out simulations of capillary
condensation and evaporation in nanopores with a long lateral dimension
$L_{z}=30\sigma$ to allow for the sampling of symmetry breaking configurations
containing bubbles and liquid bridges. Periodic boundary conditions are
applied along this lateral direction $z$. We also calculate explicitly the
interactions between Argon atoms up to a distance of $5\sigma$ and neglect the
fluid-fluid interactions beyond that cutoff distance. Finally, in the rest of
this work, we use the conventional set of reduced units Allen and Tildesley
(1987), with respect to the Lennard-Jones parameters of the fluid.

##  III Simulation method

To span the entire pathway underlying the capillary condensation and
evaporation processes, we carry out a series of $\mu VT-S$ simulations and
gradually vary the entropy of the system to drive the phase transition within
the confined fluid. We briefly outline here the principles of the $\mu VT-S$
method (more details may be found in previous work Desgranges and Delhommelle
(2016a, b)). $\mu VT-S$ simulations are carried out in the grand-canonical
ensemble $(\mu,V,T)$, i.e. at constant chemical potential $\mu$, temperature
$T$ and volume $V$. This statistical ensemble is especially well suited to
study adsorption phenomena, since the number of atoms $N$ adsorbed in the
nanopore is allowed to vary in the $(\mu,V,T)$ ensemble. Thus, this ensemble
mimics what is observed experimentally during the adsorption/desorption
processes. This ensemble has also the advantage of providing a direct way of
estimating the entropy of the confined fluid. Taking the total Helmoltz free
energy $A=U-TS$ of the confined fluid to be equal to $N\mu$, we obtain the
following equation for the entropy of the adsorbed fluid $S_{loading}$

$S_{loading}={{U-\mu N}\over T}$ (2)

In this equation, $S_{loading}$ denotes the total entropy of the system and,
as such, increases with the total number of atoms $N$ adsorbed in the
nanopore. We add that the $pV$ term is omitted in Eq. 2 since it is negligible
when compared to the other terms ($pV$ actually accounts for less than $0.5$ %
of the lowest values sampled for $TS_{loading}$). This approximation is
similar to the one made by Waghe, Rasaiah and Hummer in their calculation of
the entropy of water adsorbed in carbon nanotubes Waghe et al. (2012). To
sample the transition pathway, we use $S_{loading}$ as the reaction coordinate
and drive the nucleation event through an umbrella sampling approach Torrie
and Valleau (1977). For this purpose, we define the following bias energy:

$U_{b}={1\over 2}k(S_{loading}-S_{0})^{2}$ (3)

in which $S_{0}$ is the target value for the total entropy. This bias energy
is added to the total potential energy of the system, which appears in the
Metropolis criteria used in the Monte Carlo simulations. In practice, for
capillary condensation, we perform a series of $\mu VT-S$ simulations with
increasing values for the target entropy $S_{0}$. This promotes the uptake of
additional Ar atoms since the total entropy of the system is a function of
$N$. Similarly, the pathway for the desorption process is sampled by carrying
out successive $\mu VT-S$ simulations, with decreasing values for the target
entropies $S_{0}$.

Previous work Peterson and Gubbins (1987) has shown that the conditions for
liquid-vapor equilibrium in the pore are obtained for $\mu=-10.53\epsilon$. To
observe the capillary condensation/evaporation processes, the chemical
potential must be close enough to the chemical potential at equilibrium. More
specifically, to simulate the capillary condensation process, we carry out
$\mu VT-S$ simulations for a chemical potential of $\mu=-10.48\epsilon$, while
for the evaporation process, we use a chemical potential of
$\mu=-10.54\epsilon$. These two values for $\mu$ show that our findings are
consistent with the results from Peterson and Gubbins Peterson and Gubbins
(1987). As expected, these two values of $\mu$ bracket the estimate of
Peterson and Gubbins for the chemical potential at the liquid-vapor
equilibrium, with $\mu$ for the capillary condensation process being slightly
above the chemical potential at coexistence and $\mu$ for the capillary
evaporation process being slightly below. A total of $35$ umbrella sampling
windows are performed to connect the vapor and liquid phases for the confined
fluid and to sample the entire condensation and evaporation pathways. For each
window, we first carry out an equilibration run of $1\times 10^{8}$ MC steps
to allow the system to relax towards the target value for the entropy. We then
perform a production run of $2\times 10^{8}$ MC steps, during which averages
are collected collect and structural analyses are run.

##  IV Results and Discussion

We start by discussing the results obtained for the capillary condensation
process. To sample the entropic pathway for this process, we carry out a
series of $\mu VT-S$ simulations at $\mu=-10.48\epsilon$. We gradually
increase the value for the target entropy, as we go from one umbrella sampling
window to another, and, as a result, increase $S_{loading}$ for the confined
fluid. We show in Fig. 1 how the system responds to this increase in the
entropy of the adsorbed fluid. The top panel of Fig. 1 consists of a plot of
the number of Argon atoms $N$ adsorbed in the nanopore as a function of the
reduced entropy of the adsorbed fluid (noted as $S^{*}_{loading}$). We observe
the following trends. First, we find that $N$ increases with
$S^{*}_{loading}$, starting from about $900$ Ar atoms adsorbed for
$S^{*}_{loading}=3300$, and reaches a value of about $1600$ Ar atoms for
$S^{*}_{loading}=6800$. This plot shows that the successive $\mu VT-S$
simulations, with increasing values for the total entropy of the system, allow
to simulate the phase transition process from a metastable phase of low
density to the stable phase of high density. The bottom panel of Fig. 1 shows
the corresponding variation of the interaction energy within the confined
fluid during the $\mu VT-S$ simulations. Overall, we observe a decrease in the
interaction energy as $S^{*}_{loading}$ increases. This behavior is consistent
with the increase in the amount of $Ar$ adsorbed in the cylindrical nanopore
observed when $S^{*}_{loading}$ increases.

Figure 1: Capillary condensation ($\mu=-10.48\epsilon$): variation of the
number of atoms adsorbed $N$ (top) and of the interaction energy for the
confined fluid $U^{*}_{potential}$ (bottom) as a function of the total entropy
of the adsorbed fluid $S^{*}_{loading}$. Figure 2: Free energy profile for the
capillary condensation at $\mu=-10.48\epsilon$. The origin for the energy is
set to $0$ for the starting point (metastable phase of low density - see the
local minimum on the left of the plot). After condensation, the system reaches
a free energy minimum corresponding to the stable phase of high density (see
the minimum on the right of the plot).

We plot in Fig. 2 the free energy profile obtained for the capillary
condensation process. The starting point for the $\mu VT-S$ simulations is a
metastable phase of low density adsorbed in the nanopore, with a reduced
entropy of $S^{*}_{loading}=3300$. As shown on the left of Fig. 2, the free
energy profile exhibits a local minimum for this metastable phase. For
convenience, we have chosen to assign the origin for the free energy to this
local minimum. Several regimes can be identified in this plot of the free
energy as a function of the entropy of the adsorbed phase. First, the free
energy profile exhibits a steep increase, over entropies ranging from
$S^{*}_{loading}=3300$ to $S^{*}_{loading}=4800$. This increase in free energy
corresponds to the steady increase in $N$, and accordingly to the steady
decrease of $U^{*}_{interaction}$, over the same range of entropies. Then, for
$S^{*}_{loading}=4800$ to $S^{*}_{loading}=6000$, the free energy profile
exhibits a second regime which can be characterized as an almost flat top,
that very slowly decreases with $S^{*}_{loading}$. This second regime is
associated with changes in both the slopes for the variations of $N$ and
$U^{*}_{interaction}$ in Fig. 1. Then, for $S^{*}_{loading}=6000$ to
$S^{*}_{loading}=6850$, we observe a third regime for the free energy profile
and a steady decrease in free energy. During this stage, the system completes
the phase transition towards the stable phase of high density. As expected, we
obtain a free energy for the stable phase that is below, by $11$ $k_{B}T$, the
free energy of the metastable phase of low density. This establishes that the
end point for the $\mu VT-S$ simulations is indeed the stable phase. As for
the previous two regimes in the free energy profile, we find that the decrease
for the free energy, observed during the third regime, is connected to sharp
changes in the variations of both $N$ and $U^{*}_{interaction}$ for the
adsorbed phase. This indicates that significant structural changes take place
within the fluid confined in the nanopore.

Figure 3: Capillary condensation ($\mu=-10.48$): density profiles across the
nanopore ($r=0$ denotes the center of the pore) for increasing values of
$S^{*}_{loading}$. The profile in black ($S^{*}_{loading}=3364$) is for the
metastable phase of low density, while the profile in cyan
($S^{*}_{loading}=6859$) corresponds to the stable phase of high density.

To elucidate the mechanism underlying capillary condensation, we carry out a
series of structural analyses. We begin by determining the density profiles
for the adsorbed fluid across the nanopore. The results obtained during
several of the $\mu VT-S$ simulations are shown in Fig. 3. At the start of the
capillary condensation process, i.e. for the metastable phase of low density,
the adsorbed fluid mostly consists of two layers close to the wall, as shown
by the two peaks on the density profile for radii of about $4$ $\sigma$ (fluid
layer closest to the wall) and $3$ $\sigma$ (second layer). As
$S^{*}_{loading}$ increases, $Ar$ atoms start to populate the center of the
nanopore, and a third peak for a radius around $2$ $\sigma$ develops (see the
density profile in red for $S^{*}_{loading}=4290$). We detect the formation of
a fourth and then of a fifth peak when the system reaches the top of the free
energy barrier (see density profile in green). Then, as $S^{*}_{loading}$
further increases, the height of the peaks associated with the innermost fluid
layers continue to increase, until the nanopore has completely filled and the
system reaches the stable phase of high density.

(a)  (b)  (c)

Figure 4: Mechanism underlying capillary condensation. Radial density profile
along the side ($z$) of the nanopore for increasing values of the entropy
$S^{*}_{loading}=3364$ in (a), $S^{*}_{loading}=4929$ in (b) and
$S^{*}_{loading}=6859$ in (c), together with the corresponding snapshots. In
the snapshots, $Ar$ atoms are shown as orange spheres, while the regions of
low density are highlighted in cyan. The bright spots on the left of the
density profile around $z=25$ indicate the formation of the liquid bridge in
(b), as shown by the discontinuity of the cyan region in the corresponding
snapshot.

We now more closely look at the radial density profiles of Fig. 3 to determine
why, for intermediate values of $S^{*}_{loading}$, we find lower peaks found
for the innermost fluid layers. This prompts the following question. Are these
lower peaks associated with a uniform distribution of Ar atoms along the side
of the tube (along the $z$ direction)? Alternatively, do these lower peaks
signal the formation of a liquid contact within the nanotube? To answer this
question, we need to determine how the distribution of atoms evolves along the
side of the cylindrical nanopore, i.e. along the $z$ axis. For this purpose,
we refine our structural analysis by computing the density profiles against
both the radius and the $z$ coordinate and show how the radial density varies
along the nanopore. We show in Fig. 4(a) the results obtained for the
metastable phase of low density. The plot for the density profile shows that
the atoms adsorbed are very predominantly located close to the wall. More
specifically, the fluid is organized into two layers close to the wall as
shown by the two bright lines (i.e. regions of high density) on the density
profile of Fig. 4(a). This is consistent with the profiles of Fig. 3 that
showed the formation of two layers close to the wall. The additional
information brought by Fig. 4(a) is the fact that $Ar$ atoms in these two
layers are distributed uniformly along $z$. We also show in Fig. 4(a) a
snapshot of a configuration of the system for the metastable phase of low
density and highlight in cyan the region where the void dominates (we define
void any region of space for which the reduced density is below a threshold
value). We carry out the same analysis along the rest of the capillary
condensation pathway. At the top of the free energy barrier, we obtain the
density profile and snapshot shown in Fig. 4(b). Both show that, at this
stage, the $Ar$ atoms located in the inner part of the nanopore are not
uniformly distributed along z and that a liquid-like region has developed for
a value of $z$ around $20$ $\sigma$ on the density profile. This can best be
seen on the density profile with the onset, for $z=20$, of 3 bright spots for
$r$ around $0$, $1$ and $2$, that indicate the formation of 3 partial liquid
layers. The snapshot of a configuration of the system at the top of the free
energy barrier (see Fig. 4(b)) confirms that a liquid bridge has developed
across the nanopore for this value of $z=20$. As the capillary condensation
further proceeds, the liquid bridge becomes wider and wider along the
nanopore, resulting in a recess of the void region. Finally, as the system
reaches the stable high density phase, we see that the structure of the
adsorbed fluid exhibits a uniform structure with multiple layers (see Fig.
4(c)). The mechanism observed for capillary condensation is thus consistent
with the scenario proposed by Everett and Haynes Everett and Haynes (1972) for
a macroscopic capillary and the results from gauge cell Monte Carlo
simulations from Neimark et al. Vishnyakov and Neimark (2003b), which both
pointed to a liquid bridge-mediated capillary condensation process. We add
that the free energy barriers are of the same order as those previously
obtained for the capillary condensation process by Vishnyakov and Neimark
Vishnyakov and Neimark (2003b) using the gauge cell method.

Figure 5: Free energy profile for the capillary evaporation at
$\mu=-10.54\epsilon$. The origin for the energy is set to $0$ for the starting
point (metastable phase of high density). After evaporation, the system
reaches a free energy minimum corresponding to the stable phase of low
density.

We now turn to the results obtained for the capillary evaporation process
($\mu=-10.54\epsilon$). We show in Fig. 5 the free energy profile obtained as
the system moves from a metastable phase of high density to the stable phase
of low density. In line with our results for the capillary condensation
process, we set the origin for the free energy axis to the free energy of the
metastable phase. The main difference with the case of the condensation
process, studied in the first part of the paper, is that we now carry out a
series of $\mu VT-S$ simulations with decreasing values for the total entropy
of the system. Therefore, $S^{*}_{loading}$ now varies in the opposite sense
to the one followed for the condensation process. Looking at the free energy
profile as a function of $S^{*}_{loading}$, we find three successive regimes
for the variations of the free energy of the system. First, for
$S^{*}_{loading}$ varying from $7000$ (metastable phase of high density) to
$6500$ (top of the free energy barrier), we observe a steep increase in free
energy. Then, for $S^{*}_{loading}$ between $6500$ and $4800$, the free energy
profile exhibits a much more gradual behavior, as the free energy slowly
decreases with the entropy of the system, in an almost linear fashion. The
third regime is observed for lower values of $S^{*}_{loading}$ (between $4800$
and $2900$). For this range of entropies, the free energy decreases more
rapidly with entropy, until the system reaches the stable phase of low
density. At this point, the free energy of the system is $55$ $k_{B}T$ lower
than the metastable state from which the evaporation pathway started. This
confirms that the system has reached a thermodynamically stable phase at the
end of the $\mu VT-S$ simulations.

Figure 6: Capillary evaporation ($\mu=-10.54\epsilon$): density profiles
across the nanopore ($r=0$ denotes the center of the pore) for decreasing
values of $S^{*}_{loading}$. The profile in black ($S^{*}_{loading}=6980$) is
for the metastable phase of high density, while the profile in cyan
($S^{*}_{loading}=2887$) corresponds to the stable phase of low density.

We now analyze the structure of the adsorbed fluid along the evaporation
pathway. Fig. 6 shows the density profiles across the cylindrical nanopore for
different values of the entropy. The starting point for the $\mu VT-S$
simulations displays the expected multi-peaked density profile. This density
profile is consistent with the organization of the metastable phase of high
density in multiple fluid layers within the nanopore ($S^{*}_{loading}=6980$).
Close to the top of the free energy barrier ($S^{*}_{loading}=4761$), we
observe a small, but noticeable, increase in the number of $Ar$ atoms adsorbed
in the nanopore. This indicates that the capillary evaporation process starts
with the adsorption of a few extra $Ar$ atoms, amounting to about $1-2\%$ of
the overall number of $Ar$ atoms within the system. These extra $Ar$ atoms
destabilize the organization within the nanopore and trigger the desorption
process, as shown by the decrease in the amplitude of the peaks in the density
profile that follows (see e.g. the density profile obtained for
$S^{*}_{loading}=5750$). This leads to a decrease in the free energy as
$S^{*}_{loading}$ further decreases. Then, at the end of the second regime
found in the free energy plot, we observe that the density profile exhibits
two main peaks, the third inner most layer showing a reduced fluid density
($S^{*}_{loading}=4761$). Finally, we recover the structure of the metastable
phase of low density (see the density profile for at $S^{*}_{loading}=2887$),
for which the $Ar$ atoms are very predominantly adsorbed close to the wall of
the cylindrical nanopore.

Figure 7: Bubble formation during capillary evaporation. Radial density
profile along the side ($z$) of the nanopore for $S^{*}_{loading}=5750$ and
its corresponding snapshot. The darker spots on the left of the density
profile around $z=15$ indicate the formation of the bubble, which is
highlighted in cyan in the snapshot.

We finally comment on the local structure within the fluid along the side of
the nanopore. We observe that once the organization of the high density phase
has been destabilized, capillary evaporation proceeds through the formation of
a bubble within the nanopore. This can best be seen in Fig. 7 for
$S^{*}_{loading}=5750$. The density profile of Fig. 7 shows that a dark spot
develops for values of $z$ around 15 and that the innermost layers of the
adsorbed fluid start to be depleted. Furthermore, this plot establishes that a
non-uniform distribution of $Ar$ atoms has formed along the side of the
nanopore. In this case, the dark spots on the density profile on the left of
Fig. 7, as well as the reduction in the height of the peaks for the innermost
layers in Fig. 6, are due to the formation of a bubble, as shown through the
low density region indicated on the snapshot on the right of Fig. 7.

##  V Conclusions

In this work, we carry out $\mu VT-S$ simulations to elucidate the capillary
condensation and evaporation processes. Using the total entropy for the
adsorbate, we are able to drive these processes along the entropic pathways
underlying these processes. This allows us to shed light on the structural
changes that occur within the confined fluid and give rise to the phase
transition. Considering the example of Argon adsorbed in a smooth cylindrical
nanopore typical of the MCM-41 silica adsorbent, we start by focusing on
capillary condensation. Our simulations allow us to identify a complex free
energy profile, corresponding to the successive stages in the condensation
process. Our results show that capillary condensation from a metastable phase
of low density starts with the nucleation of a liquid bridge within the
nanopore. This liquid bridge then expands, becomes wider and wider, as the
total entropy of the system increases, and finally takes over the system to
yield the stable phase of high density. Applying this approach to the
phenomenon of capillary evaporation similarly uncovers a multi-stage process.
In this case, the layered structure of the metastable phase of high density is
first destabilized by the adsorption of a few extra atoms. This is then
followed by the nucleation of a bubble within the nanopore. Overall, the
mechanism identified for capillary condensation is consistent with the
scenario of Everett and Haynes Everett and Haynes (1972) for a macroscopic
system and with gauge cell Monte Carlo simulations by Neimark et al.
Vishnyakov and Neimark (2003b). The results obtained here further establish
the key role played by the nucleation of liquid bridges in the capillary
condensation process, and of the nucleation of bubbles during capillary
evaporation, and shed light on the interplay between the entropy and the
structure of the adsorbed fluid throughout these processes.

Acknowledgements Partial funding for this research was provided by NSF through
CAREER award DMR-1052808. Acknowledgement is made to the Donors of the
American Chemical Society Petroleum Research Fund for partial support of this
research through grant 548002-ND10.

## References

  * Derjaguin (1992) B. Derjaguin, Prog. Surf. Sci. 40, 46 (1992).
  * Coasne and Pellenq (2004) B. Coasne and R.-M. Pellenq, J. Chem. Phys. (2004).
  * Saugey et al. (2005) A. Saugey, L. Bocquet, and J. Barrat, J. Phys. Chem. B 109, 6520 (2005).
  * Major et al. (2006) R. Major, J. Houston, M. McGrath, J. Siepmann, and X.-Y. Zhu, Phys. Rev. Lett. 96, 177803 (2006).
  * Casanova et al. (2007) F. Casanova, C. E. Chiang, C.-P. Li, and I. K. Schuller, Appl. Phys. Lett. 91, 243103 (2007).
  * Puibasset (2008) J. Puibasset, Langmuir 25, 903 (2008).
  * Edison and Monson (2013) J. R. Edison and P. A. Monson, J. Chem. Phys. 138, 234709 (2013).
  * Szoszkiewicz and Riedo (2005) R. Szoszkiewicz and E. Riedo, Phys. Rev. Lett. 95, 135502 (2005).
  * Bruschi et al. (2010) L. Bruschi, G. Mistura, L. Liu, W. Lee, U. Gsele, and B. Coasne, Langmuir 26, 11894 (2010).
  * Greiner et al. (2010) C. Greiner, J. R. Felts, Z. Dai, W. P. King, and R. W. Carpick, Nano Lett. 10, 4640 (2010).
  * Lin et al. (1994) M. Lin, S. Sinha, J. Drake, X.-l. Wu, P. Thiyagarajan, and H. Stanley, Phys. Rev. Lett. 72, 2207 (1994).
  * Willett et al. (2000) C. D. Willett, M. J. Adams, S. A. Johnson, and J. P. Seville, Langmuir 16, 9396 (2000).
  * Gogotsi et al. (2001) Y. Gogotsi, J. A. Libera, A. Güvenç-Yazicioglu, and C. M. Megaridis, Appl. Phys. Lett. 79, 1021 (2001).
  * He et al. (2001) M. He, A. Szuchmacher Blum, D. E. Aston, C. Buenviaje, R. M. Overney, and R. Luginbühl, J. Chem. Phys. 114, 1355 (2001).
  * Heuberger et al. (2001) M. Heuberger, M. Zäch, and N. Spencer, Science 292, 905 (2001).
  * Patel et al. (2002) N. M. Patel, M. R. Dodge, J. I. D. Alexander, L. A. Slobozhanin, P. Taylor, and C. Rosenblatt, Phys. Rev. E 65, 026306 (2002).
  * Maeda and Israelachvili (2002) N. Maeda and J. N. Israelachvili, J. Phys. Chem. B 106, 3534 (2002).
  * Jang et al. (2003) J. Jang, G. C. Schatz, and M. A. Ratner, Phys. Rev. Lett. 90, 156104 (2003).
  * Weeks et al. (2005) B. L. Weeks, M. W. Vaughn, and J. J. DeYoreo, Langmuir 21, 8096 (2005).
  * Berim and Ruckenstein (2008) G. O. Berim and E. Ruckenstein, J. Chem. Phys. 129, 014708 (2008).
  * Thommes et al. (2006) M. Thommes, B. Smarsly, M. Groenewolt, P. I. Ravikovitch, and A. V. Neimark, Langmuir 22, 756 (2006).
  * Fortini and Dijkstra (2006) A. Fortini and M. Dijkstra, J. Phys. Condens. Matt. 18, L371 (2006).
  * Horikawa et al. (2011) T. Horikawa, D. Do, and D. Nicholson, Adv. Colloid Interface Sci. 169, 40 (2011).
  * Monson (2012) P. Monson, Microporous Mater. 160, 47 (2012).
  * Gommes (2012) C. J. Gommes, Langmuir 28, 5101 (2012).
  * Mszr et al. (2013) Z. E. Mszr, G. Hantal, S. Picaud, and P. Jedlovszky, J. Phys. Chem. C 117, 6719 (2013).
  * Zeng et al. (2014) Y. Zeng, C. Fan, D. Do, and D. Nicholson, J. Phys. Chem. C 118, 3172 (2014).
  * Hiratsuka et al. (2016) T. Hiratsuka, H. Tanaka, and M. T. Miyahara, ACS Nano (2016).
  * Restagno et al. (2000) F. Restagno, L. Bocquet, and T. Biben, Phys. Rev. Lett. 84, 2433 (2000).
  * Talanquer and Oxtoby (2001) V. Talanquer and D. Oxtoby, J. Chem. Phys. 114, 2793 (2001).
  * Ustinov and Do (2005) E. Ustinov and D. Do, J. Phys. Chem. B 109, 11653 (2005).
  * Men et al. (2009) Y. Men, Q. Yan, G. Jiang, X. Zhang, and W. Wang, Phys. Rev. E 79, 051602 (2009).
  * Zhang and Chakrabarti (1994) Z. Zhang and A. Chakrabarti, Phys. Rev. E 50, R4290 (1994).
  * Gelb and Gubbins (1997) L. D. Gelb and K. Gubbins, Phys. Rev. E 56, 3185 (1997).
  * Landman (1998) U. Landman, Solid State Commun. 107, 693 (1998).
  * Leung and Luzar (2000) K. Leung and A. Luzar, J. Chem. Phys. 113, 5845 (2000).
  * Yasuoka et al. (2000) K. Yasuoka, G. Gao, and X. C. Zeng, J. Chem. Phys. 112, 4279 (2000).
  * Liu and Grest (1991) A. J. Liu and G. S. Grest, Phys. Rev. A 44, R7894 (1991).
  * Gac et al. (1994) W. Gac, A. Patrykiejew, and S. Sokołowski, Surf. Sci. 306, 434 (1994).
  * Gelb et al. (1999) L. D. Gelb, K. Gubbins, R. Radhakrishnan, and M. Sliwinska-Bartkowiak, Rep. Prog. Phys. 62, 1573 (1999).
  * Bock and Schoen (1999) H. Bock and M. Schoen, Phys. Rev. E 59, 4122 (1999).
  * Bolhuis and Chandler (2000) P. G. Bolhuis and D. Chandler, J. Chem. Phys. 113, 8154 (2000).
  * Stroud et al. (2001) W. J. Stroud, J. E. Curry, and J. H. Cushman, Langmuir 17, 688 (2001).
  * Liu and Monson (2006) J.-C. Liu and P. Monson, Ind. Eng. Chem. Res. 45, 5649 (2006).
  * Mota and Esteves (2007) J. P. Mota and I. A. Esteves, Adsorption 13, 21 (2007).
  * Winkler et al. (2010) A. Winkler, D. Wilms, P. Virnau, and K. Binder, J. Chem. Phys. 133, 164702 (2010).
  * Nguyen et al. (2011) V. T. Nguyen, D. Do, and D. Nicholson, J. Phys. Chem. B 115, 7862 (2011).
  * Gor et al. (2012) G. Y. Gor, C. J. Rasmussen, and A. V. Neimark, Langmuir 28, 12100 (2012).
  * Siderius and Shen (2013) D. W. Siderius and V. K. Shen, J. Phys. Chem. C 117, 5861 (2013).
  * van Meel et al. (2015) J. van Meel, Y. Liu, and D. Frenkel, Mol. Phys. 113, 2742 (2015).
  * Kornev et al. (2002) K. G. Kornev, I. K. Shingareva, and A. V. Neimark, Adv. Colloid Interface Sci. 96, 143 (2002).
  * Neimark et al. (2003) A. V. Neimark, P. I. Ravikovitch, and A. Vishnyakov, J. Phys. Condens. Matt. 15, 347 (2003).
  * Vishnyakov and Neimark (2003a) A. Vishnyakov and A. V. Neimark, Langmuir 19, 3240 (2003a).
  * Vishnyakov and Neimark (2003b) A. Vishnyakov and A. V. Neimark, J. Chem. Phys. 119, 9755 (2003b).
  * Neimark and Vishnyakov (2005a) A. V. Neimark and A. Vishnyakov, J. Chem. Phys. 122, 174508 (2005a).
  * Neimark and Vishnyakov (2005b) A. V. Neimark and A. Vishnyakov, J. Chem. Phys. 122, 054707 (2005b).
  * Neimark and Vishnyakov (2005c) A. V. Neimark and A. Vishnyakov, J. Phys. Chem. B 109, 5962 (2005c).
  * Everett and Haynes (1972) D. Everett and J. Haynes, J. Colloid Interface Sci. 38, 125 (1972).
  * Oxtoby and Evans (1988) D. W. Oxtoby and R. Evans, J. Chem. Phys. 89, 7521 (1988).
  * ten Wolde and Frenkel (1998) P. R. ten Wolde and D. Frenkel, J. Chem. Phys. 109, 9919 (1998).
  * Ten Wolde et al. (1999) P. R. Ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J. Chem. Phys. 110, 1591 (1999).
  * Shen and Debenedetti (1999) V. K. Shen and P. G. Debenedetti, J. Chem. Phys. 111, 3581 (1999).
  * Desgranges and Delhommelle (2011) C. Desgranges and J. Delhommelle, J. Am. Chem. Soc. 133, 2872 (2011).
  * Patel et al. (2011) A. J. Patel, P. Varilly, D. Chandler, and S. Garde, J. Stat. Phys. 145, 265 (2011).
  * Desgranges and Delhommelle (2014) C. Desgranges and J. Delhommelle, J. Am. Chem. Soc. 136, 8145 (2014).
  * Lauricella et al. (2015) M. Lauricella, S. Meloni, S. Liang, N. J. English, P. G. Kusalik, and G. Ciccotti, J. Chem. Phys. 142, 244503 (2015).
  * Desgranges and Delhommelle (2016a) C. Desgranges and J. Delhommelle, J. Chem. Phys. 145, 204112 (2016a).
  * Desgranges and Delhommelle (2016b) C. Desgranges and J. Delhommelle, J. Chem. Phys. 145, 234505 (2016b).
  * Torrie and Valleau (1977) G. M. Torrie and J. P. Valleau, J. Comput. Phys. 23, 187 (1977).
  * Ravikovitch et al. (1997) P. Ravikovitch, D. Wei, W. Chueh, G. Haller, and A. Neimark, J. Phys. Chem. B 101, 3671 (1997).
  * Allen and Tildesley (1987) M. P. Allen and D. J. Tildesley, _Computer Simulation of Liquids_ (Clarendon, Oxford, 1987).
  * Waghe et al. (2012) A. Waghe, J. C. R. Rasaiah, and G. Hummer, J. Chem. Phys. 137, 044709 (2012).
  * Peterson and Gubbins (1987) B. K. Peterson and K. E. Gubbins, Mol. Phys. 62, 215 (1987).


# Exploring the link between coffee matrix microstructure and flow properties
using combined X-ray microtomography and smoothed particle hydrodynamics
simulations

Chaojie<EMAIL_ADDRESS>Basque Center for Applied Mathematics (BCAM),
Alameda de Mazarredo 14, 48009 Bilbao, Spain  Aircraft and Propulsion
Laboratory, Ningbo Institute of Technology, Beihang University, Ningbo 315100,
P. R. China  Richard Johnston  Faculty of Science and Engineering, Swansea
University, Swansea, SA1 8EN, UK  Luciano Navarini  Illycaffè S.p.A, Via
Flavia 110, Trieste 34147, Italy  Furio Suggi Liverani  Illycaffè S.p.A, Via
Flavia 110, Trieste 34147, Italy  Marco Ellero  Basque Center for Applied
Mathematics (BCAM), Alameda de Mazarredo 14, 48009 Bilbao, Spain  Zienkiewicz
Centre for Computational Engineering (ZCCE), Swansea University, Bay Campus,
Swansea SA1 8EN, UK  IKERBASQUE, Basque Foundation for Science, Calle de María
Díaz de Haro 3, 48013 Bilbao, Spain

## Abstract

Coffee extraction involves many complex physical and transport processes
extremely difficult to model. Among the many factors that will affect the
final quality of coffee, the microstructure of the coffee matrix is one of the
most critical ones. In this article, we use X-ray micro-computed (microCT)
technique to capture the microscopic details of coffee matrices at particle-
level and perform fluid dynamics simulation based on the smoothed particle
hydrodynamics method (SPH) with the 3D reconstructured data. Information like
flow permeability and tortuosity of the matrices can be therefore obtained
from our simulation. We found that inertial effects can be quite significant
at the normal pressure gradient conditions typical for espresso brewing, and
can provide an explanation for the inconsistency of permeability measurements
seen in the literature. Several types of coffee powder are further examined,
revealing their distinct microscopic details and resulting flow features. By
comparing the microCT images of pre- and post-extraction coffee matrices, it
is found that a decreasing porosity profile (from the bottom-outlet to the
top-inlet) always develops after extraction. This counterintuitive phenomenon
can be explained using a pressure-dependent erosion model proposed in our
prior work. Our results reveal not only some important hydrodynamic mechanisms
of coffee extraction, but also show that microCT scan can provide useful
microscopic details for coffee extraction modelling. MicroCT scan establishes
the basis for a data-driven numerical framework to explore the link between
coffee powder microstructure and extraction dynamics, which is the
prerequisite to study the time evolution of both volatile and non-volatile
organic compounds and then the flavour profile of coffee brews.

##  1 Introduction

A roasted and ground coffee matrix is a highly complex porous medium
characterized by multiscale features. A typical coffee matrix used for double-
shot espresso brew is about 15 ml in volume and consists of millions of coffee
particles. These coffee particles are produced by grinding roasted coffee
beans, and have an approximately bimodal size distribution[1]. The first peak
of the distribution profile is always at around $30\sim 40\ \mathrm{\mu m}$,
representing the so-called “fines”, i.e., inner walls cellular fragments (to
be very rigorous, this first peak is preceded by a small submicron fraction).
The second peak represents the so-called “coarses” and their mean size depends
on the grinder, with values ranging between $200\sim 1000\ \mathrm{\mu m}$.
Moreover, the coffee particles are not entirely solid but porous material
themselves, characterized by many cell-pockets of dimensions $30\sim 60\
\mathrm{\mu m}$ inside the particles[1]. Finally, the cell-walls are also
porous at nanoscale level.

During coffee extraction, as water flows through this complex porous
structure, many chemical compounds including among many others caffeine,
trigonelline, polyphenols and polysaccharides, are diffused and convected out
of the solid matrix to finally enter into the cup and determine the quality
(e.g. flavour and mouthfeel) of the final beverage.

This coffee extraction process is especially difficult to model. This is not
simply because of the complex topological features related to the multiscale
nature of coffee bed geometry as explained above, there are also many other
physical processes concomitant to the simple diffusion and convection. For
example, the porous medium undergoes geometric changes due to the fine
migration[2], the particle erosion[3, 4], and particle swelling[3, 5]; there
is also solubilization of many hydrophilic substances[6], the emulsification
of insoluble coffee oils[7], suspension of solid coffee cell-wall fragments
(fines)[7], and $\mathrm{CO}_{2}$ degassing and supersaturation[8]. These
physical processes are all interconnected, rendering the modelling of coffee
extraction very challenging.

In order to predict the coffee quality and propose better brewing protocol a
priori, many works have tried to model the extraction process using a variety
of simulation techniques. We mention early attempts using cellular automata
based model for percolation simulation[2, 9]. More recently, Lattice Bolztmann
Method (LBM) has been also used to simulate the extraction kinetics[10] and to
consider the swelling and erosion effects[4]. Multiscale models have also been
proposed. In the work of Moroney et al[11] an upscaling procedure is applied
to relate the conservation equations at different scales. The developed
multiscale model is able to quantitatively reproduce the experimental
extraction profile after being parameterized by experimentally measured coffee
bed parameters. In the work of Melrose et al[12] and Cameron et al[13],
another multiscale model, which can incorporate the bimodal characteristics of
the coffee particle, is proposed. In this model the macroscopic mass transport
is modelled using a 1-dimensional convection-diffusion equation, while the
microscopic extraction is modelled through diffusion equation inside
representative coffee particles in different layers of the coffee matrix. The
equations at different scale are coupled through source terms and boundary
conditions. This model is also found to be able to capture the extraction
kinetics. Besides these models, mesoscopic model based on smoothed dissipative
particle dynamics (SDPD) has been also recently proposed to incorporate the
fine migration[14], particle erosion[15] and swelling effect[16] during coffee
extraction.

However, in most of the aforementioned modelling and simulations, the fluid
dynamics through the coffee matrix is either characterized by average porosity
and permeability[10, 11, 12, 13] or is resolved by flows through artificial
packing of many spheres/disks[14, 4, 15, 16]. These simulations may not be
able to capture some important features of coffee extraction. For example, the
coffee particles size is usually of smooth bimodal distribution[1] and the
particles are not spheric[3], therefore many details are lost if the simple
Carman-Kozeny relation and Darcy’s law are used to characterize the flow, or
if some synthetic porous medium formed by packing ideal spheres/disks are used
for percolation simulations. In fact, it is known that the particle size
distribution is quite influential to the final quality of coffee[17].
Therefore, it is important to accurately model the whole spectrum of particle
size in coffee extraction simulation. A straightforward high-resolution
mapping - capturing all the microscopic details of a coffee matrix - can be
obtained by using X-ray microtomography (microCT)/X-ray microscopy (XRM) to
scan the coffee matrix, and build a porous medium model using the
3D-reconstructured inner topology[18, 19]. In this way, not only the full
particles size distribution can be considered, but the spatial permeability
fluctuation in the matrix can also be captured. The spatial permeability
fluctuation could correlate to the “partially clogged” phenomenon which may
cause very fine coffee powder to have less extraction rate[13].

With X-ray microCT, it is also possible to scan both pre-extraction (dry)
coffee matrix and post-extraction (wet) coffee matrix. The comparison could
reveal some information about how the porous structure changes during
extraction process. This information could provide important guidance to model
different physical processes that lead to geometric change and matrix
restructuring during extraction. As far as we are able to determine, microCT
reconstruction of coffee matrix microstructure has not been used for direct
coffee extraction simulation. This will be the subject of the current article.

In this article, we will use microCT reconstructured data to build a static
porous medium model, and simulate the percolation process using the smoothed
particle hydrodynamics method (SPH). We determine the simulated medium
permeability at different pressure gradient and compare them to experimental
data in literatures. We show that inertial effects are quite influential at
the pressure gradients typically considered in espresso brewing. As a result,
the inconsistency of many results of coffee matrix permeability found in
literature can be explained by the above-mentioned inertial effect. Several
types of coffee powder are also examined, and their permeability and
tortuosity are measured systematically using SPH simulation. This allows to
establish a link between coffee powder granulometry and other microstructural
features with extraction efficiency. Finally, microCT images of post-
extraction coffee matrices are also analyzed, it is revealed that a decreasing
porosity profile (from the bottom-oulet to the top-inlet) always develops
after extraction. This counterintuitive phenomenon can be explained using a
specific pressure-dependent erosion model proposed in our previous work[15].

The article is organized as follows: In section 2 we introduce the coffee
powders we study (section 2.1), the particle size analysis method (section
2.2), the microtomography techinique (section 2.3), and the SPH-based
percolation simulation technique (section 2.4). In section 3, the powder
granulometry is first presented (section 3.1). Then we show a brief analysis
of the microCT images in section 3.2. Discussion of percolation and tortuosity
of the coffee matrix is done in section 3.3. The analysis of the post-
extraction microCT images is reported in section 3.4. Finally conclusions are
drawn in section 4.

##  2 Materials and methods

###  2.1 Coffee powders

Four different types of ground coffee obtained by the same medium roasted 100%
Coffea arabica blend were used: Espresso (type E) characterized by the typical
size distribution for traditional espresso brew preparation, IPSO (type H)
typical for capsule espresso brew preparation , Refilly Moka (type M) typical
for moka-pot preparation, and PillowPack (type F) typical for drip-filter
preparation. In the case of type H dark roasting degree was also used. The
coffee powders have been packaged in single dose capsule for espresso brew
(Iperespresso®). The capsule (single shot espresso) is approximately
cylindrical around 14.6 mm in height by 32.5 mm in width. Each capsule
contains 6.7$\pm$0.1 gram of coffee powder. The capsules used in our
experiment are sent by illycaffè S.p.A.

###  2.2 Particle size analysis

The powders are examined with laser diffraction to obtain the particle size
characteristics. A Mastersizer 3000 (Malvern Instrument, UK) is used.

###  2.3 X-ray microtomography

X-ray 3D imaging is used to obtain the internal microscopic structure of the
coffee matrix. A Zeiss Xradia Versa 520 (Carl Zeiss XRM, Pleasanton, CA, USA)
was used to carry out high-resolution XRM non-destructive imaging; this was
achieved using a CCD (charge coupled device) detector system with
scintillator-coupled visible light optics and a tungsten transmission target.
3D scans were performed with an X-ray tube voltage of 80 kV, a tube current of
87 mA, an exposure of 3000 ms. A total of 3201 projections were collected. An
objective lens giving an optical magnification of 0.4x was selected with
binning set to 1, producing an isotropic voxel (three-dimensional pixel) size
of around $16.8\mathrm{\mu m}\times 16.8\mathrm{\mu m}\times 16.8\mathrm{\mu
m}$. The tomograms were reconstructed from 2D projections using a Zeiss
Microscopy commercial software package (XMReconstructor) and an automatically
generated cone-beam reconstruction algorithm based on the filtered back-
projection. The powder are in Illy Hyper Espresso capsule (see detailes
above).

Segmentation is acquired by choosing a global greyscale threshold to binarise
the microCT imaging. However, the resolution of the microCT scan is of the
same order as the smaller fines, and some of these fines are occupying the air
space, therefore the greyscale of the air space is actually a function of the
partial volume of air and fines. Furthermore, the coffee particles also have
many internal pores of the same scale as the microCT resolution. So the
grayscale of the coffee is also a function of the partial volume of solid and
air inside the coffee particles. It is therefore quite challenging to
determine the proper greyscale threshold distinguishing the coffee and air
space. Moreover, the greyscale distribution of the microCT imaging of the
coffee bed is unimodal and of Gaussian shape, and conventional segmentation
techniques[20, 21] typically fails on these images. Here we determine the
greyscale threshold by matching the bed porosity from the microCT imaging and
that from Pycnometer data. The procedure is as follows:

  1. 1.

The microCT images are inverted. In the original image brighter pixels
represent solid. We first invert the values of the pixels.

  2. 2.

Assuming the greyscale threshold is known as $T_{g}$, the microCT imaging are
binarized using this threshold.

  3. 3.

Then a connected-region labeling algorithm[22] is used to identify the
connected and isolated voids. These clusters of void will consist of one major
cluster with very large volume and many minor clusters with small volume.
Apparently, the major cluster represents the inter-granular pores. The
isolated minor clusters are closed pores inside the grains. We can get the bed
porosity, closed porosity, and the total porosity from the volume of these
clusters.

  4. 4.

To avoid dealing with the cylindrical surface of the coffee cake, we extract a
cuboid of $600\times 600\times 600\mathrm{voxels}$ ($1\times 1\times
1\mathrm{cm}$) from each of the microCT imaging of the coffee cakes. Varying
$T_{g}$, the porosity of this cuboid sample at different threshold can be
determined. An example obtained by processing the microCT imaging of type E
powder is shown in Fig. 1. As can be seen from the figure, the total porosity
is almost the same as the bed porosity, the total volume of the intergranular
pores is also orders higher than the intragranular closed pores, indicating
that most of the closed pores are not captured in our microCT imaging. This is
reasonable since the resolution of the microCT imaging is $16.759\mathrm{\mu
m^{3}/voxel}$ and the intragranular pores are of the same scale. The
resolution is not high enough to accurately capture the intragranular pores.
Therefore the porosity consists overwhelmingly of the bed porosity of the
coffee cake. By matching the bed porosity obtained from binarized images with
the typical bed porosity of coffee matrix available in literatures[23, 13], we
can determine the correct grayscale threshold and use that threshold to
binarize the microCT imaging. In the example presented in Fig. 1, the porosity
to be matched is 0.17, the grayscale threshold can be determined to be 200.
Then this grayscale threshold can be used to binarize the microCT image.

Figure 1: (a) Bed porosity and total porosity of the cuboid sample at
different grayscale thresholds. The dash black line mark the typical bed
porosity for a realistic coffee matrix. (b) Volume of the intergranular pores
and the intragranular closed pores at different grayscale thresholds.

###  2.4 Percolation simulation using SPH

After the segmentation, we have a label matrix from the voxel model, with the
solid voxels labelled by 1, and void/air voxels labeled by 0. For the
simulations we use Smoothed Particle Hydrodynamics (SPH) method (see
Supplementary Material). A simulation box of specified size is filled with SPH
particles in specified number density and in some lattice form. By mapping
this simulation box to a 3-D region inside the voxel model, the label matrix
can be interpolated to classify each SPH particle as fluid particle or solid
particle. The construction of SPH particle model from microCT imaging is
finished after the classification.

An example of SPH particle model generated with the above procedure is shown
in Fig. 2. In this example a cube of length $L=708\mathrm{\mu m}$ is extracted
from the microCT imaging and used for the above procedure. Buffer fluid
regions at both the top and bottom have been added. The SPH particle density
is specified to be the same as the number density of the voxels, and the
particles are at face-centered cubic lattice points.

Figure 2: An example of SPH particle model with $L_{x}=L_{z}=708\mathrm{\mu
m}$, $L_{y}=991\mathrm{\mu m}$. Blue particles: solid grains, red particles:
fluid.

With the SPH particle model ready, we can set periodic boundary condition at
all boundaries, and applied body force on the $y$ direction to drive fluid to
percolate through the porous medium. The readers are referred to the
Supplementary Material for the details of the SPH formalism and the SPH
parameters used in our simulations.

##  3 Results and discussion

###  3.1 Coffee powder granulometry

Figure 3: (a) Density distribution of the particle size for different types of
ground coffee. (b) Cumulative distribution of the particle size.

The particle distributions for the coffee powders measured with laser
diffraction are shown in Fig. 3. The density distributions are typically
bimodal, with the first peak representing the fine particles and the second
peak the coarse particles. The first peak is always around $25\sim
50\mathrm{\mu m}$, because the fines are cell fragments and the cell size is
in the range $25\sim 50\mathrm{\mu m}$. While the position of the coarse peak
will shift right with coarser grind. More results from the laser diffraction
measurement, including the uniformity, specific surface area, Sauter diameter
(d[3,2]), De Brouckere mean diameter (d[4,3]), and the volume fraction of
particles below $100\mathrm{\mu m}$, are summarised in Tab. 1. It can be seen
from both Fig. 3 and Tab. 1 that from type E to type F, the volume fraction of
fines is decreasing. If we assume that all particles smaller than
$\mathrm{100\mu m}$ are fines then there are 29.21% of fines in type E powder
but only 9.7% of fines in type F powder. The mean diameter also increases from
type E to type F, and the uniformity and specific surface area decreases
concomitantly.

Table 1: Results from laser diffraction on different types of coffee powder Type | Uniformity | Specific Surface Area ($\mathrm{m^{2}/kg}$) | d[3,2] $\mathrm{\mu m}$ | d[4,3] $\mathrm{\mu m}$ | Volume below 100 $\mathrm{\mu m}$
---|---|---|---|---|---
E | 0.646 | 74.09 | 81.0 | 308 | 29.21%
H | 0.600 | 66.60 | 90.1 | 341.6 | 27.52%
M | 0.508 | 51.66 | 116 | 372 | 20.26%
F | 0.473 | 28.05 | 214 | 673 | 9.7%

###  3.2 Xray microCT analysis

Figure 4: A vertically-cut section of the microCT imaging for the type H
coffee powder.

A vertical digital section of the microCT/XRM imaging data of type H powder is
shown in Fig. 4 (Enlarged cutting section images for all the powders are
provided in the Supplementary Materials). Both the coffee grains and the
plastic/metal structure skeletons of the capsule are clearly visible in the
image.

The microCT imaging is first binarised following the segmentation protocol
introduced above. We show some cropped sections of the microCT imaging before
and after segmentation in Figure 5. Note that in this figure the images before
segmentation have already had their pixels inverted so that brighter regions
represent voids. As can be seen from the images, the coffee grains are
becoming coarser and coarser from type E to type F. This is in agreement with
our particle size measurement.

Figure 5: Cropped horizontally-cut sections of the microCT imaging before and
after segmentation. (a)-(c) Type E powder; (d)-(f) Type H powder; (g)-(i) Type
M powder; (j)-(l) Type F powder. The first column shows the original microCT
images. The second column shows the histogram of the grayscale. The red lines
mark the grayscale thresholds determined by matching the bed porosity. The
third column shows the thresholded images.

Then the binarised images are used to analyse the spatial distribution of
porosity of the coffee cakes. Because the coffee cakes are of cylinder shape,
we divide the cakes into many annular section with the same radial thickness
and calculate the volume fraction of the void in each vertical layer, so that
the variation curves of porosity against the vertical coordinate at different
radial distance to the centre are acquired as shown in Fig. 6. It can be seen
that with the decreasing volume fraction of fines, from (a) to (d) the curves
fluctuates more significantly, but in each capsule the porosity is relatively
homogenous across the matrix. In most cases, there is no significant variation
trend in the radial direction. Except that for type H capsule, there is an
apparent trend of higher porosity near the top center. This could be a result
of different transportation condition for the capsule. In all cases, there is
no significant variation trend in the vertical direction.

Figure 6: Local porosity of the coffee cakes at different position. $y$ is the
vertical position, $r$ is the radial distance to the center of the cake. (a)
Type E; (b) Type H; (c) Type M; (d) Type F.

###  3.3 Permeability and tortuosity

Figure 7: Flow rate at different pressure gradients and the fit using Darcy’s
law. The pressure gradient is relatively low here. (Type H powder)

As a start, we investigate the permeability of type H powder in detail. The
porous flow is simulated with SPH. A cropped cubic sample of length
$L=708\mathrm{\mu m}$ is first examined. The SPH particle model has been shown
in Fig. 2. Fluid buffer regions with a total height
$L_{\mathrm{buf}}=283\mathrm{\mu m}$ are added in the $y$ direction. When a
body force $g$ is applied in the $-y$ direction, the effective pressure
gradient in the sample region will be $\rho g(1+L_{\mathrm{buf}}/L)$. Figure 7
shows the relation of the effective pressure gradient and the resultant flow
rate in $-y$ direction. The relation is linear and the fit using Darcy’s law:

$\nabla P=-\frac{\mu}{k}q$ (1)

gives a permeability $k=5.51\times 10^{-13}\ \mathrm{m^{2}}$. This value is
smaller than the measurements by Gianino[24] ($k=2.3\times 10^{-12}\
\mathrm{m}^{2}$) and King[25]($k=8.98\times 10^{-13}\ \mathrm{m}^{2}$). But it
is larger than the estimation given by Navarini et al[26] ($k=0.7\sim 4\times
10^{-13}\ \mathrm{m^{2}}$) and Corrochano et al[27] ($k=2.59\times
10^{-14}\sim 3.36\times 10^{-13}\ \mathrm{m^{2}}$). The disagreement is
probably caused by the different overpressure ($\Delta P$ across the coffee
bed), as in Gianino’s[24] and King’s[25] experiments the overpressure is only
in the order of $0.01\ \mathrm{bar}$, while in the experiments of Navarini et
al[26] and Corrochano et al[27] the overpressure is in the order of $1\
\mathrm{bar}$. Here in our SPH model (Fig. 2), the pressure gradient $8\sim
40\ \mathrm{bar/m}$ represents an overpressure between $0.16\ \mathrm{bar}$
and $0.8\ \mathrm{bar}$ on a coffee bed of $2\ \mathrm{cm}$ height. So the
overpressure is intermidiate here, resulting in a permeability higher than
that of Navarini et al[26] and Corrochano et al[27]. When the pressure
gradient is significantly increased as shown in Fig. 8 (a) the permeability
obtained using Darcy’s law decreases with the increasing pressure head,
achieving a better agreement with Navarini’s[26] and Corrochano’s[27] results.
In Fig. 8 (a) the overpressure is in $1\sim 8\ \mathrm{bar}$ when the pressure
gradient increase from $50\ \mathrm{bar/m}$ to $400\ \mathrm{bar/m}$ if a
coffee bed of $2\ \mathrm{cm}$ height is considered. With the increase of
pressure head the permeability decreases from about $5.2\times 10^{-13}$ to
$2.5\times 10^{-13}\ \mathrm{m^{2}}$. This kind of pressure dependence was
also reported in Navarini’s experimental investigation[26]. In the experiment
the dependence of permeability on pressure gradient might be attributed to the
deformation of the bed under pressure. But in our simulation the porous coffee
bed is static and not deformable, and the dependence is still observed, so the
dependence can not be all attributed to the deformation of the porous bed
under pressure. Here, we believe that the permeability decreases with higher
pressure as a result of the flow inertial effects. In Fig. 8 (b), the flow
rate at different pressure gradients and the fit using Darcy’s law with
Forchheimer’s modification[28] (Eq. 2) is presented.

$\nabla P=-\frac{\mu}{k}q-\frac{\rho}{k_{1}}|q|q$ (2)

The fit is very good and gives a permeability ($k=8.0\times 10^{-13}\
\mathrm{m^{2}}$) and an additional inertial permeability $k_{1}=2.17\times
10^{-13}\ \mathrm{m^{2}}$. The inertial effects cause the flow rate to
increase in a rate slower than linear relation, so when the Darcy’s law is
used it gives a permeability decreasing with an increasing pressure gradient.
If we use the $d_{4,3}=341.6\ \mathrm{\mu m}$ as the characteristic length,
the Reynolds number ($\mathrm{Re}=d_{4,3}\rho q/\mu$) in Fig. 8 is between
0.84 to 3.86, indicating that the inertia is indeed important at high pressure
gradient.

Figure 8: (a) When the pressure gradient is relatively high, the permeability
obtained from Darcy’s law is dependent on the pressure gradient. (b) Flow rate
at different pressure gradients and the fit using Darcy’s law with
Forchheimer’s modification. (Type H powder)

Assuming that there is no recirculation, the tortuosity of the excerpt coffee
bed sample can be determined by tracing the trajectories of the fluid
particles:

$\tau=\frac{\sum_{i}^{N}\int_{\mathbf{x}_{i}\in\Omega_{c}}||\mathrm{d}\mathbf{x}_{i}||}{\left|\sum_{i}^{N}\int_{\mathbf{x}_{i}\in\Omega_{c}}\mathrm{d}\mathbf{x}^{\beta}_{i}\right|},$
(3)

where $\tau$ is the tortuosity, $N$ is the total fluid particles, $\Omega_{c}$
is the coffee sample domain, $\mathbf{x}_{i}$ is the trajectory of particle
$i$, $\beta$ is the direction of the filtration.

For one of the simulations in Fig. 8, the trajectories of 3706 labeled
particles whose initial positions are at the top of the upper buffer fluid
region, are shown in Fig. 9. Because all boundaries are periodic and the
coordinates have been unwrapped, the trajectories span outside the simulation
box. Using Eq. 3 the tortuosity can be determined to be $\tau=1.80$. The power
law relationship has been frequently used to describe the tortuosity-porosity
correlation[29]:

$\tau=\left(\frac{1}{\epsilon}\right)^{n},$ (4)

where $n$ is a parameter dependent on the packing properties. The porosity of
the sample is known to be 0.17, using Eq. 4, it can be determined that
$n=-\log\tau/\log\epsilon=0.33$.

Figure 9: Trajectories of fluid particles through the coffee bed sample. The
arrow marks the direction of the pressure gradient. The black frame is the
simulation box with dimension $L_{x}=L_{z}=708\mathrm{\mu m}$,
$L_{y}=991\mathrm{\mu m}$, all boundaries are periodic, all coordinates have
been unwrapped. (Type H powder)

The results shown above are only obtained from a small sample of the whole
coffee bed. The bed is however quite inhomogenous, and the porosity can vary
significantly across the whole bed in both radial and axial direction as shown
in Fig. 6. Therefore, to properly characterize the property of the coffee bed
we extract more samples from the bed and perform simulations on them. Results
and the locations of the samples are summarized in Tab. 2 and Fig. 10. In Fig.
10, the discharge - pressure gradient relation for different samples have been
fitted using both Darcy’s law and Forchheimer’s modification, in low $\nabla
P$ regime and relatively higher $\nabla P$ regime, respectively. All the
samples show significant deviation from Darcy’s law in the relatively high
$\nabla P$ regime. But the Forchheimer’s formula accurately captures the
$q-\nabla P$ relationship in the high $\nabla P$ regime. Note that in Fig. 10
the transitions from Darcy to Forchheimer do not correspond to one single
critical Reynolds number. This phenomenon has also been reported by other
researchers in the field of hydrology: coarser particles lead to higher
critical Re [30]. The permeabilities obtained using Darcy’s law in the low
$\nabla P$ regime are presented in Tab. 2 as $k^{\mathrm{D}}$, the
permeabilities obtained using Forchheimer’s modification in the relatively
high $\nabla P$ regime are presented in Tab. 2 as $k^{\mathrm{F}}$ and
$k^{\mathrm{F}}_{1}$ for inertia permeability. The porosity and tortuosity for
each sample are also presented in Tab. 2. It can be seen that both the
porosity and tortuosity fluctuate slightly, and there is a strong correlation
between the porosity and the measured $k^{\mathrm{D}}$. Higher porosity tends
to lead to larger $k^{\mathrm{D}}$.

Table 2: Summary of the properties of the coffee samples (type H powder). The parameter $k^{\mathrm{D}}$ is the permeability obtained using Darcy’s law in the low $\nabla P$ regime, $k^{\mathrm{F}}$ and $k_{1}^{\mathrm{F}}$ are the permeability obtained using Forchheimer’s formula in the high $\nabla P$ regime. Sample | Location (mm) | Porosity | Tortuosity |  $k^{\mathrm{D}}$ ($10^{-13}\mathrm{m}^{2}$) |  $k^{\mathrm{F}}$ ($10^{-13}\mathrm{m}^{2}$) |  $k^{\mathrm{F}}_{1}$ ($10^{-9}\mathrm{m}^{2}$)
---|---|---|---|---|---|---
1 | (1.68, 20.35, 10.06) | 0.170 | 1.80 | 5.59 | 7.95 | 4.21
2 | (1.68, 16.76, 0) | 0.183 | 1.91 | 2.58 | 9.15 | 0.41
3 | (0, 14.45, 3.35) | 0.226 | 1.56 | 35.1 | 67.1 | 68.9
4 | (-3.35, 13.64, 6.73) | 0.162 | 1.91 | 7.99 | 14.9 | 5.77
5 | (0, 16.99, 0) | 0.180 | 1.65 | 14.3 | 24.9 | 17.5
6 | (-1.68, 18.67, 1.68) | 0.180 | 1.78 | 18.6 | 29.4 | 39.3
Figure 10: Discharge at different pressure gradient of different excerpt
coffee samples. The dash lines are the fittings using Darcy’s law, the solid
lines are the fittings with Forchheimer’s formula. (Type H powder)

Percolation simulations are also run on samples from type E, M and F powder.
The results are summarized in Tab. 3-5. The permeability $k^{D}$ for different
types of powder at low pressure gradient is also shown in Fig. 11. It can be
seen that, with the coffee particle becoming coarser the average permeability
increases, and the variance of the permeability also increases.

Table 3: Summary of the properties of the coffee samples (type E powder). The parameter $k^{\mathrm{D}}$ is the permeability obtained using Darcy’s law in the low $\nabla P$ regime, $k^{\mathrm{F}}$ and $k_{1}^{\mathrm{F}}$ are the permeability obtained using Forchheimer’s formula in the high $\nabla P$ regime. Type E powder is the finest powder, therefore in some regions the flow is nearly clogged. Some parameters cannot be obtained and we mark them with dash in the table. Sample | Location (mm) | Porosity | Tortuosity |  $k^{\mathrm{D}}$ ($10^{-13}\mathrm{m}^{2}$) |  $k^{\mathrm{F}}$ ($10^{-13}\mathrm{m}^{2}$) |  $k^{\mathrm{F}}_{1}$ ($10^{-9}\mathrm{m}^{2}$)
---|---|---|---|---|---|---
1 | (0, 12.89, 0) | 0.141 | 2.61 | 1.34 | 2.29 | 0.45
2 | (3.35, 14.56, -5.02) | 0.196 | 1.55 | 27.6 | 69.8 | 53.1
3 | (-10.04, 16.24, 3.35) | 0.204 | 1.60 | 13.4 | 28.8 | 20.7
4 | (6.70, 17.91, -10.04) | 0.132 | 2.71 | 0.27 | - | -
5 | (-1.67, 19.58, 3.35) | 0.131 | - | 0 | - | -
6 | (-8.34, 20.42, 8.34) | 0.164 | 1.76 | 4.59 | 7.40 | 2.17
Table 4: Summary of the properties of the coffee samples (type M powder). The parameter $k^{\mathrm{D}}$ is the permeability obtained using Darcy’s law in the low $\nabla P$ regime, $k^{\mathrm{F}}$ and $k_{1}^{\mathrm{F}}$ are the permeability obtained using Forchheimer’s formula in the high $\nabla P$ regime. Sample | Location (mm) | Porosity | Tortuosity |  $k^{\mathrm{D}}$ ($10^{-13}\mathrm{m}^{2}$) |  $k^{\mathrm{F}}$ ($10^{-13}\mathrm{m}^{2}$) |  $k^{\mathrm{F}}_{1}$ ($10^{-9}\mathrm{m}^{2}$)
---|---|---|---|---|---|---
1 | (5.02, 21.18, -6.70) | 0.193 | 1.66 | 38.5 | 93.5 | 79.4
2 | (0, 19.51, 0) | 0.149 | 2.36 | 1.74 | 3.17 | 0.46
3 | (10.06, 17.83, -10.06) | 0.236 | 1.40 | 29.4 | 121.0 | 45.1
4 | (1.68, 16.99, 3.35) | 0.169 | 1.64 | 10.3 | 17.6 | 8.91
5 | (-5.03, 15.32, -3.35) | 0.158 | 1.94 | 7.81 | 13.9 | 8.19
6 | (10.06, 13.64, 3.35) | 0.202 | 1.59 | 13.0 | 23.4 | 16.1
Table 5: Summary of the properties of the coffee samples (type F powder). The parameter $k^{\mathrm{D}}$ is the permeability obtained using Darcy’s law in the low $\nabla P$ regime, $k^{\mathrm{F}}$ and $k_{1}^{\mathrm{F}}$ are the permeability obtained using Forchheimer’s formula in the high $\nabla P$ regime. Sample | Location (mm) | Porosity | Tortuosity |  $k^{\mathrm{D}}$ ($10^{-13}\mathrm{m}^{2}$) |  $k^{\mathrm{F}}$ ($10^{-13}\mathrm{m}^{2}$) |  $k^{\mathrm{F}}_{1}$ ($10^{-9}\mathrm{m}^{2}$)
---|---|---|---|---|---|---
1 | (1.68, 20.35, 10.06) | 0.193 | 1.41 | 66.6 | 213.6 | 161.4
2 | (0, 20.35, 0) | 0.183 | 1.84 | 27.4 | 94.8 | 41.0
3 | (10.06, 18.67, -10.06) | 0.185 | 2.60 | 1.38 | 1.92 | 0.527
4 | (-10.06, 15.32, 10.06) | 0.211 | 1.62 | 20.7 | 66.2 | 24.1
5 | (0, 14.48, 3.35) | 0.182 | 1.56 | 32.1 | 60.3 | 75.1
6 | (5.03, 12.80, -6.70) | 0.171 | 2.41 | 1.19 | 1.98 | 0.43
Figure 11: Permeability of different powder.

From the results presented above we can see that the variance of permeability
across the bed is quite large. Despite this, the trend of the average values
of the permeability (see Fig. 11) is in agreement with our qualitative
prediction based on the coarseness of the coffee particle: the more fines
there are, and the larger the specific surface area, the smaller the
permeability. Even though only a small number of samples is used for analysis,
different powders are still able to show qualitative differences. Given enough
computation resource quantitative comparison of extraction dynamics for
different powders will be possible.

###  3.4 Analysis of Post-extraction images

In order to study the link between changing microstructure and extraction, two
additional illy iperespresso capsules are further extracted, and the post-
extraction microCT imaging are analysed. These capsules are extracted using an
X7.1 Iperespresso - Capsules Coffee Machine at water temperature $97\pm
3^{\circ}$C and pressure $13\pm 1$ bar. The capsules are scanned with microCT
before and after the extraction so that a comparison can be made to
investigate the morphological change caused by the extraction. Because we do
not know the porosity of the extracted coffee cakes a priori, the segmentation
threshold cannot be determined by matching the porosity. However, since after
wetting the distinction between coffee grains and intergranular pores are much
more significant, we can manually select the grayscale threshold to capture
the phase boundary. An example of microCT imaging of the extracted dark roast
capsule is shown in Fig. 12 (a). As can be seen, the boundary between the
grains and the void is much clearer here than that in Fig. 5 (a). Fig. 12 (b)
also shows that the distribution of the grayscale in this microCT image is
bimodal. Therefore we can select the threshold manually to binarise the
microCT image: we choose the threshold to be the horizontal center of the two
distribution peaks. A threshold selected in this way makes the segmentation
least sensitive to the threshold value. The segmentation result of Fig. 12 (a)
is shown in Fig. 12 (c).

Figure 12: (a) Original microCT image of type H dark roast powder after
extraction. (b) Histogram of grayscale of the pixels. The red line marks the
manually-selected grayscale threshold. (c) Thresholded image.

A capsule with type H medium roast powder and a capsule with type H dark roast
powder are scanned before and after extraction. The comparison of the
binarized microCT images before and after extraction are shown in Fig. 13 and
15. The comparison of the porosity before and after extraction are presented
in Fig. 14 and 16. It can be seen that for all the capsules, after extraction
the coffee cakes become less porous, the porosity decreases significantly.
That is because after wetting the coffee grains swell[31, 3] and stick to
their neighbours, liquid fills the intragranular pores and also lots of
intergranular pores. The microCT cannot properly distinguish the liquid phase
from the solid phase, resulting in an underestimation of the porosity.
However, we focus on the relative spatial porosity variation in the coffee
cakes rather than the absolute porosity values. In Fig. 14 the porosity of the
type H medium roast powder before extraction and 3 days after extraction are
shown. Before extraction the porosity is relatively homogenous. On the
opposite, 3 days after the extraction, the porosity at the top center becomes
much larger than in other regions of the matrix. This is probably due to the
drying of the matrix, i.e. the drying is faster at the top center. However,
for the regions away from the center the porosity is always larger on the
bottom (outlet) than at the top (inlet). This phenomenon is counterintuitive
since we would expect the gravity will drive the residual water and other
fragments to the bottom making the lower part of the matrix less porous.
During extraction the applied external pressure should also have consolidated
the matrix more efficiently on the bottom. But surprisingly, after extraction
we observe larger porosity in the lower part (outlet) than in the upper part
(inlet). To exclude the effect of drying, in Fig. 16 the porosity of the type
H dark roast powder before extraction and right after extraction are shown.
Without the further complexity induced by the drying process, we see clearly
that the porosity is always higher on the bottom after the extraction. More
experiments have been performed with the post-extraction matrix being scanned
1 hr or 24 hrs after the extraction, we confirm that the porosity profiles
always have an apparent decreasing trend from the bottom to the top of the
matrix.

Besides the drying process and consolidation, other physical processes that
could also affect the geometry of the matrix are fragment migration[2, 14] and
grain swelling[3, 5]. In our prior works, we have modelled both the fine
migration[14] and particle swelling[16] using SPH. But the fine migration will
tend to lead to lower porosity near the bottom as the filter on the bottom
prevents the fines from leaving the capsule. On the other hand, if swelling
process is initiated earlier at the top it will cause the matrix to be less
porous at the top. This is indeed the case, since the water is infused from
the top to the bottom. However, the difference of porosity caused by this
effect is expected to be minor after extraction, as the infusion process is
fast ($\approx 1\mathrm{s}$), and the swelling process of all grains will
eventually saturate in about 4min[3], while we still see the porosity
difference 1hr after the extraction.

The most probable cause for the observed spatial porosity distribution is the
heterogenous erosion process. Scanning electron microscopy observation has
shown that immersed coffee grain has gas entrapped[31]. Coffee particles are
not fully wetted during extraction due to the abundant gas entrapped in the
intra-grain pores. Therefore, some part of the coffee particles remain dry and
brittle in the extraction process. And for brittle material, the failure
criterion is usually pressure-dependent. Note that the pressure is also quite
non-uniform in the coffee matrix as a result of the filtration process. The
heterogeneity of the post-extraction porosity is probably correlated to the
pressure drop across the matrix through the pressure-dependent failure
criterion of eroding material. As a matter of fact, in our previous numerical
simulation study[15], we propose a bottom-up mesoscopic erosion model in the
framework of SPH that incorporates the Mohr-Coulomb yield criterion[32, 33]
and the simple shear-erosion model[34, 4]. Simulations using this model show
that heterogeneity in filtration direction does occur as a result of the
pressure-dependent erosion: an initially relatively uniform porosity
distribution becomes highly non-uniform after filtration with higher porosity
on the bottom (the flow outlet) and lower porosity at the top (the flow
inlet). The non-uniformity is mainly controlled by the internal friction angle
of the eroding material. The post-extraction porosity profiles reported here
are in qualitative agreement with the prediction of our erosion model[15].
This suggests that our proposed model[15] is suitable to capture the special
characteristics of coffee cake erosion.

Figure 13: (a) Thresholded microCT image of type H medium roast powder, (a)
before extraction, (b) 3 days after extraction. Figure 14: Porosity of the
type H medium roast powder, (a) before extraction, (b) 3 days after
extraction. Figure 15: (a) Thresholded microCT image of type H dark roast
powder, (a) before extraction, (b) right after extraction. Figure 16: Porosity
of the type H dark roast powder, (a) before extraction, (b) right after
extraction.

##  4 Summary and conclusions

To summarize, we use high-resolution X-ray microCT technique to capture the
microscopic details of coffee matrices at particle-level. The 3D
reconstructured data are then used to build an SPH-based digital twin for the
percolation process.

We explore the link between microstructure and flow properties from two
perspectives: 1) The permeability and tortuosity. Using 3D microCT
reconstruction data combined with SPH percolation simulations, both the
permeability and tortuorsity of the matrices can be determined. We found that
the inertial effects are actually quite significant at normal pressure
gradient of espresso brewing. There have been some reports in literatures
showing that larger pressure gradient seems to lead to lower coffee matrix
permeability. We demonstrate here that this can be explained by considering
the inertial effect. We also examine several types of coffee powder, the
permeability and tortuosity are measured systematically using SPH simulation.
Our results show that the variance of the permeability at different location
of the same matrix is quite large, and the coarser the powder the larger the
variance. Despite this, the average of the permeability steadily increases as
the powder becomes coarser. 2) The post-extraction microstructure. We analyzed
the microCT images of post-extraction coffee matrices, and found that a
decreasing porosity profile (from the bottom outlet to the top inlet) always
develops after extraction. This counterintuitive phenomenon can be explained
using a pressure-dependent erosion model proposed in our prior paper[15]. It
is suggested that the erosion of coffee powder is to some extent pressure-
dependent. This pressure-dependence probably originates from the gas entrapped
inside the coffee particles: grains not fully wetted will remain sufficiently
brittle, and the failure of brittle material is usually pressure dependent.
Another physical scenario consistent with this pressure-dependent erosion
model is connected to the presence of individual fines adhering
electrostatically to the coarse particles, and being removed by the flow when
their superficial attraction force is overcome by hydrodynamic normal-
tangential forces.

To summarize, using the microCT and SPH simulations we have obtained two key
findings including the identification of the Darcy-Forchheimer transition and
the heterogeneous porosity profile in the filtration direction after
extraction. The Darcy-Forchheimer transition suggests that the decrease of
permeability under higher extraction pressure found in many experiments can be
attributed to inertial effects. It helps to clarify the discrepancy in
experiments from different research groups and provides guidance for new
measurement experiments. The post-extraction heterogeneity in the filtration
direction indicates that the erosion of coffee particles can be pressure-
dependent. It highlights the role played by the erosion mechanism of coffee
grains under hydrodynamic forces.

We mainly explored the link between the microstructure obtained from the
microCT and the flow properties in a qualitative way. Constrained by the
computational resource, a 3D model of a full realistic size can not be
simulated to establish quantitative links. Nevertheless, our results show that
microCT scan can provide us many microscopic details of a coffee matrix. These
information is critical if we want to accurately model the percolation
process. Moreover, by comparing the pre- and post-extraction scans we can
obtain some crucial information about how the geometry changes during
extraction process. This can in turn reveal some details of how the
interconnected complex physical processes develop during coffee extraction,
and helps us to propose better modelling methods to account for these
processes.

## Acknowledgement

This research is supported by the illycaffè S.p.A. through the project
“Modelling and simulating espresso coffee extraction at the mesoscopic
scales”. Financial support from the BERC 2022-2025 program and by the Spanish
State Research Agency through BCAM Severo Ochoa Excellence Accreditation
CEX2021-001142-S/MICIN/AEI/10.13039/501100011033 and through the project
PID2020-117080RB-C55 (“Microscopic foundations of soft matter experiments:
computational nano-hydrodynamics” - acronym “Compu-Nano-Hydro”) are also
acknowledged.

The microCT work was supported by the Advanced Imaging Materials (AIM) core
facility (EPSCR Grant No. EP/M028267/1), and the European Social Fund (ESF)
through the European Union Convergence programme administered by the Welsh
Government (80708).

## Competing interests statement

The co-authors affiliated to illycaffè S.p.A. declare that they have no known
competing financial interests or personal relationships that could have
appeared to influence the work reported in this paper. All other authors also
declare no competing interests.

## Data availability

The datasets used and/or analysed during the current study available from the
corresponding author on reasonable request.

## Bibliography

  * [1] Megan J Povey, Melvin J Holmes, Sameera Rafiq, Elena Simone, Michael Rappolt, and Mat Francis, editors.  Physics in Food Manufacturing.  2053-2563. IOP Publishing, 2020.
  * [2] M Petracco and F Suggi Liverani.  Espresso coffee brewing dynamics: development of mathematical and computational models.  In Colloque Scientifique International sur le Cafe, volume 15, pages 702–702. ASIC ASSOCIATION SCIENTIFIQUE INTERNATIONALE, 1993.
  * [3] Verena Bernadette Hargarten, Michael Kuhn, and Heiko Briesen.  Swelling properties of roasted coffee particles.  Journal of the Science of Food and Agriculture, 100(10):3960–3970, 2020.
  * [4] André F.V. Matias, Rodrigo C.V. Coelho, José S. Andrade, and Nuno A.M. Araújo.  Flow through time–evolving porous media: Swelling and erosion.  Journal of Computational Science, 53(Jan 2021), 2021.
  * [5] Matthew J. Maille, Kyle Sala, David M. Scott, and Hannah Zukswert.  Critical examination of particle swelling during wetting of ground coffee.  Journal of Food Engineering, 295:110420, 2021.
  * [6] A Voilley and D Simatos.  Modeling the solubilization process during coffee brewing.  Journal of Food Process Engineering, 3(4):185–198, 1980.
  * [7] Andrea Illy and Rinantonio Viani.  Espresso coffee: the science of quality.  Academic Press, 2005.
  * [8] Ernesto Illy and Luciano Navarini.  Neglected food bubbles: the espresso coffee foam.  Food biophysics, 6(3):335–348, 2011.
  * [9] R. Cappuccio, G. Cattaneo, G. Erbacci, and U. Jocher.  A parallel implementation of a cellular automata based model for coffee percolation.  Parallel Computing, 27(5):685–717, 2001.
  * [10] R. H. Rosa, G. V. Von Atzingen, V. Belandria, A. L. Oliveira, S. Bostyn, and J. A. Rabi.  Lattice Boltzmann simulation of cafestol and kahweol extraction from green coffee beans in high-pressure system.  Journal of Food Engineering, 176:88–96, 2016.
  * [11] K. M. Moroney, W. T. Lee, S. B.G. O’Brien, F. Suijver, and J. Marra.  Modelling of coffee extraction during brewing using multiscale methods: An experimentally validated model.  Chemical Engineering Science, 137(December):216–234, 2015.
  * [12] John Melrose, Borja Roman-Corrochano, Marcela Montoya-Guerra, and Serafim Bakalis.  Toward a New Brewing Control Chart for the 21st Century.  Journal of Agricultural and Food Chemistry, 66(21):5301–5309, 2018\.
  * [13] Michael I. Cameron, Dechen Morisco, Daniel Hofstetter, Erol Uman, Justin Wilkinson, Zachary C. Kennedy, Sean A. Fontenot, William T. Lee, Christopher H. Hendon, and Jamie M. Foster.  Systematically improving espresso: Insights from mathematical modeling and experiment.  Matter, 2(3):631–648, 2020.
  * [14] M. Ellero and L. Navarini.  Mesoscopic modelling and simulation of espresso coffee extraction.  Journal of Food Engineering, 263(May):181–194, 2019.
  * [15] Chaojie Mo, Richard Johnston, Luciano Navarini, and Marco Ellero.  Modeling the effect of flow-induced mechanical erosion during coffee filtration.  Physics of Fluids, 33(9), 2021.
  * [16] Chaojie Mo, Luciano Navarini, Furio Suggi Liverani, and Marco Ellero.  Modeling swelling effects during coffee extraction with smoothed particle hydrodynamics.  Physics of Fluids, 34(4):043104, 2022.
  * [17] Michael Kuhn, Sandra Lang, Franziska Bezold, Mirjana Minceva, and Heiko Briesen.  Time-resolved extraction of caffeine and trigonelline from finely-ground espresso coffee with varying particle sizes and tamping pressures.  Journal of Food Engineering, 206:37–47, 2017.
  * [18] C. D’Onofrio, R. van Loon, S. Rolland, R. Johnston, L. North, S. Brown, R. Phillips, and J. Sienz.  Three-dimensional computational model of a blood oxygenator reconstructed from micro-CT scans.  Medical Engineering and Physics, 47:190–197, 2017.
  * [19] H. C. Burridge, G. Wu, T. Reynolds, D. U. Shah, R. Johnston, O. A. Scherman, M. H. Ramage, and P. F. Linden.  The transport of liquids in softwood: timber as a model porous medium.  Scientific Reports, 9(1):1–14, 2019.
  * [20] M.-O. Baradez, C.P. McGuckin, N. Forraz, R. Pettengell, and A. Hoppe.  Robust and automated unimodal histogram thresholding and potential applications.  Pattern Recognition, 37(6):1131 – 1148, 2004.
  * [21] Nicolas Coudray, Jean Luc Buessler, and Jean Philippe Urban.  Robust threshold estimation for images with unimodal histograms.  Pattern Recognition Letters, 31(9):1010–1019, 2010.
  * [22] Itseez.  Open source computer vision library, 2015.
  * [23] Borja Roman-Corrochano.  Advancing the Engineering Understanding of Coffee Extraction.  PhD thesis, The University of Birmingham, 2015.
  * [24] Concetto Gianino.  Experimental analysis of the Italian coffee pot “moka”.  American Journal of Physics, 75(1):43–47, January 2007.
  * [25] Warren D King.  The physics of a stove-top espresso machine.  American Journal of Physics, 76(6):558–565, 2008.
  * [26] L. Navarini, E. Nobile, F. Pinto, A. Scheri, and F. Suggi-Liverani.  Experimental investigation of steam pressure coffee extraction in a stove-top coffee maker.  Applied Thermal Engineering, 29(5-6):998–1004, 2009.
  * [27] B. R. Corrochano, J. R. Melrose, A. C. Bentley, P. J. Fryer, and S. Bakalis.  A new methodology to estimate the steady-state permeability of roast and ground coffee in packed beds.  Journal of Food Engineering, 150:106–116, 2015.
  * [28] A. Bejan.  Convection Heat Transfer.  Wiley, 2013.
  * [29] Ricardo Dias, José A. Teixeira, Manuel Mota, and Alexander Yelshin.  Tortuosity variation in a low density binary particulate bed.  Separation and Purification Technology, 51(2):180–184, 2006.
  * [30] Zhongxia Li, Junwei Wan, Hongbin Zhan, Xi Cheng, Wei Chang, and Kun Huang.  Particle size distribution on Forchheimer flow and transition of flow regimes in porous media.  Journal of Hydrology, 574(April):1–11, 2019.
  * [31] Maria L. Mateus, Martine Rouvet, Jean C. Gumy, and Rémy Liardon.  Interactions of water with roasted and ground coffee in the wetting process investigated by a combination of physical determinations.  Journal of Agricultural and Food Chemistry, 55(8):2979–2984, 2007\.
  * [32] Richard M. Christensen.  The Theory of Materials Failure.  Oxford University Press, 2013.
  * [33] Joseph F Labuz and Arno Zang.  Mohr–coulomb failure criterion.  In The ISRM Suggested Methods for Rock Characterization, Testing and Monitoring: 2007-2014, pages 227–231. Springer, 2012.
  * [34] Shang Huan Chiu, M. N.J. Moore, and Bryan Quaife.  Viscous transport in eroding porous media.  Journal of Fluid Mechanics, 893, 2020.