2.1 Grain geometry

It is widely reported in the literature [Reference Illy and Viani21, Reference Moroney, Lee, O’Brien, Suijver and Marra25, Reference Corrochano35, Reference Uman, Colonna-Dashwood, Colonna-Dashwood, Perger, Klatt, Leighton, Miller, Butler, Melot, Speirs and Hendoon43] that ground coffee has a bi-modal particle size distribution. As such, grains are often divided into two groups, termed boulders and fines. Boulders (larger than 100 $\mu$ m in diameter) are shown to be porous, and their internal porosity, $\varphi _{lb}$ , varies with bean origin, pretreatment and roasting processes, typically in the range $35\%-55\%$ [Reference Schenker, Handschin, Frey, Perren and Escher37]. In our model, the porous volume in any particular boulder is assumed to be instantaneously filled with liquid when liquid arrives outside the grain. Fines (smaller than 100 $\mu$ m) are assumed to have an internal structure that cannot be accessed by water. Ground coffee particles come in a wide variety of irregular shapes but for simplicity, we will approximate both types of grains as spherical with radii $a_i^*$ , for $i=f,b$ denoting fines and boulders respectively and the superscript ‘ $*$ ’ is used to indicate dimensional quantities. We define the radial coordinate $0 \leq r_i^* \leq a_i^*$ to be the distance from the centre of a grain of type $i$ .

2.2 Bed geometry

Espresso is brewed in a cylindrical container of height $L^*$ and radius $R_0^*$ . It is assumed that the properties of the bed are homogeneous in any circular cross-section, the liquid is uniformly injected at the inlet and the drag caused by the cylinder walls is negligible. Thus, macroscopically, i.e., on the scale of the whole bed (much larger than individual grains), the model for flow and transport in the intergranular pore space is one-dimensional. Hence, it is sufficient to form one-dimensional equations on the macroscale and define the coffee bed to lie in $0\leq z^* \leq L^*$ (see Figure 1). The inlet is located at $z^*=0$ and the outlet at $z^*=L^*$ .

The entirety of the coffee bed volume comprises of spherical grains and the space between them. The ratio of the total volume of grains to the bed volume is given by $\phi _s$ . The space between the grains (intergranular pore space) is initially occupied by air and is subsequently filled with liquid and has a volume fraction of the entire bed $\phi _l$ . The grain volume consists of the two size classes of grains, namely fines and boulders and their volume fractions (measured with respect to the whole bed) are $\phi _f$ and $\phi _b$ , respectively. Since the bed is homogeneous, these volume fractions are constant in any circular cross-section of the bed. Thus,

(1) \begin{align} \phi _{l} + \phi _s &= 1, & \phi _{s} & = \phi _{f} + \phi _{b}. \end{align}

The boulders are themselves porous, and their internal pores are accessible by liquid. Before a boulder comes in contact with liquid, its pore space is filled with air which (we assume) is instantaneously replaced by liquid on wetting. The fraction of the boulders’ volume (measured with respect to the boulders’ volume) accessible by liquid is taken to be $\varphi _{lb}$ , introduced above. We note that grains may swell when in contact with the liquid and increase their volume; however, there is disagreement in the literature over whether swelling is negligible [Reference Maille, Sala, Scott and Zukswert23] or significant ( $\sim 15 \%$ ) [Reference Hargarten, Kuhn and Briesen20]. Thus, due to the lack of consensus and in the interest of simplicity, we will not consider this phenomenon. Therefore, the total porosity $\phi _T$ of the bed is the sum of the intergranular pore space and the pore space within the boulders, that is

(2) \begin{equation} \phi _{T} = \phi _{l} + \phi _{b}\varphi _{lb}. \end{equation}

We model the transfer of solubles from grains to liquid as taking place on the spherical outer surface of the grains, i.e., any dissolution of solubles that occurs within the pore space of infiltrated boulders is not included as grain-to-liquid transfer. Equations for the rate of transport within grains (in the case of boulders this includes both the internal solid and liquid phases) and across the grain/liquid interface (the spherical outer surface of the grains) shall be formulated shortly. We therefore introduce the reactive spherical outer surface area of the grains (on which solubles can be transferred from grains into liquid) as the Brunauer-Emmett-Teller surface area, $b_i^*$ , defined to be the reactive surface area of grains of species $i$ per unit volume of bed. These are given by

(3) \begin{align} b_i^* = 4 n_i \pi {a_i^*}^2 \qquad \text{and} \qquad \phi _i = n_i \frac {4}{3} \pi {a_i^*}^3, \end{align}

where $n_i$ is the number of grains of type $i$ per unit volume of the bed. Hence, on eliminating $n_i$ between the two expressions in (3) we arrive at

(4) \begin{align} b_i^* = \frac {3 \phi _i}{a_i^*}. \end{align}

Figure 1.

A cartoon of the bed geometry. The liquid is injected at the inlet, located at $z^* = 0$ and exits the bed through the outlet at $z^* = L^*$ . During infiltration, the wetting front is located at $z^* = s_w^*(t^*)$ . The local volume fractions of the liquid, fines and boulders are given by $\phi _l, \phi _f$ and $\phi _b$ respectively. The volume fraction of pore space within a boulder is given by $\varphi _{lb}$ .

2.3 Liquid flow

We model the flow in the intergranular pore space using Darcy’s law with the incompressibility condition [Reference Allaire1, Reference Polishevsky33]. While the density of the liquid varies only slightly with soluble coffee content, the viscosity may be significantly higher than that of pure water at higher concentrations [Reference Petracco31]. Nonetheless, in the interests of simplicity, we elect to follow Moroney et al. [Reference Moroney, Lee, O’Brien, Suijver and Marra25] and adopt a constant density and viscosity approximation. The choice of boundary condition at the inlet depends upon the espresso machine of interest and in particular how its pump operates; most commonly this can apply either a fixed flow rate or a fixed pressure. Both situations are readily modelled, but here we choose to focus primarily on the fixed flow case for convenience. Readers interested in the fixed pressure case are referred to the supplementary material where this problem is presented.

On denoting the volumetric flux/Darcy flux (measured in units of m/s) by $q^*$ , incompressibility amounts to the statement that $\partial q^* / \partial z^* =0$ . On supplementing this with a fixed flow rate boundary condition at the inlet, we immediately learn that throughout the wet region we have

