Multiphase modeling of precipitation-induced membrane formation

TL;DR

This paper develops a multiphase model for membrane formation driven by precipitation reactions in flow, capturing the dynamic growth and fixed positioning of membranes without prior assumptions, applicable to various scientific contexts.

Contribution

It introduces a novel multiphase framework that models the dynamic, precipitation-induced membrane growth and fixed positioning in flowing fluids, verified through an exact solution.

Findings

Model exhibits no-slip behavior on the membrane.

Membrane position emerges dynamically without prior assumptions.

Applicable to hydrothermal vents and fluid-structure interactions.

Abstract

We formulate a model for the dynamic growth of a membrane developing in a flow as the result of a precipitation reaction, a situation inspired by recent microfluidic experiments. The precipitating solid introduces additional forces on the fluid and eventually forms a membrane that is fixed in the flow due to adhesion with a substrate. A key challenge is that the location of the immobile membrane is unknown $\textit{a priori}$. To model this situation, we use a multiphase framework with fluid and membrane phases; the aqueous chemicals exist as scalar fields that react within the fluid to induce phase change. To verify that the model exhibits desired fluid-structure behaviors, we make a few simplifying assumptions to obtain a reduced form of the equations that is amenable to exact solution. This analysis demonstrates no-slip behavior on the developing membrane without $\textit{a priori}$ assumptions on its location. The model has applications towards precipitate reactions where the precipitate greatly affects the surrounding flow, a situation appearing in many laboratory and geophysical contexts including the hydrothermal vent theory for the origin of life. More generally, this model can be used to address fluid-structure interaction problems that feature the dynamic generation of structures.

Equations

$$aA\text{(aq)}+bB\text{(aq)}\to cC\text{(aq)}$$

$$C\text{(aq)}\to C\text{(s)}.$$

$${\text{Ni}}^{2+}\text{(aq)}+2{\text{(OH)}}^{-}\text{(aq)}\to {\text{Ni(OH)}}_2\text{(aq)}$$

$${\text{Ni(OH)}}_2\text{(aq)}\to {\text{Ni(OH)}}_2\text{(s)}.$$

$${\mathcal{M}}_f={\rho }_f{\theta }_s{V}_0=\left({\rho }_s+\sum M_i{\psi }_i\right){\theta }_s{V}_0$$

$${\dot{\psi }}_A^{(aq)}=-ar{\psi }_A{\psi }_B,{\dot{\psi }}_B^{(aq)}=-br{\psi }_A{\psi }_B,{\dot{\psi }}_C^{(aq)}=cr{\psi }_A{\psi }_B.$$

$${\psi }_A+\Delta {\psi }_A=\frac{n_A}{({\theta }_s+\Delta {\theta }_s)V_0}.$$

$${\psi }_A+\Delta {\psi }_A=\frac{n_A}{{\theta }_s{V}_0}\left(1-\frac{\Delta {\theta }_s}{{\theta }_s}\right)={\psi }_A\left(1-\frac{\Delta {\theta }_s}{{\theta }_s}\right).$$

$${\psi }_C+\Delta {\psi }_C=\frac{n_C+\Delta n_C}{({\theta }_s+\Delta {\theta }_s)V_0}.$$

$${\dot{\psi }}_A^{(p)}=-{\psi }_A{\dot{\theta }}_s/{\theta }_s,{\dot{\psi }}_B^{(p)}=-{\psi }_B{\dot{\theta }}_s/{\theta }_s,{\dot{\psi }}_C^{(p)}=-{\psi }_C{\dot{\theta }}_s/{\theta }_s+\alpha {\dot{\theta }}_s/{\theta }_s.$$

$${\dot{\psi }}_A$$

$${\dot{\psi }}_B$$

$${\dot{\psi }}_C$$

$$\Delta {\mathcal{M}}_f=V_0{\rho }_s\Delta {\theta }_s+V_0{\theta }_s\sum _i{M}_i\Delta {\psi }_i+V_0\Delta {\theta }_s\sum _i{M}_i{\psi }_i.$$

$$\Delta {\mathcal{M}}_f=V_0{\rho }_s\Delta {\theta }_s+V_0{\theta }_s{\psi }_A{\psi }_B(c{M}_c-a{M}_A-b{M}_B)+V_0{M}_C\Delta n_C.$$

$$a{M}_A+b{M}_B=c{M}_C.$$

$$\Delta n_C=\left(\frac{{\rho }_m-{\rho }_s}{M_C}\right)\Delta {\theta }_s.$$

$$\frac{\partial \varphi }{\partial t}+\nabla \cdot \mathbf{J}=\Gamma$$

$$\displaystyle\frac{\partial(M_{A}\psi_{A}\theta_{s})}{\partial t}+\nabla\cdot\big{(}M_{A}\psi_{A}\theta_{s}\mathbf{v}_{s}-\kappa_{A}\nabla(M_{A}\psi_{A})\big{)}$$

$$\displaystyle\frac{\partial(M_{B}\psi_{B}\theta_{s})}{\partial t}+\nabla\cdot\big{(}M_{B}\psi_{B}\theta_{s}\mathbf{v}_{s}-\kappa_{B}\nabla(M_{B}\psi_{B})\big{)}$$

$$\displaystyle\frac{\partial(M_{C}\psi_{C}\theta_{s})}{\partial t}+\nabla\cdot\big{(}M_{C}\psi_{C}\theta_{s}\mathbf{v}_{s}-\kappa_{C}\nabla(M_{C}\psi_{C})\big{)}$$

$${\Gamma }_A=-ar{M}_A{\theta }_s{\psi }_A{\psi }_B$$

$${\Gamma }_B=-br{M}_B{\theta }_s{\psi }_A{\psi }_B$$

$${\Gamma }_C=cr{M}_C{\theta }_s{\psi }_A{\psi }_B+\alpha M_C{\dot{\theta }}_s$$

$${\theta }_s+{\theta }_m=1$$

$$\frac{\partial ({\rho }_s{\theta }_s)}{\partial t}+\nabla \cdot ({\rho }_s{\theta }_s{\mathbf{v}}_s)=R_s$$

