3D Pore-Scale Mixing Interface Evolution

All work presented here is based on solute transport data obtained by the simulations of Sole-Mari et al. [34]. We provide a brief description here but direct the interested reader to that paper for full details. In order to gain a better understanding of experimental observations of mixing such as those of Gramling et al. [16], Sole-Mari et al. [34] built a digital random porous medium made of spherical solid grains of uniform diameter $d_{0}$, resembling a sand-filled column, and simulated flow and transport using the OpenFOAM suite [28]. The medium was created by letting grains settle by gravity using the software Blender [6]. The interstitial space was meshed with a cubic regular mesh with a cell size of $d_{0}/60$ which was then transformed into an unstructured one to capture grain-fluid interface geometries accurately. Figure 1(a) shows the full column dimensions, an example portion of it, and a finite volume mesh. This mesh was then used to solve the steady-state, incompressible Navier-Stokes equations,

to obtain the velocity field at the pore scale. Here, $\vec{u}$ denotes the velocity vector, $\rho$ is the fluid density, $\nu$ is the kinematic viscosity, and $\nabla p$ is the pressure gradient. No-slip boundary conditions are imposed at the fluid-grain boundaries, while a full-slip boundary condition is imposed on the column lateral boundaries so as to minimize their influence on flow and ultimately transport. The advection-diffusion equation was then solved to simulate the transport of a solute,

where the solute concentration is denoted by $C$, and the diffusion coefficient by $D$.

The function $C(t=0)=1-H(x-x_{0})$ was used as the initial condition for the transport simulation to mimic an initial sharp flat interface between two solutions near the inlet, where $H(x)$ is the Heaviside step function. No-flux boundary conditions were imposed at the fluid-grain boundaries and at the column lateral walls. The simulation was run for a sufficiently long time for the invading fluid to arrive at the column exit, whereupon it is terminated. Six Péclet number cases ($Pe={10,32,100,316,1000,3160}$) were simulated by modifying the molecular diffusion coefficient, where $Pe=\overline{u}d_{0}/D$, where $\overline{u}$ is the average fluid velocity and $D$ the diffusion coefficient. A full description of the numerical simulation setup is given in [34].

The interface between two solutions at the pore scale is not expected to be sharply defined as diffusion leads to a transition zone rather than an abrupt boundary. The concentration gradually shifts from $1$ to $0$. We identify the mixing interface as the iso-concentration contour for the midpoint (i.e., $0.5$) concentration value. A representative example of such an isosurface is shown in Figure 2.

The 3D isosurfaces are generated by interpolating the concentration field using Matlab’s isosurface function which is based on the Marching Cubes algorithm ([24]). These surfaces are described by a fine triangular mesh that is able to resolve structures at smaller scales than the pore size (Figure 2(c)).

The complex filamentary structures that arise from this type of analysis (Figure 2(c)) are typically referred to as lamellae ([20],[37]), and will be called as such from now on in this paper. Their overall geometry at different times for Péclet numbers 10 and 1000 is depicted in Figure 3(a). As expected, for the higher Péclet number, advective stretching is stronger relative to diffusion and this leads to more elongated lamellae.

We quantify the deformation of the iso-concentration interface by the growth of its total surface area,

where $A_{0}$ is the initial flat interface area at time $t=0$.

By computing this quantity over time and across the range of Péclet simulations available, we obtain the evolution profiles presented in Figure 3(b), where time is normalized by the advection time, defined as $t_{a}=\frac{d_{0}}{\bar{u}}$. The main features that stand out from Figure 3 can be summarized as

1. 1.

   The interface area grows with time but eventually approaches a stagnation plateau.

2. 2.

   At early times, the interface area growth exhibits a faster than linear behavior.

3. 3.

   The plateau, i.e., the late time stagnation area value, increases with Péclet number, but not linearly.

In order to understand this behavior, rather than focusing directly on the system at hand, in the following section we will begin by exploring the mechanisms at play in the idealized context of a single parabolic fluid velocity profile.

In Figure 4(a), a lamella is stretched by the velocity profile $v(x)$, since the middle travels at a higher velocity than the borders. The total stretch rate can be obtained by integrating the velocity gradient $\partial v/\partial x$ over the half-width of the lamella,

