Modeling Multiphase Flow Within and Around Deformable Porous Materials: A Darcy‐Brinkman‐Biot Approach

We present a new computational fluid dynamics approach for simulating two‐phase flow in hybrid systems containing solid‐free regions and deformable porous matrices. Our approach is based on the derivation of a unique set of volume‐averaged partial differential equations that asymptotically approach the Navier‐Stokes Volume‐of‐Fluid equations in solid‐free regions and multiphase Biot Theory in porous regions. The resulting equations extend our recently developed Darcy‐Brinkman‐Biot framework to multiphase flow. Through careful consideration of interfacial dynamics (relative permeability and capillary effects) and extensive benchmarking, we show that the resulting model accurately captures the strong two‐way coupling that is often exhibited between multiple fluids and deformable porous media. Thus, it can be used to represent flow‐induced material deformation (swelling, compression) and failure (cracking, fracturing). The model's open‐source numerical implementation, hybridBiotInterFoam, effectively marks the extension of computational fluid mechanics into modeling multiscale multiphase flow in deformable porous systems. The versatility of the solver is illustrated through applications related to material failure in poroelastic coastal barriers and surface deformation due to fluid injection in poro‐visco‐plastic systems.

larger flow paths with length scales on the order of ∼10 −2 m) (Bottero et al., 2010). Similarly, the propagation of flow-driven fractures in porous materials and the propagation of waves in coastal barriers involve feedbacks between flow and mechanics in porous and solid-free domains. In the present study, we develop a framework capable of representing multiphase flow through and around deformable porous systems, as required to simulate many of the aforementioned phenomena.
The starting point for our study is based on the present ample understanding of multiphase flow dynamics within and around static porous materials, from viscous and capillary fingering (Ferer et al., 2004;Lenormand & Zarcone, 1989;Lenormand et al., 1988) to temperature and surface tension driven flows (Shih & Megaridis, 1996), all the way to turbulent multiphase flows (Colombo & Fairweather, 2015;Soulaine & Quintard, 2014). This knowledge, in conjunction with numerical techniques such as the Lattice Boltzmann Method, the Finite Volume Method, Homogenization Theory, and Averaging Theory, forms the basis of fast and accurate models that are routinely applied to help design and improve hydrocarbon production (Burrus et al., 1991;Mehmani & Tchelepi, 2019), CO 2 sequestration (Hassan & Jiang, 2012), and even nuclear reactors (Tentner et al., 2008). However, the study of multiphase flow across different scales remains limited as shown by the absence of well-established approaches to describe how bubbles or waves propagate into an unsaturated porous medium or how a multiphase fluid mixture is pushed out of a porous medium into open space. Improved understanding of such processes would have a direct and immediate impact in the design of batteries, coastal barriers, natural gas extraction from shales, biochemical gas production, and many other areas.
A similar situation pertains with regard to the coupling between fluid flow and solid mechanics. Theoretical and numerical approaches based on Biot's Theory of poroelasticity (Biot, 1941), Terzaghi's effective stress principle (Terzaghi, 1943), and Mixture Theory (Siddique et al., 2017) have been successful at modeling systems with flow in deformable porous media including arteries, biofilms, boreholes, hydrocarbon reservoirs, seismic systems, membranes, soils, swelling clays, and fractures (Auton & MacMinn, 2017;Barry et al., 1997;Jha & Juanes, 2014;Lo et al., 2002Lo et al., , 2005MacMinn et al., 2016;Mathias et al., 2017;Santillán et al., 2017). However, as mentioned above, we still have very little understanding of how flow-induced deformation of these solid materials affects the macroscopic flow around them (and thus their boundary conditions) or how fluid-fluid interfaces behave when pushed against a soft porous medium and vice versa. These type of systems are particularly challenging to model because of the multiscale nature of the governing physics, where the curvature of a particular interface and the magnitude of the capillary and/or viscous forces can differ by several orders of magnitude depending on the presence or absence of a porous solid.
Three major approaches have been proposed to resolve the challenge posed by fluid flow in porous media containing both solid-free regions and porous domains (hereafter referred to as multiscale systems). The most straightforward of these involves performing direct numerical simulations (DNS) throughout the entire multiscale domain, both within and outside the porous medium (Breugem & Boersma, 2005;Hahn et al., 2002;Krafczyk et al., 2015). Although rigorous, this technique is impractical in situations with a large difference in length scales between the largest and smallest pores, where it requires exceedingly fine grids and tremendous computational resources.
To save time and resources, other studies have relied on hybrid DNS-Darcy approaches, where fluid and solid mechanics within a porous medium are modeled as averaged quantities through Darcy's law, pore-network models, or Biot's theory of poroelasticity (Ehrhardt, 2010;Weishaupt et al., 2019). One such approach relies on the use of the Beavers-Joseph (BJ) boundary condition to couple fluid flow in solid-free domains (simulated using the Navier-Stokes Equations) and in porous domains (simulated using Darcy's law) for single phase flow and static porous media (Beavers & Joseph, 1967;Fetzer et al., 2016). Recent studies have extended this BJ approach to allow multiphase flow in the solid-free domain (Baber et al., 2016) or to include the effects of poroelasticity within the porous medium (Lacis et al., 2017;Zampogna et al., 2019). However, to the best of our knowledge, no BJ-based technique has yet been developed to couple solid mechanics with multiphase flow simultaneously within the solid-free and porous domains.
The Darcy-Brinkman (DB) approach presents a well-known alternative to the BJ interface matching technique. The crux of the DB approach is the use of a spatially dependent penalization term within the Navier-Stokes fluid momentum equation. This term effectively creates an equation that approximates Navier-Stokes within solid-free domains and Darcy's law within porous domains. Although initially implemented as an empirical approach (Brinkman, 1947), this technique has since been formalized and rigorously derived from first principles through volume averaging theory (Soulaine et al., 2016;Whitaker, 2013). The resulting so-called "microcontinuum" approach has been extensively used to solve single phase flow through static multiscale porous media, such as flow in biofilms (Kapellos et al., 2007) and in rocks containing unresolved porosity (Guo et al., 2018;Kang et al., 2019;Singh, 2019). The approach has proved highly flexible as illustrated by its uses to represent embedded solid boundaries in low-permeability media (Khadra et al., 2000) and the evolution of solid grain morphologies caused by mineral dissolution (Soulaine et al., 2017(Soulaine et al., , 2019. Recently, a study by Carrillo and Bourg (2019) introduced a Darcy-Brinkman-Biot (DBB) approach capable of accurately representing single phase flow in multiscale deformable media including elastic porous membranes and plastic swelling clays. Simultaneously, studies by Soulaine et al. (2019) and  extensively benchmarked and released an open-source extension of the microcontinuum framework for multiphase flow in static multiscale porous media. This allowed accurate modeling of complex systems such as multiphase flow in a fractured porous medium, methane extraction from tight porous media, and wave absorption in coastal barriers. In the present paper, we build upon these previous studies to create the first model representing coupled fluid and solid mechanics during multiphase flow in multiscale deformable porous media: the multiphase DBB model (Figure 1). This study is organized as follows. Section 2 introduces the concept of volume averaging and describes the derivation of the governing equations for coupled fluid and solid mechanics. Section 3 explains the numerical implementation and algorithm development for the coupled mass and momentum equations and introduces the resulting open-source solver "hybridBiotInterFoam." Section 4 presents five test cases that verify the implementation of different coupling terms within the model, with an emphasis on fracturing mechanics. Section 5 then presents two alternative applications that illustrate the versatility of the model: wave absorption in poroelastic coastal barriers and surface deformation due to fluid injection in poro-visco-plastic geologic formations. Lastly, Section 6 concludes with a summary of the paper and a discussion on future work.

Volume Averaging
In this section, we introduce the concept of volume averaging. This technique forms the basis of the microcontinuum equations, as it allows the classical mass and momentum conservation equations to account for the coexistence of solid (s), wetting fluid (w), and nonwetting fluid (n) within a given control volume. It is well suited for use in conjunction with the Finite Volume Method (FVM) (Patankar, 1980), as the numerical grid elements used in the FVM provide an intuitive and straightforward numerical interpretation of what we will define as the averaging volume (V). In keeping with standard volume averaging theory, we start by defining the volume averaging operator where β i is a function defined in each phase's respective volume    , , i V i w n s , and V = sum(V i ). We also define the phase averaging operator CARRILLO AND BOURG 10.1029/2020WR028734 3 of 27 The model considers wetting properties, interface mechanics, and irreducible saturations when averaging over a REV. Note that the stated relation between the averaging volume's length scale L V and the porous length scale L P is required for the creation of a REV, and thus, for the application of this model.
where μk −1 is the drag coefficient (a function of the fluid viscosities and permeability k), s U is the averaged U is a solid-fluid momentum exchange term that accounts for a moving porous medium in an Eulerian frame of reference, and F c,i represents the forces emanating from fluid-fluid and fluid-solid capillary interactions. As shown in , where p c is the average capillary pressure within a given averaging volume, γ is the fluid-fluid interfacial tension, M i = k 0 k i,r /μ i is the mobility of each fluid (a function of absolute permeability k 0 and relative permeability k i,r ), and M = M w + M n is the single-field fluid mobility. Lastly, n w,n is the unit normal direction of the fluid-fluid interface as calculated by the Continuum Surface Force (CSF) formulation (Brackbill et al., 1992). The equations presented above tend toward the standard Navier-Stokes Volume-of-Fluid approach in solid-free regions (where the drag term becomes negligible) and toward the multiphase Darcy equations in porous regions. The latter can be explained by the fact that the viscous stress tensor   S becomes negligible under the scale-separation assumption, inertial terms become negligible under the assumption of low Reynold's number flow in the porous medium, and the F c , terms are set to fit standard multiphase Darcy's law Whitaker, 1986): in solid free regions Eqn.14 in porous regions

