1. Introduction
A wide range of flow phenomena involves a moving interface, where one fluid displaces another in quasi-two-dimensional geometries. Such interfaces can become unstable to small fluctuations, leading to the development of complex patterns. Shear-dominated flows, such as those occurring in Hele-Shaw cells or porous media, can undergo viscous fingering or Saffman–Taylor instabilities in which characteristic fingering patterns emerge (Saffman & Taylor Reference Saffman and Taylor1958; Paterson Reference Paterson1981). By contrast, extensionally dominated flows, where interfacial traction is weak or absent, can give rise to tearing patterns (Sayag & Worster Reference Sayag and Worster2019a , Reference Sayag and Worsterb ; Sayag Reference Sayag2019; Ball & Balmforth Reference Ball and Balmforth2021). These tearing patterns arise in various natural or industrial phenomena. For example, a paste extruded axisymmetrically between parallel disks can develop tearing patterns at the unconfined edge (e.g. Mascia et al. Reference Mascia, Patel, Rough, Martin and Wilson2006). On a geophysical scale, unconfined ice shelves spreading into the ocean can form tearing patterns (Hughes Reference Hughes1983; Fricker et al. Reference Fricker, Young, Coleman, Bassis and Minster2005; Rignot et al. Reference Rignot, Jacobs, Mouginot and Scheuchl2013; Wearing, Kingslake & Worster Reference Wearing, Kingslake and Worster2020) that evolve into ice rifts, which enhance ice ablation and sea-level rise through calving and melting (Rignot et al. Reference Rignot, Jacobs, Mouginot and Scheuchl2013; Walker et al. Reference Walker, Bassis, Fricker and Czerwinski2013; Borstad, McGrath & Pope Reference Borstad, McGrath and Pope2017; Lipovsky Reference Lipovsky2018).
Numerical studies of extensional flows are typically focused on modelling the large-scale dynamics of spreading ice shelves and ice streams in Greenland or Antarctica. These include, for example, simulating the flow of an ice stream into a shelf (MacAyeal Reference MacAyeal1989), the evolution of an ice shelf over millennia (Wolff et al. Reference Wolff2025), predictions of whether a confined ice shelf can recover following a breakup (Akesson et al. Reference Akesson, Morlighem, Nilsson, Stranne and Jakobsson2022) and how strongly glaciers are buttressed by ice shelves (Greene et al. Reference Greene, Gardner, Schlegel and Fraser2022; Gudmundsson et al. Reference Gudmundsson, Barnes, Goldberg and Morlighem2023). While these models provide valuable insights into the evolution of specific ice shelves, the complex geometry and multiple interacting processes that they simulate make their results difficult to validate with laboratory experiments or theoretical analyses, which limits the development of physical insights. Furthermore, these numerical simulations are not designed to accurately resolve front propagation, hindering the exploration of rift formation and evolution. Resolving rift formation may be improved through recent advancements in numerical techniques that incorporate mechanical weaknesses in ice (e.g. Clayton et al. Reference Clayton, Duddu, Siegert and Martínez-Pañeda2022; Kachuck et al. Reference Kachuck, Whitcomb, Bassis, Martin and Price2022).
Radially symmetric flows provide more idealised configurations that are more amenable to theoretical and experimental methods. The patterns formed in such flows vary significantly depending on the mechanical properties of the spreading fluid. For Newtonian fluids, experiments have shown that the advancing front remains axisymmetric (Pegler & Worster Reference Pegler and Worster2012; Sayag, Pegler & Worster Reference Sayag, Pegler and Worster2012; Pegler & Worster Reference Pegler and Worster2013). In these studies, Newtonian fluids were discharged axisymmetrically into a denser, relatively inviscid ambient fluid, resulting in a flow with two distinct regions: an inner, shear-dominated circular region and an outer, extensionally dominated annular region. Both boundaries of the annular region acted as free contact lines, evolving continuously while maintaining axisymmetry. The inner contact line (grounding line) propagated along the base of the tank, while the outer contact line (front) floated. A thin-film model of this flow predicts similarity solutions, where the outer front advances as
$t^{1/2}$
at early times and transitions to linear growth at later times (Pegler & Worster Reference Pegler and Worster2012).
In contrast to Newtonian fluids, strain-rate softening polymeric solutions in similar settings develop complex tearing patterns at the spreading front (Sayag et al. Reference Sayag, Pegler and Worster2012; Sayag & Worster Reference Sayag and Worster2019a
; Ball & Balmforth Reference Ball and Balmforth2021). When the grounding line is pinned, these patterns consist of rectangular floating tongues separated by rifts, exhibiting
$k$
-fold rotational symmetries, where the integer wavenumber
$k$
decreases over time (Sayag & Worster Reference Sayag and Worster2019a
). A linear stability analysis reveals that extensional flows of uniform thickness become unstable for sufficiently strain-rate softening power-law fluids, while Newtonian and other less softening fluids remain stable (Sayag & Worster Reference Sayag and Worster2019b
). The instability arises when a circular fluid front under azimuthal extension is perturbed geometrically, causing localised flow convergence and outward radial pressure gradients for positive perturbations, and divergence with inward pressure gradients for negative perturbations. These radial forces amplify the initial perturbation amplitudes, thereby driving the instability (Sayag & Worster Reference Sayag and Worster2019b
; Sayag Reference Sayag2019). Similar tear-like patterns emerge also in viscoplastic fluids under extensional flow in similar settings (Ball & Balmforth Reference Ball and Balmforth2021). Specifically, considering a Herschel–Bulkley fluid and accounting for radial thinning during spreading, non-axisymmetric patterns can develop if the fluid is sufficiently strain-rate softening or has a yield stress. Within a linear stability framework, yield stress stabilises the flow at early times but amplifies the instability at later stages.
While the process that drives the instability is independent of the nature of the fluid, to linear order, the magnitude of the stabilising processes varies with the fluid’s rheology. Specifically, stabilising secondary flow vortices become increasingly weaker as the fluid becomes more strain-rate softening, resulting in the rise of instability when a certain fluid exponent is reached (Sayag & Worster Reference Sayag and Worster2019b ). However, it is not known whether the nonlinearity of the flow could contribute to the breakup of symmetry and the rise of instability in linearly stable flows, and particularly, Newtonian flows. Our aim in this work is to examine the nonlinear regime of Newtonian extensional flows. More specifically, we explore whether instabilities, which are suppressed in the linear stability theory for Newtonian fluids, could be driven by large-amplitude perturbations in the nonlinear regime and give rise to patterns reminiscent of rifts. Resolving the evolution in the nonlinear regime and the impact of large-amplitude perturbations requires solving the full nonlinear model numerically in two horizontal dimensions. For this, we use a quasi-three-dimensional (quasi-3-D) numerical simulation for radially spreading extensional flows in circular geometry, which solves the full nonlinear model for a fluid layer with non-uniform thickness and an evolving front.
We begin by introducing the theoretical model (§ 2). Then, we describe the numerical simulation and verify its consistency with known axisymmetric solutions (§ 3). Next, we examine the evolution of the front and stability of the flow due to small geometrical perturbations, which we contrast with linear stability theory, and due to large-amplitude perturbations to explore the nonlinear regime (§ 4). Finally, we analyse the stabilising and destabilising forces, and the potential impact of a shear-thinning rheology in the nonlinear regime (§ 5).
2. Theoretical model
We consider a buoyancy-driven viscous flow of an annular layer of a Newtonian fluid in the Stokes regime, where inertia is negligible (low Reynolds number). The fluid of density
$\rho _i$
is released axisymmetrically and at constant flux
$Q$
from the inner boundary of constant radius
$r_g$
, and spreads under gravity into the surrounding immiscible fluid layer of density
$\rho _w\gt \rho _i$
. Consequently, a fluid–fluid interface is formed whose front at radius
$r_N$
can vary with the position azimuth
$\theta$
and with time
$t$
(figure 1). We assume zero friction (free slip) along the entire fluid–fluid interface, implying that the flow of the intruding fluid is extensional. In the lubrication limit, we use the shallow-shelf approximation (SSA) to model the stress balance in the fluid layer with negligible basal friction, as in modelling of spreading ice shelves (Morland Reference Morland1987; MacAyeal Reference MacAyeal1989; Weis, Greve & Hutter Reference Weis, Greve and Hutter1999). Within the SSA, vertical shear stresses are negligible compared with horizontal shear and normal stresses. Therefore, depth integration of Stokes momentum balance, of the fluid layer with free top and bottom surfaces, leads to
\begin{equation} \boldsymbol{\nabla }(2\mu H \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u})+\boldsymbol{\nabla }\boldsymbol{\cdot }(2\mu H \kern-2pt \boldsymbol{E})=\rho _i g' H \boldsymbol{\nabla }H, \end{equation}
where
$H$
is the fluid thickness field,
$\boldsymbol{u}$
is the fluid horizontal velocity field,
$\mu$
is the constant dynamic viscosity,
$\boldsymbol{\nabla}$
is the horizontal (planar) gradient operator,
$g'\equiv g\,(1-\rho _i/\rho _w)$
is the reduced gravity with
$g$
the acceleration due to gravity and
$\boldsymbol{E}$
is the symmetric strain-rate tensor, whose components in cylindrical coordinates are
\begin{align} E_{\textit{rr}}=\frac {\partial u_r}{\partial r},\quad E_{r\theta }=\frac {1}{2} \left ( \frac {1}{r} \frac {\partial u_r}{\partial \theta } + \frac {\partial u_\theta }{\partial r} - \frac {u_\theta }{r} \right ),\quad E_{\theta \theta }=\frac {1}{r} \left ( \frac {\partial u_\theta }{\partial \theta } + u_r \right ), \end{align}
where
$u_r,u_\theta$
are the radial and azimuthal velocity components, respectively. The evolution equation for the fluid thickness, obtained by depth integration of the mass conservation equation, is
\begin{equation} \frac {\partial H}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(H \boldsymbol{u})=0. \end{equation}
We impose Dirichlet boundary conditions for the velocity and thickness at the inner boundary
\begin{equation} \boldsymbol{u}=\frac {Q}{2\pi r_gH_0 }\boldsymbol{\hat {r}},\quad H=H_0\quad \text{on}\ r=r_g, \end{equation}
where
$H_0$
is constant. Along the interface between the advancing and ambient fluids, the relatively inviscid ambient fluid applies zero tangential stress on the more viscous advancing front, while the normal stress along the front is set by the ambient fluid hydrostatic pressure. Therefore,
\begin{equation} \boldsymbol{\sigma } \boldsymbol{\cdot }\hat {\boldsymbol{n}} = -\rho _w g z \,\hat {\boldsymbol{n}} \quad \textrm{on}\quad r_N(\theta ,t), \end{equation}
where
$\boldsymbol{\sigma }=-p\boldsymbol{I}+2\mu \boldsymbol{E}$
is the full stress tensor of the viscous fluid,
$\boldsymbol{I}$
is the unit tensor and
\begin{align} \hat {\boldsymbol{n}}=\frac {\hat {\boldsymbol{r}}-\frac {r_N'}{r_N}\hat {\boldsymbol{\theta }}}{\sqrt {1+\left (\frac {r_N'}{r_N}\right )^2}} \end{align}
is the local normal to the instantaneous fluid–fluid interface, where prime represents a partial derivative with respect to
$\theta$
. The pressure distribution is obtained by integration of the vertical component of the momentum balance, which in the SSA has the form
\begin{equation} \frac {\partial \sigma _{zz}}{\partial z}=\rho _i g \end{equation}
between
$z$
and the top free surface
$s$
, together with the normal stress boundary condition along the viscous fluid free surface
$\sigma _{zz}(s)=0$
. This leads to
\begin{equation} p(z)=-(\tau _{\textit{rr}}+\tau _{\theta \theta })+\rho _ig(s-z), \end{equation}
where we used the fluid incompressibility to replace
$E_{zz}=-(E_{\textit{rr}}+E_{\theta \theta })$
, and where
$\tau _{\textit{rr}},\tau _{\theta \theta }$
are the deviatoric stress
$\boldsymbol{\tau }=2\mu \boldsymbol{E}$
normal components in the
$\hat {\boldsymbol{r}}\hat {\boldsymbol{r}}$
and
$\hat {\boldsymbol{\theta }}\hat {\boldsymbol{\theta }}$
directions, respectively, with the latter also familiar as the hoop stress (e.g. Pegler & Worster Reference Pegler and Worster2012). Substituting (2.8) into (2.5) and integrating over the whole fluid thickness results in
\begin{equation} (\tau _{\textit{rr}} + \tau _{\theta \theta }) \hat {\boldsymbol{n}} + 2\mu \boldsymbol{E} \boldsymbol{\cdot }\hat {\boldsymbol{n}} = \frac {1}{2} \rho _i g' H\hat {\boldsymbol{n}} \quad \textrm {on}\quad r_N(\theta ,t). \end{equation}
In the numerical simulation for the above-mentioned model, the propagation of the front
$r_N(\theta ,t)$
is implicit. That is, the bulk flow equations are solved in the entire numerical domain, which is larger than the domain containing the discharged volume, and the front emerges and evolves naturally. Therefore, we do not specify any particular condition to evolve the front (e.g.
$\dot {\boldsymbol{r}}_{\boldsymbol{N}}\boldsymbol{\cdot }\hat {\boldsymbol{n}}=\boldsymbol{u}(r_N)\boldsymbol{\cdot }\hat {\boldsymbol{n}}$
). The dynamic condition (2.9) is imposed along the fixed edge of the numerical domain, yet it is always satisfied to good accuracy along the curve
$r_N(\theta )$
that we identify as the front. We elaborate on these numerical aspects in § 3.1 and in Appendix A.
Figure 1. Diagram of the flow geometry, showing the viscous fluid (green patch) spreading into an ambient inviscid fluid (grey patch) in plan view (an
$r{-}\theta$
cross-section) and side view (an
$r{-}z$
cross-section). Gravity is in the
$-z$
direction.
3. Quasi-three-dimensional numerical simulation
We solve the SSA model described by (2.1), (2.3), (2.4), (2.9) numerically using Úa, a finite-element flow model (Gudmundsson Reference Gudmundsson2020). This computational framework includes two horizontal spatial dimensions and a variable thickness field, and we use a linear triangular unstructured mesh. The variational form of the governing equations is solved using a streamline-upwind Petrov–Galerkin approach, where the weak form of the SSA equations is integrated over each element, and the resulting system of equations is solved for the nodal values of the velocity and pressure. The mesh can be refined using adaptive mesh refinement strategies that increase resolution in areas of interest, such as regions with complex stress patterns, grounding lines or regions undergoing significant deformation. The time evolution of the fluid thickness and velocity is computed using an implicit time-stepping scheme. This approach ensures numerical stability, especially when solving nonlinear systems. The implicit scheme solves for the velocity and thickness at the next time step simultaneously, ensuring a consistent coupling between the two fields. Specifically, the time-stepping is based on a Crank–Nicholson method. At each time step, a nonlinear system of equations is formed due to the dependence of viscosity on the strain rate. This system is solved iteratively using a Newton–Raphson method, with convergence criteria based on the equation residual. The computational domain we use (figure 2) is in the range
$r_g \leqslant r \leqslant r_d$
and
$0\leqslant \theta \leqslant 2\pi$
, where
$r_g=1.5\,\mathrm{cm}$
and
$r_d=30\,\mathrm{cm}$
. We set the depth of the ambient fluid to
$100\,\mathrm{cm}$
uniformly throughout the domain, which is deep enough to ensure flotation of the viscous fluid over the ambient fluid with free bottom and top surfaces. The ambient fluid is essentially inviscid, having kinematic viscosity of
$10^{-4}\nu$
, where
$\nu \simeq 500\,\mathrm{cm^2\,s}^{-1}$
is the kinematic viscosity of the viscous fluid, which results in an extensionally dominated flow. Each simulation is initiated with a specified fluid–fluid interface at
$r_g\lt r_N(t=0,\theta )\ll r_d$
, where the initial thickness distribution of the viscous fluid
\begin{equation} H(t=0,r,\theta )= \begin{cases}H_0, & r_g\leqslant r \leqslant r_N, \\ 10^{-4}H_0, & r_N \lt r \leqslant r_d,\end{cases} \end{equation}
fills the entire domain. Assuming that
$H\ne 0$
everywhere at any time, we solve the equations on the entire domain without explicitly treating the front as a fluid interface. On the inner boundary, we specify constant thickness and flux of the viscous fluid, (2.4). The dynamic boundary conditions (2.9) are satisfied by construction of the numerical scheme on the outer boundary of the numerical domain,
$r_d$
. This implies that no boundary conditions are explicitly imposed on the fluid–fluid interface
$r_N$
, which evolves in the interior of the computational domain. This inconsistency is merely apparent because the space between the actual fluid front and the computational boundaries contains a film of fluid that is
$10^4$
times thinner, (3.1). Consequently, the impact of the dynamic boundary conditions imposed on the computational boundaries is communicated through the thin film with negligible modification. We verify this by finding no significant difference between simulations with identical initial conditions, but with various numerical domain geometries and sizes, as well as variable thin-film thicknesses. Moreover, we evaluate the explicit stress field at the moving fluid front and find that it satisfies the dynamic conditions (2.9) with negligible discrepancy (Appendix A).
Figure 2. (a) Snapshot of a part of the numerical domain, showing the inner boundary
$r_g$
near the origin (black), the outer boundary
$r_d$
(black), the front of the fluid layer (blue) and the mesh configuration (grey), with high element density in the vicinity of the front. (b) Radial distribution of the mesh spatial resolution (inverse area of elements) at
$t=6,15,24$
(solid) and the corresponding positions of the fluid fronts (dash). (c) Radial distribution of the normalised fluid thickness at
$t=4,22$
(solid turquoise, left axis) and the simultaneous time derivative of the thickness (solid grey, right axis) used to identify the instantaneous position of the fluid front (dashed, light blue).
We use an adaptive mesh that can be spatially non-uniform, keeping high resolution at the vicinity of the moving front (figure 2
a,b). The mesh is generated with linear triangulation, based on the Mesh2D, with the Delfront structure (Engwirda Reference Engwirda2014). A typical simulation consists of
${\sim}10^4$
nodes and elements initially, and grows up to
${\sim}10^6$
elements.
3.1. Resolving the evolving fluid front
Exploring the evolution of the fluid front, and particularly its potential to become unstable and develop complex patterns, requires an ability to trace it accurately. We accomplish this by defining the fluid front at each instant as the curve where the fluid thickness exhibits the largest rate of growth (figure 2
c). This definition coincides with the curve where the instantaneous magnitude of the horizontal flux divergence
$\boldsymbol{\nabla }\boldsymbol{\cdot }(H\boldsymbol{u})$
is maximal, (2.3). We expect that this maximum coincides with the front position, since both the fluid thickness and the flow normal to the front are finite, yet vary abruptly at the front, implying that both
$u\boldsymbol{\cdot }\boldsymbol{\nabla }H$
and
$H\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}$
become very large at the front. The capacity to accurately capture the position of the fluid front depends crucially on the spatial resolution. Particularly, low spatial resolution can lead to substantial smearing of the fluid front, thereby introducing greater discrepancy in the instantaneous position of the interface and consequently in its future evolution. To determine the necessary mesh resolution needed to accurately resolve the front, we increased the spatial resolution gradually until the front shape and position became insensitive to the mesh resolution. We chose to use a non-uniform mesh, consisting of lower spatial resolution away from the front and 3–4 orders of magnitude higher spatial resolution in the vicinity of the front (figure 2
b), to optimise the use of our computational resource. Maintaining such a mesh distribution requires a continuous tracing of the fluid front and simultaneous re-meshing of the domain such that the front remains at the centre of the high-resolution mesh regime (figure 2
b).
3.2. Validation with axisymmetric solutions
To validate our quasi-3-D numerical scheme and evaluate its accuracy, we examine its performance with respect to known axisymmetric solutions. In particular, we consider a Newtonian fluid layer that spreads axisymetrically from a vertical line source at the origin (
$r_g\rightarrow 0$
) over a denser, inviscid fluid. This system is described by the axisymmetric form of (2.1), (2.3), which have the dimensionless form
\begin{align} \frac {\partial }{\partial r}\left [H\left (2\frac {\partial u}{\partial r}+\frac {u}{r}\right )\right ]+H\frac {\partial }{\partial r}\left (\frac {u}{r}\right )&=\frac {H}{2}\frac {\partial H}{\partial r}, \\[-12pt] \nonumber \end{align}
\begin{align} \frac {\partial H}{\partial t}+\frac {1}{r}\frac {\partial }{\partial r}(rHu)&=0, \\[9pt] \nonumber \end{align}
together with the boundary conditions (2.4), (2.5),
\begin{align} {ru}=1,\quad H=1 &\quad \text{on}\ r\rightarrow 0, \\[-12pt] \nonumber \end{align}
\begin{align} 2\frac {\partial u}{\partial r}+\frac {u}{r}=\frac {H}{4} &\quad \text{on}\ r=r_N, \\[9pt] \nonumber \end{align}
in which
$r,t,H$
and
$u$
were respectively non-dimensionalised by
\begin{align} R,\quad T,\quad H_0\quad \text{and}\quad R/T, \end{align}
where the natural time and horizontal length scales are
\begin{align} T\equiv \frac {\nu }{g'H_0},\qquad R\equiv \left (\frac {\nu Q}{2\pi g'H_0^2}\right )^{1/2}. \end{align}
Such flow exhibits self-similar and steady-state solutions (Pegler & Worster Reference Pegler and Worster2012) that we describe briefly here. In the early stage of the flow (
$r_N\ll 1$
), the impact of buoyancy forces is substantially weaker compared with viscous forces. Neglecting buoyancy, the flow admits a similarity solution of the form
\begin{equation} \eta \equiv t^{-1/2} r, \quad H \equiv f(\eta ), \quad u \equiv t^{-1/2}s(\eta ), \end{equation}
where
$\eta$
is the early-time similarity coordinate,
$f(\eta )$
and
$s(\eta )$
are respectively the self-similar thickness and radial velocity. This form of solution implies that the early-time evolution of the front position is
$r_N\propto t^{1/2}$
(Pegler & Worster Reference Pegler and Worster2012).
At steady state, the axisymmetric SSA model (3.2) takes the form
\begin{align} -\textit{HH}^{\prime \prime }+\left (H^{\prime }\right )^2+\frac {\textit{HH}^{\prime }}{2 r}&=\frac {1}{4} r H^3 H^{\prime}, \\[-12pt] \nonumber \end{align}
\begin{align} \textit{ruH}&=1, \\[9pt] \nonumber \end{align}
where prime denotes a derivative with respect to
$r$
. The solution for the steady-state thickness far from the origin
\begin{align} H\sim \frac {\sqrt {6}}{r} \end{align}
describes the envelope that bounds the transient unsteady thickness solutions (Pegler & Worster Reference Pegler and Worster2012).
In the late stage of the flow, buoyancy forces become of the same order as viscous forces, leading to a different similarity solution of the form
\begin{equation} \xi \equiv t^{-1} r, \quad H \equiv t^{-1} f(\xi ), \quad u \equiv u(\xi ), \end{equation}
where
$\xi$
is the similarity coordinate, and
$f(\xi )$
and
$u(\xi )$
are respectively the dimensionless thickness and radial-velocity solutions in the similarity space. This form of solution implies that the late-time evolution of the front position is
$r_N\propto t$
(Pegler & Worster Reference Pegler and Worster2012).
To compare these theoretical axisymmetric solutions with our quasi-3-D simulation, we initiate the simulation with an axisymmetric ring of inner radius
$r_g$
, uniform thickness
$H_0$
and radial width
$\delta _0$
(figure 3
a). We find that the resulted flow, including the velocity and thickness fields and the front pattern, evolves highly axisymetrically (figure 3
b,c). Averaging the simulation results over
$\theta$
, we compare the radial distribution of each field with the axisymmetric prediction and find good agreement between them for both the steady and similarity solutions (figure 4).
Figure 3. Simulation results for an axisymmetric initial condition (base state), showing snapshots of the fluid velocity (arrows) and thickness (colour) distribution at times
$t=0,25,50$
. (a) Initial condition consisting of an axisymmetric ring of uniform thickness
$H_0$
and radial width
$\delta _0=0.3$
, and a slim fluid layer of thickness
$10^{-4} H_0$
(pale colour) that fills the rest of the numerical domain between
$r_N=r_g(1+\delta _0)$
and
$r_d=20 r_g$
. (b,c) Thickness, velocity and fluid front remaining axisymmetric throughout the evolution.
Figure 4. Validation of the quasi-3-D simulation with axisymmetric solutions. (a) Radial cross-section of the evolving fluid layer, showing results of the quasi-3-D simulation (solid, green) at
$t=1,10,20,30$
, compared with the axisymmetric solutions at early time (dash-dotted, red), with axisymmetric solutions at late times (dashed, orange) and with the steady-state solution (dashed, blue). The horizontal line (solid, dark) represents sea level. (b) Simulation results at
$t=20,30,40,53$
(solid, green) mapped to the similarity space, showing convergence towards the late-time similarity solution (solid, black). (c) Simulation results of the evolution of the fluid front
$r_N$
(solid) compared with the theoretically predicted similarity solutions at early time (dash-dotted, red) and late time (dashed, orange). (d) Simulation results of the radial velocity at
$t = 20, 30, 40, 53$
(solid, green) mapped to the similarity space, showing the process of convergence towards the late-time similarity solution (solid, black).
We note that in the axisymmetric model, the inner boundary is specified at the origin
$r_g=0$
to represent a line source, while in our simulation, the source is cylindrical with finite radius
$r_g\gt 0$
. This difference leads to the discrepancy in the vicinity of the origin in the thickness and in the early-time evolution of the front (figure 4
a,c). However, the details near the origin have no impact on the late-time evolution of the flow, which is consistent with the late-time similarity solution of the axisymmetric model and with the steady-state envelope towards which the late-time similarity solution converges (figure 4
a,c).
Another feature to note is the small kinks that form in each front at
$z=0$
(figure 4
a,b). These kinks are artefacts of the numerical method that we use, in which the viscous fluid fills the entire numerical domain. That is, the thickness of the viscous fluid varies smoothly from
${\sim} H_0$
in the shelf region to an order of
$10^{-4}H_0$
in the no-shelf region. Consequently, the free surfaces at the top and the base of the viscous fluid layer change continuously from the thick shelf region, across the shelf front, to the thin film at the no-shelf region, where they are significantly closer to sea level
$z=0$
. Therefore, looking at the vertical cross-section of the viscous fluid, which consists of its top and base free surfaces (figure 4
a,b), the viscous film at the no-shelf region appears as a thread of fluid emerging from the shelf front through a kink at sea level. Larger kinks in the late-time profiles are due to coarser spatial resolution chosen to sustain computationally long-time simulations and thereby resolve long shelves that approach the self-similar asymptotic profile. Due to limited computational resources, we restricted the domain size to
$r_d=20\, r_g$
at most, which appears insufficient to reach a late-enough time when a complete convergence to the similarity solution can be observed for the thickness and velocity. Nevertheless, their trends appear to be on the path of convergence to the late-time similarity solution (figure 4
b,d).
Our quasi-3-D simulations have demonstrated consistency with the predicted axisymmetric solutions for a range of dimensional parameters, including viscosity, flux, inner boundary radius and thickness. Moreover, our quasi-3-D results demonstrate axisymmetric patterns consistently with the predicted solutions and throughout the entire evolution, even though it is not constrained to evolve axisymmetrically (figure 3). Therefore, these results enhance the confidence in our numerical scheme and particularly in its capacity to accurately resolve the evolution of the leading front.
4. Flow evolution and stability under geometric perturbations
Having validated the performance of the numerical scheme with axisymmetric solutions, we now investigate the evolution of the flow following geometric perturbations of the front. In particular, we explore the stability of the flow, and the evolving patterns of the front and in the secondary flow.
4.1. Perturbing the axisymmetric base state
Figure 5. Perturbation parameters of the quasi-3-D numerical simulations that we performed, represented (a) in a table form and (b) in the
$k_0{-}A_0$
space, where each shape represents a different
$\delta _0$
.
We consider the axisymmetric solution presented in § 3.2 as the base-state solution and perturb it by adding a harmonic pattern to the axisymmetric front. Normalising the radius of the perturbed front by
$r_g$
, it has the dimensionless form
\begin{equation} r_N(\theta ,t=0)=1+\delta _0+A_0\cos (k_0\theta ), \end{equation}
where
$\delta _0$
is the dimensionless radial width of the initial (base-state) annular fluid layer,
$k_0$
is the perturbation azimuthal wavenumber and
$A_0=a_0\delta _0$
is the perturbation amplitude, which is a fraction
$a_0$
of
$\delta _0$
. We performed several simulations for a range of
$\delta _0,k_0$
and
$A_0$
that are
$10\,\%$
to
$75\,\%$
of
$\delta _0$
(figure 5). We traced the resulted thickness and flow distributions and the front shape (figure 6), and compared them with the patterns of the corresponding axisymmetric base-state. In all the simulations, we considered a Newtonian fluid injected through the inner boundary at
$r_g=1.5\,\mathrm{cm}$
, with constant mass flux
$Q=6.4\,\mathrm{cm}^3\,\mathrm{s}^{-1}$
and with initial uniform thickness
$H_0=2\,\mathrm{cm}$
.
Figure 6. Plan-view snapshots of the fluid thickness field (colour) at dimensionless times (left to right)
$t=0,10,20$
of two simulations with initially perturbed fronts, having perturbation amplitudes and wavenumbers (a)
$A_0=\delta _0/10$
,
$k_0=5$
, and (b)
$A_0=\delta _0/2$
,
$k_0=14$
.
4.2. Evolution of the perturbation amplitude and volume
A common feature among all simulations is that perturbations do not grow, though the
$k_0$
-fold rotational symmetry of the initial perturbation persists over time. In particular, two major evolution patterns emerge. When the perturbation amplitude and wavenumber are relatively small, the amplitude declines relatively fast (figure 6
a), whereas a perturbation with relatively large amplitude and wavenumber relaxes more slowly, and evolves into a pattern with pointy edges having a
$k_0$
-fold symmetry (figure 6
b).
We evaluate the evolution of the perturbation amplitude through the the ratio
$a(t)=A(t)/\delta (t)$
, where the instantaneous amplitude
$A(t)$
is the average radial distance between troughs and crests of the fluid front, and
$\delta (t)$
is the instantaneous width of a reference base-state simulation. We find that
$a(t)$
declines in a power-law form
\begin{align} a(t)\propto (t/T_p)^{-7/10} \end{align}
in all simulations (figure 7 a), where
\begin{align} T_p\equiv \frac {2\pi r_g^2H_0A_0(1+\delta _0)}{Q} \end{align}
is a time scale representing the ratio of the volume of a ring of fluid having a radius of the initial base-state front
$r_N(0)$
and a width of the perturbation amplitude
$A_0$
, to the flux
$Q$
. This decline implies that the perturbations relax and that the simulated flows are stable.
Figure 7. Evolution of (a) the perturbation relative amplitude
$a=A/\delta$
, and of (b) the volume ratio
$V_{{crest}}/V$
, showing power-law relaxation of the front perturbations. Regression analyses of the power-law relaxation intervals (black, dashed) yield the exponents −7/10 and −8/10 for the relative amplitude and volume ratio, respectively.
The relative amplitude of the perturbation
$a(t)$
reflects only part of the perturbation evolution because it does not account for variations in the fluid thickness. A more indicative quantity for the evolution of the perturbation is the fluid volume in a perturbation crest, given by
\begin{equation} V_{\textit{crests}}(t)=\int _{\theta _{\textit{crests}}}^{}\int _{1+\delta (t)}^{r_N(\theta ,t)} H(r,\theta ,t)\,r \,{\rm d}r \,{\rm d}\theta , \end{equation}
compared with the fluid volume of the corresponding base state
$V(t)$
. We find that
\begin{align} \frac {V_{\textit{crests}}}{V}\propto (t/T_p)^{-8/10} \end{align}
in all simulations (figure 7 b), which reinforces that the fluid front is stable and implies that the relaxation of the perturbation is faster when accounting also for the thickness variations.
4.3. Evolution of the kinetic energy
Although the amplitude and volume perturbations in the front region diminish over time (§ 4.2), a more thorough stability analysis should account for the evolution of perturbations in the entire fluid volume. This can be achieved by examining the evolution of the total kinetic energy
$E_k(t)$
, given by
\begin{equation} E_k(t)=\frac {\rho _i}{2}\int _{0}^{2\pi }\int _{r_g}^{r_N(\theta ,t)} \boldsymbol{u}^2(\theta ,r,t)H(r,\theta ,t) \,r \,{\rm d}r\,{\rm d}\theta , \end{equation}
compared with the simultaneous base-state kinetic energy
$E_{k}^{b}(t)$
, which has the same form as definition (4.6) only with the base-state quantities. Though
$E_{k}^{b}(t)$
grows in time, as more energy is pumped into the system through the entry flux
$Q$
than dissipated through viscous deformation, we expect that if
$\lim _{t\to \infty } E_k(t)/E_{k}^b(t)= 1$
, then perturbations decay and the flow is stable. Indeed, computing the ratio
$E_k(t)/E_{k}^b(t)$
, we find that it converges towards 1 for each simulation shown in figure 5, including in particular simulations with large perturbation amplitudes (figure 8). This implies that the energy of the perturbations decays over time in the entire fluid volume, reinforcing more substantially that the flow is globally stable.
Figure 8. Evolution of the kinetic energy ratio
$E_k/E_k^b$
for simulations with a wide range of perturbations (figure 5), indicating global stability.
4.4. Generation of super- and sub-harmonics
Even though the flow remains stable independently of the perturbation amplitude and wavenumber, the perturbed front can develop non-axisymmetric patterns. The resulted patterns vary with the perturbation amplitude and wavenumber through the emergence of super- and sub-harmonics of the fundamental perturbation mode, which we investigate by analysing the power spectrum of the evolving fluid front (figure 9).
Figure 9. (a) Front shape of a simulation with a perturbation wavenumber
$k_0=9$
and amplitude
$A_0=\delta /2$
at times
$t=0,15,75$
, and (b) the corresponding power spectrum at each instant normalised by the global maximum over the entire simulation time.
Each perturbation consists of a single mode
$k_0$
, which dominates the initial evolution of the front. As time progresses, additional harmonics of
$k_0$
gradually emerge. This includes super-harmonics of the mode
$k_0$
(integer multiplications of
$k_0$
), whose number at any given time grows with
$k_0 A_0$
, while their relative power diminishes monotonically with the wavenumber (figures 10 and 11
a). Consequently, as time progresses, the initial sinusoidal pattern of the front sharpens with
$k_0 A_0$
, forming cusp-like patterns while maintaining the initial
$k_0$
-fold rotational symmetry (e.g. figure 9
a). This implies that although the power of each super-harmonic grows with time, that of the initial mode remains sufficiently dominant throughout the evolution in determining the leading order shape of the front. The emergence of cusp-like patterns when large-amplitude perturbations are introduced could be a consequence of the SSA equations becoming hyperbolic. Specifically, the SSA force balance is elliptic in the velocity, while the thickness evolution equation is parabolic in the thickness. The large spatial gradients introduced by a large-amplitude perturbation may modify this normally elliptic–parabolic system, making it partially or fully hyperbolic, consequently admitting solutions with discontinuities in the derivatives such as cusps.
Figure 10. Evolution of the instantaneous power spectrum of the front shape, for simulations with varying perturbation wavenumbers and perturbation amplitude (specific values for each panel (bottom-left to top-right) are:
$k_0=4,4,5,8,8,9,18,18,18$
and
$A_0=0.04,0.1,0.15,0.05,0.1,0.15,0.04,0.075,0.15$
). Wavenumbers are normalised by
$k_0$
and the power spectrum amplitude is normalised by the global maximal amplitude.
Figure 11. (a) Number of super-harmonics at
$t=50,100,200$
(marker’s size grows with time) as a function of perturbation parameters
$A_0k_0$
. (b) Power spectrum evolution of a simulation with
$A_0=0.04$
,
$k_0=4$
showing the faster growth of a sub-harmonic mode (
$k/k_0=0.3$
) compared with the fundamental (
$k/k_0=1$
) and the first super-harmonic (
$k/k_0=2$
) modes.
On a longer time scale (dimensionless time greater than 300), we observe the generation of sub-harmonics of the perturbation wavenumber
$k_0$
(figure 10, left column). The power of these sub-harmonics grows in time faster than that of the fundamental mode and its super-harmonics (figure 11
b), which suggests that over time, the
$k_0$
-fold pattern of the front may undergo coarsening into a lower
$k$
-fold symmetry. Our present simulations do not reach the stage when coarsening becomes observable, though such coarsening occurs experimentally in similar flows with strain-rate softening fluids (Sayag & Worster Reference Sayag and Worster2019a
), and was implied theoretically, though inconsistently, in linear stability analysis (Sayag & Worster Reference Sayag and Worster2019b
).
4.5. The secondary flow field
The persistence of a
$k_0$
-fold rotational symmetry of the front pattern suggests that the secondary flow has a similar symmetry. In the case of a uniform fluid thickness (Sayag & Worster Reference Sayag and Worster2019b
), the secondary flow,
$\boldsymbol{u_1}$
, derived by a linear perturbation analysis, has a
$k_0$
-fold symmetry consisting of two types of streamline patterns, whose dominance depends on
$K(t)\equiv k_0 \delta (t)$
(figure 12
a). The dominating pattern at early time (lower
$K$
values) consists of open streamlines that are primarily azimuthal, diverging from the centres of frontal troughs and converging towards frontal crests (figure 12
a,
$K=1$
). As time progresses,
$\delta$
grows and a vortex pattern emerges near the inner boundary, having a similar
$k_0$
-fold symmetry as the early-time pattern, with their centres along the boundaries between frontal crests and troughs (figure 12
a,
$K=2.5$
). At later times (higher
$K$
values), the vortex pattern dominates the flow, while maintaining the same
$k_0$
-fold symmetry (figure 12
a,
$K=7$
). This vortex expansion reduces the perturbation outflow near the front crests and inflow near the front troughs, thereby enhancing the stability of the flow (Sayag & Worster Reference Sayag and Worster2019b
). The corresponding perturbation vorticity
$\boldsymbol{\nabla }\times \boldsymbol{u}_1$
is more intense at the front near the trough–crest boundaries and it also has a
$k_0$
-fold symmetry in phase with the streamline pattern (figure 12
a).
Figure 12. Comparison of the secondary flow streamline and vorticity (colour) patterns between the linear perturbation field
$\boldsymbol{u}_1$
of the uniform-thickness theory (a) calculated analytically (Sayag & Worster Reference Sayag and Worster2019b
), and the simulation secondary flow
$\boldsymbol{u}_{\textit{II}}$
of a non-uniform fluid thickness (b). All panels have
$k_0=3$
and the same colour scale, and columns correspond to
$K=1,2.5$
and
$7$
. The perturbation amplitude used in the simulation is
$A_0k_0=0.09$
and
$a_0=10\,\%$
.
Even though the simulations we consider involve non-uniform fluid thickness, we find similarities in the secondary flow with the linearised uniform-thickness theory. For the simulations, we define the secondary flow as the difference between the full field and the base flow:
\begin{equation} \boldsymbol{u}_{\textit{II}}\equiv \boldsymbol{u}-\boldsymbol{u_b}, \end{equation}
and examine the resulted streamline and vorticity patterns at relatively low and high perturbation amplitudes.
To compare the secondary flow fields with the linear theory, we considered the simulation with the same
$k_0$
and with relatively small perturbation amplitude of
$a_0=10\,\%$
(
$A_0k_0=0.09$
), at times corresponding to
$K=1,2.5$
and 7 (figure 12
b). To leading order, we find a remarkable qualitative agreement between the linear theory of uniform sheet thickness and the fully nonlinear secondary flow in the simulation, particularly at large
$K$
values (figure 12,
$K=7$
column). Similar to the linear case, vortices emerge from the inner boundary and gradually dominate the flow field as
$K$
(and therefore time) grows. This similarity includes both the front rotational symmetry and the flow orientation. Discrepancies between the two fields are more apparent at early time (figure 12,
$K=1$
column), where a streamline pattern similar to the early-time pattern of the linearised uniform-thickness theory appears to occupy only a thin margin close to the fluid front. In addition, the vortex pattern seems to emerge and dominate the flow earlier than in the linear calculation (figure 12,
$K=1,2.5$
). Though a full vortex is not clearly visual in those snapshots, we believe that this is a graphical artefact, as in fact the secondary velocities
$\boldsymbol{u}_{\textit{II}}$
are identically zero along the inner boundary (
$r=1$
). The vorticity field
$\boldsymbol{\nabla }\times \boldsymbol{u}_{\textit{II}}$
in the simulation also appears consistent with the linear prediction, to leading order. Consequently, we find that the linear analysis of a uniform-thickness layer (Sayag & Worster Reference Sayag and Worster2019b
) predicts surprisingly well the leading order of the full secondary flow of the simulations with non-uniform thickness.
Increasing the perturbation amplitude, new patterns emerge in the secondary flow that are not predicted by the linearised uniform-thickness theory. Specifically, at low
$A_0k_0$
, a
$k_0$
-fold pattern of vortices fills the entire space between the inner boundary and the front, and the intensity of those vortices declines with time (figure 13
a). However, as the magnitude of
$A_0k_0$
grows, the vortex structure in the inner region collapses and streamlines becomes more distorted and random, even though a
$k_0$
-fold vortex pattern is still maintained in a smaller annular region adjacent to the front (figure 13
b). The corresponding vorticity also grows quite dramatically with
$A_0k_0$
(figure 13
b).
Figure 13. Secondary-flow streamline and vorticity (colour) fields of the simulation, showing (a) an early (left) versus a late (right) instant of a low amplitude simulation (
$A_0k_0=0.09$
), and (b) a medium amplitude simulation (
$A_0k_0=0.36$
, left) versus a higher amplitude simulation (
$A_0k_0=1.35$
, right), at nearly the same instant and having the same wavenumber
$k_0=9$
. Vorticity signature corresponds to clockwise and counter-clockwise flows (blue and red colour, respectively).
5. Physical mechanism and the likelihood of instability
The successful prediction of the secondary flow in a fluid layer of non-uniform thickness by the linear theory for uniform thickness (figure 12) suggests that the stability mechanism predicted by that theory may also apply to the non-uniform flow in the simulations, particularly when the nonlinearity
$k_0A_0$
is weak. The linear stability theory for the case of uniform thickness (Sayag & Worster Reference Sayag and Worster2019b
) considers a power-law fluid having viscosity field of the form
\begin{align} \mu \propto E_{\textit{II}}^{1/n-1}, \end{align}
where
$E_{\textit{II}}=\sqrt {\boldsymbol{E:E}}$
is the second invariant of the rate-of-strain tensor, and the exponent
$n$
represents a Newtonian fluid (
$n=1$
) and shear-thinning or -thickening fluids (
$n\gt 1$
or
$n\lt 1$
, respectively). In that model, the perturbation radial force
\begin{align} -\frac {\partial p_1}{\partial r} = -\frac {1}{r} \frac {\partial \tau _{1r\theta }}{\partial \theta } -\left (\frac {\partial \tau _{\textit{1rr}}}{\partial r} + \frac {\tau _{\textit{1rr}} - \tau _{1\theta \theta }}{r}\right ), \end{align}
where
$p_1$
is the perturbation pressure field and the subscript 1 denotes a perturbation or secondary flow field, provides insights into the forces that stabilise or destabilise the front and to the potential impact of non-Newtonian rheology on rift formation.
The radial force due to the shear-stress divergence (first term on the right-hand side of (5.2)) enhances the perturbation radial force
$-\partial p_1/\partial r$
at the front, making it positive along perturbation crests and negative along perturbation troughs. This is evident when expressing that term using the boundary condition of no tangential stress along the base-flow front
$R_0$
\begin{align} \tau _{1r\theta } = \frac {2\tau _{0\theta \theta }}{R_0} \frac {\partial R_1}{\partial \theta } \qquad \text{on } R_0, \end{align}
where
$R_1(\theta )\propto e^{ik\theta }$
is the perturbation amplitude of the front, to get
\begin{align} -\frac {1}{r}\frac {\partial \tau _{1r\theta }}{\partial \theta } = -\frac {2\tau _{0\theta \theta }}{R_0^2}\frac {\partial ^2 R_1}{\partial \theta ^2} = \frac {2k^2}{R_0^2}\tau _{0\theta \theta }R_1 = \frac {4k^2}{R_0^4}R_1 \qquad \text{on } R_0, \end{align}
since the base-flow hoop stress
$\tau _{0\theta \theta }=2\mu _0E_{0\theta \theta }=2/R_0^2$
is extensional (subscript 0 denotes a base-flow quantity). Consequently, the perturbation radial force grows proportionally to
$R_1(\theta )$
due to this term, which drives the instability.
The radial force due to the normal-stress divergence (second term on the right-hand side of (5.2)) acts to attenuate the perturbation along the front. To see that, we substitute
$\tau _{\theta \theta }=-\tau _{\textit{rr}}$
, due to incompressibility in this approximation, which implies that
\begin{align} -\left (\frac {\partial \tau _{\textit{1rr}}}{\partial r}+\frac {\tau _{\textit{1rr}}-\tau _{1\theta \theta }}{r}\right ) = -\left (\frac {\partial \tau _{\textit{1rr}}}{\partial r}+\frac {2\tau _{\textit{1rr}}}{r}\right ) = -\frac {1}{n}\left (\frac {\partial }{\partial r}(2\mu _0E_{\textit{1rr}})+\frac {4\mu _0E_{\textit{1rr}}}{r}\right ). \end{align}
In the Newtonian limit
$n= 1$
, the base-state viscosity is the constant
$\mu _0=1$
and the perturbation problem has an exact closed-form solution (Sayag & Worster Reference Sayag and Worster2019b
), which suggests (Appendix B) that to leading order the normal-stress divergence equals to
\begin{align} -\left (\frac {\partial \tau _{\textit{1rr}}}{\partial r}+\frac {2\tau _{\textit{1rr}}}{r}\right ) \approx -\frac {\partial \tau _{\textit{1rr}}}{\partial r} \approx -\frac {8k^2}{R_0^4}R_1 \qquad \text{on } R_0, \end{align}
which grows proportionally to
$-R_1(\theta )$
and, consequently, stabilises the perturbation radial force.
Combining (5.4) and (5.6) into (5.2), we find that for a Newtonian fluid,
\begin{align} -\frac {\partial p_1}{\partial r} \approx -\frac {4k^2}{R_0^4}R_1 \end{align}
to leading order, making Newtonian flows stable regardless of
$k$
and
$R_0$
. The decline of the perturbation pressure gradient with the base-state radius
$R_0$
suggests that the destabilising force gradually becomes comparable to the stabilising force as time progresses (
$R_0=\sqrt {1+2t}$
). This implies that the flow may become marginally stable at late time. We find this approximation consistent with the exact theoretical predictions in the uniform thickness limit (figure 14).
Figure 14. Perturbation radial force
$-{\rm d}p_1/{\rm d}r$
at the front (normalised by
$k_0^2$
) variation with
$R_0$
and
$k_0$
– linear-stability theory versus simulation. Exact theoretical predictions shown for base-state radius
$1\leqslant R_0\leqslant 20$
and perturbation wavenumbers
$k_0=2,4,8,20,40$
(solid, coloured), and the approximated theoretical value by (5.7) (dashed, black). Simulation were initialised with
$\delta _0=0.3, a=10\,\%$
and
$k_0=3$
(
$A_0k_0=0.09$
, blue marker) and
$k_0=7$
(
$A_0k_0=0.21$
, green marker).
The theoretical solutions based on linear-stability analysis (Sayag & Worster Reference Sayag and Worster2019b
) provide a good prediction for the simulation results in the nonlinear regime when the perturbation amplitude is low. In particular, the perturbation pressure gradient in the simulations is stabilising along the advancing front and its magnitude declines as time progresses towards neutral stability (e.g.
$A_0k_0=0.09,0.21$
, figure 14). This consistency arises despite the varying fluid thickness in the simulation, while it is assumed uniform in the theoretical calculation. Moreover, the role of the stabilising and destabilising contributions to the perturbation radial force
$-{\rm d}p_1/{\rm d}r$
at the front, as described in (5.4), (5.6) and (5.2), persists into the interior of the fluid layer (figure 15). Specifically, the pattern of the radial force due to the shear-stress divergence is proportional to the perturbation pattern
$R_1$
(figure 15
a) and opposite to the contribution of the normal stress divergence (figure 15
b). The latter is dominant in magnitude, making the total perturbation radial force proportional to
$-R_1$
(figure 15
c), which stabilises the flow. The qualitative impact of the destabilising and stabilising contributions on the radial force remains consistent with the theoretical predictions even for relatively large perturbation amplitudes (e.g.
$A_0k_0=1.35$
, figure 16, top row), with correspondingly larger magnitudes. In particular, compared with a low perturbation simulation
$A_0k_0=0.21$
having the same base state (figure 16, bottom row), the full field pressure gradient becomes weaker when the perturbation amplitude is smaller. This implies that the stabilising force dominates under large-amplitude perturbation, while it is more closely balanced by the destabilising force at lower amplitude perturbations. Therefore, stronger nonlinearity tends to make the radial force more stabilising.
Figure 15. Simulation results of low perturbation amplitude (top row) compared with theoretical predictions (bottom row), showing the secondary (perturbation) fields of the (a) destabilising and (b) stabilising components of the radial force, (c) total radial force, and (d) of the second invariant of the strain-rate tensor. The magnitude of each field is shown in colour and the geometrical perturbation of the fronts in solid black. The simulation parameters are
$\delta _0=0.3,\, a_0=10\,\%,\, k_0=7$
and
$A_0k_0=0.21$
.
Figure 16. Simulation results of high perturbation amplitude
$a_0=50\,\%$
(top), along with lower perturbation amplitude
$a_0=10\,\%$
(bottom), showing the full fields of the (a) destabilising and (b) stabilising components of the radial force, (c) total radial force, and (d) the second invariant of the strain-rate tensor at the same time. The magnitude of each field is shown in colour and the geometrical perturbation of the fronts in solid black. The simulation parameters are
$\delta _0=0.3,\, a_0=50\,\%,\, k_0=9$
and
$A_0k_0=1.35$
(top), and
$\delta _0=0.3, a_0=10\,\%, k_0=7$
and
$A_0k_0=0.21$
(bottom).
An insight into the possible emergence of instability and the formation of rifts is provided by the second invariant of the strain rate tensor
$E_{\textit{II}}$
(figures 15
d and 16
d). Specifically, in both low and high perturbation amplitudes, the strain-rate second invariant is larger at perturbation troughs than at perturbation crests. Though insignificant for Newtonian fluids, this pattern suggests that shear-thinning fluids (
$n\gt 1$
) may have a smaller viscosity inside troughs and larger viscosity along crests (5.1). In addition, the effective stress, which is proportional to
$E_{\textit{II}}^{1/n}$
, is larger in troughs than it is in crests, particularly under large perturbations (figure 16
d). This combination of more mobile fluid together with larger stress field inside troughs increases the likelihood of rift formation in shear-thinning fluids. Moreover, the stabilising contribution to the radial force (5.5) would potentially become weaker as the fluid exponent
$n$
gets larger (i.e. becoming more shear thinning) due to the decline of the
$1/n$
coefficient on one hand and the bounded stress field in that limit on the other. This weakening may result in the stabilising contribution becoming less dominant as
$n$
grows and possibly smaller than the destabilising shear-stress divergence. Consequently, we expect that shear-thinning fluids having sufficiently large
$n\gt 1$
would develop radial force
$-{\rm d}p_1/{\rm d}r\propto R_1$
that would drive an instability and, potentially, the development of rifts.
6. Conclusions
Though there is no experimental evidence for instability of Newtonian flows spreading in circular geometry under frictionless conditions, predictions made by linear-stability analysis (Sayag & Worster Reference Sayag and Worster2019b ) suggest that the instability mechanism is active in such flows, only suppressed by stabilising processes. Therefore, our exploration here addresses the possibility that large-amplitude perturbations may lead to nonlinear effects that would drive an instability. We use a quasi-3-D numerical simulation of the shallow-shelf equations with non-uniform fluid thickness and an evolving fluid font, which we validate with theoretical predictions of an axisymmetric model. Thereafter, our exploration focuses on the evolution of non-axisymmetric initial conditions of harmonic structure and varying amplitude. In particular, we evaluate the stability of the flow and the evolving patterns, and explore the nonlinear effects on the physical mechanism.
Surprisingly, we find that such flows can be non-axisymmetric, even though they are globally stable, as indicated by the global kinetic energy. In particular, we reveal the emergence of cusp-like patterns in the fluid front and complex secondary flow patterns. The emerging patterns depend on the magnitude of the nonlinearity, described by the dimensionless quantity
$A_0k_0$
. Specifically, when
$A_0k_0$
is relatively small, the perturbations relax relatively fast in dimensionless time and the resulted patterns are similar to predictions made by linear stability analysis on similar flows with uniform thickness. These patterns include laminar eddies in the secondary flow and coarsening of the front shape through the growth of sub-harmonics. As the nonlinearity
$A_0k_0$
grows, the resulted patterns are strongly modified through the growth of cusps along the front, having
$k_0$
-fold rotational symmetry. The secondary flow in the large-amplitude regime involves intense and distorted eddies along the front, while the interior flow consists of disordered streamlines. Though there is no experimental evidence for non-axisymmetric patterns in extensional flows of Newtonian fluids, our numerical prediction of the emergence of cusp-like patterns may motivate conducting such experiments with fixed grounding lines and with large-amplitude disturbances induced at the front.
Despite the new non-axisymmetric patterns in the nonlinear flow regime, large-amplitude perturbations are insufficient to give rise to rift-like patterns in Newtonian fluids. Consistently with linear stability theory, our analysis demonstrates that stability persists primarily by the radial force due to normal-stress divergence, which dominates over the destabilising radial force due to shear-stress divergence. We expect the stabilising force to decline the more shear thinning the fluid is, thereby facilitating the emergence of instabilities and the likelihood of rift formation. This argument is strengthened by the distribution of the strain-rate and stress invariants, which suggest that a shear-thinning fluid could be more mobile and have a more intense stress field in frontal troughs than in frontal crests. This implies that the formation of ice rifts may strongly depend on the nonlinear rheology of ice. In Part 2 of this study, we focus on the impact of the fluid nonlinearity on the flow stability in the nonlinear regime, and on the formation and evolution of rift-like patterns.
Acknowledgements
We thank S. Pegler for useful discussions on the validation part and Y. Ashkenazy for assisting with computational resources.
Funding
L.S. is grateful to the Center for Integration in Science in the Ministry of Aliyah for financial support. This research was supported by the Israel Science Foundation (grant no. 3006/21).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Dynamic conditions along the moving fluid front
The dynamical boundary conditions (2.9) are imposed on the boundaries of the numerical domain. To explore how accurately these conditions are satisfied along the fluid front, we estimate the stress field along the resolved fluid front of a base-state simulation. Due to the perfectly circular front evolution
$r_N=r_N(t)$
, the stress balance along it should satisfy (2.9),
\begin{align} 2\tau _{\textit{rr}}+\tau _{\theta \theta }&=\frac {1}{2}\rho _i g'H, \\[-12pt] \nonumber \end{align}
\begin{align} \tau _{r\theta }&=0. \\[9pt] \nonumber \end{align}
We find that the mean of the simulation tangential and normal stress conditions at
$r_N$
are in a good agreement with the conditions (A1), with rather small standard deviation (figure 17). Therefore, even though the boundary conditions are not explicitly imposed on the fluid front, they are transmitted in good accuracy to the fluid front regardless the shape and distance of the computational boundary.
Appendix B. Analytic solution from linear stability analysis of a Newtonian fluid
In the Newtonian limit (
$n = 1$
), the governing equations for the perturbations resulting from linear stability analysis (Sayag & Worster Reference Sayag and Worster2019b
) are
\begin{align} 2\frac {\partial E_{\textit{1rr}}}{\partial r}+\frac {4}{r}E_{\textit{1rr}}+\frac {2}{r}\frac {\partial E_{1r\theta }}{\partial \theta }&=\frac {\partial p_1}{\partial r}, \\[-12pt] \nonumber \end{align}
\begin{align} \frac {2}{r}\frac {\partial E_{1\theta \theta }}{\partial \theta }+\frac {4}{r}E_{1r\theta }+2\frac {\partial E_{1r\theta }}{\partial r}&=\frac {1}{r}\frac {\partial p_1}{\partial \theta }, \\[-12pt] \nonumber \end{align}
\begin{align} E_{\textit{1rr}}+E_{1\theta \theta }&=0, \\[9pt] \nonumber\end{align}
representing respectively the force balance in the radial (5.2) and azimuthal directions, and continuity. In those equations, we substituted the base-state viscosity
$\mu _0=1$
and the perturbation strain-rate components given by (2.2) in terms of the perturbation velocity. The corresponding boundary conditions are
\begin{align} u_1=\frac {\partial u_1}{\partial r}=0&\qquad \text{ on } r=1, \\[-12pt] \nonumber \end{align}
\begin{align} E_{1r\theta }=\frac {2}{R_0^3}\frac {\partial R_1}{\partial \theta },\quad p_1=2E_{\textit{1rr}}+\frac {4}{R_0^3}R_1&\qquad \text{ on } r=R_0. \\[9pt] \nonumber\end{align}
Equation (B1) can be joined into a single fourth-order ordinary differential equation for the perturbation of the radial velocity component
\begin{align} \frac {1}{2k^2} r^4 u_1^{(4)} + \frac {3}{k^2} r^3 u_1^{(3)} + \left ( \! \frac {5}{2k^2} - 1 \! \right ) r^2 u_1^{(2)} - \left (\! 1 + \frac {1}{2k^2} \! \right ) r u_1^{(1)} + \left ( \! \frac {k^2}{2} + \frac {1}{2k^2} - 1 \! \right ) u_1 = 0. \end{align}
The general solution of this ordinary differential equation is
\begin{align} u_1(r) = c_1 r^{1+k} + c_2 r^{1-k} + c_3 r^{-1+k} + c_4 r^{-1-k}, \end{align}
with the four coefficients determined by the boundary conditions (B2), resulting in (Sayag & Worster Reference Sayag and Worster2019b )
\begin{align} c_1 &= \frac {k \left ( R_0^{2k} + k \left ( R_0^2 - 1 \right ) + 1 \right )} {k^2 R_0^{k+1} - 2 \left ( k^2 - 1 \right ) R_0^{k+3} + k^2 R_0^{k+5} + R_0^{3-k} + R_0^{3k+3}} R_1, \nonumber \\ c_2 &= \frac {k \left ( - \left ( k + 1 \right ) R_0^{2k} + k R_0^{2k+2} - 1 \right )} {k^2 R_0^{k+1} - 2 \left ( k^2 - 1 \right ) R_0^{k+3} + k^2 R_0^{k+5} + R_0^{3-k} + R_0^{3k+3}} R_1, \nonumber \\ c_3 &= -\frac {k \left ( k \left ( R_0^2 - 1 \right ) + R_0^2 \left ( R_0^{2k} + 1 \right ) \right )} {k^2 R_0^{k+1} - 2 \left ( k^2 - 1 \right ) R_0^{k+3} + k^2 R_0^{k+5} + R_0^{3-k} + R_0^{3k+3}} R_1, \nonumber \\ c_4 &= \frac {k \left ( k R_0^{2k} - \left ( k - 1 \right ) R_0^{2k+2} + R_0^2 \right )} {k^2 R_0^{k+1} - 2 \left ( k^2 - 1 \right ) R_0^{k+3} + k^2 R_0^{k+5} + R_0^{3-k} + R_0^{3k+3}} R_1. \end{align}
The stabilising term evaluated at the base-flow front
$R_0$
is
\begin{align} & -\left (\frac {\partial \tau _{\textit{1rr}}}{\partial r}+\frac {2\tau _{\textit{1rr}}}{r}\right ) \nonumber \\ & = -4k\frac {k^3 R_0^{2 k-2}-2 k \left (k^2-1\right ) R_0^{2 k}+k \left (k^2+2\right ) R_0^{2 k+2}+(2 k+1) R_0^{4 k}+2 k-1}{k^2 R_0^{2 k+2}-2 \left (k^2-1\right ) R_0^{2 k+4}+k^2 R_0^{2 k+6}+R_0^{4 k+4}+R_0^4}R_1 \end{align}
where
\begin{align} \tau _{\textit{1rr}}=2\mu _0E_{\textit{1rr}}. \end{align}
In the limit of
$k\gt 1$
and since
$R_0\gt 1$
, dominant balance implies
\begin{align} -\left (\frac {\partial \tau _{\textit{1rr}}}{\partial r}+\frac {2\tau _{\textit{1rr}}}{r}\right ) &\approx -\frac {8k^2R_0^{4k}}{R_0^{4k+4}}R_1 = -\frac {8k^2}{R_0^4}R_1 \quad \text{on } R_0, \end{align}
consistently with (5.6).
Open access