where $c_{A}$ is a dimensionless stretching proportionality constant, representing the geometry and spread of a velocity profile acting on the lamella contour. For the two-dimensional lamella representation proposed above, it is $3/2$. This changes, for instance, with dimension and boundary geometry, such that, for a 3D lamella with a circular paraboloidal profile, we have $c_{A}=2$. One should therefore select the suitable value of $c_{A}$ depending on the dimensionality and geometry of the system at hand. Aside from that, the stretching rate is only dependent on the average flow velocity.

In a porous medium, we may assume that lamellae are generated by a similar stretching mechanism to the above, as near-parabolic velocity profiles occur between nearby solid walls, albeit with the added distortion and heterogeneity brought about by the random geometry [10]. However, the parallel plates model’s assumption of constant velocity along streamlines is a crucially unrealistic representation of an actual porous medium. Within the latter, velocity correlation along the flow direction can be statistically characterized by a typical distance $\lambda$. A hypothetical long lamella such that $s\gg\lambda$ would tend to experience a significantly diminished stretching rate. Indeed, given two points of the same lamella separated by a sufficiently large longitudinal distance, their two-point advective stretching rate (their mutual separation velocity) would appear random with near-zero mean due to decorrelation.

To account for this, we introduce a penalty to the stretching rate, which is determined by the ratio of lamella length ($s$) to velocity correlation length ($\lambda$), such that

When lamellae are short ($s<\lambda$) their limited spatial extent experiences a highly correlated velocity profile, resulting in no significant penalty to the stretching rate. By contrast, for longer lamellae ($s>\lambda$) the penalty becomes substantial as the less-correlated velocity field loses its stretching ability. The exact format of the proposed correction of the stretching rate (Equation 7) will be discussed and justified based on the observed data.

Associated with the idealized parabolic lamella profile, we can construct an also idealized 2D concentration field $C$, defined by Equation 8, which is a superposition of two complementary error functions (erfc()) perpendicular to the lamella direction (Figure 4(c)). Such functions were chosen because it is a typical solution for the diffusion equation. They are centered at the lamella contour line $x_{L}$, which is the inverse function of the proposed parabolic profile Equation 4 ($x_{L}(y)=y_{L}^{-1}(x)$). The parameter $Dt_{0}$ sets a non-zero initial horizontal spread to create a continuous concentration field (Figure 4) that smoothes out with time.

Molecular diffusion ($D$) has an impact on the lamella contour line by changing concentration field $C$. We assume that only transverse diffusion contributes significantly to this effect, and this assumption is corroborated by numerical observations in a parallel plate setup. We use Equation 8 to evaluate the change in the lamella contour ($C(x,y_{L},t)=0.5$) after a small change in time $dt$. Figure 5 shows how diffusion changes the lamella contour line at its tip after one and two time steps ($dt$).

By numerically evaluating the change in $s$ after a time increment $dt$ under diffusion, and how this rate is affected by geometric parameters $s$ and $h$, we find the following differential equation:

where $c_{D}$ is a dimensionless shrinking proportionality constant. Since this mechanism generates a negative change in the length, we call it shrinking rate. This scaling was found by dimensional analysis, and in fact, this general diffusive shrinking behavior, other than the particular value $c_{D}$, does not depend on the parabolic shape assumption or lateral boundary conditions. More details on how Equation 9 was obtained can be found in Appendix A.

From the previously derived relationships we may write a differential equation for the lamella length:

The equilibrium state ($\frac{ds}{dt}=0$) is attained when stretching and shrinking rates become equal in magnitude, resulting in a stabilization of both lamella length and contour length, which corresponds to our plateau regime. This equilibrium involves a balance between the mechanisms of stretching driven by heterogeneous advection and shrinking driven by molecular diffusion.

Defining Peclet number as $Pe=\frac{\bar{v}h}{D}$, $s^{*}=s/h$, $\lambda^{*}=\lambda/h$, $c=\frac{c_{A}}{c_{D}}$, and setting Equation 10 to zero, we obtain a dimensionless lamella length at equilibrium as