$$\frac{\partial ({\rho }_m{\theta }_m)}{\partial t}=R_m$$

$$\mathcal{M}(V_0)={\int }_{V_0}{\rho }_s{\theta }_s+{\rho }_m{\theta }_m+\sum M_i{\psi }_i{\theta }_{s}dV$$

$$\frac{d}{dt}\mathcal{M}(V_0)+{\int }_{\partial V_0}\underset{\text{boundary flux}}{\underbrace{\left({\rho }_s{\theta }_s{\mathbf{v}}_s+\sum {\mathbf{J}}_i\right)\cdot \hat{\mathbf{n}}}}dS={\int }_{V_0}\underset{\text{transfer \& reaction}}{\underbrace{R_s+R_m+\sum {\Gamma }_i}}dV.$$

$$R_s=-\frac{{\rho }_s}{{\rho }_m}{R}_m$$

$$\nabla \cdot {\mathbf{q}}_s=0.$$

$$\frac{d\mathcal{M}}{dt}=-\int_{\partial\Omega}\sum M_{i}\big{(}\mathbf{q}_{s}\psi_{i}-\kappa_{i}\nabla\psi_{i}\big{)}\cdot\hat{\bf{n}}\,\,dS\,.$$

$$R_m=\beta {\psi }_C{\theta }_s\mathcal{H}({\psi }_C-{\psi }_C^{\ast })$$

$$\nabla \cdot \mathbf{T}-{\theta }_s\nabla P-\xi {\mathbf{v}}_s=0$$

Full text

I Introduction

One hypothesis for the “origin of life” is that the first biomolecules were formed in undersea hydrothermal vents. In this theory, passive, anisotropic diffusion across a membrane supports the transmembrane gradients necessary for the first biochemical molecules (Martin et al., 2008). An experimental approach to study this theory examines simpler systems in microfluidic reactors which allow for the controlled study of the prebiotic chemistry in hydrothermal vent chimneys (Barge and White, 2017).

Microfluidic devices have become an important tool in modern chemistry and biomedical analytics (Nge et al., 2013). One application is the possibility of a “lab on a chip”, i.e. the miniaturization of chemical separation and analysis procedures onto a disposable device as small as a few square centimeters. The devices are typically made from etched glass or lithographically-processed elastomers and the fluid flow is usually controlled mechanically by external pumps or electrically via electro-osmotic flows (Mark et al., 2010).

Recent studies have used microfluidic methods to form inorganic membranes within Y-shaped devices (Batista and Steinbock, 2015; Wang et al., 2017). The membranes result from chemical reactions between two different solutions that are injected separately but later merge in a long reaction channel that brings the reactants into direct contact. This merging is usually performed under low Reynolds number (Re) conditions and for miscible liquids, such as aqueous solutions of NaOH, and a dissolved metal salt, such as NiCl2. At the reactive interface between these liquids, a precipitate, such as Ni(OH)2, swiftly forms a thin porous membrane (Figure 1). This precipitate reaction typically involves the formation of microscopically small colloid particles and their aggregation or addition to the membrane. This phenomenon is related to so-called “chemical gardens” which consist of thin cylindrical precipitate membranes separating a metal salt solution from silicate or hydroxide solutions (Roszol and Steinbock, 2011; Makki et al., 2009).

Ding et al. (2016) showed experimentally that membrane thickness increases with the square root of time, indicating diffusion-controlled growth. The membrane thickening occurs only in the direction of the metal salt solution (e.g. NiCl2) and not in the direction of the anionic precipitation partner (e.g. OH-), indicating that the membrane is more permeable to anions than cations. This phenomenon has been qualitatively explained by the charged nature of the membrane that suppresses the transmembrane transport of the positive metal ion (Batista and Steinbock, 2015; Wang et al., 2017).

The modeling challenges presented by this experiment involve a confluence of topics that have been studied before, namely ionic reactions (Sircar et al., 2013; Yang et al., 2019), precipitation (Zhang and Klapper, 2010; Agarwal and Peters, 2014), passive diffusion through a membrane (Mouritsen, 2005; Cogan and Chellam, 2009; Ding et al., 2016), and fluid-structure interaction (Vogel, 1996; Childress et al., 2012; Ristroph et al., 2012; Strychalski et al., 2015; Strychalski and Guy, 2016; Moore, 2017; Quaife and Moore, 2018). The particular combination of these aspects provides the opportunity for a new model that captures them all. One key challenge is that the “structure” in this problem is generated dynamically according to equations governing the chemistry. We choose to model the fluid-structure combination as a single multiphase material: one component “fluid” or solvent and one component “structure” or precipitate membrane. Such multiphase models have proven useful in a variety of complex-fluid applications, such as bacterial biofilms (Cogan and Keener, 2004; Cogan and Guy, 2010), tumor-growth (Byrne and Preziosi, 2003; Preziosi and Tosin, 2009; Frieboes et al., 2010; Sorribes et al., 2019), and biological membranes (Magi and Keener, 2017); their formulation is based on averaging momentum and stresses in separated, multi-component fluids (Drew, 1983; Drew and Passman, 2006).

The multiphase framework developed here builds on previous ones (Cogan and Guy, 2010; Zhang and Klapper, 2010; Leiderman and Fogelson, 2011), but with some keys differences that are guided by a combination of physical principles, model simplicity, and the micro-fluidic experiments mentioned above. First, our formulation conserves the total mass of the components — solvent, dissolved species, and precipitate membrane — throughout evolution. In particular, the model accounts for changes in solute concentrations that result from the formation of new membrane and the associated exclusion of solvent volume. This effect is neglected in previous models that treat reaction chemicals as scalar fields distinct from the multiphase material, but is essential for overall mass conservation. The treatment of reaction chemicals as additional components of a multiphase material has been successfully modeled by many (Nunziato and Walsh, 1980; Yang et al., 2019, see) but greatly complicates the analysis, interpretation, and simulation of the governing equations. Second, by making certain choices in the averaging procedure for the multicomponent stress, our formulation becomes equivalent to an incompressible Brinkman system with variable permeability. This equivalence is important for a few reasons. First, it guarantees that the model reduces to the Stokes equations in the fluid limit and to Darcy’s equation in the porous-medium limit (Brinkman, 1949; Durlofsky and Brady, 1987; Hill et al., 2001). In particular, it guarantees that when membrane is fully developed, the interface behaves as an impermeable surface with a no-slip condition on fluid velocity. As demonstrated in section III.2, many existing multiphase models fail to exhibit this behavior (Breward et al., 2003; Cogan and Keener, 2005; Cogan and Guy, 2010; Sorribes et al., 2019), as they were developed primarily for highly permeable systems. Second, the equivalence to Brinkman significantly simplifies the overall structure of the partial differential equation (PDE) system by eliminating certain cross-terms in the stress divergence that arise in other models. This simplification is one key that will allow reduction to a non-trivial case where chemical and phase dynamics can be solved exactly.