(5) \begin{align} q^*|_{z^*=0} = q^*_{\text{app}}, \qquad \text{and thus,} \qquad q^* = q^*_{\text{app}}. \end{align}

Provided that the shot is brewed properly (no channelling or clogging), the interface between the wet and dry parts of the bed appears to be sharp, as can be seen in experiments when espresso is brewed in a transparent portafilter [17]. The physics of the wetting process involves a complex interplay of the applied pressure, surface tension, capillary forces, swelling of the coffee cellular matrix, release of trapped gases (mainly carbon dioxide [Reference Clarke, Coffee and Macrae7, Reference Smrke, Wellinger, Suzuki, Balsiger, Opitz and Yeretzian39]) and potentially the nature of the dissolved solids in the water phase [Reference Ellero and Navarini9]. Here, we opt to avoid having to treat these phenomena and instead will make the assumption that the infiltration of the water into the intergranular pore space is described by Darcy’s law, as above, and that infiltration into the intragranular pore space is instantaneous on contact with water. Once water has entered the intragranular pore space, we will assume that it can no longer flow. On denoting the location of the wetting interface by $s_w^*(t^*)$ we have

(6) \begin{align} \phi _T \frac {ds_w^*}{dt^*} = q^*|_{z^*=s_w^*(t^*)} \end{align}

which states that the volumetric flux must be equal to the volume of pore space swept through by the wetting interface per unit time. We supplement this ordinary differential equation (ODE) for $s_w^*$ with the initial condition

(7) \begin{align} s_w^*|_{t^*=0} = 0 \end{align}

and note that when the waterfront arrives at the outlet of the machine $z^* = L^*$ , there is no longer a need to track the position of the interface. Solving (6) with (5) and (7) yields

(8) \begin{align} s_w^*(t^*) = \begin{cases} \frac {q_{\text{app}}^*}{\phi _{T}}t^* \quad &\text{for} \quad t^*\lt t_w^* \\ L^* \quad &\text{for} \quad t^* \geq t_w^* \end{cases}, \end{align}

such that at a time $t^*=t_{w}^*= \phi _T L^* / q_{\text{app}}^*$ the whole bed is saturated with water.

2.4 Transport in the intergranular pore space

We will model the change in concentration of coffee solubles in the three distinct volumes, namely the intergranular pore space and the space occupied by fines and boulders. It has been assumed that the macroscopic variables are one-dimensional, depending only on depth in the bed ( $z^*$ ). In the subsequent subsection, we will further assume that the microscale transport ( i.e., that of solubles in the interior of the grains) is also one-dimensional and depends only on radial distance from the centre of a grain $r^*$ . The model that we are in process of forming can aptly be termed a pseudo 2D model, similar to the Newman model [Reference Richardson, Denuault and Please34] which is the standard electrochemical model used to describe the operation of lithium-ion batteries. Thus, the model variable denoting the concentration of solubles in the liquid, $c_l^*=c_l^*(z^*,t^*)$ , is defined in the wet part of the bed $(0\leq z^*\leq s_w^*)$ and measured as the mass of solubles per unit volume of wet intergranular pore space. The concentration within the boulders is denoted by $c_b^*=c_b^*(r_b^*,z^*,t^*)$ and that in fines by $c_f^*=c_f^*(r_f^*,z^*,t^*)$ ; both of which are defined within the respective grains $(0\leq r_i^*\leq a_i^*) \times (0\leq z^*\leq L^*)$ for $i=b,f$ and measured per unit volume of grain. In the case of boulders, and consistent with previous assumptions, this grain volume is taken to include both the internal solid and liquid phases. We note that although the concentrations in the fines and boulders depend on two spatial variables, $r^*$ and $z^*$ , it is apposite to say that the microscopic equations are one-dimensional because (as we will see shortly) only derivatives with respect to $r^*$ appear in the governing PDEs, see (13)-(14).

We now formulate the equations governing the transport in the intergranular pore space, and those governing the transport in the grains will be treated in Section 2.5. The coffee solubles migrate through the liquid via a combination of diffusion and advection by the flow therefore, at macroscopic length scales, we have

(9) \begin{align} \phi _l \frac {\partial c_l^*}{\partial t^*} + \frac {\partial {\mathcal F}_l^*}{\partial z^*} = b_{f}^* G_f^* + b_{b}^* G_b^*, \quad {\mathcal F}_l^* = - D_{\text{eff}}^* \frac {\partial c_l^*}{\partial z^*} + q^* c_l^*, \end{align}

where $c_l^*$ is the mass concentration of solubles in the liquid, ${\mathcal F}_l^*$ is the flux of these solubles, $D_{\text{eff}}^*$ is the effective macroscopic diffusivity and $G_i^*$ are the mass fluxes of solubles exiting the grains and entering the water, measured on the surfaces of the fines ( $i=f$ ) and boulders ( $i=b$ ). The effective macroscopic diffusivity $D_{\text{eff}}^*$ accounts for both diffusion and also mechanical (hydraulic) dispersion, which arises from solute mixing due to microscale variations in the flow, at macroscopic scales. Ideally, this effective transport property should be computed via a suitable homogenisation procedure, see e.g., [Reference Foster, Gully, Liu, Krachkovskiy, Wu, Schougaard, Jiang, Goward, Botton and Protas16]. Alternatively one could approximate its value using, e.g., Bruggeman’s relation [Reference Bruggeman5].

The physically correct condition to apply at the inlet is that of continuity of flux. However, we do not wish to expand the model to involve the flow of water upstream of the inlet. To this end, and with foresight that parameters estimates will soon tell us that advective effects are dominant over diffusive ones, see Section 3, we choose to say impose that the liquid entering the bed is solubles-free and therefore we impose

(10) \begin{align} c_l^* |_{z^*=0} =0. \end{align}

An alternative, that also avoids needing to include a model for the flow upstream of the inlet, would be to require that the flux of solubles at the inlet is zero, i.e., ${\mathcal F}_l^*|_{z^*=0}=0$ ; however in the relevant scenario, where advective effects dominate diffusive ones, this collapses to (10). At the wetting front, we enforce conservation of solubles accounting for the movement of the front and allowing solubles to move from the liquid in the intergranular pore space and into the pore space inside the boulders. To formulate this condition, we apply a ‘Gaussian pillbox’ [Reference Griffiths19] argument; details are given in the supplementary material. The model doesn’t capture the flow after it leaves the bed, therefore we require a boundary condition at the outlet $z^* = L^*$ . We choose to impose zero diffusive flux at the outlet; however, this boundary condition is passive as it has no effect on the asymptotic analysis that follows. The set of boundary conditions supplementing (9) is as follows:

(11) \begin{align} -\frac {d s_w^*}{dt^*} \phi _l c_l^*|_{z^*=s_w(t)} + {\mathcal F}_l^*|_{z^*=s_w^*(t^*)} = \frac {d s_w^*}{dt^*} \phi _b \varphi _{lb} c_l^* |_{z = s_w^*(t^*)} \, \, \text{for} \, \, t^*\lt t_w^*, \end{align}

(12) \begin{align} \frac {\partial c_l^*}{\partial z^*}|_{z^*=L^*} = 0 \, \, \text{for} \, \, t^*\geq t_w^*. \end{align}

2.5 Transport in the intragranular pore space (coffee grains)

We follow previous studies [Reference Cameron, Morisco, Hofstetter, Uman, Wilkinson, Kennedy, Fontenot, Lee, Hendon and Foster6, Reference Moroney, O’Connell, Meikle-Janney, O’Brien, Walker and Lee26] and model the transport of solubles in the grains by Fickian diffusion. Each grain contains a concentration of soluble coffee $c_f^*$ or $c_b^*$ for fines and boulders, respectively. Since boulders have an internal pore fraction $\varphi _{lb}$ , the concentration $c_b^*$ should be understood as the mass of solubles per unit volume within the boulder including both the solid and liquid phases, by contrast, that in fines, $c_f^*$ , is the mass of solubles per unit volume of solid (the fines are assumed to have zero internal porosity). When the grains become wet, the solubles migrate from the centre towards the surface of the grain and once they are at the surface, they are transported to the surrounding liquid at a rate $G_i^*$ , to be defined shortly. The multiscale nature of the model requires that, at each location in macroscopic space $z^*$ , there is a microscopic coordinate $r_i^*$ denoting the distance from the centre of the grains. The concentration of coffee solubles is assumed to be spherically symmetric and that solubles become immediately mobile after wetting. The latter assumes that dissolution of solid coffee within the grains takes place on a much faster timescale than diffusion of dissolved coffee out of the grains. Thus

(13) \begin{align} \frac {\partial c_{i}^*}{\partial t^*} + \frac {1}{r_i^{*2}} \frac {\partial }{\partial r_i^*} \left ( r_i^{*2} {\mathcal F}_i^* \right )=0, \quad {\mathcal F}_i^* = - D_{si}^* \frac {\partial c_i^*}{\partial r_i^*}, \quad \text{for} \quad i=b,f, \end{align}

with boundary conditions

(14) \begin{align} {\mathcal F}_i^* \Big {|}_{r_i^*=0} = 0, \qquad {\mathcal F}_i^* \Big {|}_{r_i^*=a_i^*} = G_i^*, \end{align}