From the role of $Pe$ in the above equation, we note that the higher the advection forces in comparison to diffusion forces, the longer a lamella can get at equililbrium, aligning with intuitive expectations. Additionally, we can identify two limiting cases:

When the two rates are not equal, we are in a transient regime of the lamella evolution. At early times when $s^{*}=0$, the lamella is expected to grow (Figure 3) because the stretching rate overcomes the shrinking rate. Defining non-dimensional time $t^{*}=t\,\bar{v}/h$ we may rewrite Equation 10 as:

which can be integrated to obtain the transient $s^{*}(t)$ for any $Pe$. While we cannot obtain a general analytical solution, we can do so for specific limiting cases (similar to the two equilibrium limiting cases shown above):

In the first limit (Equation 15), lamellae are shorter than the velocity correlation length, which occurs at early times and for lower Peclet numbers. Conversely, the second limit (Equation 16) characterizes the evolution of elongated lamellas, observed at later times and for higher Peclet numbers. The complete solution (obtained numerically) offers a more comprehensive description of the transition between these two regimes.

As noted in Section section 2, the measured mixing interface is an area in a 3D medium and a length in a 2D medium. In order to compare such measurements to our model predictions, one needs to convert lamella length, $s$, to either interface length, $L$, or interface area, $A$.

Based on the parabolic shape described by Equation 4, in two dimensions we have:

where $G$ is the relative interface length growth. By exploiting radial symmetry, the two-dimensional (2D) parabolic lamella can be converted into a three-dimensional (3D) paraboloid. In this case we have

where $G$ becomes the relative interface area growth, and $A_{0}=\frac{\pi h^{2}}{4}$ is the flat initial interface area.

From both equations (17) or (18), two scaling limits can be identified by using the leading terms of series expansions: when the lamella is much shorter than its width ($s^{*}\ll 1$) the interface grows quadratically with $s^{*}$, and when the lamella is much longer than its width ($s^{*}\gg 1$) the interface grows linearly with $s^{*}$.

By combining Equation 14 with Equation 17 or Equation 18, we can estimate interface growth as a function of time ($G(t)$) for any given $Pe$. By combining Equation 11 with Equation 17 or Equation 18, we can estimate the equilibrium interface area as a function of $Pe$, $G_{equil}(Pe)$.

Based on the SPLM model described, we can anticipate the expected scaling regimes for the equilibrium interface area with Péclet (Figure 6(a)). And in the same way, predict regimes for the evolution of the interface with time (Figure 6(b)). Below we describe the expected regimes for $\lambda^{*}>1$, as is the case for our porous medium data presented in section 2. The regimes for the equilibrium interface area are:

1. 1.

   For low $Pe$ cases, characterized by $cPe\ll 1$, $G_{equil}\sim Pe^{2}$.

2. 2.

   For intermediate $Pe$ cases, characterized by $1\ll cPe\ll\lambda^{*}$, $G_{equil}\sim Pe^{1}$.

3. 3.

   For high $Pe$ cases, characterized by $cPe\gg\lambda^{*}$, $G_{equil}\sim Pe^{1/2}$.

1. 1.

   At very early times, characterized by $s^{*}\ll 1$, $G\sim t^{2}$.

2. 2.

   During intermediate times, when $1\ll s^{*}\ll\lambda^{*}$, $G\sim t$.

3. 3.

   At later times, with $s^{*}\gg\lambda^{*}$, $G\sim t^{1/2}$.

4. 4.

   The interface ultimately stagnates, reaching a plateau regime ($s^{*}=s^{*}_{equil},G=G_{equil}$).

However, not all these regimes may manifest clearly under specific conditions. For instance, if the equilibrium is reached in earlier stages ($s^{*}_{equil}<\lambda^{*}$ or $s^{*}_{equil}<1$), as shown in the parallel plates case in Figure 6(a) or if velocity correlation length is not much larger than the width ($\lambda^{*}\sim 1$), as shown in the porous medium case in Figure 6(a). Likewise, Figure 6(b) illustrates the different stages of lamella growth, which follows a vertical line from $G=0$ to $G=G_{equil}$. Depending on the value of $\lambda^{*}$ and on that of $G_{equil}$ (which depends on $cPe$ and $\lambda^{*}$), some of the temporal growth regimes may not occur.