To demonstrate that the new framework possesses the desired properties listed above, we consider a simplified system in which the incoming reactant concentrations are held fixed via chemostat (Rubinow, 1975). By assuming parallel flow and neglecting solute diffusion, the governing equations reduce to a planar system of ordinary differential equations (ODEs). This nonlinear system can be linearized around the fixed points, and eigenvalue analysis provides an estimate for the rate at which new membrane forms. Moreover, we find that the equation for the aqueous product is a second-order nonlinear ODE known as the Ricatti equation (Riccati, 1724; Tenenbaum and Pollard, 1985). Exact solutions to the Ricatti equation give explicit formulas for the time dependence of the chemical product and, consequently, the formation of new membrane. Once the membrane dynamics are known exactly, the flow profile can be obtained through the numerical solution of a simple boundary value problem (BVP). Access to the resulting flow profile allows a careful comparison between variants of the multiphase framework. In particular, we demonstrate that the framework developed here properly captures the transition from one-channel to two-channel flow as membrane develops.

The paper proceeds as follows: in section II we develop the governing equations for both the reacting chemicals and multiphase material such that the total mass is conserved. Section III contains analysis and results based on simplifying assumptions. These assumptions generate a reduced form for an idealized scenario which can be solved with a combination of analytic and simple numerical methods. Finally in section IV the predictive power of the model and further applications are discussed.

II Mathematical Model

The model requires the accurate description of several aspects of the experiment, including the flow transport of the two ionic species and their reaction to form a product, the precipitation of the product out of solution, and finally the response of the bulk fluid motion to the dynamically-generated precipitate membrane. Advection-Diffusion-Reaction (ADR) equations are derived for the aqueous chemical concentrations, while the fluid and membrane dynamics are described by multiphase mass and momentum balance equations. In many multiphase models, either constituent can be viscous, viscoelastic, poroelastic, or otherwise. Here, since the membrane adheres to the substrate, it can be treated as an immobile solid, leading to considerable simplifications.

We assume that aqueous reactants and product contribute mass, but not volume, to the fluid phase. The solvent and membrane each have their own distinct mass densities, and any arbitrary control volume can be divided into solvent and membrane volume fractions. The formation of new membrane involves the precipitation of product out of solution and the sequestration of solvent. A key modeling assumption is that the volume of fluid sequestered equals the volume of the resulting membrane. As shown in section II.2, this assumption ensures incompressibility of the phase-averaged velocity field, i.e. the so-called Darcy velocity.

We now detail the model equations. First we derive evolution equations for the reaction of aqueous ionic species, then we list mass balances for all chemical species as well as solvent and membrane phases, and finally we describe the momentum equation for the fluid. In the end we obtain a closed, coupled PDE system governing the chemistry and physics of the system, where total mass is conserved throughout aqueous reactions and phase transitions.

II.1 Model for Chemical Reactions

In this section we derive equations for the chemical reactions. We follow the “nucleation and growth” model of precipitation (Matsue et al., 2018) and separate the reaction into two sequential parts: in the first, two reactants come together to form an aqueous product, and the second describes the aggregation of the aqueous product into a solid precipitate. While many aqueous chemical reactions do not alter the solution volume significantly, the formation of a membrane excludes fluid volume and therefore can alter the local concentration of the dissolved species. Accordingly, our model neglects the volume occupied by the aqueous species but does account for changes in species concentration that are due to the precipitated solid excluding fluid volume. This effect introduces additional terms in the aqueous reaction equations that are required for mass conservation. To our knowledge these additional terms are not accounted for in the multiphase precipitation literature that treats aqueous chemicals as scalar fields distinct from of the multiphase material.

The aqueous reaction is written as a generic net ionic equation

[TABLE]

where $A\text{(aq)}$, $B\text{(aq)}$, and $C\text{(aq)}$ are chemicals in the aqueous phase and $a$, $b$, and $c$ are their respective stoichiometric coefficients; the precipitation reaction is written simply as

[TABLE]

As a concrete example consider the reaction described in the introduction,

[TABLE]

Then $A=\text{Ni}^{2+}$, $B=\text{(OH)}^{-}$ and $C=\text{Ni(OH)}_{2}$ and $a=c=1$, $b=2$.

The aqueous chemicals will be measured with a variable for the number of chemicals per unit solvent volume, i.e. molarity, which we will call $\psi_{i}$ for chemical species $i\in\{A,B,C\}$. Reaction rates depend on a reactants’ molarity, and molarity can change due to two independent factors: either the number of molecules changes due to the aqueous reaction, or the solvent volume changes due to precipitation. Because either one can occur in a precipitation reaction, these two competing effects must be carefully considered when formulating the reaction equations.

We begin by deriving equations for how the aqueous reaction proceeds in a spatially homogeneous environment; later in section II.2 the effects of advective and diffusive spatial fluxes will be added. Suppose the chemicals exist in some aqueous solution of fixed control volume $V_{0}$. The chemicals undergo both the aqueous and precipitation reactions which results in fluid mass and solvent volume being converted to membrane mass and volume (see figure 2). The fluid component has mass

[TABLE]