where $D_{si}^*$ is the diffusivity of solubles within the grains for $i = b,f$ . In light of our previous assumptions, these material properties should strictly be thought of as effective diffusivities that characterise the rate of transport through the grains‘ internal multiphase structure. Roasted coffee beans can be considered to consist of porosity made up of larger spherical pores (cell lumina with radii of approximately $2\times 10^{-5}$ m, connected by nanoscale pores in the cell walls (the walls have thicknesses of approximately $10^{-5}$ m. The internal structure of the two grain types being different, e.g., in contrast to the fines the boulders have much more significant internal pore space, therefore likely means that $D_{sb}^*$ and $D_{sf}^*$ differ.Footnote ¹ Despite this, and in order to simplify the subsequent analysis, we will later assume they may be taken to be equal as a first approximation.

The initial concentration of solubles coffee per unit volume of grain, in fines and boulders (prior to wetting) is

(15) \begin{align} c_i^*|_{t^*=0} = c_{i,\text{init}}^*. \end{align}

In accordance with our previous assumptions, in the case of boulders, this is the be understood as the concentration per unit volume of both internal solid and pore space. Thus, the maximum mass of solubles that can be extracted from a given grain is $4\pi a_i^3/3$ for $i=b,f$ . When the wetting front moves through the bed, the liquid instantaneously fills the internal pores of the boulders. The infiltrating liquid may well be carrying solubles and therefore, immediately after the boulder is wet, the concentration of solubles in the boulder (measured with respect to both solid and now wet pore space) is

(16) \begin{align} c_b^*|_{z^*=s_w^*(t^*)}= c_{b,\text{init}}^* + c_l^*|_{z^*=s_w^*(t^*)}\varphi _{lb}. \end{align}

Thus, it is conceivable that shortly after wetting the concentration in the boulders is higher than the saturation concentration, i.e., $c_b^* \gt c_{\text{sat}}^*$ . Should this occur, in reality we would expect that solubles in the solid phase of the boulders would remain in solid immobile form until the local concentration within the pore space of the boulder drops sufficiently to allow internal dissolution to take place. We have adopted an effective/homogenised description of the internals of the boulders and assumed that dissolution of the solid phase solubles is always instantaneous; therefore, our model is incapable of capturing this phenomena systematically. Alleviation of this shortcoming is left open for future investigation but interested readers are directed to [Reference Moroney and Vynnycky30] which may provide a starting point for the remedy.

2.6 Mass transfer between grains and liquid in the wet region

Selecting a suitable extraction rate is non-trivial; however, we will assume that a sensible rate should have the following properties: (i) there is no transfer when the liquid immediately outside the grain is saturated, or when the concentrations of the liquid immediately outside is equal to that on its surface; (ii) mass transfer from regions with low concentration to regions with high concentrations is higher than regions with similar concentrations. With these in mind, we define an extraction rate similar to Cameron et al. [Reference Cameron, Morisco, Hofstetter, Uman, Wilkinson, Kennedy, Fontenot, Lee, Hendon and Foster6] and say that

(17) \begin{align} G_i^{*} = k^*\begin{cases} c^*_{\text{sat}}-c^*_{l} \quad \text{if} \quad c_{i}^*|_{r^*_i=a^*_i} \geq c^*_{\text{sat}} \\ c^*_{i}|_{r^*_i=a^*_i}-c^*_{l} \quad \text{if} \quad c^*_{i}|_{r^*_i=a^*_i} \lt c^*_{\text{sat}} \end{cases} \end{align}

where $k^*$ is the mass transfer rate constant that has units $ \text{m s}^{-1}$ and $c^*_{\text{sat}}$ is the saturation concentration of solubles in the liquid. We emphasise that more sophisticated choices for the rate of mass transfer are certainly possible. For example, if data were to become available that shed light on the equilibrium relationship between the concentration of coffee in the grains and that in the liquid one could include then include and select an appropriate value for a partition coefficient [Reference Bear and Cheng2, Reference Melrose, Roman-Corrochano, Montoya-Guerra and Bakalis24, Reference Corrochano35]. In fact, there are good reasons to expect that this may be warranted. First, the concentrations in the different phases (boulders, fines and liquid) are measured with respect to different reference volumes: in the grains, it is measured per unit volume of grain including both intragranular pores and the internal solid structure whereas in the liquid it is per unit volume of liquid only. Thus, these may differ from the true local concentrations within the grain’s internal pores, which is what we might expect to be in equilibrium with the concentration in the liquid. Second, chemical binding and interactions of different coffee molecules with the intragranular matrix may also result in alterations to the equilibrium concentrations thereby warranting the inclusion of a partition coefficient. This being said, owing to the lack of requisite data, we follow [Reference Melrose, Roman-Corrochano, Montoya-Guerra and Bakalis24] and assume a partition coefficient of unity for both grain species. Before moving on, we point out that the inclusion of a partition coefficient makes little difference to the overall structure of ensuing analysis and therefore does not hamper our central aim of elucidating the mathematical structure of the problem. More details on the partition coefficient and attempts to measure it can be found in [Reference Corrochano35].

2.7 Does saturation occur?

As we will demonstrate, it is possible to distinguish two varieties of espresso recipe, namely those that give rise to extractions in which the liquid that initially leaves the bed is so concentrated that it is fully saturated vs. those in which it is not. We have conducted experiments in which the first drips of liquid from a café-style espresso recipe are segregated from the remainder of the shot. The concentration of solubles in these initial drips appears to be fixed, independent of the details of the recipe, therefore suggesting that we have reached a maximum obtainable concentration which we expect to be the saturation limit of the water. More details on these experiments will be given in a forthcoming article. Whether saturation occurs is primarily determined by a competition between the rate of extraction from the fines (significantly faster than that from boulders owing to the size disparity between the two grain types) and the flow rate. When extraction from fines is fast compared to the flow rate, a saturated layer will form immediately above the wetting front because solubles are extracted from fines faster than they are transported away from the wetting front. This being said, fast transport within fines alone is not sufficient to form this saturated layer. In addition, fines must contain enough solubles (per unit volume of bed, $c_{f,\text{init}} \phi _f$ ) to saturate the surrounding liquid. Conversely, when extraction from fines is relatively slower and/or there are less solubles in the fines then saturation does not occur. Later in this work, in particular see (74) and the discussion shortly thereafter, we will be able to be more precise about demarcating these two regimes, however, for the time being, it is sufficient to note that for relevant parameter values, experimental evidence suggests that a saturated layer does form. As such, we shall, in the body of this paper present the analysis of regime in which a saturated layer forms; likewise is done for the case of no saturation in the supplementary material.

3.1 Parameter estimates

The asymptotic analysis that follows requires estimates of the sizes of the dimensionless parameters. Currently, as far as we are aware, there is no complete set of parameters for espresso brewing. In an attempt to formulate such a set, we combine the estimates presented in [Reference Cameron, Morisco, Hofstetter, Uman, Wilkinson, Kennedy, Fontenot, Lee, Hendon and Foster6, Reference Moroney, Lee, O’Brien, Suijver and Marra25, Reference Moroney, O’Connell, Meikle-Janney, O’Brien, Walker and Lee26, Reference Petracco32, Reference Uman, Colonna-Dashwood, Colonna-Dashwood, Perger, Klatt, Leighton, Miller, Butler, Melot, Speirs and Hendon42].

The parameters $Q_i$ , $c_{i,\text{init}}$ , $b_i$ depend simply on the particle size distribution and bed packing. They are calculated based on estimates of $a_i^*$ in [Reference Moroney, Lee, O’Brien, Suijver and Marra25, Reference Moroney, O’Connell, Meikle-Janney, O’Brien, Walker and Lee26] listed in Table 1, and would change slightly based on bed packing and the grind setting determining the size distributions of grains. The remaining parameters $D_{\text{eff}}$ , $D_{si}$ and $\epsilon$ depend on the diffusivities as well as the fluid velocity $q_{\text{app}}^*$ and the bed height $L^*$ and are comparable to the estimates given in [Reference Cameron, Morisco, Hofstetter, Uman, Wilkinson, Kennedy, Fontenot, Lee, Hendon and Foster6]. In the absence of detailed experimental data, we shall now assume that the diffusivities in the boulders and fines are equal despite there being reasons to believe that this may not be the case; see the discussion under (14). Before moving on to the asymptotic analysis we point out that despite the naivety of assuming that the diffusivities in the boulders and fines are equal, there is little to no impact on the reduced model that we ultimately derive. In what follows we will examine the asymptotic limit that $D_{sf} \to \infty$ . As an upshot, the exact value of the diffusivity in the fines does not appear at leading order, and therefore our reduced model is insensitive to its value. Thus, the ensuing analysis holds provided future evidence does not indicate that the diffusivities change markedly from the estimates that we use here.

Table 1.

Typical model parameter values. Those marked with superscript 1 are fitted to data, those with superscript 2 are available or can be estimated from the literature [Reference Cameron, Morisco, Hofstetter, Uman, Wilkinson, Kennedy, Fontenot, Lee, Hendon and Foster6, Reference Moroney, Lee, O’Brien, Suijver and Marra25, Reference Moroney, O’Connell, Meikle-Janney, O’Brien, Walker and Lee26, Reference Petracco32, Reference Uman, Colonna-Dashwood, Colonna-Dashwood, Perger, Klatt, Leighton, Miller, Butler, Melot, Speirs and Hendon42], and those marked with 3 are directly available from the experiment. Parameters marked with both superscripts 1 and 2 are fitted but their values are comparable to those available in the literature. We note that the diffusivities $D_{sf}^*,D_{sb}^*$ in the two grain species may differ due to their internal structure; however in the absence of data on those values, we use the same estimate

4.1 Outline of the asymptotic analysis

For $t=O(1)$ the asymptotic structure of the solution, in the limit that $\epsilon \rightarrow 0$ , separates into three regions as shown in Figure 2. In Section 4.2, we analyse region (iii), adjacent to the wetting front, and in which the quasi-equilibrium for a saturated liquid pertains. Thus, in region (iii) there is little extraction and the solution structure is particularly simple. In Section 4.3, we examine region (ii), a narrow region separating regions (i) and (iii), in which the fines undergo rapid extraction, thereby allowing the concentration in the liquid to change markedly. Region (i), adjacent to the inlet, is the subject of Section 4.4. Here, the concentrations at the surfaces of the grains (of both types) and in the liquid are in a quasi-equilibrium. The boulders and fines undergo slow extraction leading to a concentration profile in the liquid that slowly increases with depth into the bed. In Section 4.5, we carry out the asymptotic matching required to tie all three regions together thereby closing the problem. A noteworthy feature of the analysis that follows is that the position of region (ii) cannot be determined in advance; instead must be determined as part of the solution to the problem. Thus, region (ii) can aptly be termed an intrinsic moving internal boundary, and the mathematics required to resolve this is interesting in and of itself.

For short times of $t=O(\epsilon )$ , before the saturated layer has formed the asymptotic structure is made up of only a single narrow spatial region. This single region is of width $O(\epsilon )$ and is located adjacent to the inlet. Here, the fines rapidly dissolve in order to boost the concentration in the incoming water to saturation. We shall not tackle this short time problem because our primary interest is predicting the concentration at the exit of the bed. The error engendered in omitting to analyse and include this short time problem (which would alter the ‘initial conditions’ for the $t=O(1)$ problem that is the focus of the remainder of this section) is small in the sense that it makes only an $O(\epsilon )$ difference to the prediction of the outlet concentration.

Figure 2.

Schematics to illustrate the asymptotic structure of the problem. The left panel shows a snapshot of the different regions of the bed arising in the asymptotic analysis along with the moving boundaries $s_d(t)$ and $s_{w}(t)$ . The right panel shows the location of these moving boundaries evolving in the $z$ – $t$ plane. The wetting front, $s_w(t)$ , is represented by the solid black line. Region (iii) is located adjacent to the wetting front, region (i) is adjacent to the inlet. The saturation interface represented by the solid blue line and located at $z = s_d(t)$ , is surrounded by a narrow region (ii).

4.2 Region (iii): saturated liquid region

This region is located immediately behind the wetting front at $z=s_w(t)$ . The liquid is close to its saturation concentration and hence dissolution rates are almost zero. Grains in this region have not yet been in contact with unsaturated water and therefore will remain close to their initial concentration, with the boulder’s internal solubles concentration having been boosted by the infiltration of the saturated liquid. We seek an asymptotic solution in which all dependent variables are expanded in a regular asymptotic series in small $\epsilon$ as follows

(32) \begin{align} c_{l} = c_{l,0}^{\text{(iii)}}+ {o(1)}, \quad {c_{f}} = {c^{\text{(iii)}}_{f,0}}+ {o(1)}, \quad c_{b} = c_{b,0}^{\text{(iii)}} + {o(1)}, \end{align}

(33) \begin{align} {\mathcal F}_{l} = {\mathcal F}_{l,0}^{\text{(iii)}} + {o(1)}, \quad {\mathcal F}_{f} = {\mathcal F}_{f,0}^{\text{(iii)}} + {o(1)}, \quad {\mathcal F}_{b} = {\mathcal F}_{b,0}^{\text{(iii)}} + {o(1)}, \end{align}

(34) \begin{align} {G_{f}= G_{f,0}^{\text{(iii)}}+ {o(1)}, \quad G_{b}= G_{b,0}^{\text{(iii)}}+ {o(1)}}. \end{align}

In what follows we seek only the lowest-order terms, and thus, the size of the corrections is of secondary interest. Inserting these expansions in the definition of the reaction rates (27) and balancing leading order terms asserts that

(35) \begin{align} c_{l,0}^{\text{(iii)}} = 1, \end{align}

provided that $c_{f,0}^{\text{(iii)}}|_{r_f=1} \gt 1$ and $c_{b,0}^{\text{(iii)}}|_{r_b=1} \gt 1$ , which as we will show is typical in this region for our parameter set. The other possible balances $c_{f,0}^{\text{(iii)}}|_{r_f=1}=c_{l,0}^{\text{(iii)}}$ and $c_{b,0}^{\text{(iii)}}|_{r_b=1}=c_{l,0}^{\text{(iii)}}$ are rejected on the basis that the fines within region (iii) are yet to come into contact with liquid that is not saturated and therefore internal concentration remains at its initial value and boulders are topped up with saturated liquid. This result (35) also follows from the solution of (25) after inserting the expansions (32)-(34) and balancing leading order terms,

(36) \begin{align} c_{f,0}^{\text{(iii)}}= c_{f,\text{init}}, \quad {\mathcal F}_{f,0}^{\text{(iii)}} = 0. \end{align}

It follows from (28) that the boulders which have just been infiltrated with water, at $z=s_w(t)$ , have an internal concentration $c_{b,0}^{\text{(iii)}}|_{z=s_w(t)} = c_{b,\text{init}} + \varphi _{lb}\gt 1$ for a typical choice of parameters. Thus, since the boulders within region (iii) have also yet to come into contact with unsaturated water, and their internal concentration is larger than the saturation concentration of water, no mass transfer occurs at their surfaces and hence

(37) \begin{align} c_{b,0}^{\text{(iii)}} = c_{b,\text{init}} + \varphi _{lb}, \quad {\mathcal F}_{b,0}^{\text{(iii)}} = 0. \end{align}

It is now straightforward to determine the leading order flux of solubles in the liquid by back substituting (35) into (23) and balancing at leading order. This yields

(38) \begin{align} {\mathcal F}_{l,0}^{\text{(iii)}} = q. \end{align}

Balancing leading order terms in (27) directly yields the source terms in the region:

(39) \begin{align} G_{f,0}^{\text{(iii)}}=0 \quad \text{and} \quad G_{b,0}^{\text{(iii)}}=0. \end{align}

4.3 Region (ii): local to the saturation front

Region (ii) is a slender region adjacent to region (iii), in which the fines rapidly undergo dissolution thereby boosting the concentration in the intergranular liquid to saturation. In order to appreciably change $c_l$ across this narrow region, and match together regions (i) and (iii), we require a large mass transfer rate from one of the sets of particles. It is only the smaller fines that are able to supply the requisite size flux through their surfaces. We use the method of dominant balance to determine the local variable transform in region (ii), for which the dissolution terms are in balance with diffusion in the fines and with advection and diffusion in the liquid terms. The solution in this region is found by making the following change of variables

(40) \begin{align} \epsilon w =z-s_d, \end{align}

under which equations (23)–(25) become:

(41) \begin{align} \frac {\phi _l}{\phi _T}\left ( \frac {\partial c_l}{\partial t} - \frac {\dot {s_d}}{\epsilon }\frac {\partial c_l}{\partial w} \right ) + \frac {1}{\epsilon }\frac {\partial {\mathcal F}_l}{\partial w} = { G_f + G_b}, \quad \quad {\mathcal F}_{l} = -\hat {D}_{\text{eff}} \frac {\partial c_{l}}{\partial w} + q c_l, \end{align}

(42) \begin{align} \frac {\partial c_{f}}{\partial t} - \frac {\dot {s_d}}{\epsilon }\frac {\partial c_f}{\partial w} = \frac {1}{r_f^{2}} \frac {\partial }{\partial r_f} \left ( r_f^{2} \frac {\hat {D}_{sf}}{\epsilon } \frac {\partial c_f}{\partial r_f} \right ),\end{align}

(43) \begin{align} \mbox{with} \quad - \frac {\hat {D}_{sf}}{\epsilon } \frac {\partial c_f}{\partial r_f} \Big {|}_{r_f=0} = 0 \quad \mbox{and} \quad -\frac {\hat {D}_{sf}}{\epsilon } \frac {\partial c_f}{\partial r_f} \Big {|}_{r_f=1} = Q_fG_f.\end{align}

(44) \begin{align} \frac {\partial c_{b}}{\partial t} - \frac {\dot {s_d}}{\epsilon }\frac {\partial c_b}{\partial w} = \frac {1}{r_b^{2}} \frac {\partial }{\partial r_b} \left ( r_b^{2} D_{sb} \frac {\partial c_b}{\partial r_b} \right ),\end{align}

(45) \begin{align} \mbox{with} \quad - D_{sb} \frac {\partial c_b}{\partial r_b} \Big {|}_{r_b=0} = 0 \quad \mbox{and} \quad -D_{sb} \frac {\partial c_b}{\partial r_b} \Big {|}_{r_b=1} = Q_bG_b.\end{align}

These equations need to match the solutions in region (iii) as $w \to \infty$ and the solution in region (i) as $w \to -\infty$ . The matching conditions can be written as

(46) \begin{align} \left .\begin{array} c_l \to 1 + \cdots \\ c_{f} \to c_{f,\text{init}} + \cdots \\ c_b \to c_{b,\text{init}} + \varphi _{lb} \end{array}\right \} \mbox{as} \quad w \to +\infty, \end{align}

(47) \begin{align} \left .\begin{array} c_l \to c_{l,0}^{\text{(i)}}|_{z=s_d(t)} + \cdots \\ c_{f}|_{r_f=1} \to c_{l,0}^{\text{(i)}} |_{z=s_d(t)} + \cdots \\ c_b|_{r_b=1} \to c_{l,0}^{\text{(i)}} |_{z=s_d(t)} + \cdots \end{array}\right \} \mbox{as} \quad w \to -\infty . \end{align}