Let us consider a continuous tracer injection into a fixed pressure-driven flow between two parallel plates (2D Poiseuille Flow). The mixing interface (0.5 concentration isoline) develops into a single lamella, as shown in Figure 7(a), which matches the idealized model, yet without the parabolic shape. This setup is a traditional benchmark problem ([10, 29]) which has both numerical and semi-analytical solutions that enable accurate calculations of interface size evolution over time and at equilibrium.

In this configuration, the constant velocities along streamlines imply that we are in the limit case of $\lambda^{*}\rightarrow\infty$. Therefore, the model for the equilibrium interface length is defined by Equations 12 and 17, and the model for its temporal evolution is defined by Equations 15 and 17. From the Stokes flow solution of the velocity profile, we have $c_{A}=3/2$, thus the sole parameter in the model that remains to be determined is the proportionality constants ratio $c=\frac{c_{A}}{c_{D}}$.

We compare the solution of the SPLM with a semi-analytical solution [15] for a wide range of $Pe$, from $1$ to $10^{5}$. Using the equilibrium interface length ($G_{equil}$) relation with Péclet , we determine the proportionality constants ratio value $c=1/32$ to adjust the model to the reference solution Figure 7(b).

The entire transient regime of the interface growth behaves as predicted, with a near-perfect agreement between the model and the semi-analytical solution Figure 7(c) for the entire range of $Pe$ tested.

Both the equilibrium and transient behaviors (Figures 7(b) and 7(c), respectively) show scaling regimes as described in subsection 3.5 and Figure 6, with the exception (in both cases) of regime 3, which cannot be reached since $\lambda^{*}\rightarrow\infty$. For low $Pe$ cases ($Pe<32$), the transient regime 2 in which $G\sim t$ is skipped because $s^{*}_{equil}<1$ and therefore the condition $1\ll s^{*}\ll\lambda^{*}$ never occurs.

It has been verified (not shown here for brevity) that the same type of results and model matching are obtained from a three-dimensional parallel plates (square pipe) setup, after adjusting $c$ adequately, confirming that the applicability of the model does not depend on the system’s dimensions.

In this scenario, we use the numerical simulation data described in section section 2. The SPLM model is applied, using the grain size diameter $d_{0}$ as the representative single lamella width $h$, and assuming a circular Poiseuille flow, which gives paraboloid velocity profile with $c_{A}=2$. Dealing in this case with a 3D problem, we compute the interface area growth by Equation 18 to compare it to the data.

For the velocity correlation length $\lambda$, we estimate it based on the velocity variogram along streamlines, which has been found to be about one grain diameter ($1.2\pm 0.6d_{0}$). This represents the typical distance along a streamline over which velocity retains a certain degree of autocorrelation, and such an estimate is in line with our expectations that it should be on the same order of magnitude as the grain diameter (e.g., [7]).

Figure 8 shows the comparison between SPLM and the porous medium measurements for several values of $Pe$. Due to small fluctuations in the data at late times, we use the average of the last 10 time points as representative of the equilibrium interface area for all cases. In the same way that we did for the parallel plates configuration, choosing an adequate $c$ to fit the equilibrium data ($c=1/18$), the model captures reasonably well the behavior seen in the data (Figure 8(b)). Although significant efforts were undertaken to minimize such effects, since the data is based on a numerical simulation, it can be impacted by numerical dispersion, particularly for our largest Peclet number. An assessment under idealized uniform flow conditions, has allowed us to estimate numerical dispersion for the largest Péclet case to be at least $D_{num}\approx 25\%D_{m}$. This artifact may be interpreted as an effectively simulated Péclet number which is lower than the intended value. This is illustrated in Figure 8(b) by a leftward shift of the respective marker and in Figure 8(a) by the dashed line representing the SPLM model for a Péclet number after correction. For other Péclet numbers ($Pe\leq 10^{3}$), this impact was negligible.