where $\rho_{s}$ is the constant solvent mass density (without any reactants or products present), $\theta_{s}$ is the solvent volume fraction, and $M_{i}$ is the molar mass of chemical species $i\in\{A,B,C\}$. The summation represents the contribution of the chemical species to fluid mass, so that the fluid mass density $\rho_{f}$ is not constant. The membrane component has mass $\mathcal{M}_{m}=\rho_{m}\theta_{m}V_{0}$ where $\rho_{m}$ is the constant membrane mass density and $\theta_{m}$ is the membrane volume fraction. Physically, membrane mass is composed of both precipitated chemical $C$ and sequestered solvent mass.

The change in $\psi_{i}$ purely due to aqueous reaction, i.e. no precipitation, can be modeled as a second-order kinetics reaction

[TABLE]

where the dot indicates a derivative with respect to time, and $r$ is the rate of aqueous reaction per chemical concentration. More general power laws are sometimes used to model chemical kinetics, but here we use purely second-order kinetics for simplicity (see Chang and Goldsby, 2013, pp. 573 - 575). None of the analysis, however, depends specifically on this choice and the results could be carried forward for other kinetics.

To derive equations for the change in $\psi_{i}$ purely due to precipitation we appeal to ideas from continuum mechanics. The concentration of ions $A$ in the control volume is written as $\psi_{A}=n_{A}/(\theta_{s}V_{0})$ where $n_{A}$ is the number of $A$ ions in $V_{0}$. Note that this formulation makes explicit the dependence of $\psi_{A}$ on both $n_{A}$ and $\theta_{s}$. Consider the change in a small increment of time $\Delta t$. Then the time-dependent variables are updated so that

[TABLE]

Recall that $n_{A}$ is constant during precipitation as only $C\text{(aq)}$ precipitates. Approximating for small $\Delta\theta_{s}$ and neglecting higher-order terms gives

[TABLE]

Then, cancelling the $\psi_{A}$, dividing both sides by $\Delta t$, and letting $\Delta t\to 0$ gives the change in $\psi_{A}$ purely due to precipitate reaction as $\dot{\psi}^{(p)}_{A}=-\psi_{A}\dot{\theta}_{s}/\theta_{s}$. By symmetry, a similar formula holds for $\dot{\psi}^{(p)}_{B}$. Note that both of these are essentially applications of the product rule for $\partial_{t}(\psi_{i}\theta_{s})=0$, which physically means that the total number of ions of $i\in\{A,B\}$ in the control volume does not change in time due to precipitation.

A similar procedure can be followed for $\psi_{C}$, except now the number of aqueous chemicals $n_{C}$ changes as $C\text{(aq)}$ precipitates into membrane,

[TABLE]

Above, both $n_{C}$ and $\theta_{s}$ change in time. Expanding both expressions while linearizing for small $\Delta\theta_{s}$, dividing by $\Delta t$, and taking the limit as $\Delta t\to 0$ one obtains $\dot{\psi}^{(p)}_{C}=-\psi_{C}\dot{\theta}_{s}/\theta_{s}+\dot{n}_{C}/\theta_{s}$. The first term in this expression is analogous to those obtained for reactants $A$ and $B$, and simply describes the effect on concentration when solvent volume is changing. The second term, however, is new and describes the effect on $\psi_{C}$ as aqueous $C$ molecules are converted into membrane. We write $\dot{n}_{C}=\alpha\dot{\theta}_{s}$ where the specific value of $\alpha$ will be found shortly to guarantee conservation of mass throughout the entire reaction. The expressions for the rate of change of aqueous chemical concentrations due to precipitation are thus:

[TABLE]

Assuming that the aqueous and precipitate reactions act independently, $\dot{\psi}_{i}=\dot{\psi}^{(aq)}_{i}+\dot{\psi}^{(p)}_{i}$, gives

[TABLE]

These equations describe the dynamics of aqueous species concentrations in the absence of spatial fluxes.

To obtain the value of $\alpha$ that guarantees conservation of mass, we again apply a continuum mechanics argument. The change in membrane mass after a small time step is $\Delta\mathcal{M}_{m}=\rho_{m}\Delta\theta_{m}V_{0}$. To simplify the expression for change in fluid mass, we expand $\Delta\mathcal{M}_{f}$ while neglecting second order terms to get

[TABLE]

Replacing the $\Delta\psi_{i}$ with their respective differential terms in equations (11) and performing some algebraic manipulation produces

[TABLE]

Conservation of mass during the aqueous reaction (1) implies

[TABLE]

Thus, the term in parenthesis in (13) vanishes. Meanwhile, conservation of mass of the entire system implies $\Delta\mathcal{M}_{m}=-\Delta\mathcal{M}_{f}$, i.e. the mass lost by the fluid equals the mass gained by membrane. Additionally, the assumption that fluid volume is converted perfectly to membrane volume implies $\Delta\theta_{m}=-\Delta\theta_{s}$. Using the respective definitions of $\mathcal{M}_{i}$ and solving for $\Delta n_{C}$ gives

[TABLE]

Diving by $\Delta t$ and taking the limit $\Delta t\to 0$ gives $\dot{n}_{C}=\alpha\dot{\theta}_{s}$ where $\alpha=(\rho_{m}-\rho_{s})/M_{C}$. Physically, this value of $\alpha$ corresponds to the concentration of $C\text{(aq)}$ that must leave the fluid phase during precipitation in order for mass to be conserved.

The reaction equations derived in this section, along with the specific $\alpha$ term, will be used to provide reaction terms for the chemistry mass balance equations, as described in the next section.

II.2 Mass Balance Equations

In the experiments, the aqueous reaction occurs within the flow of a microfluidic device and therefore spatial fluxes must be considered. To describe these fluxes, consider the general conservation law for the chemical mass per unit control volume $\phi=M_{i}\psi_{i}\theta_{s}$,

[TABLE]

where $\mathbf{J}$ is the flux of $\phi$ and $\Gamma$ is a transfer term for the rate that $\phi$ enters the system.

We choose $\mathbf{J}$ to account for advection and diffusion of the chemical concentrations,

[TABLE]

where $\mathbf{v}_{s}$ is the (tracer) velocity of the solvent and $\kappa_{i}$ are diffusion coefficients which possibly depend on the solvent volume fraction. Note that the diffusive flux used above transports mass according to gradients in molarity $\psi_{i}$, not gradients in $\phi_{i}$. This choice produces the physically realistic steady state of uniform molarity in a quiescent, non-reacting fluid that has inhomogeneous volume fraction.