Derivation of the Solid Mechanics Equations
We proceed with the derivation of the microcontinuum solid mechanics equations by starting from the equations presented in Carrillo and Bourg (2019) for solid mass and momentum conservation in systems with a single incompressible solid phase.
where σ is the volume-averaged solid elastic (or plastic) stress tensor and    We now seek closure of these two coupling terms. By conservation of momentum, we know that any drag-induced momentum lost by the fluid must be gained by the solid. Therefore, we can use the drag term used in Equation 14 to obtain (Carrillo & Bourg, 2019)   Closure of the capillarity-induced interaction term B cap is obtained by combining the solid and fluid momentum equations within the porous medium at low Reynold numbers and low permeability, which yields In multiphase porous systems with incompressible grains and no swelling pressure (i.e.,    p s τ ), Biot s s f f and p c is the capillary pressure (Jha & Juanes, 2014;Kim et al., 2013). This expression is satisfied by the previous equation in the absence of capillary forces, where F c,1 , F c,2 , B cap , and p c equal zero (Carrillo & Bourg, 2019). In the presence of capillary forces, however, it imposes the following equality Given that F c,1 = −p c ∇α w in the porous domains , the previous equation can be rearranged to obtain Equation 24 gives closure to the last coupling parameter and marks the end of this derivation. The result is a solid conservation equation that tends toward Biot Theory in porous regions and toward an infinitely deformable solid with no momentum sources in solid-free regions.