The difference between the lines SPLM($\lambda\rightarrow\infty$) and SPLM($\lambda=1.2$) shows the impact of velocity correlation length and how it is crucial to capture the data. With regards to the transient behavior our model captures the general trend, but broadly does not reflect the early time behaviors as well as one would hope, missing the spread observed among different $Pe$ at very early times (Figure 8(a)). We suspect that there is a sizable impact of pore-space heterogeneity in the early-time advective stretching, such that in the actual porous medium a variety of lamellae are locally generated at various rates and with uneven geometry. In contrast, the model is based on single effective values of lamella width $h$ and velocity $\bar{v}$, thus disregarding heterogeneity, which explains the much lesser impact that $Pe$ has on early-time interface growth compared to the data. This effect becomes less important at later times, when the various lamellae have a statistically similar history and can thus be modeled by a single set of effective parameters. The fact that the plateau is reached in all cases within a few advective times suggests that transient effects wash out fairly quickly and are only important within a few grains of the initial sharp interface.

Finally, we can back up our proposed stretching rate formulation (Equation 7) based on the observed porous media data scalings in Figure 8(b). We note that $G_{equil}\sim\sqrt{Pe}$ for large $Pe$, which means $s_{equil}\sim\sqrt{Pe}$ for $s^{*}\gg 1$. Given that $\frac{\partial s}{\partial t}|_{shrink}\sim s^{1}$ (Equation 9), it must be that $\frac{\partial s}{\partial t}|_{stretch}\sim s^{-1}$ for sufficiently large $Pe$, such that at equilibrium $s_{equil}\sim\sqrt{Pe}$ be satisfied. As expected, the larger the value of Péclet, the longer the lamellae become and a more significant limitation is imposed by velocity decorrelation effects (separation between orange and black lines in Figure 8(b)).

Regarding observed regimes (see Section 3.5 and Figure 6), it is worth noting that regime 2 for the equilibrium ($G_{equil}\sim Pe^{1}$) is skipped or very short. This is due to the fact that the velocity correlation length is very similar to the grain size, or the lamella width, i.e., $\lambda^{*}\sim 1$, hence the condition $1\ll cPe\ll\lambda^{*}$ is very short-lived. The same idea can explain why regime 2 ($G\sim t$ for $1\ll s^{*}\ll\lambda^{*}$) does not appear clearly in the transient behavior. Regime 1, both equilibrium ($G_{equil}\sim Pe^{2}$) and transient ($G\sim t^{2}$), are also difficult to observe in the data due to the numerical limitations for very low $Pe$ and for very short times, respectively.

As noted, the early time predictions from our current model do not capture the spread of behaviors with $Pe$ seen in the numerical data, which we primarily attribute to a heterogeneous distribution of parameters in the actual system as opposed to the simple model. To this end we delve a little deeper to better try and capture early time effects. Applying series expansions on the transient lamella length equations 15 and 16, we notice that the stretching proportionality constant $c_{A}$ controls the lamella growth at early times in our model:

The spread seen at early times in the porous medium data across different Péclet numbers (Figure 8(a)) suggests a dependency of $c_{A}$ on Péclet. Executing a fitting procedure on $c_{A}$ we observe $c_{A}\sim log(Pe)$ (Figure 9(a)). Note that to keep the equilibrium solution, which matches the data this would mean $c_{D}$ must have the same scaling with Péclet ($c_{D}\sim log(Pe)$).

Assuming this relationship comes from the fact that we are modeling an ensemble of lamellae instead of a single one, we can think of this as a filter applied on the original model (Equation 14) that arises due to averaging to represent the porous medium behavior also at early times. In other words

where $<>$ denotes the averaging operator. When accounting for this a reasonable match for all times, including early ones, is obtained as shown in Figure 9(b).

Despite our attempts, at this stage we do not have a mathematical demonstration of how this factor would arise. However, we do note that similar scalings depending on $log(Pe)$ arise naturally in classical studies of upscaling dispersion in porous media such as that off Saffman ([31]) when deriving residence times in a pore-network approximation of an ensemble of particles. We believe that there may be an analogy here.