We seek an asymptotic solution using the following expansions

(48) \begin{align} c_{l} = c_{l,0}^{\text{(ii)}}+ {o(1)}, \quad {c_{f}} = {c_{f,0}}^{\text{(ii)}}+ {o(1)}, \quad c_{b} = c_{b,0}^{\text{(ii)}} + {o(1)}, \end{align}

(49) \begin{align} {\mathcal F}_l = {\mathcal F}_{l,0}^{\text{(ii)}} + {o(1)}, \quad {\mathcal F}_f = {\mathcal F}_{f,0}^{\text{(ii)}} + {o(1)}, \quad {\mathcal F}_b = {\mathcal F}_{b,0}^{\text{(ii)}} + {o(1)}, \end{align}

(50) \begin{align} G_f = \frac {1}{\epsilon } G_{f,-1}^{\text{(ii)}}+ {o(1)}, \quad G_{b}= o \left ( \frac {1}{\epsilon } \right ).\end{align}

In what follows we seek only the lowest-order terms. Substituting the expansions (48)–(50) and balancing leading order terms in (41) yields

(51) \begin{align} - \dot {s_d} \frac {\phi _l}{\phi _T} \frac {\partial c_{l,0}^{\text{(ii)}}}{\partial w} + \frac {\partial {\mathcal F}_{l,0}^{\text{(ii)}}}{\partial w} = { G_{f,-1}^{\text{(ii)}} }\quad \mbox{where} \quad {\mathcal F}_{l,0}^{\text{(ii)}} = -\hat {D}_{\text{eff}} \frac {\partial c_{l,0}^{\text{(ii)}}}{\partial w} + q c_{l,0}^{\text{(ii)}}. \end{align}

