1 Introduction
In his lucid 2015 book Questions About Elastic Waves Engelbrecht (2015), Engelbrecht asks “What is a wave?” and answers “As surprising as it may sound, there is no simple answer to this question.” Indeed, the definition of ‘wave’ depends on the physical context at hand Christov (2014). Although most wave phenomena in classical continuum mechanics are described by hyperbolic (wave) equations, one of the surprises of 20^{th} century research into nonlinear partial differential equations (PDEs) is that certain parabolic (diffusion) equations also yield structures with finite speed of propagation. Two examples are (i) a linear diffusion equation with a nonlinear reaction term Kolmogorov et al (1991); Fisher (1937), and (ii) a diffusion equation that is nonlinear due to a concentrationdependent diffusivity Barenblatt (1952).^{1}^{1}1More specifically, Barenblatt (1952) (see also (Ostriker et al, 1992, p. 13)) credits the observation of finitespeed of propagation in a nonlinear diffusion equation to a difficulttofind 1950 paper by Zeldovich and Kompaneets. Indeed, it is known that certain aspects of wave phenomena can be reduced to a problem of solving a parabolic PDE, as gracefully illustrated by Engelbrecht (Engelbrecht, 1997, Ch. 6) through a series of selected case studies; further examples include, but are not limited to: electromagnetic waves propagating along the earth’s surface Vlasov and Talanov (1995), seismic waves Jerzak et al (2002), underwater acoustics Jensen et al (2011), and the classical theory of nerve pulses (Engelbrecht, 1997, §6.4.2), which nowadays has been updated by Engelbrecht et al (2018a) to a nonlinear hyperbolic (wave) model in the spirit of the Boussinesq paradigm Engelbrecht et al (2018b); Christov et al (2007).
Of special interest to the present discussion are physical problems that are modeled by nonlinear parabolic PDEs. These nonlinear problems lack general, allencompassing solution methodologies. Instead, finding a solution often involves methods that are specific to the nature of the governing equation or the physical problem that it describes (Evans, 2010, Ch. 4) (see also the discussion in Christov (2014) in the context of heat conduction). The classical examples of nonlinear parabolic PDEs admitting traveling wave solutions come from heat conduction (Zel’dovich and Raizer, 2002, Ch. X) (see also Straughan (2011)) and thermoelasticity Straughan (2011); Berezovski and Ván (2017). The sense in which these nonlinear parabolic PDEs admit travelingwave and ‘wavefront’ solutions now rests upon solid mathematical foundations Gilding and Kersner (2004); Vázquez (2007), including the case of gradientdependent nonlinearity Tedeev and Vespri (2015) (e.g., the last case in table 1 to be discussed below).
A classical example of a nonlinear parabolic PDE governing the finitespeed wavelike motion of a substance arises in the study of an ideal gas spreading in a uniform porous medium Barenblatt (1952). A similar nonlinear parabolic equation was derived for the interface between a viscous fluid spreading horizontally underneath another fluid of lower density ( between the fluids) Huppert (1982). The motion of the denser fluid is dictated by a balance of buoyancy and viscous forces at a low Reynolds number (viscous forces dominate inertial forces). Such viscous gravity current flows are characterized by ‘slender’ fluid profiles i.e., they have small aspect ratios (, where and are typical vertical and horizontal length scales, respectively). Therefore, these flows can be modeled by lubrication theory (Leal, 2007, Ch. 6). Generically, one obtains a nonlinear parabolic equation for the gravity current’s shape as a function of the flowwise coordinate and time . The case of the spreading of a fixed mass of Newtonian fluid was originally explored contemporaneously by Didden and Maxworthy (1982) and Huppert (1982).
Being governed by a parabolic (irreversible) equation, these currents ‘forget’ their initial conditions after some time has elapsed; this is Barenblatt’s concept of intermediate asymptotics Barenblatt and Zel’dovich (1972); Barenblatt (1996). Moreover, the PDE (1
) can be reduced to an ordinary differential equation (ODE) through a selfsimilarity transformation. If the similarity variable can obtained by a scaling (dimensional) analysis, then the solution is termed a selfsimilar solution of the
first kind (Barenblatt, 1996, Ch. 3). Specifically, the transformation is ( in units^{2}^{2}2Throughout the chapter, we use SI units for all dimensional quantities. of meters), where is the similarity variable (dimensionless), is the selfsimilar profile to be determined by solving an ODE, and and are dimensional consistency constants. The exponents and are obtained through scaling (dimensional analysis) of the governing PDE. As a representative example, consider the onedimensional (1D) spreading of a fixed mass of fluid having an arbitrary ‘wavy’ initial shape, as shown in figure 3. Suppose the fluid’s shape is governed by the linear diffusion equation (taking m^{2}/s in this example without loss of generality) subject to . The initial condition (IC) is quickly ‘forgotten,’ and the ultimate asymptotic state (here, flat) is achieved after passing through an intermediate asymptotic regime. It is straightforward to determine the selfsimilarity transformation: and (see, e.g., Barenblatt (1996) and the Appendix). Here, depends on the initial condition and . The convergence of the rescaled profiles towards can be clearly observed in figure 3(b). The IC is forgotten, and the profile converges onto the Gaussian intermediate asymptotic shape Barenblatt (1996). The profile is termed ‘universal’ because it is independent of .Having illustrated the notion of firstkind selfsimilarity as intermediate asymptotics, let us summarize its use in studying the gravitational spreading of Newtonian viscous fluids in a variety of physical scenarios. For example, gravity currents arise in geophysical applications associated with flows through porous rocks Woods (2015) such as in ground water extraction Bear (1988), during oil recovery Huppert (2000); Felisa et al (2018), and during CO_{2} sequestration Huppert and Neufeld (2014). In these examples, represents an interface between two immiscible fluids in the limit of large Bond number (gravity dominates surface tension). There is now an extensive literature featuring a wealth of exact and approximate analytical selfsimilar solutions for gravity currents in porous media, e.g., Barenblatt (1952); Huppert and Woods (1995); Anderson et al (2003); Lyle et al (2005); Vella and Huppert (2006); Hesse et al (2007); Anderson et al (2010); De Loubens and Ramakrishnan (2011); Ciriello et al (2013); Huppert et al (2013); Zheng et al (2014) amongst many others.
In this chapter, we focus on the propagation of nonNewtonian gravity currents, specifically ones for which the denser fluid obeys a powerlaw rheology. This tractable model of nonNewtonian rheological response is also known as the Oswald–de Weale fluid Bird et al (1987). In unidirectional flow, the powerlaw model simply dictates that fluid’s viscosity depends upon a power of the velocity gradient. Di Federico et al (2006) generalized Huppert’s problem Huppert (1982) to powerlaw fluids, although Gratton et al Gratton et al (1999); Perazzo and Gratton (2003) had also considered some related problems. Even earlier, Kondic et al (1996, 1998) derived the governing equations for powerlaw fluids under confinement (i.e., in HeleShaw cells) using the lubrication approximation. These works have contributed to the use of a modified Darcy law to model the flow of nonNewtonian fluids in porous media using the analogy to flow in HeleShaw cells. Aronsson and Janfalk (1992) were perhaps the first to combine a Darcy law for a powerlaw fluid with the continuity equation to obtain a single PDE, of the kind studied herein, governing the gravity current’s shape. Recently, Lauriola et al (2018) highlighted the versatility of this approach by reviewing the existing literature and extending it to twodimensional axisymmetric spreading in media with uniform porosity but variable permeability. All these flows are of interest because exact analytical selfsimilar solutions in closed form have been derived previously Gratton et al (1999); Perazzo and Gratton (2005); Di Federico et al (2012, 2017); Ciriello et al (2016). Specifically, the solution of Ciriello et al (2016) will be used in §4.1 below to verify the truncation error of the proposed numerical method.
For a selfsimilar solution to exist, both the governing PDE and its boundary conditions (BCs) must properly transform into an ODE in with suitable BCs. A number of studies have specifically shown that the volume of fluid within the domain can be transient, varying as a power law in time, (), and a selfsimilar solution still exists (see, e.g., Barenblatt (1952); Lyle et al (2005); Hesse et al (2007); Di Federico et al (2012); Zheng et al (2014); Di Federico et al (2017) and the references therein). However, the nonlinear ODE in often cannot be integrated exactly in terms of known function, except for . In §2.2 below, we discuss how a constraint of the form can be implemented numerically through flux BCs at the computational domain’s ends.
With increasing complexity of the flow physics incorporated in the model, finding a selfsimilarity transformation may no longer be possible simply by scaling (dimensional) arguments. Gratton and Minotti (1990)classified a number of such situations, including the socalled ‘focusing’ flows involving fluid axisymmetrically flowing towards the origin on a flat planar surface. Further examples involving confined currents in channels with variable width, and/or in porous media whose permeability and porosity are functions of , were proposed by Zheng et al (2014), as illustrated in figure 4. These gravity currents do enter a selfsimilar regime, even though a selfsimilar transformation cannot be obtained by scaling arguments alone. The exponents and in the transformation are unknown a priori, hence this situation represents a selfsimilarity of the second kind (Barenblatt, 1996, Ch. 4)
. The governing equation can be transformed to an ODE, following which a nonlinear eigenvalue problem must be solved for
and through a phaseplane analysis Gratton and Minotti (1990); Gratton (1991). Alternatively, experiments or numerical simulations are necessary to determine and . For example, early numerical simulations were performed to this end by Diez et al (1992). However, a ‘prewetting film’ ahead of the current’s sharp wavefront ( where ) was required to avoid numerical instabilities. The scheme therein was also firstorder accurate in time only. In this chapter, we propose a modern, highorderaccurate implicit numerical method for use in such problems.Specifically, we develop and benchmark a strongly implicit and conservative numerical scheme for 1D nonlinear parabolic PDEs arising in the study of gravity currents. We show how the proposed scheme can be used to simulate (with high accuracy and at low computational expense) the spreading of 1D nonNewtonian viscous gravity currents in variable geometries (specifically, HeleShaw cells with widths varying as a power law in ). To this end, we build upon the work of Zheng et al (2014), which introduced this type of finitedifference scheme for simulating the spreading of a finite mass of Newtonian fluid in a variablewidth HeleShaw cell. Owing to its accuracy and stability, this finitedifference scheme has been recently applied by Alhashim and Koch (2018) to study hydraulic fracturing of lowpermeability rock.
This chapter is organized as follows. In §2, we briefly summarize existing models describing certain flows of viscous gravity currents. Then, we introduce a convenient general notation for such nonlinear parabolic PDEs. In §3.1, we introduce the 1D equispaced but staggered grid upon which the proposed finitedifference scheme is to be implemented. The derivation of the BCs for the PDE, from the mass conservation constraint, is discussed in §2.2. Then, we construct the nonlinear Crank–Nicolson scheme in §3.2 and discuss the discretized form of the nonlinear flux BCs in §3.4. Continuing, in §4.1, the scheme’s accuracy is justified by comparing the numerical solution provided by the finitedifference scheme (up to a specified physical time) against an analytical solution obtained through a selfsimilar transformation of the PDE. Specifically, this approach involves three validation cases: (i) a symmetric (about ) lump of fixed fluid mass spreading in two directions (convergence is independent of BCs), (ii) a fixed fluid mass spreading away from the origin () (requires only no flux BCs), and (iii) a variable fluid mass injected at the origin spreading away from it (requires careful implementation of the nonlinear BCs). In all three cases, the scheme is shown to be capable of accurately computing the evolution of gravity current’s shape. In §4.2, we analyze the scheme’s conservation properties by verifying numerically that it respects the mass constraint . We consider two validation cases: (i) release of a fixed fluid mass (), and (ii) fluid mass injection into the domain (). In both cases, we specifically focus on the challenging case of a nonNewtonian (powerlaw) displacing fluid in a variablewidth channel. As a benchmark, we use previously derived firstkind selfsimilar solutions from the literature, which are discussed in the Appendix.
2 Preliminaries
In this section, we summarize the mathematical model for viscous gravity currents in a selected set of applications involving Newtonian and nonNewtonian fluids. We study their spreading in a fixed or variablewidth channel geometry (also known as a “HeleShaw cell”), as well as flows in heterogeneous porous media with independently variable permeability and porosity. Our goal is to highlight the fact that all these models can be concisely summarized by a single nonlinear parabolic PDE supplemented with a set of nonlinear Neumann (flux) BCs.
2.1 Fluid Domain and Flow Characteristics
The flow domain is assumed to be long and thin. For example, it can be a channel existing in the gap between two impermeable plates, i.e., a HeleShaw (HS) cell, which may or may not have variable transverse (to the flow) width as shown in figure 4(a); or, it can be slab of uniformthickness heterogeneously porous material, as shown in figure 4(c). The viscous gravity current consists of one fluid displacing another immiscible fluid. Therefore, a sharp interface separates the two fluids at all times. The present study considers the limit of negligible surface (interfacial) tension (compared to gravitational forces). The density difference between the two fluids is large compared to the density of the lighter fluid, and the denser fluid flows along the bottom of the cell, which is a horizontal impermeable surface. In doing so, the denser fluid displaces the lighter fluid out of its way. Here, the geometry is considered to be vertically unconfined so that the details of the flow of the upper (lighter) fluid can be neglected.
We are interested in the evolution of the interface between the two fluids. Owing to the vertically unconfined, long and thin geometry of the flow passage, the denser fluid has a slender profile (small aspect ratio), and the fluid flow can be described by lubrication theory. The lubrication approximation also requires that viscous forces dominate inertial forces; this is the limit of small Reynolds number. In this regime of small Reynolds number but large Bond number, the flow is governed by a balance of viscous forces and gravity. Furthermore, the lubrication approximation allows for (at the leading order in the aspect ratio) the variation of quantities across the transverse direction, as well as the vertical velocities of the fluids to be neglected.
As shown in figure 4(a), for the flow in a HS cell, we allow the cell’s width to vary as a powerlaw of the streamwise coordinate , i.e., , where is a dimensionless exponent, and is a dimensional consistency constant having units m. Since the cell has a variable width, it originates from a cell ‘origin,’ which is always taken to be such that . As discussed in Zheng et al (2014), in such a flow geometry, the lubrication approximation may fail when is an increasing function of i.e., . In such quicklywidening cells, the transverse variations of properties become significant. We ensure the validity of the lubrication approximation, and models derived on the basis of it, by only considering such that remains a decreasing function of .
The porosity can also be varied by filling the HS cell with beads of fixed diameter, as illustrated in figure 4(b). We also consider a gravity current spreading horizontally in a porous slab of constant transverse width () with heterogeneous porosity and permeability , as shown in figure 4(c). Here, are dimensionless exponents and are dimensional constants needed for consistency with the definitions of porosity and permeability, respectively; specifically has units of units of m, and has units of m. These variations are illustrated by the streamwise changes of bead radii in figure 4(c). Now, the point at which the porosity and permeability vanish is the origin of the cell. Another interesting case, that of a medium with vertically heterogeneous porosity, has been explored by Ciriello et al (2016). In this chapter, we limit our discussion to flow in a completely porous (i.e., unobstructed, ) HS cell of variable width as in figure 4(a). However, the numerical scheme developed herein can readily treat any of these cases, taking the appropriate parameter definitions from table 1 in §2.2.
We allow the denser fluid to be nonNewtonian. Specifically, it obeys the powerlaw rheology. In unidirectional flow, the one unique nontrivial shear stress component is given by , where the dynamic viscosity depends on the shear rate as . Here, is the flow consistency index (units of Pas), and () is the fluid’s rheological index. Fluids having are termed shearthinning (e.g., blood), and fluids with are termed shearthickening (e.g., dense particulate suspensions). In the special case , the powerlaw model reduces to the Newtonian fluid. As stated above, the flow of the displaced fluid is immaterial to the dynamics of the gravity current, as long as the viscosity and density contrasts are large. This condition is satisfied, e.g., by assuming (for the purposes of this chapter) the displaced fluid is air.
Finally, the volume of the fluid in the cell itself may be either fixed (constant mass) or vary with time (injection). Consistent with the literature, we consider the instantaneous volume of fluid in the cell to increase as a power law in : , where is the initial volume of fluid in the HS cell (measured in m^{3}), is a dimensionless exponent, and is an injection pseudorate (in units ms), becoming precisely the injection rate for . Next, we discuss how this assumption leads to BCs for the physical problem and for the numerical scheme.
2.2 Governing Equation, Initial and Boundary Conditions
The propagation of a viscous gravity current is described by a diffusion equation for the interface , which is the shape of profile of the denser fluid. The models are derived either from porous medium flow under Darcy’s law and the Dupoit approximation (Bear, 1988, Ch. 8) or using lubrication theory with noslip along the bottom of the cell and zero shear stress at the fluid–fluid interface (Leal, 2007, Ch. 6C). The resulting velocity field is combined with a depthaveraged continuity equation to derive the nonlinear parabolic PDE for . We propose to summarize all gravity current propagation along horizontal surfaces through a single ‘thinfilm’ Oron et al (1997) equation:
(1) 
According to Engelbrecht (2015, Ch. 5), eq. (1) can be classified as an ‘evolution equation.’ The term in the parentheses on the righthand side of eq. (1), roughly, represents a fluid flux balanced by the change in height on the lefthand side. The multiplicative factor arises due to (i) geometric variations of the flow passage in the flowwise direction, (ii) porosity variations in the flowwise direction, or (iii) from the choice of coordinate system in the absence of (i) or (ii). Here, is dimensional constant depending on the flow geometry, the domain, and the fluid properties. Additionally, and are dimensionless exponents that depend on the flow geometry and fluid rheology. The quantity denoted by represents specifically the nonlinearity in these PDEs. Thus, it is necessarily a function of , and possibly for a nonNewtonian fluid (as in the third and fifth rows of table 1).^{3}^{3}3Interestingly, an ‘Laplacian’ PDE, similar to eq. (1) for a powerlaw fluid in a HS cell (third row of table 1), arises during fluid–structure interaction between a powerlaw fluid and an enclosing slender elastic tube Boyko et al (2017). This PDE can also be tackled by the proposed finitedifference scheme.
As stated in §1, several versions of eq. (1) will be explored herein, incorporating geometric variations, porosity variations, nonNewtonian behavior. The pertinent physical scenarios that will be tackled herein (using the proposed numerical scheme) are presented in table 1, which lists expressions for , , and . From a dimensional analysis of the PDE (1), it follows that the constant must have units of ms, as long as the nonlinearity has units of length (as is the case for all the models summarized in Table 1). It is worth noting that in the case of 1D, linear diffusion ( and ), becomes the ‘diffusivity’ in units of m/s.
The PDE (1) is solved on the finite spacetime interval . Here, and represent the initial and final times of the numerical simulation’s run, respectively. An initial condition (IC) is specified at , so that is known. Meanwhile, is a small positive value (close to 0). Boundary conditions (BCs) are specified at and . These involve some combination of and . The reason for taking becomes clear below.
Thus, let us now discuss such a suitable set of BCs. The BCs are based on the imposed mass conservation/growth constraint. Consider the case of a viscous gravity current in a porous slab with variable porosity , and transverse width Then, the conservation of mass constraint (see Zheng et al (2014)) takes the form
(2) 
where . In the parallel case of a HS cell with variable width and porosity , which can either be set to unity or absorbed into via , the mass constraint becomes
(3) 
Taking a time derivative of eq. (3) and employing eq. (1), we obtain
(4) 
Here, in this case of interest, as described in table 1, and Thus, we have obtained conditions relating at and to . These conditions, if satisfied, automatically take into account the imposed volume constraint from eq. (3). The calculation starting with eq. (2) is omitted as it is identical, subject to proper choice of .
For the case of propagation away from the cell’s origin (i.e., any injection of mass must occur near , specifically at ), to satisfy eq. (4), we can require that
(5a)  
(5b) 
where . Recall, the case of represents mass injection. Although eq. (3) and eqs. (5) are equivalent, the imposition of the nonlinear BC in eq. (5a) must be approached with care. It should be clear that to impose a flux near the origin (at ), we need to be finite. Then, as . On the spatial domain , such an asymptotic behavior is possible for . However, in a variablewidth cell (), the local profile and slope as blow up if they are to satisfy as . To avoid this uncomputable singularity issue, we defined the computational domain to be , where is ‘small’ but . The BC from eq. (5a) at can then be rewritten as
(6) 
It may also be of interest to consider the case of a gravity current released a finite distance away from the origin and then spreading towards . In this case, an additional length scale arises in the problem: the initial distance of the current’s edge from the origin, say . The existence of this extra length scale complicates the selfsimilarity analysis, leading to solutions of the secondkind (Barenblatt, 1996, Ch. 4), as discussed in §1. However, the numerical scheme can handle this case just as well; in fact, it requires no special consideration, unlike spreading away from the origin. Now, we may simply take and consider spreading on the domain subject to the following BCs:
(7a)  
(7b) 
which together allow us to satisfy eq. (4) and, thus, eq. (3) for all .
The most significant advantage of defining nonlinear flux BCs, such as those in eqs. (5) or (7), is that a nonlinear nonlocal (integral) constraint, such as that in eq. (2) or (3), no longer has to be applied onto the solution . Furthermore, if we start with compact initial conditions, i.e., there exists a nose location such that , then the finitespeed of propagation property of the nonlinear PDE (1) Gilding and Kersner (2004); Vázquez (2007) ensures that this nose exists for all and as well. The proposed fullyimplicit scheme inherits this property of the PDE. Therefore, we can solve the PDE on the fixed domain , without any difficulty, instead of attempting to rescale to a moving domain on which is one of the endpoints with as the BC applied there. The latter approach proposed by Bonnecaze et al (1993) (and used in more recent works Acton et al (2001) as well) leads to a number of additional variablecoefficient terms arising in the PDE (1), due to the nonGalilean transformation onto a shrinking/expanding domain. From a numerical methods point of view, having to discretize these additional terms is not generally desirable.
Having defined a suitable set of BCs, the last remaining piece of information required to close the statement of the mathematical problem at hand is the selection of pertinent initial conditions (ICs). For the case of the release of a finite fluid mass (), an arbitrary polynomial IC may be selected, as long as it has zero slope at the origin (), leading to satisfaction of the noflux boundary condition (5a). To this end, let the IC be given by
(8) 
where is a ‘releasegate’ location defining the initial position of the current’s nose, i.e., and . The constant is an arbitrary dimensionless exponent. Finally, (units of m) is set by normalizing such that the initial volume of fluid corresponds to the selected initial fluid volume, , via eq. (3).
The case of the release of a finite mass of fluid is particularly forgiving in how we set the IC, and its slope at . In fact, we could even take in eq. (8) and the scheme will provide an initial flux of fluid at , with thereafter. On the other hand, the case of mass injection () governed by the nonlinear BCs is not as forgiving. By virtue of the ‘pointsource’ mass injection at
, the slope at the origin rises sharply from the moment of mass injection. This very sharp rise has a tendency to introduce unphysical oscillations in the current profile when starting from the IC in eq. (
8). To avoid this, we must select a ‘better’ IC, which has a shape more similar to the actual solution’s singularity near . Having tested a few different options, we found that an exponential function works well:(9) 
Here, (dimensionless) and (units of m) are positive constants, ensures that the IC has no negative values and a sharp wavefront, and (units of m) is set by normalizing to the selected intial volume via eq. (3), as above.
Finally, it should be noted that the IC from eq. (8) is not used in the convergence studies for finite initial mass (§4.1 and §4.1). Rather, the IC is taken to be the exact selfsimilar solution of Ciriello et al (2016) for a powerlaw fluid in a uniformwidth () HS cell (see also the Appendix). The reasoning behind this particular choice is further expounded upon in §4.1.
Case/Variable  [ms]  [–]  [–]  [m]  