Interfacial Conditions Between Solid-Free Regions and Porous Regions
One of the most important features within the framework presented above is the existence of an interface between solid-free and porous domains. Although the creation of a rigorous unaveraged description of this interface is still an open question, we approximate a solution to it by guaranteeing its necessary components within our fluid and solid averaged equations.
An accurate description of fluid behavior at the interface requires three components: (1) mass conservation across the interface, (2) continuity of stresses across the interface, and (3) an interfacial wettability condition. Components 1 and 2 are intrinsically fulfilled by our solver due to its single-field formulation for velocity and pressure within the fluid conservation equations (Equations 8 and 14). As shown in Neale and Nader (1974) and Carrillo and Bourg (2019), these two components are necessary and sufficient to model single-phase flow within a multiscale system. Furthermore, these conditions have also been used for closure when modeling multiphase flow in moving porous media Lacis et al., 2017;Zampogna et al., 2019). The required wettability condition at the porous interface (Component 3) is included in our model through the implementation of a penalized contact angle condition (Equation 33) following the steps outlined in Horgue et al. (2014) and .
The complementary solid conditions at the porous interface are very similar: (1) solid mass conservation across the interface, (2) continuity of fluid-induced stresses across the interface, and (3) a discontinuity of solid stresses at the interface. Just as before, the first two conditions are intrinsically fulfilled through the use of a single set of mass and momentum conservation equations across both domains and have also been used as closure conditions in previous studies (Lacis et al., 2017;Zampogna et al., 2019). The third condition is enforced by the use of volume-averaged solid rheology models that tend toward infinitely deformable materials in solid-free regions, as shown in Carrillo and Bourg (2019). When volume-averaged, the behavior of the solid's stress tensor is domain dependent (i.e., solid fraction dependent). Thus, in solid regions, the elasticity and viscosity of the porous medium is determined by standard averaged rheological properties (the elastic and viscoplastic moduli). Contrastingly, in solid-free regions, the solid fraction tends to zero and, as such, said properties do as well. The result is a stress-free "ghost" solid that does not apply resistance to the porous region, creating the required stress discontinuity at the porous interface.
Although necessary, these conditions represent but an approximation to the complete description of fluid and solid mechanics at the porous interface. However, to the best of our knowledge, there does not exist an alternative set of interfacial conditions that can or have been used to model multiphase flow in multiscale porous media.

Model Summary
The final set of equations in our proposed multiphase DBB framework now follows. The combination of these solid and fluid conservation equations leads to a model that tends toward multiphase Navier-Stokes in solid-free regions and toward Biot Theory in porous regions, as described in Figure 1.
All that is left is stating the closed-form expressions of the multiscale parameters μk −1 , F c,i , and U r , which are defined differently in each region. A full derivation and discussion of these parameters can be found in .
where C α is an interface compression parameter (traditionally set to values between 1 and 4 in the Volume-of-Fluid method), k 0 is the absolute permeability, k r,i and M i = k 0 k i,r /μ i are the relative permeability and mobility of each fluid, and M = M w + M n . Lastly, θ is the imposed contact angle at the porous wall, and n wall and t wall are the normal and tangential directions relative to said wall, respectively.
Finally, closure of the system of equations requires appropriate constitutive models describing the averaged behavior of the different phases within the porous regions. For the purpose of validating our multiphase DBB approach, in the present study, we use the following well-established constitutive models: absolute permeability is modeled as isotropic and porosity-dependent through the well-known Kozeny-Carman (1 ) f f k k ; relative permeabilities and average capillary pressures within the porous domains are represented using the Van Genutchen (van Genuchten, 1980) and Brooks-Corey (Brooks & Corey, 1964) models (Appendix A); plasticity is described through the Herschel-Bulkley model, were the solid viscously deforms only after local stresses become higher than the material yield stress (Appendix B1); the solid's yield stress and plastic viscosity are modeled as solid fraction-dependent based on the Quemada fractal model (Quemada, 1977;Spearman, 2017) (Appendix B2); finally, elastic solids are modeled as averaged linear-elastic materials, such that their averaged elastic coefficients scale linearly with respect to the solid fraction (Appendix B3). The last three choices imply that solid rheological properties are modeled as isotropic and independent of saturation, a significant simplification that is sufficient for the purpose of testing and validating the present framework. For the reader's convenience, a full implementation of this framework and its related models are included in the accompanying simulation "toolbox." If necessary, more complex constitutive models, such as the saturation-depended solid rheology models presented in Wan et al. (2014), Oldecop and Alonso (2003), Buscarnera and Einav (2012), and Di Donato et al. (2003) can be readily implemented into our code by virtue of its open-source implementation.

Numerical Platform
The implementation of the multiphase DBB model was done in OpenFOAM®, a free, open-source, parallelizable, and widely used computational fluid mechanics platform. This C++ code uses the Finite Volume Method to discretize and solve partial differential equations in complex 3-D grids. Its object-oriented structure and multitude of supporting libraries allows the user to easily customize each simulation's setup with different numerical discretization schemes, time-stepping procedures, matrix-solution algorithms, and supporting physical models. The implementation described below represents the natural extension of the multiphase microcontinuum toolkit "hybridInterFoam"  to systems with deformable solids. In particular, its solution algorithm stems directly from that used by "hybridInterFoam" and its precursor "interFoam."