Integrating the former of the two equations in (51) across region (ii) and imposing the matching conditions (46) and (47) leads to

(52) \begin{align} \left ( - \dot {s_d} \frac {\phi _l}{\phi _T} + q \right ) \left ( 1 - c_{l,0}^{\text{(i)}}|_{z=s_d(t)} \right ) = \int _{-\infty }^{\infty } {G_{f,-1}^{\text{(ii)}} }dw. \end{align}

Substituting the expansions (48)–(50) and balancing leading order terms in (42)–(43) leads to

(53) \begin{align} -\dot {s_d}\frac {\partial c_{f,0}^{\text{(ii)}}}{\partial w} = \frac {1}{r_f^{2}} \frac {\partial }{\partial r_f} \left ( r_f^{2} \hat {D}_{sf} \frac {\partial c_{f,0}^{\text{(ii)}}}{\partial r_f} \right ), \end{align}

(54) \begin{align} \mbox{with} \quad - \hat {D}_{sf} \frac {\partial c_{f,0}^{\text{(ii)}}}{\partial r_f} \Big {|}_{r_f=0} = 0 \quad \mbox{and} \quad -\hat {D}_{sf} \frac {\partial c_{f,0}^{\text{(ii)}}}{\partial r_f} \Big {|}_{r_f=1} = Q_fG_{f,-1}^{\text{(ii)}}. \end{align}