0  0  







3 The Numerical Scheme
The proposed numerical method is a finitedifference scheme using the Crank–Nicolson approach toward implicit timestepping. Our presentation follows recent literature, specifically the construction in (Zheng et al, 2014, Appendix B). The proposed scheme’s truncation error is formally of second order in both space and time, and we expect the scheme to be unconditionally stable. Furthermore, the scheme is conservative in the sense that it maintains the imposed timedependency of the fluid volume with high accuracy via a specific set of nonlinear BCs. This section is devoted to discussing all these topics one by one.
3.1 Notation: Grids, Time Steps, and Grid Functions
The PDE (1) is solved on an equispaced 1D grid of nodes with grid spacing . The solution values are kept on a staggered grid of cellcenters, which are offset by with respect to the equispaced grid points. As a result, there is a node lying a halfgridspacing beyond each domain boundary. It follows that the location of the th grid point on the staggered grid is , where . A representative grid with 12 nodes is shown in figure 5. The use of a staggered grid affords additional stability to the scheme and allows us to evaluate derivatives with secondorder accuracy via central differences, by default, using only two cellcentered values.
As stated in §2.2, the PDE (1) is solved over a time period , such that , where both the initial time and the final time of the simulation are user defined. The scheme thus performs discrete time steps each of size . The th time step advances the solution to , where . Finally, we define the discrete analog (‘grid function’) to the continuous gravity current shape, which we actually solve for, as .
3.2 The Nonlinear Crank–Nicolson Scheme
Let us denote by the continuous spatial operator acting on on the righthand side of eq. (1), i.e.,
(10) 
Since is a secondorder spatial operator and, thus, eq. (1) is a diffusion equation, we are inclined to implement a secondorderaccurate timestepping by the Crank–Nicolson scheme Crank and Nicolson (1947). The Crank–Nicolson scheme is fully implicit, which avoids the stringent restriction () suffered by explicit time discretizations of diffusion equations (Strikwerda, 2004, Ch. 6). Then, the timediscrete version of eq. (1) is
(11) 
where is the discrete analog to the continuous spatial operator defined in eq. (10). Based on the approach of Christov and Homsy (2009), the discrete spatial operator is constructed via fluxconservative central differencing using two cellface values, while staggering the nonlinear terms:
(12a)  
(12b) 
where is the slope of the gravity current’s shape. Note that the nonlinear terms, denoted by , have been evaluated the same way, i.e., at the midtimestep , for both and .
Substituting eqs. (12) into eq. (11) results in a system of nonlinear algebraic equations because is evaluated at midtimestep and, thus, depends on both (known) and
(unknown). This system must be solved for the vector
(), i.e., the approximation to the gravity current’s shape at the next time step. Solving a large set of nonlinear algebraic equations can be tedious and computationally expensive. A simple and robust approach to obtaining a solution of the nonlinear algebraic system is through fixedpoint iterations, or ‘the method of internal iterations’ (Yanenko, 1971). Specifically, we can iteratively compute approximations to , the grid function at the new time step, by replacing it in eq. (11) with , where . Then, the proposed numerical scheme takes the form:(13)  
The key idea in the method of internal iterations is to evaluate the nonlinear terms from information known at iteration and the previous time step , while keeping the linear slopes from the next time step at iteration . This manipulation linearizes the algebraic system, at the cost of requiring iteration over . Upon convergence of the internal iterations, is simply the last iterate . Before we can further discuss the iterations themselves or their convergence, we must define our discrete approximations for and .
The operator is essentially a second derivative, so we take inspiration from the standard way of constructing the threepoint central finitedifference formula for the second derivative Strikwerda (2004). Therefore, can be discretized using a twopoint centraldifference approximation on the staggered grid. For example, at any time step:
(14) 
Next, following Zheng et al (2014); Christov and Deng (2002), we evaluate at by averaging the known values at and or and , respectively. Likewise, to approximate , we average the known values: at and at the previous internal iteration. In other words, our approximation of the nonlinear terms is
(15a)  
(15b) 
Equations (15) afford improved stability for nonlinear PDEs, while preserving the conservative nature of the scheme (as will be shown in §4.2), as discussed by Von Rosenberg (1975) who credits the idea of averaging nonlinear terms across time stages and staggered grid points to the seminal work of Douglas Jr. et al (1959, 1962). The scheme thus described is depicted by the stencil diagram in Figure 6.
Here, it is worthwhile noting that, the classical Crank–Nicolson Crank and Nicolson (1947) scheme is only provably unconditionally stable Strikwerda (2004) when applied to a linear diffusion equation. It was suggested by Christov and Homsy (2009) that the current approach provides additional stability to this nonlinear scheme for large time steps. But, since our problem is nonlinear, some care should be taken in evaluating how large of a time step could be taken. Nevertheless, it is still expected that the largest stable will be independent of .
A complication arising in the present context is that we focus on the case of a powerlaw nonNewtonian viscous gravity current spreading in a variablewidth cell. As a result, recalling table 1, this model features in , unlike the Newtonian case. While the temporal accuracy of the scheme is ensured through the robust implementation of the nonlinear Crank–Nicolson timestepping, the spatial accuracy is contingent upon the discretization of in . A further consequence is that, once we discretize , the discretization of becomes nonlocal (i.e., it requires information beyond the th grid point). Nevertheless, the overall scheme still only requires a threepoint stencil for . In particular, for interior grid points, we use a centraldifference formula, giving rise to the expression (at any time step):
(16) 
This choice of approximation ensures secondorder accuracy at all interior grid nodes. However, at the second () and the penultimate () nodes, the secondorder accurate approximation to in as defined in eqs. (15) requires the unknown values and , respectively. To resolve this difficulty, we use ‘biased’ (backward or forward) threepoint difference approximations:
(17a)  
(17b) 
Finally, substituting the discretization for from eq. (14) into eq. (13), it is possible to rearrange the scheme into a tridiagonal matrix equation:
(18) 
for the interior grid points . In eq. (18), the righthand side and the variable coefficients in brackets on the lefthand side are both known, based on , at any given internal iteration . Then, each internal iteration involves the inversion of a tridiagonal matrix to solve for the grid function . The inversion of this tridiagonal matrix can be performed efficiently with, e.g., ‘backslash’ in Matlab. Subsequently, the coefficient matrix must be recalculated for each internal iteration because of the dependency of on arising from eqs. (15), (16) and (17)
The iterations in eq. (18) are initialized with () and continue until an iteration is reached at which a relative error tolerance is met. Specifically,
(19) 
Only a small number (typically, less than a dozen) of internal iterations are required at each time step, making the scheme quite efficient overall.
A detail remains, however. The algebraic system defined in eq. (18) applies to all interior nodes, i.e., . To complete the system, we must define rows and , which arise from the discretization of the nonlinear BCs, which comes in §3.4. Upon completing the latter task successfully, becomes the grid function at the next time step upon the completion of the internal iterations, and the time stepping proceeds.
3.3 The Special Case of Linear Diffusion
A noteworthy special case of the proposed finitedifference scheme arises from setting the dimensionless exponents (i.e., no spatial variation of the diffusivity) and (linear diffusion). Then, eq. (18) can be simplified and rearranged in the form ():
(20) 
If the grid function represents the temperature field along a 1D rigid conductor situated on , eq. (20) is then the original secondorder (in space and time) numerical scheme proposed by Crank and Nicolson (1947) to solve a linear (thermal) diffusion equation (Strikwerda, 2004, §6.3). As such, this simplification helps illustrate the mathematical roots of the current scheme, and how we have generalized the classical work.
3.4 Implementation of the Nonlinear Boundary Conditions
As discussed in §2.2, the boundary conditions are a manifestation of the global mass conservation constraint, eq. (2) or (3), imposed on eq. (1). The BCs described in eqs. (5) and (7) are defined at the ‘real’ boundaries of the domain, i.e., at and . The numerical scheme is implemented over a staggered grid. This allows for derivatives at and to be conveniently approximated using central difference formulas using two nearby staggered grid points. In this manner, the BC discretization maintains the scheme’s second order accuracy in space and time. Accordingly, for the case of a current spreading away from the cell’s origin, eqs. (5) are discretized in a ‘fullyimplicit’ sense (to further endow numerical stability and accuracy to the scheme Christov and Deng (2002)) as follows:
(21)  
(22) 
Within the internal iterations, however, is known independently of and . Hence, we can express the first () and last () equations, which defined the respective rows in the tridiagonal matrix stemming from eq. (18), as
(23a)  
(23b) 
Similarly, we can derive the discretized BCs for spreading towards the origin, upon its release a finite distance away from the origin, from eqs. (7). Then, the first () and last () equations, which defined the respective rows in the tridiagonal matrix, as
(24a)  
(24b) 
4 Convergence and Conservation Properties of the Scheme
At this point, the numerical scheme and boundary conditions defined in eqs. (18) and (23) or (24) form a complete description of the numerical solution to the parabolic PDE from eq. (1), for a gravity current propagating away from the origin. We have claimed that the finitedifference scheme is conservative (i.e., it accurately maintains the imposed timedependency of the fluid volume set by eq. (3)) and has secondorder convergence. These aspects of the scheme will be substantiated in §4.1 and §4.2, respectively. The computational domain’s dimensions, which are set by and , and the properties of fluid being simulated are summarized in table 2. For definiteness, in this chapter we select the fluid properties to be those of a 95% glycerolwater mixture in air at 20°C (see Cheng (2008); Volk and Kähler (2018)).
4.1 Estimated Order of Convergence
First, we seek to justify the formal accuracy (order of convergence) of the proposed scheme through carefully chosen numerical examples. To do so, we pursue a series of benchmarks that are successively ‘more complicated’ (from a numerical perspective). First, we simulate the case of a centrally released fixed mass of fluid propagating in two directions (§4.1). Second, we simulate the unidirectional spreading of a fixed mass of fluid (§4.1). Last, we simulate the unidirectional spreading of a variable fluid mass (§4.1) by taking into account injection of fluid at the boundary.
Parameter  Value  Units 

