Growth of a viscoplastic blister underneath an elastic sheet

Inspired by the formation of geological structures as earth's crust deforms by magmatic intrusions, we investigate the elastohydrodynamic growth of a viscoplastic blister under an elastic sheet. By combining experiments, scaling analysis and numerical simulations we reveal a new regime for the growth of the blister's height $\sim t^{5/9}$ and radius $\sim t^{2/9}$. A plug like flow inside the blister dictates its dynamics, whereas the blister takes a quasi-static self-similar shape given by a balance in the pressure gradient induced by bending of the elastic sheet and the fluid's yield stress.

The intrusion of a liquid underneath an elastic sheet resting on a non-deformable solid has been deployed as a generic model to understand the formation of geological structures such as magmatic sills and laccoliths [1][2][3], and ice sheet relaxation [4]. In models of e.g., sill and laccolith emplacement, the intruding magma is commonly assumed to be Newtonian, and the solid is assumed to deform by elastic bending and tensile fracturing. However, there is growing evidence that some magmas exhibit highly non-Newtonian behaviour [5][6][7], such that the magma can transport as plug flow [8]. Conversely, rock layers along which magma is emplaced can exhibit complex viscoplastic rheological response to magma flow [9,10]. To which extent these non-linear phenomena affect the dynamics of these geophysical systems is unclear.
In this context, the dual nature of a viscoplastic fluid, i.e. when the applied stress is above the yield stress τ 0 it flows like a liquid and below τ 0 it is a solid, is particularly interesting to capture the consequences of the fluid-and-solid nature of the flowing matter. We describe in here the spatiotemporal growth of a viscoplastic blister underneath an elastic sheet.
When an elastic sheet is resting on a solid pre-wetted by a fluid film, it has been shown that the height and radius of the formed blister is determined by local effects at the intrusion tip, which couple to the interior quasi-static solution of the shape of the sheet [11]. Other effects have also been demonstrated to affect the blister's dynamics, e.g. vapour at the intrusion tip [12], adhesion to the supporting solid [13][14][15], gravitational flow [16], the height of the pre-coated film layer [17] as well as stochastic effects [18]. The influence of a deformable elastic sheet on the resulting flow also has been shown in other settings such as the suppression of viscous fingering in a Hele-Shaw geometry [19] and the effect of the porous-base [4,20] .
Fluid models that account for viscoplasticity have been used to provide a more realistic representation of mud, glacier and magma flows on the Earth's surface [5,21,22]. Recently, viscoplastic fluids have also attracted much in- * acarlson@math.uio.no terest in interfacial flows, demonstrated to affect coating, drop spreading and coalescence [23][24][25].
FIG. 1: a) Schematic drawing of the studied system. An elastic sheet, of Young's modulus E = 0.126 MPa, thickness d = 10 mm, and Poisson's ratio ν = 0.5, rests on top of pre-wetted layer of yield stress fluid Carbopol, of thickness h 0 . Solutions of the same Carbopol as that of pre-wetted layer is injected at constant Q beneath the elastic sheet through a small tube of radius r t = 2 mm. The resulting blister is measured by tracking the laser-line and extracting the height profile, h(r, t). Experimental profiles h(r, t) of the sheet at t = [20, 340, 660] s when injecting a fluid with b) low (τ 0 = 2.6 P a) and c) high (τ 0 = 31.8 P a) yield stress.
We combine here the effects of having an elastic interface that can bend and look at the formation and growth of a blister as we inject a viscoplastic fluid. Intuitively, we might suspect that the large stress at the peeling front, as predicted from the Newtonian case [11], will locally fluidize the material, but how is the peeling front dynamics affected by the yield stress threshold? To answer this, we make experiments adapting the design from previous Newtonian elastohydrodynamics studies [11,26]. Fig.  1a shows a schematic representation of the experimental arXiv:2203.04193v1 [physics.flu-dyn] 8 Mar 2022 system, where we use Carbopol as a model yield-stress fluid as it has well characterised rheological properties and small thixotropic effects [27]. We prepare four different Carbopol solutions, giving a range in yield stresses τ 0 ∩ [2.6 − 47.2] Pa when fitting the Herschel-Bulkley model to the rheological measurements, see Supplemental Material (SM) for details.
A solid plate is coated by a thin layer of Carbopol with height h 0 ≈ [0.15 − 0.35] mm. The elastic sheet initially rests upon the Carbopol layer, before injecting Carbopol with a controlled flux Q∩[1.7, 3.3, 6.7]·10 −7 m 3 /s through the tube at the center of the plate. The sheet is made of a silicone based elastomer (Zhermack, Elite Double) with a Young's modulus E = 0.126 MPa. A laser line is used to visualise the sheet's deflections h(r, t) as the blister grows in time and space by using side view imaging.
The spatiotemporal growth of the blister for two Carbopol solutions, τ 0 = 2.6 Pa and τ 0 = 31.8 Pa, are shown in Fig. 1b-c, respectively. There are some noteworthy features to highlight in Fig. 1b. A distinct effect of τ 0 on the growth of the blister's radius R(t) and its height at the center h(0, t) is observed. The ratio between the deflection and the radius remains small throughout the experiments (h(0, t)/R(t) ≤ 0.1), and the flow is highly viscous with a small Reynolds number. Our measurements of the blister profiles show a high degree of axissymmetry around the tube where Carbopol is injected in the experiments.
To give a theoretical description of the elastohydrodynamics of the viscoplastic intrusion observed in Fig. 1b, we use the small aspect ratio and the viscous flow allowing us adopt the axisymmetric lubrication theory. There is a minor effect from gravity (illustrated in the SM), which we have neglected here for clarity, an assumption strictly valid for lengths below the elastogravity length L eg = B ρg 1/4 ≈ 0.035 m, with ρ the liquid density and g gravitational acceleration. The bending modulus of the sheet is fixed B = Ed 3 12(1−ν 2 ) = 0.014 N·m, with d the sheet thickness and ν = 0.5 the Poisson's ratio. We use the fact that there are only small deformations of the sheet and we can then neglect in-plane tension [28]. Thus the bending pressure is described as p(r, t) = B∇ 4 h(r, t) for small slopes. To capture the principal viscoplastic effect we use the Bingham model, which can be combined with the integrated lubrication equations to derive the expression for the vertical yielding line Y (r, t) of the fluid, where the shear stress τ rz equals the yield stress τ 0 [5,21,29], We expect the shear profile τ rz to be symmetric, and will therefore have another corresponding yielding position along z(r, t) = h(r, t) − Y (r, t). Between the two yielding lines we have a plug flow, which flows with a constant speed and zero shear rate [30]. We follow the standard procedure to derive the thin film equation, where an additional term including Y (r, t) appear from the viscoplastic Bingham model [5,21,29,31]: rY (r, t) 2 (3h(r, t) − 2Y (r, t)) ∂p(r, t) ∂r +w(r). (2) We have here introduced w(r) = 2Q for r ≤ r t , which is the Poisseuille flow through the tube with radius r t allowing us to prescribe the influx Q. Note that a no-slip condition has been applied at both surfaces and that µ is the effective fluid viscosity. Equation (2) has two particularly interesting limits. If Y (r, t) = h(r, t) Newtonian flow is recovered and predicts the blister to grow with a height ∼ t 8/22 [11]. As Y (r, t) → 0 we expect a plug flow to form and the flow to be dominated by the effects from τ 0 . Before moving on to the experiments and the numerical simulations of equation (2), we derive a scaling law for the spatiotemporal blister growth. If we scale the height h(r, t) with h(0, t) and the radial direction r with R(t), while assuming unidirectional flow from r = 0 and outwards, in the limit Y (r, t) → 0, equation (1)  and r max = R(t)T 2/9 Lt 2/9 . h(r,t)∂p(r,t) 2∂r ∼ Bh 2 (0,t)
We make equation (1) and equation (2) nondimensional by scaling h(r, t) and r with L and time t with T , giving only one non-dimensional number appearing in front of the first term on the right-hand-side of equation (2) i.e. Π = B 6µQ , which gives the ratio between elastic and viscous forces. We assume the effective viscosity in our experiments to be set by the shear rate in the regions where the fluid is yielded, because the flows here are determining our viscous timescale. We estimate the shear rate by assuming a parabolic velocity profile for z > h(r, t) − Y (r, t) and z < Y (r, t), and obtainγ max ≈ ∂u ∂z | z=0,z=h ∼ Q R(t)Y (r,t) 2 . Setting Q ≈ 10 −7 m 3 /s, R(t) ≈ 10 −1 m, and Y (r, t) ≈ 10 −4 m givesγ max ≈ 100 s −1 . This is a conservative estimate for the shear rate and we therefore assumeγ ≈ We use the blister's maximal height h(0, t) as a metric to characterise the dynamics to test the predictions from equation (3) and equation (4). The bending stiffness B is fixed, while we vary the yield stress τ 0 and the flux Q in the experiments shown in Fig. 2a. It is clear that h(0, t) is affected by both parameters by steepening of the blister profile. With these experiments at hand, the experimental data can be rescaled by h(0, t)/L and t/T , as shown in Fig. 2b, and collapse onto a single curve following the proposed scaling law of equation (3). At late times in the experiments, gravity is starting to affect the dynamics and we enter a transitional regime where bending and gravitational forces become comparable.
The experiments are complemented by numerical solutions of equation (2), which is discretised by linear finite elements and the solution method is described in [17]. We use an adaptive time-stepping routine, with a time step limit of ∆t/T = 9 · 10 −6 , and a discretisation in space ∆r/L ∈ [0.0005 − 0.005]. We have defined Y (r, t) = max , h(r, t)/2 − 1 | ∂ ∂r ∇ 4 h(r,t)| , with a regularisation parameter = 10 −6 similar to [24], and we tested that the results are unaffected by this choice of . The simulations are started with a pre-wetted layer h 0 /L = 0.0069, and with Y (r, t) = 0, which is gradually introduced until t/T = 10 −4 . The results from the numerical simulations are shown with the experiments in Fig. 2b, they are fairly insensitive to the non-dimensional number Π, and recover the scaling law in equation (3). The small time shift separating the experiments and the simulations (a factor of 2.5 in time) may stem from the minimal rheological model we use to describe Carbopol [27], wall slip or the simplified influx condition. Regardless, the simulations overall represent well the elastohydrodynamic growth of the viscoplastic blister.
Next, we want to understand the interface dynamics of the blister and if it adopts a self-similar shape. Fig.  3a shows the blister shape in time from the experiments. To better understand the viscoplastic flow in the blister, we plot in Fig. 3b the non-dimensional height profile and the yielded lines for two time instances from a numerical simulation. It is clear that a major part of the area is encompassed between the yielding lines (dashed lines) in Fig. 3b with a significant plug flow, whereas a flow with a viscous gradient is limited to thin areas near the two bounding solid surfaces. This further supports our argument allowing us to derive equation (3) and equation (4). To go one step further, we now rescale the height profiles from the experiments and the numerical simulations data (eleven profiles are shown in-between the two timepoints in Fig. 3b, using the predictions from equation (3) and equation (4), where we normalise each profile with: Lt 5/9 and radial value, r max = R(t)T 2/9 Lt 2/9 , see SM. The experimental and numerical data fall onto a profile for a universal blister shape, see Fig. 3c.
The self-similar blister shape can be understood if we look back to equation (1) in the limit Y (r, t) → 0 that reduces to a non-linear ordinary differential equation (ODE): h(r)/d∇ 4 h(r)/dr = 2τ 0 sgn(h 0 −h) [24]. By scaling h(r) and r with L, equation (1) becomes parameter free. To determine the blister shape we use a numerical profile from our simulation, Π = 1000 and t/T = 0.75, as an initial guess. We need five boundary conditions to define our problem. We choose three conditions at the centerline, where we pin the height to the maximum height of the simulation profile, h max , h(0) = h max , and set two symmetry conditions, ∂h(0) ∂r = 0 and ∂ 3 h(0) ∂r 3 = 0. We supplement with two boundary conditions atR = 1.3 × R, where we set h(R) = h 0 and ∂ 2 h(R) ∂r 2 = 0. The solver Solve bvp retrieved from the scipy package in Python is used to solve the ODE. In Fig. 3c, the solution for the height profile is shown together with the experiments and the numerical simulations of equation (2), which gives a perfect agreement for the predicted shapes.
Intrusion of flowing matter with complex rheological properties underneath an elastic layer or inside an elastic matrix appear in a myriad of geological processes. Both non-linear magma plug flow [8] and viscoplastic deformation of magma intrusion's host rocks [10] have been documented to play a key role in magma emplacement in the Earth's crust. Our study quantifies the dynamics of a viscoplastic blister growth underneath an elastic sheet, and evidences a viscoplastic regime for the growth of its height and radial expansion. The blister is shown to adopt a self-similar quasi-static shape in time that is set by the balance between the pressure gradient induced by the bending of the elastic sheet and the yield stress of the intruding fluid. These results suggest that the growth of geological structures such as sills and laccoliths exhibit distinct dynamics than those predicted from a Newtonian model if the intruding fluid can support a critical stress before starting to flow. Our work gives a first look at the elastohydrodynamics of a growing viscoplastic blister, but highlight that there are many more questions in these configurations that remain to be understood such as additional non-Newtonian effects in the bulk and the elastic sheet [32] and to combine the yield-stress intrusion model with an elastic fracture toughness at the tip.