Solution Algorithm
The solution of the governing equations is done in a sequential manner, starting with the fluid mechanics equations and following with the solid mechanics equations for every time step. Of particular importance is the handling and modification of the velocity-pressure coupling required for modeling incompressible fluids in conjunction with a moving solid matrix. For this step, we based our solution algorithm on the Pressure Implicit Splitting-Operator (Issa, 1986). First, we explicitly solve the fluid saturation equation (Equation 9) for  1 t w through the Multidimensional Universal Limiter of Explicit Solution algorithm (Márquez & Fich, 2013). This allows for stable numerical advection of the saturation field by the application of Flux Corrected Transport Theory (Rudman, 1997). Then, we update the boundary values of U f and U r in addition to the cell-centered values of the permeability k t+1 , density  1 where   These equations can be recast into a single coupled equation which is then used to implicitly solve for pressure. This step can be done through several generalized matrix solvers that are standard in OpenFOAM®.
After solving for pressure p*, velocity can be recalculated from Equation 36. This semi-implicit pressure-velocity correction step is repeated until the desired convergence is reached. It has been shown that at least two pressure-velocity correction loops are required to ensure mass conservation (Issa, 1986).
At this point, 1 t f U and p t+1 are set and used as input values for updating the drag and pressure source terms present in the solid mechanics momentum equation (Equation 29). Then, in the case of visco-poro-plasticity, said equation is discretized in a similar way as the fluid momentum equation (Equation 27) and used to implicitly solve for 1 t s U . In the case of poroelasticity, the solid mechanics equation is solved through the algorithm presented in Jasak and Weller (2000). Here, the solid's elastic equation (Equation 29) is discretized and segregated into implicit and explicit components, after which it is iteratively solved until convergence is reached. This segregated method not only guarantees fast convergence but also memory efficiency. Finally, the updated solid velocity is used to "advect" the solid fraction field ϕ s by solving the mass conservation equation (Equation 28). At this point, the algorithm advances in time according to the imposed Courant-Friedrichs-Lewy number. A flow-chart of the complete algorithm can be found in the Supporting Information (SI), and further discussion regarding the discretization techniques and matrix-solution procedures can be found in ), Jasak (1996, and Jasak and Weller (2000).

Open-Source Implementation
The complete set of governing equations and solution algorithms, along with the necessary rheology, relative permeability, and capillary pressure models (Appendices A and B) were implemented into a single solver: "hybridBiotInterFoam." This solver, along with its representative tutorial cases, automated compilation and running procedures, and all the simulated cases presented in this paper were incorporated into an open-source CFD package of the same name. OpenFOAM® and our code are free to use under the GNU general public license and can be found at https://openfoam.org/ and https://github.com/Franjcf , respectively.

Model Validation
Most of the underlying components of the approach described above have been previously tested and verified. Carrillo and Bourg (2019) validated the momentum exchange terms as an effective coupling mechanism between a single fluid phase and a deformable plastic or elastic porous medium. The effects of confining and swelling pressures on porous media were also examined in said study. Then,  extensively validated the extension of the Darcy-Brinkman equation into multiphase flow within and around static porous media by comparison with reference test cases in a wide range of flow, permeability, capillarity, and wettability conditions. Therefore, the only thing left to validate is the ability of the multiphase DBB model to accurately predict the behavior of multiscale systems that exhibit coupling effects between multiple fluids and a deformable porous matrix.
To that point, we begin with two validation cases relating to multiphase poroelasticity and the coupling between solid deformation and fluid pressure. Then, we proceed with two poro-visco-plastic cases that validate this framework for multiscale plastic systems. Finally, we conclude with two additional cases that verify the implementation of the capillary force interaction terms. All of these can be found in the accompanying CFD simulation package. We note that all validation cases presented below were tested for grid-independence. A sample of these tests pertaining to Figures 4-7 can be found in the SI.

Terzaghi Consolidation Problem
The Terzaghi uniaxial compaction test has been extensively used as a benchmark for the validation of numerical codes relating to poroelasticity (Terzaghi et al., 1996). Its main utility is to test the accuracy of the solid-fluid couplings that relate fluid pressure to solid deformation and vice versa. The problem consists of a constrained saturated elastic porous medium that is abruptly compressed from its upper boundary by a constant uniaxial load (Figure 2). This creates a sudden increase in pore pressure, which is then dissipated by flow through the upper boundary (all other boundaries have impermeable boundary conditions). In the case of a one-dimensional porous medium, the resulting temporal and spatial evolution in fluid pressure can be described by the following simplified analytical solution (Verruijt, 2013).
where c v = (k 0 E (ν − 1))/(ρ f g (2ν 2 + ν − 1)) is the consolidation coefficient, k 0 is permeability, E is Young's modulus, ν is Poisson's ratio, h is the column height, and z is the vertical coordinate. Our equivalent numerical setup is shown in Figure 2. for all tested conditions. Further verification of these terms for an oscillating linear-elastic solid with pressure boundary conditions (as opposed to stress boundary conditions) can be found in the SI.

Capillary Pressure Effects in a Poroelastic Column
Having verified the two-way coupling between solid deformation and fluid pressure, we now verify the implementation of the capillary pressure terms within the solid mechanics equation. To do so, we simulate a poroelastic column (1 m tall, 1,500 Cells, ϕ f = 0.5) bounded by two nonwetting fluid reservoirs at its upper and lower boundaries. The column is initialized with a linear saturation profile spanning from α w = 0 to 1 (see Figure 3). Fluid saturation is kept fixed by not solving Equation 26, and the mobilities of both fluids are set to very high values (M i = 1 × 10 10 m 3 /kg.s) to minimize drag-related effects. Under these conditions, the solid's effective vertical stress is exclusively controlled by capillary effects and is described by the following analytical solution: We used the Van Genutchen capillary pressure model with m = 0.6 or 0.8 and p c,0 = 50-2,000 Pa to calculate the solutions to said problem. The resulting agreement between the numerical and analytical solutions, shown in Figure 3, confirms the accuracy of the fluid-solid capillary pressure coupling implemented in our model. Furthermore, the transitional behavior of the effective stress at the macroscopic solid-fluid interface confirms the applicability of the interfacial condition described in Section 2.4: as expected, solid stresses are dictated by standard elasticity theory in the porous region and become negligible in solid-free regions.
Given that the fluid-solid couplings in a poroelastic solid are now verified, we proceed to verify said terms for poro-visco-plastic materials.