Multiplying through by $r_f^2$ and integrating over the radial domain and imposing the boundary conditions (54) gives

(55) \begin{align} -\dot {s}_d \int _0^1 r_f^2 \frac {\partial c_{f,0}^{\text{(ii)}}}{\partial w} dr_f = -Q_fG_{f,-1}^{\text{(ii)}}. \end{align}

This can then be integrated across the width of region (ii) yielding

(56) \begin{equation} -\dot {s}_d \int _{-\infty }^{+\infty } \int _0^1 r_f^2 \frac {\partial c_{f,0}^{\text{(ii)}}}{\partial w} dr_f dw = \int _{-\infty }^{+\infty }-Q_fG_{f,-1}^{\text{(ii)}} dw. \end{equation}

We can change the order of integration on the left hand side because $r_f$ and $w$ , i.e., the micro- and macroscopic coordinates respectively, are orthogonal (a foundational assumption of multiple scales homogenisation on which our model is based). Evaluating the integrals that result leads us to

(57) \begin{align} -\frac {\dot {s}_d}{3} \left ( c_{f,\text{init}} - c_{l,0}^{\text{(i)}}|_{z=s_d(t)} \right ) = -Q_f \int _{-\infty }^{+\infty } G_{f,-1}^{\text{(ii)}} dw. \end{align}

Eliminating the integral in (57), representing the fine dissolution rate, using (52) and the matching conditions (46)–(47), gives

(58) \begin{align} {\frac {\dot {s}_d}{3Q_f} \left ( c_{f,\text{init}} - c_{l,0}^{\text{(i)}}|_{z=s_d(t)} \right ) = \left ( - \dot {s_d} \frac {\phi _l}{\phi _T} + q \right ) \left ( 1 - c_{l,0}^{\text{(i)}}|_{z=s_d(t)} \right ).} \end{align}

The result (58) is one of the key equations that will form the reduced model, see Section 4.5. Equation (58) is a concise statement of the conservation of solubles across region (ii) and serves as an ODE that (when coupled with the remainder of the reduced model) will allow us to determine the position $z=s_d(t)$ , separating regions (i) and (iii). For consistency with the initially dry coffee bed, we will require that $s_d|_{t=0}=0$ , see (74), and (58) pertains until such a time that $s_d(t)=1$ , i.e., when region (ii) reaches the bottom of the basket, after which (58) is no longer necessary because region (i) occupies the entirety of the basket and hence there is no longer a demarcation between regions (i) and (iii). Finally, we note that the leading terms in equation (44) that governs the concentration in the boulders throughout region (ii) asserts that

(59) \begin{align} \frac {\partial c_{b,0}^{\text{(ii)}}}{\partial w} = 0 \end{align}

whose solution, along with the matching condition (46) is

(60) \begin{align} c^{\text{(ii)}}_{b,0} = c_{b,\text{init}}+\varphi _{lb}. \end{align}

The remaining equations for the leading order concentrations in the fines and liquid could be solved if necessary but, as we shall show, the analysis presented here gives sufficient information to evolve the position of the saturation front.

4.4 Region (i): slowly reacting region

For $t=O(1)$ , once the fines near the inlet have been dissolved to such an extent that they are close to being in equilibrium with the surrounding water ( $c_f|_{r_f=1} \approx c_l$ ), region (i) is formed in which the dissolution continues controlled by the rate at which solubles can be transported from the interior of the boulders to their surface. We seek an asymptotic solution in region (i) of the following form

(61) \begin{align} c_{l} = c_{l,0}^{\text{(i)}}+ {o(1)}, \quad {c_{f}} = {c^{\text{(i)}}_{f,0}}+ {o(1)}, \quad c_{b} = c_{b,0}^{\text{(i)}} + {o(1)},\end{align}

(62) \begin{align} {\mathcal F}_l = {\mathcal F}_{l,0}^{\text{(i)}} + {o(1)}, \quad {\mathcal F}_f = {\mathcal F}_{f,0}^{\text{(i)}} + {o(1)}, \quad {\mathcal F}_b = {\mathcal F}_{b,0}^{\text{(i)}} + {o(1)}, \end{align}

(63) \begin{align} G_f = G_{f,0}^{\text{(i)}} + {o(1)}, \quad G_{b}= G_{b,0}^{\text{(i)}} + {o(1)}.\end{align}

As in the previous regions, in what follows we seek only the lowest-order terms. Substituting the expansions (61)–(63) and balancing leading order terms in (25) reveals that diffusion through the fines is sufficiently fast that at leading order their internal concentration is independent of $r_f$ . Hence, on multiplying (25) through by $r_f^2$ , integrating across $r_f \in (0,1)$ , and imposing the boundary conditions (26) at leading order, we discern that

(64) \begin{align} \frac {\partial c_{f,0}^{\text{(i)}}}{\partial t} = - 3 Q_f G_{f,0}^{\text{(i)}}. \end{align}

Inserting the expansion (61)–(63) into the dissolution rates, balancing at leading order, and recalling that the concentration in the fines is uniform in $r_f$ reveals that

(65) \begin{align} c_{b,0}^{\text{(i)}}(z,r_b,t)|_{r_b=1}=c_{f,0}^{\text{(i)}}(z,t)=c_{l,0}^{\text{(i)}}(z,t) \end{align}

Balancing leading order terms in (23) leads to

(66) \begin{align} \frac {\phi _l}{\phi _T} \frac {\partial c_{l,0}^{\text{(i)}}}{\partial t} + \frac {\partial {\mathcal F}_{l,0}^{\text{(i)}}}{\partial z} = {{G}_{b,0}^{\text{(i)}} + G_{f,0}^{\text{(i)}}}, \qquad \quad {\mathcal F}^{\text{(i)}}_{l,0} = q c^{\text{(i)}}_{l,0}. \end{align}

Eliminating ${\mathcal F}_{l,0}$ and $G_{f,0}$ from (66a) using (66b) and (64) leads us to an equation for conservation of solubles in the region $0\le z \le s_d(t)$