Assuming that the reactions and spatial fluxes act independently, the $\Gamma_{i}$ correspond to the rates given in equations (11). Rearranging and multiplying each equation by its respective molar mass $M_{i}$ gives

[TABLE]

where $\alpha=(\rho_{m}-\rho_{s})/M_{C}$. Now that mass balance equations for the chemistry are established, mass balance equations for the multiphase solvent-membrane system are needed.

A simple but necessary assumption is that our volume is occupied by only solvent and membrane, i.e. there are no “voids”. This no-void assumption implies

[TABLE]

everywhere. Mass balances for the solvent and membrane phases provide

[TABLE]

where $R_{i}$ denotes the rate of mass added to phase $i$. Equation (23) has no advective term since the membrane is assumed to be immobile.

To ensure conservation of total mass, the rates $R_{m}$ and $R_{s}$ must be related. To derive this relationship, let $V_{0}$ be an arbitrary control volume. The total mass (of all components) in $V_{0}$ is

[TABLE]

Summing the five mass balance equations, (17)–(19) and (22)–(23), integrating over $V_{0}$, and applying the divergence theorem gives

[TABLE]

where $\hat{\bf{n}}$ is the outward unit normal vector. Summing equations (20) and applying (14) gives $\sum\Gamma_{i}=\alpha M_{C}\dot{\theta}_{s}$. For the sake of obtaining a relationship between $R_{s}$ and $R_{m}$, briefly consider the case of zero boundary flux. Then to conserve total mass within any control volume, we must have $R_{s}+R_{m}+\alpha M_{C}\dot{\theta}_{s}=0$ holding point-wise. Substituting the value of $\alpha$, using the no-void assumption (21) and the definition of $R_{m}$ in (23) gives

[TABLE]

This relationship is required for overall mass conservation.

Now consider the so-called Darcy velocity field $\mathbf{q}_{s}=\theta_{s}\mathbf{v}_{s}$. Substituting the value of $R_{s}$ from equation (26) into equation (22), using a consequence of the no-void assumption ($\dot{\theta}_{s}=-\dot{\theta}_{m}$) and substituting $R_{m}$ by its value in equation (22) implies that the Darcy velocity field is in fact incompressible

[TABLE]

Now return to (25) and consider the entire domain $\Omega$ with total mass $\mathcal{M}=\mathcal{M}(\Omega)$. From the above argument, the right-hand-side of this equation vanishes as a necessary condition on mass conservation. Then applying incompressibility (27), gives the total mass balance

[TABLE]

As expressed in this equation, the total mass of the system is conserved as long as the chemical flux at the boundary vanishes. More generally, the total mass of the system can change according to how much chemical mass is being injected or removed via the boundary flux terms.

We now specify our choice for the form for the precipitation term $R_{m}$. Although complicated models of precipitation exist (Matsue et al., 2018; Ostapienko et al., 2018), we employ a simple model in which the rate of membrane mass growth is proportional to the amount of product, provided that the product concentration exceeds some precipitation threshold, i.e.

[TABLE]

where $\mathcal{H}$ is the standard Heaviside function and $\psi_{C}^{*}$ is the concentration threshold for precipitation to occur.

II.3 Momentum Balance Equations

Since inertial effects are assumed negligible ($\text{Re}\ll 1$), the solvent momentum balance can be written as

[TABLE]

where $\mathbf{T}$ is the multiphase stress tensor, $P$ is a so-called common pressure that is shared by the fluid and membrane phases (Cogan and Guy, 2010), and $\xi$ is a friction coefficient. Deviating slightly from the predominant multiphase literature, we chose the form of the stress tensor as

[TABLE]

where $\eta$ is the fluid viscosity. In particular, since $\mathbf{q}_{s}=\theta_{s}\mathbf{v}_{s}$, we have placed the fluid volume fraction $\theta_{s}$ inside the gradient, whereas most multiphase models place the $\theta_{s}$ outside of the gradient but inside the divergence (Byrne and Preziosi, 2003; Cogan and Keener, 2005; Cogan and Guy, 2010). Such a choice must be made for model closure, and neither is fully justified by first principles. We make the above choice to obtain equivalence to the Brinkman system, which offers considerable gains in model tractability.

The incompressibility of the Darcy velocity, while being notable in itself, allows the momentum equation to be transformed into something more familiar. Applying the divergence-free property to equation (30) reduces it to a Brinkman equation with variable coefficients

[TABLE]

We can now see that equation (30), equivalently equation (32), does not have any cross-derivative terms that appear in Cogan and Guy (2010); Keener et al. (2011). Because the membrane is assumed immobile ($\mathbf{v}_{m}\equiv\mathbf{0}$), no momentum equation is needed for it.

The friction coefficient $\xi$ should be chosen in such a way that, at high membrane volume fraction, friction becomes the dominant effect in equation (32). The choice made here, and mentioned briefly in Leiderman and Fogelson (2013), is to use the Kozeny-Carman (KC) formula for permeability as it depends on porosity (Dullien, 2012). In the present notation, the KC relationship gives the friction coefficient as

[TABLE]

where $h$ is an arbitrary constant. This friction coefficient will provide the desired no-slip behavior in the membrane limit $\theta_{m}\to 1$. Angot (1999) discusses the implications of a similar singular friction term in a Brinkman system, although their model does not include the solvent volume fraction term in front of the pressure gradient and only applies to domains with spatially discontinuous volume fractions; the current fluid-membrane model generalizes this notion by being able to account for smooth spatial and temporal gradients of the volume fraction.

As an alternative to the singular $\xi_{KC}$, the friction coefficient can be chosen to be a (non-singular) Hill function, as used in Leiderman and Fogelson (2011, 2013),

[TABLE]

The use of Hill functions is largely empirical, although it has significant advantages in that it is finite in the membrane limit and therefore more numerically stable. Additionally, $K$ determines the half-saturation point and $n$ indicates the qualitative manner at which this saturation is achieved. These parameters allow for fine-tuning to specific experimental observations.