Fluid Invasion and Fracturing in a Hele-Shaw Cell
The third verification case (and the first poro-visco-plastic case) consists in the qualitative replication of a set of fracturing experiments that examined the injection of aqueous glycerin into dry sand within a 30 by 30 by 2.5 cm Hele-Shaw cell (Huang et al., 2012a(Huang et al., , 2012b   in that the characteristic length scale of fractures in this system (∼cm) is orders of magnitude larger than that of pores within the porous matrix (∼100 μm). They are also multiphysics, as they clearly exemplify the drag-controlled transition from Darcy flow within the porous medium to Stokes flow in the open fractures and the coupling between the hydrodynamics of fluid flow and the mechanics of fracture propagation (Figure 4).
The experimental setup involved the injection of aqueous glycerin at various flow rates, q, between 5 and 50 ml/min while also varying the fluid's viscosity, μ gly , between 5 and 176 cp for different experiments. Our numerical simulations were parameterized using measured values of the glycerin-air surface tension (γ = 0.063 kg/s 2 ), the density of pure glycerin (ρ gly = 1,250 kg/m 3 ), the density of air (ρ air = 1 kg/m 3 ), the viscosity of air (μ air = 0.017 cp), and the average radius and density of sand grains (100 μm and 2,650 kg/m 3 , respectively). To mimic the sand's experimental configuration and permeability, the simulated solid fraction field was set to a random initial normal distribution such that ϕ s = 0.64 ± 0.05 and the permeability was modeled as a function of the solid fraction through the Kozeny-Carman relation with k 0 = 6.7 × 10 −12 m 2 . Relative permeabilities were calculated through the Van Genutchen model with the Van Genuchten coefficient m set to 0.99 (see Appendix A), while capillary pressures were deemed negligible (as 2γr −1 ≪ μk −1 U f L).  with yield stress τ 0 = 16.02 m 2 /s 2 (Quemada, 1977). Plasticity was used as the preferred mode of solid rheology due to its ability to account for the compressive and irreversible effects caused by fracturing within these experiments (Ahmed et al., 2007;van Dam et al., 2002).
Numerically speaking, the simulations were carried out in a 30 by 30 cm 2-D grid (500 by 500 cells) with constant velocity and zero-gradient pressure boundary conditions at the inlet, zero-gradient velocity and zero pressure boundary conditions at the boundary walls, and a solid velocity tangential slip condition at all boundaries (i.e., the solid cannot flow across the boundaries, but the fluids can). Each simulation was run in parallel for approximately 2.5 hours on a single 28-core node or until the injected glycerin reached the outer boundary. Lastly, to enable a closer comparison between our 2-D simulation and the 3-D experiment we added an additional drag term to the fluid momentum equation equal to 12μa −2 U f , which accounts for viscous dissipation through friction with the walls in a Hele-Shaw cell with aperture a (Ferrari et al., 2015).
As shown in Figure 4, a dramatic transition in the mode of fluid invasion is observed with increasing fluid injection velocity and viscosity. At low flow rates and low viscosity (q = 5 ml/min, μ = 5 cp), there is no discernible solid deformation and the main mode of fluid flow is through uniform invasion of the porous medium (Figure 4a). At intermediate flow rates and low viscosity (q = 25 ml/min to 30 ml/min, μ = 5 cp), we still observe a uniform invasion front, but small fractures begin to appear (Figures 4b and 4c). At high viscosity (μ = 176 cp), we see clear fracturing patterns preceded by a nonuniform fluid invasion front (Figures 4h and 4i). Figure 4 shows that our simulation predictions are qualitatively consistent with the experiments presented in Huang et al. (2012a) with regard to both the stability of the capillary displacement front and the observed fracturing transition behavior. As suggested above, accurate prediction of this transition requires not only proper handling of fluid-fluid interactions (surface tension and relative permeability effects), but also accurate descriptions of their relationship with solid mechanics (drag) and the proper implementation of a solid rheological model that can replicate irreversible and unstable fracturing processes. We note that in our simulations, fracture initialization and propagation are predicted based on continuum-scale equations for the rheology and mechanics of the bulk porous solid, with no specific treatment of grain-scale mechanics. Grid-level instabilities are brought about by the normally distributed porosity and permeability fields, as shown in Appendix C. The microstructural differences between the experiments and our simulations (most clear in Figures 4c, 4f, 4h, and 4k) likely arise at least in part from the fact that the solid is modeled as a continuum rather than a granular material.
This section demonstrates that the multiphase DBB model can be used to replicate and predict the main mode of fluid flow and solid deformation within fracturing systems. A comprehensive study of the controlling parameters for multiphase fracturing in the presence of both viscous and capillary stresses will be the focus of an adjacent study.

Modeling Fracturing Wellbore Pressure
Having shown that our model can qualitatively predict fracturing behavior, we now aim to determine whether it can do so in a quantitative matter. As depicted in Figure 5, fluid-induced fracturing of low-permeability rocks proceeds through the following well-established series of events: First, fluid pressure increases linearly as fracturing fluid is injected into the wellbore. Second, as wellbore pressure increases and approaches the leak-off pressure, a small amount of pressure is propagated by fluid leakage into the rock. Third, fluid pressure continues to increase until it reaches the breakdown pressure, at which point it is high enough to fracture the rock. Fourth, a fracture is initiated and propagates; the wellbore pressure slowly decreases. Fifth, injection stops, fracture propagation stops, and wellbore pressure rapidly dissipates (Abass et al., 2007;Ahmed et al., 2007;Huang et al., 2012a;Papanastasiou, 2000;Santillán et al., 2017). In this section, we aim to numerically replicate the time-dependent fracturing wellbore pressure during fracture propagation (i.e., the fourth stage outlined above) as described by an analytical solution presented in Barros-Galvis et al. (2017).
where t is the time elapsed since fracture initialization, q is the fluid injection rate, p well is the wellbore pressure, p 0 is the minimum pressure required for starting a fracture (a function of the solid's yield stress τ 0 ), h is the formation thickness, and r well is the wellbore radius. The remaining variables follow the same definitions described earlier. The general numerical setup is almost identical to the one presented in the previous section. The key difference is that we now inject aqueous glycerin into a strongly non wetting (and thus almost impermeable) porous material. This is done to ensure an accurate replication of the analytical solution and its related assumptions, where fracturing is the main mode of fluid flow and there is virtually no fluid invasion into the porous matrix. The exact simulation parameters are q = 46-110 ml/min, τ 0 = 0.2 or 2 m 2 /s 2 , k 0 = 6.7 × 10 −11 or 6.7 × 10 −12 m 2 , μ gly = 5 cp, and m = 0.05. Note that low values of m indicate that the porous formation is strongly nonwetting to the injected fluid (see Figure S4 in the SI for the resulting relative permeability curve). All other parameters are as in the previous section.  Figure 4, and p max is the maximum analytically predicted pressure in each simulation.
Lastly, as hinted at before, a notable characteristic of our model is that different normally-distributed solid fraction field initializations give different fracturing results (Appendix C). For this reason, we performed four simulations for each parameter set. In Figure 6, we present the average predicted wellbore pressure evolution with errors bar representing the 95% confidence interval. Figure 6 shows that our model can accurately and reliably predict the pressure and deformation behavior of a variety of fracturing systems, as all curves exhibit excellent agreement with their respective analytical solution. Note that the length of each curve relates inversely to the injection speed. This is because fractures at higher injection rates consistently reach the system's boundary faster than their counterparts, at which point there is a sharp decrease in pressure and the analytical solution no longer applies. Therefore, each curve's cutoff point represents the time at which the fracture effectively becomes an open channel between the wellbore and the outer boundary, normalized to the average value of that time for the slowest-moving fracture (i.e., t = t max ).
The successful replication of the analytical pressure profiles in this section verifies the model components pertaining to the pressure-velocity-deformation coupling and the two-way momentum transfer between the fluid and solid phases (drag). Therefore, the only model component left to verify is the implementation of the capillary force terms during fracturing of a plastic solid.
CARRILLO AND BOURG 10.1029/2020WR028734 16 of 27 Figure 7. Effect of capillary entry pressure on fracturing wellbore pressure. (a and b) Wellbore pressure as a function of time and entry pressure for low and high permeability systems, respectively. In (b), curves at increasingly high pressures were cut off for illustrative purposes and the solid line represents a fitted reference logarithmic pressure descent curve. (c-h) Time evolution of fractured system with a 1 kPa capillary entry pressure and high permeability. (c) Initial fluid invasion (t/t max < 0): At early times, the wellbore pressure rises rapidly and becomes larger than the entry capillary pressure. The fluid invades the porous formation symmetrically. (d) Fracture initiation (t/t max = 0): The wellbore pressure continues to rise until it is larger than the breakdown pressure, at which point small fractures start to form. Fluid invasion continues. (e and f) Fracture propagation (t/t max > 0|p well > p c,0 ): The wellbore pressure drops as fractures propagate. Fluid invasion continues asymmetrically around said fractures. (g) Fluid invasion stops (t/ t max > 0|p well ∼p c,0 ): As the wellbore pressure keeps dropping, the entry capillary pressure condition at the porous interface ensures that that wellbore pressure never goes below p c,0 , at which point fluid invasion stops. (h) Fracture reaches the simulation boundary (t/t max = 1). The color convention in Figures (c-h) is the same as in Figure 4.

Capillary Effects on Fracturing Wellbore Pressure
Our fifth verification systematically varies the capillary entry pressure within nonwetting fracturing systems to quantify its effects on wellbore pressure. For this, we consider two different complementary cases: one where capillary forces are comparable to their viscous counterparts, and another where they are significantly larger than them. All parameters are the same as in the previous experiments (Section 4.4) unless otherwise specified.
The first set of experiments expands the previous analysis (Section 4.4) into strongly nonwetting systems with the addition of a constant capillary pressure jump at the fracture interface imposed by a flat capillary pressure curve (p c = p c,0 = 0-2 kPa, τ 0 = 2 m 2 /s 2 , k 0 = 6.7 × 10 −12 m 2 , m = 0.05, and q = 78 ml/min).  Figure 7a, where we present the updated analytical results in conjunction with our equivalent numerical results, demonstrating excellent agreement between them. Note that the predicted linear relationship between wellbore pressure and capillary entry pressure is not explicitly imposed in the numerical model. On the contrary, it arises naturally from the balance of viscous, capillary, and structural forces in Equations 25-29.
The second set of experiments modifies the previous experiments by making the porous medium significantly more permeable, while still maintaining a constant capillary pressure jump at the fracture interface (p c = p c,0 = 1-3 kPa, τ 0 = 0.2 m 2 /s 2 , k 0 = 6.7 × 10 −11 m 2 , m = 0.99, and q = 78 ml/min). This results in a set of cases where the wellbore pressure is increasingly controlled by the capillary pressure drop rather than by the viscous pressure drop across the fracture and porous formation.  (Figures 7b-7h). The resulting pressure drop cannot be modeled by the previously presented analytical solution (as it violates the no leak-off assumption), but still follows a logarithm-type curve that is characteristic of flow in fracturing systems. With increasing fracture propagation, the viscous pressure drop decreases until the wellbore pressure equals the entry pressure, which is, by definition, the minimum pressure drop required for fluid flow in highly permeable nonwetting systems. Finally, we note that in cases where capillary entry pressure is high relative to the pressure required to fracture the solid (i.e., at (p c,0 > 2.250 Pa in the conditions simulated in Figure 7b), fracturing begins before the wellbore pressure can exceed p c,0 . This prevents essentially all flow into the porous formation, and the wellbore pressure is immediately stabilized at ∼p c,0 . For all cases, fractures continue to propagate until they reach the system boundary, at which point the pressure drops rapidly as noted in Section 4.4.
In this section, we reduced the inherent complexity of the model's capillary force terms F c,i (Equations 31 and 32) into a simple set of intuitive verifications. The quantitative agreement between these two analytical cases and their corresponding numerical simulations validate the implementation of the impact of capillary pressure effects on the mechanics of a ductile porous solid within our model.

Illustrative Applications
Having verified and tested the model, we now proceed with two illustrations that demonstrate how hybrid-BiotInterFoam enables the simulation of relatively complex coupled multiphase multiscale systems. The following cases serve as illustrative examples of our model's features and capabilities as well as tutorial cases within the accompanying toolbox.

Elastic Failure in Coastal Barriers
Coastal barriers are ubiquitous features in coastal infrastructure development. When designed appropriately, these structures can be very effective in regulating water levels and protecting against inclement weather (Morton, 2002). However, accurate prediction of the coupled fluid-solid mechanics of these structures (which can lead to barrier failure) is inherently challenging as it requires modeling largescale features (waves) while also considering small-scale viscous and capillary interactions within the barrier.
The following case represents the continuation of the three-dimensional coastal barrier illustration presented in  with the addition of linear-elastic poromechanics. As such, the simulation was created by initializing a heterogeneous porosity field (with k 0 = 2 × 10 −8 m 2 and ϕ f = 0.5) in the shape of a barrier within a 8.3 by 2.7 by 0.25 m rectangular grid with over 43 million cells (1600 by 540 by 50 cells). The relevant solid mechanics parameters were E = 5 MPa, ν = 0.45, and ρ s = 2350 kg/m 3 . Relative permeabilities and capillary pressures were evaluated through the Van Genuchten model with m = 0.8 and p c,0 = 1 kPa. Before the start of the simulation, the water level was set to partially cover the barrier and then allowed to equilibrate. A single wave was then initialized at t = 0. This results in a simulation that exhibits a clear wave absorption cycle that gradually dissipates in time, as seen in Figure 8. Detailed discussion on the fluid mechanics of this problem can be found in . In total, this 3-D simulation lasted for 15 simulated seconds, which took approximately 30 hours to run on 16 computational nodes with 28-Broadwell Xeon cores each.
Here, however, we are interested in evaluating the barrier's propensity to failure. We do this by applying the Von Mises yield criterion, which is commonly used to predict material failure in elastic systems. It states that if the second invariant of the solid's deviatoric stress (the Von Mises stress) is greater than a critical value (the yield strength) the material will begin to deform nonelastically (Von Mises, 1913). Although we do not specify said critical value within our simulations, we can map the time-evolution of Von Misses stresses within the coastal barrier as a result of a wave absorption cycle (Figure 8). Our results illustrate the potential utility of our simulation framework in predicting the location and time-of-formation of stress induced defects within coastal barrier as a function of wave characteristics, permeability, and barrier geometry.  Note that the largest stresses are seen during the initial wave crash and increase toward the base of the barrier due to gravitational effects.

Flow-Induced Surface Deformation
Surface deformation due to subsurface fluid flow is a common geological phenomenon occurring in strongly coupled systems and has clear implications in studies related to induced seismicity (Shapiro & Dinske, 2009), CO 2 injection in the subsurface (Morris et al., 2011), land subsidence (Booker & Carter, 1986), and the formation of dikes and volcanoes (Abdelmalak et al., 2012;Mathieu et al., 2008). In order to properly model these systems, it is necessary to be able to capture the time-evolution of surface uplift, cracks, and hydraulic fractures, as well as the effects that these features have on the overall flow field. Here, we use the terms hydraulic fracture versus crack to refer to solid failure at versus away from the injected fluid, respectively.
This illustrative case was inspired by the experiments reported by Abdelmalak et al. (2012), where the authors injected a highly viscous fluid into a dry silica powder in a Hele-Shaw cell in order to study the impact of hydraulic fractures on surface deformation, for example, during the creation of volcanic structures. The system also bears some analogy to situations involving the injection of fluids into subsurface reservoirs, for example, during geologic CO 2 sequestration (Rutqvist, 2012). The base case of our simulations consists of an impermeable rectangular container (50 by 30 cm, 500 by 300 cells) that is open to the atmosphere, is partially filled with a dry porous medium (ϕ s = 0.6 ± 0.05, ρ s = 2650 kg/m 3 , k 0 = 5 × 10 −11 m 2 ), and has an injection well at its lower boundary that injects water at q = 6.5 ml/s ( Figure 9). To account for irreversible solid deformation, the porous medium is modeled as a plastic with yield stress τ 0 = 0.22 m 2 /s 2 . The solid is represented as impermeable to the invading fluid through the use of the Van Genuchten model with m = 0.05 and p c = 0. Then, using this base case as a standard, we individually varied each of the main parameters (q, k 0 , τ 0 , m, ϕ s , μ water ) over several simulations in order to model the resulting solid deformation processes: fracturing, cracking, surface uplift, and subsidence ( Figure 9).
The resulting cases demonstrate that cracking (solid failure away from the injected fluid) is strictly dependent on the number and orientation of existing hydraulic fractures, as it only occurs when there is more than one fracture branching off from the main injection point (Figures 9b-9d, 9h, and 9I). This is likely because in cases presenting a single vertical fracture solid displacement is almost exclusively perpendicular to the fracturing direction, leading to virtually no surface deformation or cracking (Figures 9a, 9e, and 9m). Contrastingly, the creation of inclined fractures exerts vertical forces on the solid, resulting in surface uplift and crack formation. The above diagram strongly suggests that deformation is controlled by the balance between viscous and structural forces: larger fractures occur within softer solids with higher momentum transfer, and smaller fractures occur in tougher solids with lower momentum transfer. As stated above, a comprehensive examination of the parameters that control solid fracturing will be the focus of an adjacent paper.
In addition to the surface uplift presented above, subsurface subsidence is observed in the simulated system in conditions where the porous solid is rendered permeable to the invading fluid (i.e., m ≫ 0.05). This phenomenon is not primarily controlled by momentum transfer, but rather by a gravitational effect whereby the displacement of air by water within the porous medium around the advancing hydraulic fracture renders the solid-fluid mixture heavier. Once it is heavy enough to overcome the plastic yield stress, the solid subsides and compresses around the fluid source ( Figure 9m).
With these last two illustrative examples, we have shown that our modeling framework is flexible and readily applicable to a large variety of cases within elastic and plastic systems. We invite the interested reader to tune, adapt, and expand the present illustrative simulations, which are included in the accompanying CFD toolbox.

Conclusions
We derived, implemented, benchmarked, and applied a novel CFD package for simulation of multiphase flow within and around deformable porous media. This microcontinuum modeling framework is based on elementary physics and was rigorously derived through the method of volume averaging and asymptotic matching to the multiphase Volume-of-Fluid equations in solid-free regions and multiphase Biot Theory in porous regions. The result is a single set of partial differential equations that is valid in every simulated grid cell, regardless of content, which obviates the need to define different meshes, domains, or complex interfacial conditions within the simulation. The solver's numeric and algorithmic development were discussed and implemented into hybridBiotInterFoam, an open-source package accessible to any interested party.
Throughout this study and its of predecessors (Carrillo & Bourg, 2019;, we show that the Multiphase DBB model can be readily used to model a large variety of systems, from single-phase flow in static porous media, to elastic systems under compression, to viscosity-or capillarity-dominated fracturing systems, to multiscale wave propagation in poroelastic coastal barriers. We note, however, that the solver presented here cannot be liberally applied to any porous system, as it comes with the following inherent limitations. First, closure of the system of equations requires appropriate constitutive and parametric relations that describe fluid pressure, permeability, capillarity, and rheol-CARRILLO AND BOURG 10.1029/2020WR028734 20 of 27 Figure 9. Study of the impact of subsurface fluid injection on hydraulic fracturing, cracking, and surface deformation. (a-i) Representative cases showing the effects of changing permeability k 0 (purple), solid yield stress τ 0 (green), injection rate q (brown), and injected fluid viscosity μ (red) on surface deformation. The blue and yellow subsections contain the results of increasing or decreasing the controlling parameters, respectively. (j-l) Time evolution of the fracturing base case. (m) Surface subsidence example. The difference between the base case (e) and all other simulations is shown in each case's legend. Dotted white lines represent the surface height of the initial solid fraction configuration. Note that the color scheme in all simulations is the same as in Figure 4. ogy within volume-averaged porous regions. Therefore, the assumptions present in each of these models should be carefully considered. Second, volume averaging imposes important length scale restrictions in order to fulfill the scale-separation hypothesis, where the pore sizes within the averaging volume must be substantially smaller than the chosen REV, and the REV must be substantially smaller than the macroscopic length scale. Third, as currently implemented, the multiphase DBB framework only represents continuum-level elastic or plastic solid mechanics that can be described from an Eulerian frame of reference. As such, it cannot be used to model large elastic deformations or phenomena originating from sub-REV heterogeneities such as fluidization or granular mechanics (Meng et al., 2020), except insofar as they are captured in an averaged manner at the REV scale. Fourth, the use of the CSF as a representation of capillary forces within solid-free regions enforces mass conservation, but it creates a diffuse fluid-fluid interface that may generate spurious and parasitic currents.
Finally, although the modeling framework developed here opens up significant new possibilities in the simulation of coupled fluid-solid mechanics, it also creates a need for the development of constitutive relations describing the coupling between multiphase flow and poromechanics. Of particular importance is the formulation of saturation and deformation-dependent solid rheological models (both plastic and elastic), as well as the rigorous derivation of the interfacial condition between solid-free and deformable porous regions. In this study, we proposed a suitable approximation for said condition based on our single-field formulation, the implementation of a wettability interfacial condition, and the previous work done by Neale and Nader (1974) and Zampogna et al. (2019). However, the accuracy and validity of such an approximation is still an open question, one that is at the frontier of our modeling and characterization capabilities (Qin et al., 2020). The derivation and implementation of said interfacial condition, along with the addition of erosion and chemical reactions into this modeling framework, will be the focus of subsequent papers.

Acknowledgments
This work was predominantly supported by the National Science Foundation, Division of Earth Sciences, Early Career program through Award EAR-1752982. FJC acknowledges additional support from the High Meadows Environmental Institute at Princeton University through the Mary and Randall Hack '69 Research Fund. We do not report any conflicts of interest.