Channel length,  0.75  m 
Width coefficient,  0.017390  m 
Total released mass,  0.31550  kg 
Density difference,  1250.8  kg/m^{3} 
Consistency index ()  
or dynamic viscosity (),  0.62119  Pas 
In each of these three cases, there is a need for a reliable benchmark solution against which the numerical solutions on successively refined spatial grids can be compared. For the case of the release of a fixed mass of fluid, an exact selfsimilar solution is provided by Ciriello et al (2016). Specifically the solution is for a powerlaw fluid in uniform HS cell (). The derivation of the selfsimilar solution is briefly discussed in the Appendix. We use this solution as the benchmark. As mentioned in §1
, parabolic equations ‘forget’ their IC and the solution becomes selfsimilar after some time. However, for a general PDE, it is difficult (if not impossible) to estimate how long this process takes. Therefore, to ensure a proper benchmark against the exact selfsimilar solution, we start the simulation with the exact selfsimilar solution evaluated at some nonzero initial time (
). Then, we let the current propagate up to a final time , with the expectation that the current will remain in the selfsimilar regime for all . Comparing the final numerical profile with the exact selfsimilar solution at then allows for a proper benchmark.To quantify the error between a numerical solution and a benchmark solution at , we use three standard functionspace norms Evans (2010):
(25a)  
(25b)  