Finally, a friction term employed in many biofilm multiphase models is

[TABLE]

This choice has become popular in the literature (Byrne and Preziosi, 2003; Cogan and Keener, 2004, 2005; Cogan and Guy, 2010; Sorribes et al., 2019), and is often justified by the idea that friction should vanish if either phase, $\theta_{s}$ or $\theta_{m}$, is absent. While it is an intuitive notion, we will show in section III.2 that this friction coefficient does not produce the physically realistic behavior of no-slip velocity on fully developed solid surfaces. To produce this behavior, it is necessary that friction dominates, not vanishes, in the limit $\theta_{m}\to 1$.

A visual comparison of these three friction coefficients is shown in figure 3. For the sake of comparison, we have chosen the constant $h$ so that the three curves intersect at a reference porosity $\theta_{s}^{*}$, i.e. $\xi(\theta_{s}^{*})=\xi^{*}$ in each case. For $\xi_{KC}$ and $\xi_{H}$, the value $\theta_{s}^{*}$ can be interpreted as the percolation threshold, i.e. the critical porosity below which the medium essentially behaves as impermeable to flow (Golden, 1997).

II.4 Model Summary

The governing equations are now summarized for the reader. Rearranging equations (17)–(19) and applying incompressibility of $\mathbf{q}_{s}=\theta_{s}\mathbf{v}_{s}$ gives the following ADR evolution equations for aqueous chemical concentration:

[TABLE]

where $\alpha=(\rho_{m}-\rho_{s})/M_{C}$. Meanwhile, the evolution equations for the fluid and solid phases can be summarized as

[TABLE]

The momentum equation is a Brinkman equation with variable permeability

[TABLE]

where we have employed the $\xi_{KC}$ friction term. We will now consider these coupled PDEs in a simplified setting that will permit exact solutions.

III Analysis of Reduced Model

Analysis of any complicated system is aided by reduction into a form that is analytically tractable. Inspired by microfluidic experiments, we assume that some chemostat controls the influx of reactants’ molarity far upstream. All variables are kept constant along the longitudinal axis by neglecting diffusion and assuming parallel flow. This requirement on parallel flow also means that the reaction takes place everywhere along the longitudinal axis simultaneously. Finally, we assume that the precipitation threshold is negligible.

Applying these assumptions to the governing equations reduces the system considerably so that it becomes a Poiseuille analysis; these assumptions generate the following reduced system

[TABLE]

where only $\psi_{C}$, $\theta_{s}$, $\theta_{m}$ and $q_{y}$ are unknown. Note that the longitudinal axis is chosen to be $y$ such that only this component of the Darcy velocity $\mathbf{q}_{s}=q_{x}\mathbf{\hat{x}}+q_{y}\mathbf{\hat{y}}$ remains. In section III.1 we obtain an analytic estimate on the rate of membrane formation by treating $\psi_{C}$ and $\theta_{m}$ as a planar dynamical system. Then in section III.2 we solve equations (43) with a combination of analytic and numerical methods to visualize how the membrane affects the flow profile in time.

III.1 Fixed Point Analysis

Our approach is to linearize the reduced model system, then perform an eigenvalue analysis about a steady state fixed point. The benefit of this is the eigenvalue will correspond to the rate that membrane develops, a quantity that is possible to measure experimentally.

Before doing a fixed point analysis, it is helpful to understand the conditions on which the existence and stability of fixed points depends. To do so, eliminate the explicit dependence of $\dot{\psi}_{C}$ on volume fraction and replace all $\dot{\theta}_{s}$ in equation (43a) with equation (43c) to obtain a quadratic ODE of the form

[TABLE]

which is an ODE in time alone, as the $x\text{-dependence}$ of $\psi_{A}$ and $\psi_{B}$ are determined by the initial conditions. Examining the qualitative behavior of this ODE by considering $\dot{\psi}_{C}=\dot{\psi}_{C}(\psi_{C})$, it is quadratic in $\psi_{C}$, intercepts the $\dot{\psi}_{C}$ axis at $cr\psi_{A}\psi_{B}\geq 0$, is concave up, and has equilibria at

[TABLE]

whose existence depends on the sign of

[TABLE]

If $\chi>0$, equation (44) will have two fixed points, for $\chi=0$ these fixed points coalesce, and for $\chi<0$ there are no fixed points and $\dot{\psi}_{C}$ will grow without bound; see figure 4(a).

We now ask the question of whether fixed points exist in the reduced system, i.e. $\chi\geq 0$ or $\alpha^{2}\beta\geq 4cr\rho_{m}\psi_{A}\psi_{B}$? We interpret this condition based on the physical meaning of the parameters: $\alpha=(\rho_{m}-\rho_{s})/M_{C}$ has dimension of molarity and is $O(10\,\text{M})$ where M refers to molar units, $1\,\text{M}=\text{1 mol}/\text{liter}$; for example, using the reaction system in the introduction gives $\alpha\approx 30\,\text{M}$. We note that a similar analysis also justifies neglecting the precipitation threshold $\psi_{C}^{*}$, as $\psi_{C}^{*}\approx 0.001\text{ M}$. Although both $r$ and $\beta$ scale the rates of the aqueous and precipitate reactions, respectively, they have different units. $r$ has units of volume per time while $\beta$ has units of mass per time. Because experimental values of $r$ and $\beta$ are expensive to acquire, for the sake of this simplified analysis we will assume that $\beta\approx r\rho_{m}$ such that their effects don’t impact the sign $\chi$. The stoichiometric coefficient $c$ for $C\text{(aq)}$ can be assumed $O(1)$. Finally, examine the reactants $\psi_{A}$ and $\psi_{B}$. Most experiments in microfluidic chambers use molar concentrations with an upper bound of $O(1\,\text{M})$; for example, in Ding et al. (2016) the maximum concentration of reactants was $0.5\,\text{M}$. Therefore, using parameter values taken from experiments, $\chi>0$ and fixed points exist for the reduced system.

We now consider the planar dynamical system in phase space $(\psi_{C},\theta_{m})\in[0,\infty)\times[0,1]$ with fixed points $(\Psi_{C},\Theta_{m})=(\psi_{C}^{-},1)$. The dynamical system is

