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

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.


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 ∼ 40 µm, representing the socalled "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 ∼ 1000 µm.Moreover, the coffee particles are not entirely solid but porous material themselves, characterized by many cell-pockets of dimensions 30 ∼ 60 µ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 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 convectiondiffusion 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 [4,[14][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 matrixcan 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.

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±0.1 gram of coffee powder.The capsules used in our experiment are sent by illycaffè S.p.A.

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

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µm × 16.8µm × 16.8µ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.The microCT images are inverted.In the original image brighter pixels represent solid.We first invert the values of the pixels.
2. Assuming the greyscale threshold is known as T g , the microCT imaging are binarized using this threshold.
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. To avoid dealing with the cylindrical surface of the coffee cake, we extract a cuboid of 600 × 600 × 600voxels (1 × 1 × 1cm) 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µ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 [13,23], 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.

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µ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.
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

Coffee powder granulometry
The particle distributions for the coffee powders measured with laser diffraction are shown in Fig. mean diameter (d [4,3]), and the volume fraction of particles below 100µ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 100µ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.

Xray microCT analysis
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.
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.

Permeability and tortuosity
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µm is first examined.The SPH particle model has been shown in Fig. 2. Fluid buffer regions with a total height L buf = 283µ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 ρg(1 + L 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: gives a permeability k = 5.51 × 10 −13 m 2 .This value is smaller than the measurements by Gianino [24] (k = 2.3 × 10 −12 m 2 ) and King [25](k = 8.98 × 10 −13 m 2 ).But it is larger than the estimation given by Navarini et al [26] (k = 0.7 ∼ 4×10 −13 m 2 ) and Corrochano et al [27] (k = 2.59×10 −14 ∼ 3.36×10 −13 m 2 ).The disagreement is probably caused by the different overpressure (∆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 bar, while in the experiments of Navarini et al [26] and Corrochano et al [27] the overpressure is in the order of 1 bar.Here in our SPH model (Fig. 2), the pressure gradient 8 ∼ 40 bar/m represents an overpressure between 0.16 bar and 0.8 bar on a coffee bed of 2 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 ∼ 8 bar when the pressure gradient increase from 50 bar/m to 400 bar/m if a coffee bed of 2 cm height is considered.With the increase of pressure head the permeability decreases from about 5.2 × 10 −13 to 2.5 × 10 −13 m 2 .This kind of pressure dependence was also reported in Navarini's experimental investigation [26].
(a) 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.
The fit is very good and gives a permeability (k = 8.0 × 10 −13 m 2 ) and an additional inertial permeability k 1 = 2.17 × 10 −13 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 µm as the characteristic length, the Reynolds number (Re = d 4,3 ρq/µ) in Fig. 8 is between 0.84 to 3.86, indicating that the inertia is indeed important at high pressure gradient.
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:    where τ is the tortuosity, N is the total fluid particles, Ω c is the coffee sample domain, x i is the trajectory of particle i, β is the direction of the filtration.

SPH simulation Darcy
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 τ = 1.80.The power law relationship has been frequently used to describe the tortuosity-porosity correlation [29]: 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 τ / log ϵ = 0.33.9: Trajectories of fluid particles through the coffee sample.The arrow marks the direction of the pressure gradient.The black frame is the simulation box with dimension L x = L z = 708µm, L y = 991µ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 we extract more samples from the bed perform simulations on them.Results and 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 ∇P regime and relatively higher ∇P regime, respectively.All the samples show significant deviation from Darcy's law in the relatively high ∇P regime.But the Forchheimer's formula accurately captures the q−∇P relationship in the high ∇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 ∇P regime are presented in Tab. 2 as k D , the permeabilities obtained using Forchheimer's modification in the relatively high ∇P regime are presented in Tab. 2 as k F and k 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 D .Higher porosity tends to lead to larger k D .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 Table 3: Summary of the properties of the coffee samples (type E powder).The parameter k D is the permeability obtained using Darcy's law in the low ∇P regime, k F and k F 1 are the permeability obtained using Forchheimer's formula in the high ∇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 D (10 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.

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 ± 3 • C and pressure 13 ± 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 Table 5: Summary of the properties of the coffee samples (type F powder).The parameter k D is the permeability obtained using Darcy's law in the low ∇P regime, k F and k F 1 are the permeability obtained using Forchheimer's formula in the high ∇P regime.

Sample
Location (mm) Porosity Tortuosity k D (10  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   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 [3,31] 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 (≈ 1s), 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 [4,34].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.

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 (a) 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. Bibliography

Figure 1 :
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.

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

Figure 5 :
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.

Figure 6 :
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.

Figure 7 :
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)

Figure 8 :
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)

Figure 10 :
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)

Figure 12 :
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.

Figure 16 :
Figure 16: Porosity of the type H dark roast powder, (a) before extraction, (b) right after extraction.

Table 2 :
Summary of the properties of the coffee samples (type H powder).The parameter k D is the permeability obtained using Darcy's law in the low ∇P regime, k F and k F 1 are the permeability obtained using Forchheimer's formula in the high ∇P regime.D (10 −13 m 2 ) k F (10 −13 m 2 ) k F 1 (10 −9 m 2 )

Table 4 :
Summary of the properties of the coffee samples (type M powder).The parameter k D is the permeability obtained using Darcy's law in the low ∇P regime, k F and k F 1 are the permeability obtained using Forchheimer's formula in the high ∇P regime.