(67) \begin{align} \left ( \frac {\phi _l}{\phi _T} + {\frac {1}{3Q_f}} \right ) \frac {\partial c^{\text{(i)}}_{l,0}}{\partial t} + q \frac {\partial c^{\text{(i)}}_{l,0}}{\partial z} = {G_{b,0}^{\text{(i)}}} \quad \mbox{with} \quad c^{\text{(i)}}_{l,0}|_{z=0}=0. \end{align}

Dissolution from a boulder begins when the saturation front passes (i.e., when the boulder enters region (i)) and it is exposed to non-saturated water for the first time. We denote $t_0 = s_d^{-1}(z)$ to be the time at which the saturation interface passes a given location $z$ . Thus $t_0(z)$ is the first time at which the boulder at position $z$ becomes wet. The pores within the boulders entering region (i) have been filled with saturated liquid, see (37), and hence their internal concentration on entering region (i) is given by (60), in agreement with the matching condition (47). Thus, the leading order transport problem within the boulders defined for $0\le r_b \le 1$ , $t \ge t_0(z)$ is

(68) \begin{align} \frac {\partial c_{b,0}^{\text{(i)}}}{\partial t} = \frac {1}{r_b^{2}} \frac {\partial }{\partial r_b} \left ( r_b^{2} D_{sb} \frac {\partial c_{b,0}^{\text{(i)}}}{\partial r_b} \right ),\end{align}

(69) \begin{align} -D_{sb} \frac {\partial c_{b,0}^{\text{(i)}}}{\partial r_b} \Big {|}_{r_b=0} = 0, \quad c_{b,0}^{\text{(i)}}|_{r_b=1}=c_{l,0}^{\text{(i)}} \quad \mbox{and} \quad c_{b,0}^{\text{(i)}}|_{{t=t_0(z)}} = c_{b,\text{init}} + \varphi _{lb}.\end{align}

(70) \begin{align} -D_{sb} \frac {\partial c_{b,0}^{\text{(i)}}}{\partial r_b} \Big {|}_{r_b=1} = Q_bG_{b,0}^{\text{(i)}} .\end{align}

Equations (68)–(70) can be viewed as defining a relationship between the leading order concentration in the liquid $c_{l,0}$ and the leading order flux of solubles at the boulder’s surface $G_{b,0}$ . In fact, because the transport problem is linear this relationship can be determined analytically (details are given in the supplementary material); it reads as follows:

(71) \begin{align} G_{b,0}^{\text{(i)}}(z,t) &= \frac {2 D_{sb}}{Q_b} \sum _{n=1}^{\infty } e^{-n^2 \pi ^2 D_{sb}(t-t_0)}\left (c_{b,0}^{\text{(i)}}|_{t=s_d^{-1}(z)} -c_{l,0}^{\text{(i)}}(z,t_0) \right .\nonumber \\ &\left .-\int _{t_0}^{t}\dot {c_{l,0}}^{\text{(i)}}(z,\tau )e^{n^2 \pi ^2 D_{sb}\tau } \text{d}\tau \right ), \end{align}

The combination of (67) and (71) is a closed problem to be solved for the leading order concentration in the liquid, which can be written elegantly under the change of variables

(72) \begin{align} \chi = z, \qquad \tau = -\frac {1}{q} \left ( \frac {\phi _l}{\phi _T} + {\frac {1}{3Q_f}} \right ) z + t, \qquad \tilde {c}(\chi, \tau ) = c_{l,0}^{\text{(i)}}(z,t), \end{align}

as

(73) \begin{align} \frac {\partial \tilde {c}}{\partial \chi } = {\frac {\tilde {G_b}(\chi, \tau )}{q} }, \qquad \tilde {c}|_{\chi = 0} = 0, \end{align}

where $\tilde {G_b}(\chi, \tau ) = G_{b,0}^{\text{(i)}}(z,t)$ , i.e., that arrived at by applying the variable transform (72) to (71). This change of independent variables amounts to the application of the method of characteristics; the PDE (66) is reduced to an ODE, namely (73), that is satisfied along the characteristic lines. Replacing an advection-reaction PDE (66), with a relatively simple ODE leads to a significant reduction in mathematical complexity and in turn to reduced computational cost when the resulting problem is solved numerically. Furnishing solutions to this problem (72)-(73) requires a numerical method. Once $\tilde {c}$ has been found, the leading order concentrations in the grains and liquid can be found readily by backsubstitution.

4.5 The asymptotically reduced model

In the wake of the asymptotic analysis we now summarise the equations that must be solved in order to predict the quantity of primary interest, namely, the concentration of the liquid exiting the espresso bed. Rearranging (58) yields the following ODE to be solved for the position of the saturation front, in terms of the concentration in the liquid at the edge of region (i):

(74) \begin{align} \frac {ds_d}{dt} = q \left ( \frac {\phi _l}{\phi _T} + {\frac {1}{3Q_f}} \frac {c_{f,\text{init}}-c^{\text{(i)}}_{l,0}|_{z=s_d(t)}}{1-c^{\text{(i)}}_{l,0}|_{z=s_d(t)}} \right )^{-1} \quad \mbox{with} \quad s_d|_{t=0}=0. \end{align}

This must be solved in tandem with (73), an ODE for the concentration in the liquid, until such a time that $s_d(t)=1$ . Beyond this instance, the concentration at the exit is determined wholly by (73). Thus, the concentration at the exit (to leading order) for all time is given by

(75) \begin{align} c_{\text{exit}} = \begin{cases} 0 \quad \text{for} \quad 0\lt t\lt 1,\\ 1 \quad \text{for} \quad 1\lt t\lt s_d^{-1}(1),\\ \tilde {c}|_{\chi =1} \quad \text{for} \quad t\gt s_d^{-1}(1). \end{cases} \end{align}

We present a numerical method for solving this problem doing so in the subsequent section.

We return to the question raised in Section 2.7 of whether a saturated layer forms immediately above the wetting front, i.e., whether regions (ii) and (iii) form. In retrospect, we now observe that if on solving the asymptotically reduced model $s_d\gt s_w$ , at any time, then at that time the asymptotic structure of the model that we have presented breaks down and regions (ii) and (iii) along with the saturated layer cease to exist. In such a scenario, the solution throughout the wet part of the bed is unsaturated and must be found by solving the problem described in the supplementary information.