[TABLE]

where $r$, $c$, $\rho_{m}$, $\alpha$, $\beta$, $\psi_{A}$, and $\psi_{B}$ are assumed to be known and constant. The eigenvalues of the Jacobian generated by equations (47) evaluated at the fixed points provides information about the rate of growth of $\psi_{C}$ and $\theta_{m}$. This particular eigen-system is simple to interpret because the eigenvectors align with the coordinate axes and therefore the eigenvalues correspond to the rates that the physical variables $(\psi_{C},\theta_{m})$ approach their equilibria. These rates are given by:

[TABLE]

Both eigenvalues are negative, and therefore the fixed point is stable, because $\alpha\beta/\rho_{m}+\lambda_{\psi_{C}}=\alpha\beta/\rho_{m}-\sqrt{\chi}/\rho_{m}>0$ always. For this same reason, it is true that $|\lambda_{\theta_{m}}|<|\lambda_{\psi_{C}}|$, meaning the membrane volume fraction approaches its fixed point at a slower rate than the aqueous product. This agrees with our intuition, as the conversion of $C\text{(aq)}$ to $C\text{(s)}$ means we would expect $\theta_{m}$ production to lag behind $\psi_{C}$.

III.2 Poiseuille Analysis

We now solve the equations to see how the solvent velocity transitions from one-channel to two-channel flow. An exact solution for $\psi_{C}(x,t)$ can be found by solving a Ricatti equation with constant coefficients. This solution $\psi_{C}$ can be integrated using elementary functions, so $\theta_{s}(x,t)$ can be found exactly because equation 43c is separable; $\theta_{m}(x,t)$ can then be found using the no-void assumption. Finally, we solve the variable-coefficient BVP for $q_{y}$ numerically using finite differences.

For initial conditions, let $\psi_{C}(x,t)=0$ and $\theta_{s}(x,t)=1$. Fix $\psi_{A}^{0}(x)$ and $\psi_{B}^{0}(x)$ to be piecewise-constant values in $x$ such that there is only a middle region in $x\in(0,L)$ in which the reactants $A$ and $B$ are simultaneously present. Given $\psi_{A}^{0}(x)$ and $\psi_{B}^{0}(x)$, the Ricatti equation (43a) can be solved exactly for $\psi_{C}(x,t)$ (see appendix A):

[TABLE]

where

[TABLE]

Because both the antiderivative of $\psi_{C}$ can be given in terms of elementary functions (see appendix A) and equation 43c is separable, we can also solve for solvent volume fraction $\theta_{s}(x,t)$ exactly:

[TABLE]

where have implemented the initial condition $\theta_{s}(x,0)=1$. Then, the membrane volume fraction $\theta_{m}(x,t)$ can be computed easily using the no-void assumption.

Until this point, all solutions in space have been treated independently. The effect of variations in space is taken into account when solving for the longitudinal component of the Darcy velocity, $q_{y}$. We numerically solve the momentum equation (43d) for $q_{y}$ using a finite difference method. The $x$ domain is discretized into $N$ intervals of equal width $\Delta x=1/N$ such that $x_{j}=j\Delta x$, $j=1,\dots,N-1$, and use centered-difference approximations to the derivatives. Note that $q_{y}(x_{0})=q_{y}(x_{N})=0$ due to the no-slip boundary conditions. This discretization results in a tridiagonal linear system which can be solved in $O(N)$ complexity by using the Thomas algorithm (see Strikwerda, 2004, pp. 78-79). The velocity is constrained to satisfy constant flux in accordance with experiments, which mathematically is represented by

[TABLE]

This constant-flux condition allows the computation of $G(t_{n})$ at each time step. The authors use the Julia programming language to solve the BVP (Bezanson et al., 2017).

The developing membrane for the 1D reduced model geometry is shown in figure 5(a). Membrane develops within the shaded region, which in this example is 10% of the domain. Because the membrane has finite width, the constant-flux condition causes the pressure gradient $G(t)$ to increase with the developing membrane, causing the maximum speed for the two-channel flow to be slightly higher than the maximum speed for the one-channel flow. In this sense, the developing membrane splits the domain into two symmetric one-channel flows.

Figure 5(b) shows the three main variables of the reduced model as functions of time, evaluated in the middle of the reaction region ($x=L/2$). The $\psi_{C}$ variable increases immediately due to the presence of $\psi_{A}$ and $\psi_{B}$. The membrane initially has zero growth rate due to the absence of $\psi_{C}$, and grows at a slower rate than $\psi_{C}$. This ordering on the growth rates matches our expectations from the eigenvalue analysis of section III.1. The $q_{y}$ curve demonstrates the transition from one-channel to two-channel pipe flow by measuring the normalized value in the middle of the pipe as a function of time. By comparing $q_{y}$ with $\theta_{m}$, one can see the effect of the percolation threshold $\theta_{s}^{*}=0.3$. After the solvent volume fraction declines past this value, the fluid velocity begins to react more strongly to precipitating membrane.

Figure 6 demonstrates the effect of using different friction coefficients. Both figures 6(a,b) demonstrate the desired no-slip behavior in the membrane limit by using $\xi_{KC}$ and $\xi_{H}$, respectively. In figure 6(c) the effect of the friction coefficient $\xi_{B}$ is shown. While there is some effect on the flow profile, $\xi_{B}$ does not demonstrate the no-slip condition on the membrane. While the $\xi_{B}$ term was developed primarily for high permeability applications, our framework was developed to capture the transition from purely fluid behavior, to partially permeable, to a fully-developed impermeable solid. As demonstrated, this full transition requires either the $\xi_{KC}$ or $\xi_{H}$ friction coefficient.

Figure 7 shows three flows with reaction regions of various sizes, and therefore different width of membranes. The initial flow profiles of all are equivalent, as the reaction has not yet occurred and no membrane is present. However, as membrane develops, the constant-flux condition requires that for regions with thicker membranes, the flow velocity must increase in the non-reacting regions to compensate for the loss of flux in the membrane region. These results demonstrate that, once the membrane is fully developed, the flow domain treats the membrane portion as a no-slip boundary and the prescribed constant-flux conditions lead to the expected results from single-phase fluids.

