1. Introduction
Fluid-driven or hydraulic fracturing results from the injection of a pressurized fluid in low-permeability solid media. The formation and propagation of the fluid-filled tensile fracture is observed commonly in engineering and natural geophysical processes. For example, the formation of magma-driven dykes is due to density differences that generate pressure large enough to propagate a vertical fracture in the surrounding rock (Lister & Kerr Reference Lister and Kerr1991; Rubin Reference Rubin1995). The most common industrial application of hydraulic fracturing is well stimulation to facilitate the extraction of oil and gas from shale formations (Cueto-Felgueroso & Juanes Reference Cueto-Felgueroso and Juanes2013). A fluid is injected at high pressure to expand fractures initiated with small-scale explosions in unconventional reservoirs. The fractures create new flow pathways, facilitating fluid transport and storage in low-permeability and low-porosity rock formations. Other applications leverage the enhanced transport. For example, fractures connecting wells can be used to extract geothermal energy as the fluid pumped through the fracture heats up when it travels underground (Murphy et al. Reference Murphy, Tester, Grigsby and Potter1981; Caulk et al. Reference Caulk, Ghazanfari, Perdrial and Perdrial2016; Luo et al. Reference Luo, Zhu, Guo, Tan, Zhuang, Liu, Zhang, Xiang and Rohn2017). Fractures are also used for storage, including carbon sequestration (Huppert & Neufeld Reference Huppert and Neufeld2013; Jia, Tsau & Barati Reference Jia, Tsau and Barati2019) and disposal of liquid waste (Bao & Eaton Reference Bao and Eaton2016; Alessi et al. Reference Alessi, Zolfaghari, Kletke, Gehman, Allen and Goss2017).
The complex mechanics of fluid-driven fractures are controlled by the deforming boundary, fluid flow and stress singularity at the tip of the tensile fracture (Detournay Reference Detournay2016). Field testing, laboratory-scale experiments and predictive modelling evidence that the fracture propagation is characterized by multiple length and time scales. When a Newtonian fluid is injected from a point source, in an infinite, homogeneous and impermeable medium, a single fracture propagates radially. The elastic stress in the medium leads to the growth of a penny-shaped fluid-filled fracture in the direction of minimum confining stress. The propagation of such fractures has been studied extensively as it is essential to the modelling of more complex geological situations, including those involving a finite medium (Bunger & Detournay Reference Bunger and Detournay2005), complex fluids (Barbati et al. Reference Barbati, Desroches, Robisson and McKinley2016; Hormozi & Frigaard Reference Hormozi and Frigaard2017; Lai et al. Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018) and interacting fractures (O'Keeffe et al. Reference O'Keeffe, Zheng, Huppert and Linden2018b). Seminal work on the stress distribution in a penny-shaped fracture (Sneddon & Mott Reference Sneddon and Mott1946) and the injection of viscous fluids to form fractures (Khristianovic & Zheltov Reference Khristianovic and Zheltov1955; Barenblatt Reference Barenblatt1956) led to the development of self-similar solutions for fractures whose propagation is limited by the viscous dissipation in the fluid (Spence & Sharp Reference Spence and Sharp1985). Further work on the vicinity of the crack front or tip region identified two asymptotic regimes for the tip geometry and fracture propagation (Desroches et al. Reference Desroches, Detournay, Lenoach, Papanastasiou, Pearson, Thiercelin and Cheng1994; Garagash & Detournay Reference Garagash and Detournay2000, Reference Garagash and Detournay2005; Savitski & Detournay Reference Savitski and Detournay2002; Detournay & Garagash Reference Detournay and Garagash2003). In the viscous-dominated scaling, the viscous dissipation in the flow opposes the elastic stress of the deformable boundary to control the fracture evolution. In the toughness-dominated regime, the material toughness opposes the elasticity-driven propagation of the fracture and determines the system's behaviour. Laboratory-scale experiments commonly rely on clear brittle elastic gels to study the crack tip region and the penny-shaped fracture (Takada Reference Takada1990; Giuseppe et al. Reference Giuseppe, Funiciello, Corbi, Ranalli and Mojoli2009; Kavanagh, Menand & Daniels Reference Kavanagh, Menand and Daniels2013; Baumberger & Ronsin Reference Baumberger and Ronsin2020). Injections of water, glycerol and oil in gelatin and polyacrylamide have validated the existence of two propagation regimes and the corresponding scaling laws (Lai et al. Reference Lai, Zheng, Dressaire, Wexler and Stone2015, Reference Lai, Zheng, Dressaire and Stone2016; O'Keeffe, Huppert & Linden Reference O'Keeffe, Huppert and Linden2018a). The behaviour of the crack tip region was studied by injecting liquid between two plates of polymethylmethacrylate (PMMA) glued by an adhesive (Bunger & Detournay Reference Bunger and Detournay2008).
Most hydraulic fracturing processes involve multiphase flows and, in particular, displacement flows (Hormozi & Frigaard Reference Hormozi and Frigaard2017; Osiptsov Reference Osiptsov2017; Lai et al. Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018; Wang, Elsworth & Denison Reference Wang, Elsworth and Denison2018; Bessmertnykh, Dontsov & Ballarini Reference Bessmertnykh, Dontsov and Ballarini2021). During hydraulic fracturing operations, several fluids are injected, ranging from low-viscosity fluids to high-viscosity polymer solutions (Moukhtari & Lecampion Reference Moukhtari and Lecampion2018; Bessmertnykh et al. Reference Bessmertnykh, Dontsov and Ballarini2021) and proppant slurries (Hormozi & Frigaard Reference Hormozi and Frigaard2017; Wang et al. Reference Wang, Elsworth and Denison2018; Barboza, Chen & Li Reference Barboza, Chen and Li2021). This sequence of injections aims at increasing the surface area of the fracture and at keeping the fracture open during the hydrocarbon extraction. For example, carbon dioxide injection is a promising strategy to enhance oil recovery after primary production of shale oil reservoirs (Huppert & Neufeld Reference Huppert and Neufeld2013; Jia et al. Reference Jia, Tsau and Barati2019). The rapid injection of supercritical CO
$_2$ in water-filled fractures is followed by the slower permeation of CO
$_2$ into the rock and the migration of the oil into the fracture. This strategy increases the amount of oil recovered while storing carbon in the rock. Finally, enhanced geothermal systems rely on fractures in hot rocks to connect the injection and extraction wells (Murphy et al. Reference Murphy, Tester, Grigsby and Potter1981; Caulk et al. Reference Caulk, Ghazanfari, Perdrial and Perdrial2016; Parisio & Yoshioka Reference Parisio and Yoshioka2020). Fracturing fluids are injected first, and then the working fluid, commonly water or CO
$_2$, is pumped into the fracture to extract heat.
Displacement flows in porous media can give rise to complex out-of-equilibrium flow patterns when the invading fluid has a lower viscosity than the fluid that occupies the porous medium, and is referred to as the displaced or defending fluid. Practically, the patterns generated by liquid–liquid or gas–liquid displacement flows lead to preferential flow pathways in the porous medium. Extensive work has therefore been dedicated to the formation and geometry of the patterns, ranging from experimental to numerical and theoretical (Saffman & Taylor Reference Saffman and Taylor1958; Paterson Reference Paterson1981; Park & Homsy Reference Park and Homsy1984; Homsy Reference Homsy1987; Lenormand, Touboul & Zarcone Reference Lenormand, Touboul and Zarcone1988; Chen Reference Chen1989; Tanveer Reference Tanveer1993; Primkulov et al. Reference Primkulov, Pahlavan, Fu, Zhao, MacMinn and Juanes2019). Viscous and capillary forces can contribute to fluid/fluid displacement in a porous medium. A displacement flow is characterized by two dimensionless parameters: the viscosity or mobility ratio 
$M=\mu _{inv}/\mu _{def}$, and the capillary number of the invading fluid 
$Ca = {\mu _{inv} u}/{\gamma }$, where 
$\mu _{inv}$ and 
$\mu _{def}$ are the viscosities of the invading and defending fluids, respectively, 
$u$ is the characteristic velocity, and 
$\gamma$ is the surface tension of the interface between the two fluids. The influence of the two dimensionless parameters on the geometry of the invading front is summarized in Lenormand's phase diagram (Lenormand et al. Reference Lenormand, Touboul and Zarcone1988), which was revisited recently by Primkulov et al. (Reference Primkulov, Pahlavan, Fu, Zhao, MacMinn and Juanes2021) to include wettability. In summary, for large 
$Ca$, the viscous forces control the system dynamics and the interfacial forces are negligible. If the invading fluid is more viscous than the defending fluid (
$M\geq 1$), then a compact front or interface moves through the porous medium. If the invading fluid is less viscous than the defending one (
$M\ll 1$), then the Saffman–Taylor instability leads to an unstable front with the formation of a viscous fingering pattern, observed in porous media of different complexity, from Hele-Shaw cells (Saffman & Taylor Reference Saffman and Taylor1958; Tabeling, Zocchi & Libchaber Reference Tabeling, Zocchi and Libchaber1987) to intricate networks of pores and throats (Lenormand, Zarcone & Sarr Reference Lenormand, Zarcone and Sarr1983; Zhao, MacMinn & Juanes Reference Zhao, MacMinn and Juanes2016). For low 
$Ca$, the interfacial forces contribute to the system dynamics and result in more complex patterns at the interface between the two fluids, depending on the local pore geometry and wettability (Stokes et al. Reference Stokes, Weitz, Gollub, Dougherty, Robbins, Chaikin and Lindsay1986; Cottin, Bodiguel & Colin Reference Cottin, Bodiguel and Colin2010). The rich dynamics of displacement flows is reported in various model porous media, including networks of microchannels and rough fractures (Glass, Rajaram & Detwiler Reference Glass, Rajaram and Detwiler2003; Chen et al. Reference Chen, Fang, Wu and Hu2017; Yang et al. Reference Yang, Méheust, Neuweiler, Hu, Niemi and Chen2019). Yet the system geometry can delay the onset of viscous fingers and even suppress the Saffman–Taylor instability. In a converging Hele-Shaw cell, the stability of the interface depends on the mobility ratio, but also the characteristic velocity, the gradient of cell depth, and the contact angle at the interface (Al-Housseiny, Tsai & Stone Reference Al-Housseiny, Tsai and Stone2012; Al-Housseiny & Stone Reference Al-Housseiny and Stone2013; Lu et al. Reference Lu, Browne, Amchin, Nunes and Datta2019). Below a critical capillary number, a compact front is observed in a converging Hele-Shaw cell despite the unfavourable nature of the displacement. Similar results are reported in flexible cells, whose geometry depends on elastohydrodynamic interactions, such as displacement flows under elastic membranes (Pihler-Puzović et al. Reference Pihler-Puzović, Illien, Heil and Juel2012, Reference Pihler-Puzović, Périllat, Russell, Juel and Heil2013, Reference Pihler-Puzović, Juel, Peng, Lister and Heil2015; Peng et al. Reference Peng, Pihler-Puzovic, Juel, Heil and Lister2015). Two physical mechanisms contribute to the stabilization of the interface under an elastic membrane that deforms as fluid is injected. First, the local increase in cell depth leads to a depth gradient that has been shown to delay viscous fingering for rigid converging cells. Second, the increase in depth reduces the characteristic velocity or capillary number corresponding to a given injection flow rate.
The purpose of the present paper is to model axisymmetric two-phase flows in fluid-filled growing fractures. Experiments and theoretical modelling focus on immiscible two-phase flows with a mobility ratio smaller than or of order 1 and a low capillary number, ensuring the propagation of a compact front. We build on the approach of Savitski & Detournay (Reference Savitski and Detournay2002) to study the coupling between the two-phase flow and the fracture growth, in the viscous and toughness regimes (§ 2). We derive new scalings for the radius and aperture of the fracture and the position of the interface for immiscible displacement flows in elastic, brittle and impermeable media (§ 3). To test the scalings, we conduct injection experiments in gelatin, which is a common model medium (§ 4). During the two consecutive injections, we record the geometrical parameters of the fracture and compare their time dependence with scalings (§ 5). Finally, we discuss the time scales of the fracturing displacement flow (§ 6).
2. Theoretical models
In this section, we model the displacement flow that is responsible for the propagation of a crack during the successive injections of two immiscible fluids. The mathematical models presented build on the framework introduced originally by Spence & Sharp (Reference Spence and Sharp1985) and further developed in recent studies of single-fluid injection (Savitski & Detournay Reference Savitski and Detournay2002; Lai et al. Reference Lai, Zheng, Dressaire, Wexler and Stone2015; O'Keeffe et al. Reference O'Keeffe, Huppert and Linden2018a). Past work has focused on the injection of a single incompressible fluid in an elastic brittle solid through a point source (see figure 1), forming a penny-shaped crack. The fracture dynamics depends on the material properties, such as the Young's modulus 
$E$, Poisson's ratio 
$\nu$ and toughness 
$K_{IC}$, and the injection parameters, i.e. the constant flow rate 
$Q$ and liquid viscosity 
$\mu$. As the injection stops, the fracture reaches its final configuration and volume 
$V_0$.
Figure 1. Schematic of a penny-shaped fracture formed by successively injecting two liquids, first the outer fluid (dark grey) and then the inner fluid (light grey). The fracture is axisymmetric.
This study addresses the injection of an immiscible liquid in a pre-formed penny-shaped crack. The fluid is injected through the same point source at the centre of the crack. The displaced fluid fills an outer annular region of the fracture (see figure 1). In what follows, we use the subscript ‘in’ to refer to the injected liquid and ‘0’ to refer to the displaced fluid. The surface tension of the interface is denoted 
$\gamma$, and the contact angle with the solid is denoted 
$\theta$. Over the time scale of an experiment, typically a few minutes, the solid is not porous to the liquids, and the volume of the fracture is equal to the total volume of fluid injected.
We make assumptions regarding fluid flow and fracture propagation to model the system dynamics. As the injection rates are low, we assume that the fracture propagates at equilibrium, and the liquid and fracture fronts coincide at all times, with no fluid lag. The fluid injection results in linear elastic deformation of the surrounding material.
To describe the crack aperture 
$w(r,t)$, radius 
$R(t)$ and pressure 
$p(r,t)$, as well as the position of the interface 
$R_I$, we need to solve the coupled equations that describe (a) the viscous flow of the two fluid phases in the time-dependent fracture, (b) the elastic deformation of the solid material or fracture walls, (c) the stress intensity factor at the tip of the fracture, and (d) the volume conservation. These sets of equations are coupled by the net pressure in the fracture. The fluid domain is divided into two regions. The outer region is composed of the displaced fluid and bound by the liquid–liquid interface at 
$r= R_I$ and the crack tip at 
$r= R$. The injected fluid fills the inner region of the crack from the injection point to the fluid–fluid interface at 
$r= R_I$. Both regions are axisymmetric, as shown in figure 1.
2.1. Fluid flow in the crack
2.1.1. Lubrication theory
The low Reynolds number flow in the elongated fracture allows us to simplify the Navier–Stokes equation and use lubrication theory. The fluid is injected in a pre-formed crack whose aspect ratio is small:
\begin{equation} w \ll R. \end{equation}The Reynolds number of the flow through the crack is defined for the injected fluid as
\begin{equation} Re = \epsilon\,\frac{\rho Uw}{\mu} = \frac{\rho_{in} Q_{in}\,w(r=0)}{2 {\rm \pi}\mu_{in} R^2} \leq 1, \end{equation}
where 
$\epsilon$ is the aspect ratio of the crack. As a result, the flow of both fluids can be modelled with the lubrication theory, similarly to the single-phase flows in a penny-shaped fracture (Savitski & Detournay Reference Savitski and Detournay2002; Detournay Reference Detournay2004; Detournay & Peirce Reference Detournay and Peirce2014; Lai et al. Reference Lai, Zheng, Dressaire, Wexler and Stone2015; O'Keeffe et al. Reference O'Keeffe, Huppert and Linden2018a). For the two-fluid system, the lubrication equations are
$$\begin{gather} \frac{\partial w(r, t)}{\partial t}=\frac{1}{12 \mu_{in}}\,\frac{1}{r}\, \frac{\partial}{\partial r}\left(r\,w^{3}(r, t)\,\frac{\partial p}{\partial r}\right) \quad \text{for}\ 0\leq r \leq R_I, \end{gather}$$
$$\begin{gather}\frac{\partial w(r, t)}{\partial t}=\frac{1}{12 \mu_{0}}\,\frac{1}{r}\, \frac{\partial}{\partial r}\left(r\,w^{3}(r, t)\,\frac{\partial p}{\partial r}\right) \quad \text{for}\ R_I\leq r \leq R. \end{gather}$$2.1.2. Liquid interface
The interface between the two fluids moves outwards during the injection and is described by the dynamics boundary condition
\begin{align} & \boldsymbol{n}
\boldsymbol{\cdot} \left({-}p^{-} \boldsymbol{I}
+\mu_{in}\left(\boldsymbol{\nabla}
\boldsymbol{u}^{-}+ (\boldsymbol{\nabla}
\boldsymbol{u}^{-})^{\mathrm{T}}\right)\right)
\boldsymbol{\cdot} \boldsymbol{n}+\gamma \kappa_c \nonumber\\
&\quad =\boldsymbol{n} \boldsymbol{\cdot}
\left({-}p^{+}\boldsymbol{I}+\mu_{0}
\left(\boldsymbol{\nabla}
\boldsymbol{u}^{+}+(\boldsymbol{\nabla}
\boldsymbol{u}^{+})^{\mathrm{T}}\right)\right)
\boldsymbol{\cdot} \boldsymbol{n},
\end{align}
where 
$\boldsymbol {n} = \boldsymbol {e}_r$ is the vector normal to the interface, and 
$\kappa_c$ is the sum of the principal curvatures of the interface. The 
$+$ and 
$-$ exponents indicate that the value of the variable is determined at 
$r = R_I+ \delta$ and 
$r = R_I - \delta$ respectively, with 
$\delta \ll R_I$. We note that 
$\boldsymbol {{\tau }} = \mu (\boldsymbol {\nabla } \boldsymbol {u}+(\boldsymbol {\nabla } \boldsymbol {u})^{\mathrm {T}})$. We assume that the fluids are perfectly wetting the gel, and neglect the thin film deposited by the outer fluid. The normal stress balance is
\begin{align} p^{-} - p^{+} = \sigma_I &= \gamma \left(\frac{1}{R_I} + \frac{2}{w_I}\right) - {\boldsymbol{n}} \boldsymbol{\cdot} \left(\boldsymbol{\tau}^+-\boldsymbol{\tau}^-\right)\boldsymbol{\cdot} {\boldsymbol{n}} \nonumber\\ &= \gamma \left(\frac{1}{R_I} + \frac{2}{w_I}\right) - 2 \mu_{0} \left.\frac{\partial u}{\partial r}\right|_{r= R_I^{+}} - 2 \mu_{in} \left.\frac{\partial u}{\partial r}\right|_{r= R_I^{-}}, \end{align}
where 
$w_I$ is the width of the fracture at the interface. The expression can be simplified further as 
$R_I>w_{I}$, and the viscous normal stresses are negligible. Indeed, the capillary number of the invading fluid is
\begin{equation} Ca_{in} = \frac{\mu_{in} U}{\gamma} = \frac{\mu_{in} Q_{in}}{2 \gamma {\rm \pi}R\,w(r=0)} \ll 1. \end{equation}The pressure change across the interface is
\begin{equation} p^{-} - p^{+} \approx \frac{2\gamma}{w_I}. \end{equation}As the fluid–fluid interface moves, we can write the following kinematic condition using Reynolds equations:
\begin{equation} q^-{=} q^+{=} -\frac{w_I^3}{12\mu_{in}} \left.\frac{\partial p}{\partial r}\right|_{r= R_I^{-}} ={-}\frac{w_I^3}{12\mu_{0}} \left.\frac{\partial p}{\partial r}\right|_{r= R_I^{+}}. \end{equation}2.1.3. Volume conservation
Finally, through volume conservation, the volume of the crack is equal to the volume injected. We note that 
$V_0$, the volume of the pre-fracture, is equal to the volume of the outer fluid. The injection begins at 
$t=0$:
\begin{equation} V_0 + Q_{in}t= 2 {\rm \pi}\int_0^{R}r\,w(r,t)\,{\rm d} r \end{equation}and
\begin{equation} Q_{in}t= 2 {\rm \pi}\int_0^{R_I}r\,w(r,t)\,{\rm d} r. \end{equation}2.2. Fracture equations
2.2.1. Linear elasticity
Similarly to the single-fluid fracture (Savitski & Detournay Reference Savitski and Detournay2002), the linear elasticity equation is
\begin{equation} w(r, t)=\frac{8 (1-\nu^2) R}{{\rm \pi} E} \int_{r / R}^{1} \frac{\xi}{\sqrt{\xi^{2}-(r / R)^{2}}} \int_{0}^{1} \frac{x\,p(x \xi R, t)}{\sqrt{1-x^{2}}}\,\mathrm{d} x\,\mathrm{d} \xi. \end{equation}2.2.2. Fracture propagation
In a small region at the tip of the crack, the material undergoes plastic deformation. The tensile fracture propagates when the mode I stress intensity factor 
$K_I$ reaches a critical value called the toughness of the material, 
$K_{IC} = \sqrt {2 E^{'} \gamma _S}$, where 
$\gamma _S$ is the fracture surface energy of the solid material (Kanninen & Popelar Reference Kanninen and Popelar1985), and 
$E' = E/(1-\nu ^2)$. For a penny-shaped crack, the stress intensity factor near the tip is defined as (Rice Reference Rice1968)
\begin{equation} K_{{I}}=\frac{2}{\sqrt{{\rm \pi} R}} \int_{0}^{R(t)} \frac{p(r, t)}{\sqrt{R^{2}-r^{2}}}\,r\,\mathrm{d} r. \end{equation}2.3. Boundary conditions at the fracture tip and injection point
At the tip of the crack, the width 
$w(R)$ goes to zero,
\begin{equation} w=0, \quad r=R(t), \end{equation}and the flow rate also goes to zero,
\begin{equation} w^{3}(r, t)\,\frac{\partial p(r, t)}{\partial r}=0, \quad r=R(t). \end{equation}At the point source, the local flow rate is equal to the injected flow rate:
\begin{equation} 2 {\rm \pi}\lim _{r \rightarrow 0} r\,q(r, t)=Q_{in}. \end{equation}3. Scaling
The equations are non-dimensionalized by identifying the characteristic scales in both phases as listed in table 1.
Table 1. Rescaled parameters.
The characteristic radius and aperture of the fracture are 
$R$ and 
$W_0$, respectively, with 
$R$ the radius of the fracture, and 
$W_0$ the maximum aperture of the fracture at 
$r=0$. The position of the interface is 
$R_I$. The characteristic pressure in both the fluids is taken to be 
$P_{0}$. We define effective material parameters 
${\mu '}$, 
$E'$ and 
$K'$ as proposed by Savitski & Detournay (Reference Savitski and Detournay2002):
\begin{equation} \mu_{in}' = 12\mu_{in},\quad \mu_{0}' = 12\mu_{0},\quad K' = 4\left(\frac{2}{\rm \pi}\right)^{{1}/{2}}K_{IC}\quad \mbox{and}\quad E' = \frac{E}{1-\nu^2}. \end{equation}
To compare the viscosity of the two liquids, we introduce the parameter 
$M = \mu _{in}'/\mu _{0}'$. When the two fluids are present in the fracture, we used a weighted average to define the resulting viscosity 
$\mu _{e}$ and the corresponding effective value 
$\mu '_{e} = 12 \mu _{e}$.
Using the characteristic parameters summarized in table 1 and the effective material properties, we obtain the following set of equations to describe the displacement flow and fracture propagation.
(i) Lubrication theory (from (2.3) and (2.4)):
(3.2)
\begin{equation} {\frac{\mu'_e R^2}{W_0^2 P_{0}}\,\frac{\partial \hat{w}}{\partial t}= \frac{1}{\hat{r}}\,\frac{\partial}{\partial \hat{r}}\left(\hat{r} \hat{w}^{3}\, \frac{\partial \hat{p}}{\partial \hat{r}}\right)}. \end{equation}(ii) Linear elasticity (from (2.12)):
(3.3)\begin{equation} {\hat{w} = \frac{8 R}{{\rm \pi} E' W_0} \left[ P_0\int_{0}^{1} \int_{0}^{1} \frac{x\,\hat{p}(x \xi R,t)}{\sqrt{1-x^2}}\,{{\rm d}x}\,{\rm d}\xi \right]}. \end{equation}(iii) Fracture propagation (from (2.13)):
(3.4)\begin{equation} {\frac{K'}{P_{0} \sqrt{R}}= \frac{2^{7/2}}{\sqrt{\rm \pi}}\int_0^{1} \frac{\tilde{p} \tilde{r}}{\sqrt{1-\tilde{r}^2}}\,{\rm d}\tilde{r}}. \end{equation}(iv) Global mass balance (from (2.10) and (2.11)):
(3.5)and\begin{equation} \frac{Q_{in}t}{2 {\rm \pi}R_I^2 W_0}= \int_0^{1} \hat{r} \hat{w} \,{\rm d}\hat{r} \end{equation}(3.6)\begin{equation} \frac{V_0}{2 {\rm \pi}R^2 W_0} = \int_0^{1} \tilde{r} \hat{w} \,{\rm d}\tilde{r} - \left(\frac{R_I}{R}\right)^2 \int_0^{1} \hat{r} \hat{w} \,{\rm d}\hat{r}. \end{equation}
As the fluid is injected, the elastic pressure drives the propagation of the crack, which is resisted by the viscous dissipation associated with the motion of the injected and displaced fluids, and the fracture toughness. We first assume that the interfacial pressure is negligible compared to the viscous and toughness-related stresses. We then assume that one of the resisting stresses controls the propagation and balances the elastic stress. Studies on single-fluid injection have validated this approach, with experimental evidence of the two asymptotic regimes. If the viscous stresses in the fluids are larger than the fracture-opening stress, then the fracture propagation is said to be in the viscous regime. Alternatively, the propagation is in the toughness regime. For both regimes, we can derive scaling arguments from the dimensionless groups derived from (3.2)–(3.6).
3.1. Toughness regime
For the toughness scaling, we set the non-dimensional groups in (3.3)–(3.6) equal to 1. Indeed, the viscous stresses are negligible, and the crack opening is limited by the toughness of the material. The scaling relations are
$$\begin{gather} \frac{W_0E'}{{P_0}R} = 1, \end{gather}$$
$$\begin{gather}\frac{K'}{{P_0} \sqrt{R}} = 1, \end{gather}$$
$$\begin{gather}\frac{Q_{in}t}{R_I^2 W_0} = 1, \end{gather}$$
$$\begin{gather}\frac{V_0}{ W_0} = R^{2} - R_I^{2}. \end{gather}$$
We define the characteristic time scale 
$T = V_0/Q_{in}$ and the corresponding dimensionless time 
$\tilde {t} = t/T$. By combining those groups, we obtain the toughness-dominated scaling of the fracture properties:
$$\begin{gather} W_0 = \left(\frac{K'}{E'}\right)^{4/5} V_0^{1/5} \left(1+\tilde{t}\right)^{1/5}, \end{gather}$$
$$\begin{gather}R = \left(\frac{E' V_0}{K'}\right)^{2/5} \left(1+\tilde{t}\right)^{2/5}, \end{gather}$$
$$\begin{gather}R_I = \left(\frac{E' V_0}{K'}\right)^{2/5} \tilde{t}^{1/2} \left(1+\tilde{t}\right)^{{-}1/10} \end{gather}$$and
\begin{equation} {P_0} = K' \left(\frac{K'}{E' V_0}\right)^{1/5} \left(1+\tilde{t}\right)^{{-}1/5}. \end{equation}3.2. Viscous regime
For the viscous scaling, we set the non-dimensional groups in (3.2)–(3.3) and (3.5)–(3.6) equal to 1. Here, the propagation of the fluids follows the lubrication equation, and the viscous dissipation limits the propagation of the fracture:
$$\begin{gather} {\frac{W_0E'}{P_{0}R}= 1}, \end{gather}$$
$$\begin{gather}{\frac{\mu'_{e} R^2}{W_0^2 P_{0}\left(t+V_0/Q_{0}\right) } = 1}, \end{gather}$$
$$\begin{gather}{\frac{Q_{in}t}{ R_I^2 W_0} = 1}, \end{gather}$$
$$\begin{gather}{\frac{V_0}{ W_0} = R^{2} - R_I^{2}}. \end{gather}$$
We define the effective viscosity of the volume of fluid in the fracture as a weighted average. The effective viscosity depends on the viscosities of the two fluids in the fracture and their relative volumes at time 
$\tilde {t}$:
\begin{equation} \mu'_e = \frac{\mu'_{0}+\mu'_{in}\tilde{t}}{1 + \tilde{t}}. \end{equation}We obtain the viscous-dominated scaling of the variables:
$$\begin{gather} {W_0 = \left(\frac{\mu'_{e}Q_{in}}{E'}\right)^{2/9} V_0^{1/9} \left(1+\tilde{t}\right)^{1/3}\left(\alpha+\tilde{t}\right)^{{-}2/9},} \end{gather}$$
$$\begin{gather}{R = \left(\frac{E'}{\mu'_{e}Q_{in}}\right)^{1/9} V_0^{4/9} \left(1+\tilde{t}\right)^{1/3} \left(\alpha+\tilde{t}\right)^{1/9},} \end{gather}$$
$$\begin{gather}{R_I = \left(\frac{E'}{\mu'_{e}Q_{in}}\right)^{1/9} V_0^{4/9} \tilde{t}^{1/2} \left(1+\tilde{t}\right)^{{-}1/6} \left(\alpha+\tilde{t}\right)^{1/9},} \end{gather}$$
$$\begin{gather}{P_{0} = \left(\frac{\mu'_{e}Q_{in}E^{\prime 2}}{V_0}\right)^{1/3}\left(\alpha + \tilde{t}\right)^{{-}1/3}}, \end{gather}$$
with 
$\alpha = Q_{in}/{Q_{0}}$. The results are summarized in table 2.
Table 2. Scaling relations for a penny-shaped fracture driven by a displacement flow. The time evolution of the geometrical properties of the fracture depends on the dominant resisting stress, which can be viscous or toughness-related.
3.3. Discussion
We consider the single-fluid limit of the two asymptotic regimes. Indeed, for a pre-fracture volume equal to zero, 
$V_0=0$, or for large values of the injection time 
$t\gg 1$, the expressions derived above should be equal to those obtained previously for single-fluid injection. In the toughness regime, we recover the single-fluid scaling relations (Savitski & Detournay Reference Savitski and Detournay2002; Lai et al. Reference Lai, Zheng, Dressaire and Stone2016; O'Keeffe et al. Reference O'Keeffe, Huppert and Linden2018a)
$$\begin{gather} W_0 = \left(\frac{K'}{E'}\right)^{4/5} \left(Q_{in}t\right)^{1/5}, \end{gather}$$
$$\begin{gather}R = \left(\frac{E'}{K'}\right)^{2/5} \left(Q_{in}t\right)^{2/5} \end{gather}$$and
\begin{equation} P_0 = K' \left(\frac{K'}{E' Q_{in} t}\right)^{1/5}. \end{equation}In the viscous regime, we also recover the scaling relations for a single-fluid injection:
$$\begin{gather} W_0 = \left(\frac{\mu'_{in}}{E'}\right)^{2/9} Q_{in}^{1/3} t^{1/9}, \end{gather}$$
$$\begin{gather}R = \left(\frac{E'}{\mu'_{in}}\right)^{1/9} Q_{in}^{1/3} t^{4/9} \end{gather}$$and
\begin{equation} P_0 = \left(\frac{\mu'_{in}E^{\prime 2}}{t}\right)^{1/3}. \end{equation} The fluid and matrix properties determine whether the fracture propagation is in the viscous or toughness regime. Past studies in single-fluid injection have defined criteria to predict the propagation regime. The regime is defined by the largest of the two stresses that oppose the elastic stress: the toughness-related stress 
$\Delta P_m \approx {K'}/{\sqrt {R}}$, and the viscous stress 
$\Delta P_v = {\mu 'Q}/{W_0^3}$. In the viscous regime, we can use the scaling relations for 
$R$ and 
$W_0$ to estimate the ratio:
\begin{equation} \left(\frac{\Delta P_m}{\Delta P_v}\right)_v = \left(\frac{K^{\prime 9} t}{E^{\prime 13/2} Q^{3/2} \mu^{\prime 5/2}}\right)^{1/9} = \left(\frac{t}{t_{mk}}\right)^{1/9}= \kappa, \end{equation}where
\begin{equation} t_{mk} = \left(\frac{E^{\prime 13/2} Q^{3/2} \mu^{\prime 5/2}}{K^{\prime 9} }\right) \end{equation}
is the characteristic time scale of the system, and 
$\kappa$ is the dimensionless toughness. Similarly, we can define the relative magnitude of both stresses in the toughness regime using the corresponding scaling relations:
\begin{equation} \left(\frac{\Delta P_m}{\Delta P_v}\right)_m = \left(\frac{t}{t_{mk}}\right)^{2/5}= \kappa^{18/5}. \end{equation}
The viscous-dominated propagation is therefore associated with small values of 
$\kappa$, i.e. 
$\kappa \ \mbox {or}\ t/t_{mk} \lesssim 1$, and the toughness-dominated dynamics for 
$\kappa \ \mbox {or}\ t/t_{mk} \gg 1$.
Similarly, we can define the propagation regime of the fracture formed by a displacement flow. The toughness-related stress remains 
$\Delta P_m \approx K'/\sqrt {R}$, while the viscous stress becomes 
$\Delta P_v =\mu _e'Q_{in} {R^2}({R_I^2 W_0^3+ \alpha V_0 W_0^2})^{-1}$ as defined in (3.8) and (3.16).
In the viscous regime, we substitute 
$R$, 
$R_I$ and 
$W_0$ by the scaling relations defined in table 2:
\begin{equation} \left(\frac{\Delta P_m}{\Delta P_v}\right)_v = \frac{K' V_{0}^{1/9}(\alpha+\tilde{t})^{5/18}}{E^{\prime 13/18} \mu_{e}^{\prime 5/18}Q_{in}^{5/18}(1 +\tilde{t})^{1/6}}. \end{equation}Similarly, in the toughness regime, we estimate the ratio to be
\begin{equation} \left(\frac{\Delta P_m}{\Delta P_v}\right)_m = \left(\frac{\Delta P_m}{\Delta P_v}\right)_v^{18/5}. \end{equation}
The ratio of the viscous and toughness-related stresses is a function of time that increases as 
$\tilde {t}^{1/9}$ for large values of 
$\tilde {t}$. In consequence, the propagation becomes toughness-controlled for long-time or large-volume injections. This result is consistent with what is known for a single-fluid injection. Contrary to the single-fluid criterion however, there is no explicit solution for the threshold injection time in the case of displacement flow. We can determine the propagation regime at time 
$\tilde {t}$ by comparing the pressure ratio with 1.
4. Experiments
Laboratory-scale experiments commonly use hydrogels as rock analogues to study hydraulic fracturing in brittle elastic materials (O'Keeffe & Linden Reference O'Keeffe and Linden2017; Baumberger & Ronsin Reference Baumberger and Ronsin2020). In particular, gelatin is a clear gel whose elasticity can be tuned easily by varying the volume fraction of gelatin powder in water (Giuseppe et al. Reference Giuseppe, Funiciello, Corbi, Ranalli and Mojoli2009; Kavanagh et al. Reference Kavanagh, Menand and Daniels2013; Lai et al. Reference Lai, Zheng, Dressaire, Wexler and Stone2015). Because gelatin expands as it sets in the container, the material is spontaneously compressed, which furthers the analogy with soft rocks for fracture studies.
4.1. Material preparation and characterization
The properties of the gelatin and the injected fluids control the fracture dynamics and are therefore characterized systematically. The gelatin is prepared by first heating ultra-pure water to 
$60\,^{\circ }$C. While stirring the heated water, we slowly add gelatin powder (Gelatin type A; Sigma-Aldrich, USA) at a mass fraction of 10 % to 30 % to vary Young's modulus of the resulting gel. The gelatin then cools down to room temperature and sets over 24 hours prior to testing or fracturing. The Young's modulus of gelatin is measured with a custom-built displacement-controlled load frame. The sample is compressed by a ballscrew stage whose speed is set by a stepper motor Parker Compumotor OS22B controlled by a controller Parker Compumotor ZL6104. A load cell, Eaton 3108-10 (10 lb capacity) measures the force generated by the compressed sample. We record force–displacement values for cylindrical samples of gelatin of diameter and height equal to 1 in, and compute the stress–strain curves of the material. At small strain, all gelatin samples exhibit a linear elastic response to the compression. As listed in table 3, Young's modulus ranges between 
$15$ and 116 kPa with a measurement error of 
${\pm }10\,\%$. The fracture energy and Poisson's ratio of the gelatin are assumed constant in our experiments, with 
$\gamma _S \approx 1$ J m
$^{-2}$ and 
$\nu \approx 0.5$ (Menand & Tait Reference Menand and Tait2002).
Table 3. List of experiments: symbols and parameters.
To study displacement flows in fluid-filled fractures, we use immiscible Newtonian liquids. Silicone oils of different viscosity are used to form the pre-fracture. An aqueous solution composed of water, glycerol or corn syrup is then injected. The silicone oil is displaced outwards, further expanding the fracture in the gelatin. We measure the viscosity of the fluids and the surface tension at the oil/water and oil/syrup interfaces. Viscosity measurements are conducted using an MCR 92 Anton Parr rheometer with a parallel plate measuring system at 
$20\,^{\circ }$C. The values obtained have a measurement error of 
${\pm }1\,\%$ and are listed in table 3. Surface tension measurements are conducted using the pendant drop method with an Attension Theta Flex tensiometer: the surface tension between silicone oil and water is 
$\gamma _{o/w} \approx 35 \pm 2$ mN m
$^{-1}$, and between silicone oil and syrup is 
$\gamma _{o/s} \approx 50 \pm 2$ mN m
$^{-1}$. The fluid properties are selected to study the propagation of the pre-fracture and fracture in a single regime during the experiment. For example, the values of 
$\kappa$ at the end of the injection forming the pre-fracture, for experiments 1–10, vary as 
$35 \leq \kappa \leq 55$. Those experiments are expected to be in the toughness-dominated regime. Similarly, for experiments 11–16, 
$1.7 \leq \kappa \leq 1.9$ at the end of the formation of the pre-fracture: these experiments target the viscous-dominated propagation. These values are consistent with previous work on the formation of fluid-driven fractures in a block of gelatin. During the formation of the fracture through the displacement flow, we estimate the pressure ratio in the viscous regime 
$({\Delta P_m}/{\Delta P_v})_v$ at the end of the injection. The ratio varies between 24 and 58 for experiments 1–10, which are therefore expected to be in the toughness regime. The ratio varies between 1.9 and 2.5 for experiments 11–16, which are therefore expected to be in the viscous regime.
4.2. Set-up
The injection experiments are conducted in a large block of gelatin to avoid boundary effects on the propagation of the fracture. The gelatin is set in a cubic clear container (
$15\,{\rm cm} \times 15\,{\rm cm}\times 15\,{\rm cm}$) around a blunt needle as represented in figure 2(a). We use a thin needle (inner diameter 
$ID = 1$ mm) for low-viscosity injections, and a wider needle for high-viscosity injections (
$ID = 2.15$ mm). In our system, the fracture propagates horizontally, in the direction that opposes the least resistance. To avoid small tilts of the fracture that would compromise the quality of the recording, a plastic washer (outer diameter 
$OD=6.5$ mm) is placed at the tip of the needle to initiate the fracture in the horizontal plane. Two fluids are successively pumped into the gelatin. For each fluid, we use a kdScientific
Legato 200 syringe pump to set the injection flow rate at a value between 
$0.1$ and 
$20$ ml min
$^{-1}$ with accuracy 
${\pm }0.35\,\%$. Both fluid-filled syringes are connected to the same injection needle using a switch valve. First, a silicone oil of viscosity 
$\mu _0$ is injected at a constant volumetric flow rate 
$Q_0$ in the gelatin matrix to form the pre-fracture. The injection stops when a volume 
$V_0$ of silicone oil has been dispensed. The valve is then switched to inject the aqueous phase of viscosity 
$\mu _{in}$ into the pre-fracture at flow rate 
$Q_{in}$. The injection stops when the fracture tip is within 2 cm of the container walls to avoid confinement effects (Bunger, Gordeliy & Detournay Reference Bunger, Gordeliy and Detournay2013). The two fluids are dyed to allow visualizing the propagation of the fracture and interface between the two liquids in the clear gelatin: we use a blue water-based food dye for the aqueous phase, and a red oil-based food dye for the oil phase, as shown in figure 2. The propagation is recorded using a Nikon D5300 camera with a Phlox
LED panel ensuring uniform backlighting. The images are processed using a custom-made MATLAB code to determine the radius of the fracture 
$R$ and the position of the interface between the two liquids 
$R_I$.
Figure 2. (a) Schematic of the experimental set-up. (b–e) Time evolution of the fracture formed in experiment 1. The water dyed with blue food colour is injected at 
$Q_{in} = 0.5$ ml min
$^{-1}$ in a pre-fracture formed with silicone oil dyed with red food colour. The recording starts when the water injection begins, and the experimental images are taken at (b) 
$t =50$ s, (c) 
$t = 300$ s, (d) 
$t = 1300$ s, and (e) 
$t =2800$ s. Supplementary movies 1 and 2 (available at https://doi.org/10.1017/jfm.2022.954) show the complete time evolution of the pre-fracture and fracture, respectively.
4.3. Measurement of the fracture width
To measure the thickness of the fracture during the propagation, we use a light absorption method (Bunger Reference Bunger2006; Bunger & Detournay Reference Bunger and Detournay2008). This method consists in selecting a soluble dye and the corresponding optical filter. The filter should transmit light to the camera at the wavelength at which the dye absorbance 
$A$ is maximum, 
$A_{\lambda }$. The absorbance follows Beer's law:
\begin{equation} A_{\lambda} ={-} \log_{10}\left(\frac{I_{\lambda}}{I_{\lambda,0}}\right)= \epsilon_{\lambda} c h, \end{equation}
where 
$I_{\lambda,0}$ is the background intensity, and 
$I_{\lambda }$ is the intensity when light passes through a liquid layer of thickness 
$h$, with dye concentration 
$c$. The parameter 
$\epsilon _{\lambda }$ characterizes the absorbance of the dye at the selected wavelength and is obtained through calibration. To measure the thickness of both liquids in the fracture, we use two dyes, one water-soluble and one oil-soluble, and record the absorbance using a single optical filter. To get accurate measurements, we select dyes with a large absorbance at the same wavelength. In all experiments, the water-soluble dye is nigrosin (Sigma-Aldrich) at 0.05 g l
$^{-1}$. Nigrosin is a black dye that absorbs at all wavelengths. We use different dyes depending on the viscosity of the silicone oils, because of their solubility limit. We dilute Sudan red (Sigma-Aldrich) in the low-viscosity silicone oils, i.e. 10 and 20 mPa s silicone oils, at 0.05 g l
$^{-1}$, and Nile red (Sigma-Aldrich) in high-viscosity silicone oils, i.e. 10 300 and 30 000 mPa s silicone oils, at 0.2 g l
$^{-1}$. For experiments with Sudan red in the oil phase and nigrosin in the aqueous phase, we use a 520 nm optical filter. When Nile red dyes the oil phase and nigrosin the aqueous phase, we use a 632 nm optical filter.
To measure the thickness of the fracture, we first conduct calibration experiments for each dye solution. For the low-viscosity solutions (water, 10 and 20 mPa s silicone oils), we use a glass wedge with an aperture that increases linearly from 0 to 10 mm, as shown in figure 3. The wedge is placed on the LED panel and the light intensity is obtained by taking a picture of the wedge with the optical filter mounted on the camera. The background intensity corresponds to an empty wedge. The light intensity of a liquid-filled wedge decreases as the thickness of the aperture increases (see figure 3b). The grey values are used to determine the absorbance as a function of the liquid thickness, as plotted in figure 3(c). Using a linear fit, we get 
$1/{\epsilon _{\lambda } c} = 47.2$ mm for Sudan red in oil, and 
$1/{\epsilon _{\lambda } c} = 32.6$ mm for nigrosin in water, for 
$\lambda = 520$ nm. For the high-viscosity samples, we use rectangular cells that are easier to fill. The cell thickness or height ranges from 0.3 to 3 mm. For each cell, we measure the background intensity of the empty cell and the intensity of the cell filled with the viscous fluid, i.e. syrup, or 10 000 or 30 000 mPa s silicone oil. We obtain the absorbance for a set of thickness values as shown in figure 3. Using a linear fit, we get 
$1/{\epsilon _{\lambda } c} = 56.5$ mm for Nile red in oil, and 
$1/{\epsilon _{\lambda } c} = 44$ mm for nigrosin in water, for 
$\lambda = 632$ nm. The fitting parameters obtained are then used to obtain the width of the fracture using Beer's law:
\begin{equation} h = \frac{A_{\lambda}}{\epsilon_{\lambda} c} = \frac{-\log_{10}\left(\dfrac{I_{\lambda}}{I_{\lambda,0}}\right)}{\epsilon_{\lambda} c}. \end{equation}
Figure 3. Calibration. (a) Schematic of the calibration experiment for low-viscosity fluids. (b) The intensity gradient was recorded for a wedge filled with water dyed with nigrosin at 0.05 g l
$^{-1}$ through the 520 nm filter. (c) Absorbance measured for nigrosin-dyed water at 
$\lambda = 520$ nm (
$\blacktriangle$), nigrosin-dyed syrup at 
$\lambda = 632$ nm (
), Sudan-red-dyed 10 mPa s siliconeoil at 
$\lambda = 520$ nm (
$\vartriangle$), and Nile-red-dyed 10 Pa s silicone oil at 
$\lambda = 632$ nm (
). The solid lines are the best linear fit for each calibration data set.
The concentrations of dyes chosen for this study, some of which are limited by solubility, are all sufficiently low for the absorbance to vary linearly with the sample thickness or fracture aperture. The accuracy of the measurements is limited by the noise due to low absorbance values.
5. Results
In this section, we present the results of the experiments listed in table 3. We measure the radius and aperture of the pre-fracture and fracture, and compare the data with the scalings obtained for the viscous and toughness regimes.
The fracture propagates radially upon injection of the fluid. Viscous fingering is not observed in our experiments, as illustrated in figures 2(b–e). The viscosity ratio 
$M = {\mu _{in}/}{\mu _{out}}$ varies between 0.05 and 1, and the thickness of the crack is of the order of a few millimetres, so the characteristic wavelength of the instability is larger than the perimeter of the injected fluid region. We estimate the experimental error by conducting an error analysis based on the scaling laws and the measurement error of the various parameters: 
$10\,\%$ for the Young's modulus 
$E$, 
$10\,\%$ for the fracture toughness 
$K$, 
$0.35\,\%$ for the flow rate 
$Q$, 
$1\,\%$ for the viscosity 
$\mu$, and 
$5\,\%$ for the volume 
$V_0$. In the toughness regime, the experimental error on the radius is estimated to be 
${\sim }10\,\%$. In the viscous regime, the error is 
${\sim }4\,\%$.
5.1. Single-fluid injection
We prepare the pre-fracture by injecting silicone oil into the gelatin cube. After an initial pressure build-up, the oil propagates rapidly over the washer at the tip of the needle. A radial fracture then forms around the washer, propagating more slowly with the oil filling the gap between the two gelatin surfaces. The homogeneous properties of the gelatin result in axisymmetric fractures for both low- and high-viscosity fluids. The radius of the fracture is measured during the injection and compared with the scaling for the toughness- and viscous-dominated regimes for a single fluid ((3.25) and (3.28), respectively). Experiments 1–10 are in the toughness regime as low-viscosity silicone oil is injected in soft gelatin (see table 3). Experiments 11–16 are in the viscous regime as high-viscosity silicone oil is injected in harder gelatin. In figure 4(a), we plot the results of the experiments in the toughness regime. The radius increases with time, and the rescaled radius follows a 
$t^{2/5}$ power law, which is consistent with (3.25). In the log-log scale, the best fit line with slope 
$2/5$ has a prefactor 
$k = 0.7$, which is in agreement with the theoretical prefactor 0.85 derived by Savitski & Detournay (Reference Savitski and Detournay2002). These results are also consistent with previous experimental data in this regime (Lai et al. Reference Lai, Zheng, Dressaire and Stone2016; O'Keeffe et al. Reference O'Keeffe, Huppert and Linden2018a). In figure 5(a), we report the data obtained in the viscous regime. The experimental conditions for those experiments are similar, and the results demonstrate the high reproducibility of the experiments (Lecampion et al. Reference Lecampion, Desroches, Jeffrey and Bunger2017). The radius increases as a power law of time. Upon rescaling the radius with the viscous scaling parameter, the data collapse on a line of slope 
$4/9$ with prefactor 
$k=0.36$. This result differs from the theoretical prefactor 
$0.7$, yet it is comparable to the values obtained in previous experimental studies in the viscous regime (Lai et al. Reference Lai, Zheng, Dressaire, Wexler and Stone2015; O'Keeffe et al. Reference O'Keeffe, Huppert and Linden2018a). Due to the initiation transient and finite size of the container, the power-law fits span about a decade of the log-log plots, which is common for laboratory-scale experiments.
Figure 4. Dynamics of the pre-fracture for low-viscosity oils. (a) Dependence of rescaled fracture radius on time for experiments 1–10 (see table 3 for corresponding injection parameters). The radius is rescaled using (3.25) in the toughness scaling for single-fluid injection. The origin for time is set when the oil enters the gelatin. The black curve represents the best linear fit with slope 
$2/5$. Inset: dependence of the radius 
$R$ on time 
$t$. (b) Rescaled fracture thickness profiles based on (3.24) and (3.25) at 
$t = 245,370,495,620,745$ s, with time increasing from clear to dark grey. Experimental parameters: 
$E=30$ kPa, 
$\mu _0= 10$ mPa s, 
$Q_0 = 0.3$ ml min
$^{-1}$.
Figure 5. Dynamics of the pre-fracture for high-viscosity oils. (a) Dependence of rescaled fracture radius on time for experiments 11–16 (see table 3 for corresponding injection parameters). The radius is rescaled using (3.28) in the viscous scaling for single-fluid injection. The black curve represents the best linear fit with slope 
$4/9$. Inset: dependence of the radius 
$R$ on time 
$t$. (b) Rescaled fracture thickness profiles based on (3.27) and (3.28) at 
$t = 9,11.5,14,16.5,19$ s, with time increasing from clear to dark grey. Experimental parameters: 
$E=88$ kPa, 
$\mu _0= 10$ Pa s, 
$Q_0 = 10$ ml min
$^{-1}$.
For each regime, we measure the fracture aperture using a dye whose absorbance varies linearly with the aperture. The results presented show the evolution of the fracture cross-section over time for one experiment in the toughness regime (see figure 4b) and one in the viscous regime (see figure 5b). For both experiments, the curves show the aperture as a function of the radial position at different injection times. Because the needle disturbs the absorbance measurement near the centre of the fracture, the aperture is not measured near the needle, i.e. for small values of the radius. Upon integration of the thickness curve recorded when the injection is complete, we obtain the value of 
$V_0 \pm 0.5$ ml. The aperture-radius curves are rescaled using the scaling for the aperture and radius. In both regimes, the curves collapse on a self-similar profile. The two data sets presented here are representative of the pre-fracture obtained for all the experiments conducted in this study, and are similar to results reported previously (Lai et al. Reference Lai, Zheng, Dressaire and Stone2016; O'Keeffe et al. Reference O'Keeffe, Huppert and Linden2018a).
5.2. Displacement flow
Once the pre-fracture is formed, the valve is immediately switched to the second immiscible to avoid further propagation of the pre-fracture (Möri & Lecampion Reference Möri and Lecampion2021). The material properties of the gelatin contribute to the definition of the propagation regime. Low-stiffness gelatin allows the observation of the toughness regime, while the viscous regime is most commonly reached in stiffer hydrogel. Thus we characterize the displacement flow and fracture propagation in the toughness regime with the experimental systems 1–10 (see table 3), in which the pre-fracture is also formed in the toughness regime. The injected fluid is water. Similarly, we investigate the viscous regime in experiments 11–16. The displacing fluid is a syrup whose high viscosity is of the same magnitude as the silicone oil in the pre-fractures. During the displacement flow, the fracture continues its radial expansion but with a different dynamic from the one observed during the formation of the pre-fracture. To study the radial expansion of the fracture during the injection, we track the position of the interface between the two fluids, 
$R_I$, and the radius of the fracture, 
$R$, over time. In figure 6, we plot the results of the experiments in the toughness regime. The radial position of the interface increases with time, similarly to what was observed for a single-fluid injection. The annular region of displaced fluid between 
$R_I$ and 
$R$ moves outwards, and its radius 
$R-R_I$ decreases over time. Using the equations in table 2, we plot the rescaled radii with respect to the relevant dimensionless time on the log-log scale. The rescaled radial positions of the interface collapses on a line of slope 1 with prefactor 1. The rescaled radii of the fracture collapse on a line of slope 1 with prefactor 0.3. We note that for large values of time 
$\tilde {t}\gg 1$, the fracture dynamics for the displacement flow is expected to become similar to the fracture dynamics for a single fluid. Indeed, if we plot the dimensionless radius as 
$R^{5/2} ({K'}/{E'V_0})$, then the prefactor 0.3 for the displacement flow is comparable with the prefactor for a single fluid 
$0.7^{5/2} \approx 0.4$. In figure 7, we plot the results of the displacement flow experiments in the viscous regime. The results can be rescaled using equations in the corresponding column of table 2. We plot the rescaled radii with respect to the relevant functions of the dimensionless time on the log-log scale. The rescaled radial positions of the interface collapse on a line of slope 1 with prefactor 0.37. The rescaled radii of the fracture collapse on a line of slope 1 with prefactor 0.36. We note that for large values of time 
$\tilde {t}\gg 1$, the fracture dynamics for the displacement flow is expected to become similar to the fracture dynamics for a single fluid. Since at large values of 
$\tilde {t}$, 
$(1+\tilde {t})^{1/3}(\alpha +\tilde {t})^{1/9} \approx \tilde {t}^{4/9}$, the identical prefactor for the scaling of the radius for the pre-fracture and the fracture indicates that indeed the single-fluid behaviour is recovered for large injection times. We also note that for large values of the fracture radius, the experimental data fall below the trend line. This is due to the slow down of the growth due to the confinement of the gelatin block. Experimentally, this is associated with the tilting of the fracture and the formation of finger-like structures.
Figure 6. Dynamics of the fracture for the low-viscosity aqueous phase. (a) Dependence of the position of the interface on time for experiments 1–10 (see table 3 for corresponding injection parameters). (b) Rescaled interface position as a function of rescaled time, based on the toughness scaling laws in table 2. (c) Dependence of the fracture radius on time. (d) Rescaled radius as a function of rescaled time, based on the toughness scaling laws in table 2. The black curves represent the best linear fit with slope 1.
Figure 7. Dynamics of the fracture for the high-viscosity aqueous phase. (a) Dependence of the position of the interface on time for experiments 11–16 (see table 3 for corresponding injection parameters). (b) Rescaled interface position as a function of rescaled time, based on the viscous scaling laws in table 2. (c) Dependence of the fracture radius on time. (d) Rescaled radius as a function of the rescaled time, based on the viscous scaling laws in table 2. The black curves represent the best linear fit with slope 1.
For the displacement flows, we measure the fracture aperture in the two fluids simultaneously, relying on two different dyes (see figures 8 and 9). Since the two dyes absorb the light differently, the absorbance profile presents a discontinuity across the interface. The absorbance values are converted to thickness measurements using Beer's law and the calibration parameters. The experimental results plotted in figures 8(a) and 9(a) show that the aperture is a continuous function of the radial position. The aperture–radius curves are rescaled using the scalings for both aperture and radius. In the toughness regime, both the rescaled aperture and radius collapse, resulting in overlapping profiles; see figure 8(b). In the viscous regime, the rescaled profiles collapse onto a single curve; see figure 9(b). The two data sets presented here are representative of all the experiments. In summary, the experimental observations establish the existence of two regimes and validate the respective scaling relations.
Figure 8. Fracture profile for the low-viscosity aqueous phase. (a) Fracture profiles at 
$t = 150,250,350$ s, with time increasing from clear to dark grey. (b) Rescaled fracture profiles using the toughness scaling laws in table 2. Experimental parameters: 
$E=30$ kPa, 
$\mu _0= 10$ mPa s, 
$Q_0 = 0.3$ ml min
$^{-1}$, 
$V_0 = 2.5$ ml, 
$\mu _{in}= 1$ mPa s, 
$Q_{in} = 1.8$ ml min
$^{-1}$.
Figure 9. Fracture profile for the high-viscosity aqueous phase. (a) Fracture profiles at 
$t = 5,10,40$ s, with time increasing from clear to dark grey. (b) Rescaled fracture profiles using the viscous scaling laws in table 2. Experimental parameters: 
$E=88$ kPa, 
$\mu _0= 10$ Pa s, 
$Q_0 = 10$ ml min
$^{-1}$, 
$V_0 = 8$ ml, 
$\mu _{in}= 8.6$ Pa s, 
$Q_{in} = 10$ ml min
$^{-1}$.
6. Conclusion
In this study, we model the properties of fractures driven by displacement flows by revisiting the theoretical framework established for single-fluid injections. We derive scaling relationships for the position of the interface between the two fluids, the radius and the aperture of the fracture in a brittle elastic matrix in the viscous-dominated and toughness-dominated regimes. We define a dimensionless time, which is equal to the ratio of the volume of displacing fluid injected over the volume of the pre-fracture. In the toughness regime, the propagation dynamics is independent of the fluid properties. As a result, the fracture dynamics for the displacement flow is the same as the dynamics for a single fluid, with the addition of initial finite volume 
$V_0$. In the viscous regime, however, the two fluids that fill the fracture contribute to the viscous dissipation. Over time, the relative volume of the two fluids changes. To describe the viscous dissipation in the fracture, we therefore define an average viscosity that accounts for the relative volume of displaced and displacing fluid. The scalings are compared to experimental results obtained by successively injecting an oil phase and an aqueous phase in a gelatin block. The experiments confirm the existence of two regimes of fracture propagation, and are in good agreement with the derived relationships.
This study focuses on displacement flows of immiscible fluids with comparable viscosities. Other types of displacement flows are common in fractures. For example, industrial applications involve the sequential injection of miscible aqueous fluids, some of which can be complex fluids such as suspensions of particles or polymer solutions. Recent and future efforts to characterize and model multiphase flows in fractures should ultimately support efficient hydraulic fracturing operations.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2022.954.
Acknowledgements
We would like to thank Z. Zheng for helpful discussions.
Funding
The work was supported in part by the ACS PRF (grant no. 60464-ND9).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data associated with this paper can be downloaded from the repository https://doi.org/10.25349/D9W32D.