IV Conclusion

Typically, no-slip conditions on boundaries must be prescribed a priori when modeling fluid flows. However, motivated by microfluidic experiments of a precipitate reaction, there exist situations where solid materials develop dynamically and so the exact location of such a boundary is not known a priori. This situation exposes a deficiency in modeling techniques that rely on knowledge of where no-slip boundaries exist.

To formulate the model, we assumed conservation of mass for the entire chemistry-fluid-membrane system. The reacting chemical species were assumed to exist in the fluid phase as volumeless scalar fields subject to mass flux determined by a combination of advection with the fluid flow and diffusion down gradients in molarity. The careful consideration of changing solvent volume inside the domain led to extra terms being included in the reaction equations, and in particular conservation of mass for the entire system was used to determine the closure for the aqueous product’s reaction equation. A momentum equation for the fluid velocity followed the usual Cauchy stress formulation with the slight modification that all tensor quantities depend on the Darcy velocity $\mathbf{q}_{s}=\theta_{s}\mathbf{v}_{s}$, not simply the tracer velocity $\mathbf{v}_{s}$ as is usual in the single- multi-phase fluids literature. This choice, paired with incompressibility of the fluid Darcy velocity, led to a simplification of the momentum balance to a variable-coefficient Brinkman equation.

In order to demonstrate that the model reflects the expected behaviors, we first sought growth rate estimates from a linearized form of reduced equations. Then, to show that the dynamically-generated no-slip boundary corresponded to development of the membrane, we performed what was essentially a Poisseuille analysis on the reduced model to show that the no-slip behavior agreed qualitatively with expected behaviors, specifically the recovery of one-to-two channel transitions and the effect of membrane width paired with a constant-flux condition.

We examined three potential terms to be used in the variable-coefficient Brinkman friction term, both by a direct comparison as functions of fluid volume fraction and by examining their effect on fluid flow from the perspective of how well they demonstrated no-slip behavior. These results demonstrated that, qualitatively, the friction terms derived from Karman-Cozeny relationship and using Hill functions gave no-slip flow behavior on the developed membrane. Additionally, the percolation threshold $\theta_{s}^{*}$ can be chosen to reflect specific permeability properties of the structure under investigation. Both of these friction coefficients were preferable to the term commonly employed in biofilm models because $\xi_{B}$ does not recover the no-slip behavior in low-permeability regions.

To our knowledge, this is the first model that qualitatively captures the fluid-structure dynamics of a precipitate reaction in a low-Re environment where the dynamically-developing precipitate significantly affects the surrounding fluid flow. Future work will focus on numerical simulation of the full model in various geometries to be used as a predictive tool for experimentalists. In particular, exploring sufficient modeling conditions to generate the asymmetric growth in membrane in a 2D setting analogous to the experimental microfluidic domains was something that we were not able to explore in this paper, as diffusion was ignored to reduce the model equations. A model capable of accurately capturing microfluidic experiments would be valuable to researchers using these devices to study precipitate reactions at the microscale and ultimately useful in examining origin of life theories.

Acknowledgements

P.S.E. is supported by the National Science Foundation (NSF) Graduate Research Fellowship under Grant 1449440. M.N.J.M. is supported by Simons Collaboration Grant 524259. N.G.C. is supported by NSF-CBET 1510743. Q.W. and O.S. are supported by NSF Grants 1609495 and 1565734.

Appendix A Solving Ricatti’s Differential Equation

The solution technique to Ricatti’s differential equation (see Tenenbaum and Pollard, 1985, p. 97) is not very well known, so the derivation is stated here for the interested reader. A homogeneous 1st order ODE is called a Ricatti equation if it is quadratic in the unknown, i.e.

[TABLE]

If $q_{0}(x)\equiv 0$, this reduces to Bernoulli’s equation (see Tenenbaum and Pollard, 1985, pp. 95-96). In general, one can transform Ricatti’s equation to an equivalent 2nd order linear differential equation. In this appendix, we detail this transformation and explicitly solve for the case of constant coefficients $q_{0}$, $q_{1}$ and $q_{2}$.

First, define a new variable $v$ by

[TABLE]

so that equation (53) becomes

[TABLE]

where $R=q_{1}+q_{2}^{\prime}/q_{2}$ and $S=q_{0}q_{2}$. Then, introduce another variable, $u$, related to $v$ via a Cole-Hopf transform:

[TABLE]

and now the original equation, in terms of $u$, becomes

[TABLE]

For constant coefficients $R$ and $S$, equation (57) can be solved exactly. Using terminology from the present paper, the constant coefficient Ricatti equation is

[TABLE]

where $\psi_{C}=\psi_{C}(t)$ is the unknown variable and $r$, $c$, $\rho_{m}$, $\alpha$, $\beta$, $\psi_{A}$, and $\psi_{B}$ are fixed parameters. Relating this to equation (53), let $y=\psi_{C}$, $q_{0}=cr\psi_{A}\psi_{B}$, $q_{1}=-\alpha\beta/\rho_{m}$, and $q_{2}=\beta/\rho_{m}$. Therefore, in the final equation, $R=q_{1}+q_{2}^{\prime}/q_{2}=-\alpha\beta/\rho_{m}$ and $S=q_{0}q_{2}=cr\beta\psi_{A}\psi_{B}/\rho_{m}$.

The solution to equation (57) when $R$ and $S$ are constant depends on the eigenvalues of its corresponding characteristic equation. More specifically, it depends on the sign of the determinant of the root of the characteristic equation $R^{2}-4S=\chi/\rho_{m}^{2}$, where $\chi$, defined in section III.1, was determined to be positive for physically realistic parameter values.

The solution to this case – where the characteristic equation to equation (57) has two real, distinct roots – can be found in any introductory book on ODEs (see Boyce and DiPrima, 2008, pp. 137 - 143), but for completeness the solution is detailed here with homogeneous initial conditions $y(0)=0$:

[TABLE]

where

[TABLE]

and the antiderivative of this solution is given by

[TABLE]

where $C$ is an arbitrary constant of integration.