Solving Nonlinear Parabolic Equations by a Strongly Implicit Finite-Difference Scheme

by   Aditya A. Ghodgaonkar, et al.
Purdue University

We discuss the numerical solution of nonlinear parabolic partial differential equations, exhibiting finite-speed of propagation, via a strongly implicit finite-difference scheme with formal truncation error O[(Δ x)^2 + (Δ t)^2 ]. Our application of interest is the spreading of viscous gravity currents in the study of which these type of differential equations arise. Viscous gravity currents are low Reynolds number flow phenomena in which a dense, viscous fluid displaces a lighter (usually immiscible) fluid. The fluids may be confined by the sidewalls of a channel or propagate in an unconfined two-dimensional (or axisymmetric three-dimensional) geometry. Under the lubrication approximation, the mathematical description of the spreading of these fluids reduces to solving the so-called thin-film equation for the current's shape h(x,t). To solve such nonlinear parabolic equations we propose a finite-difference scheme based on the Crank--Nicolson idea. We implement the scheme for problems involving a single spatial dimension (i.e., two-dimensional, axisymmetric or spherically-symmetric three-dimensional currents) on a uniform but staggered grid. We benchmark the scheme against analytical solutions and highlight its strong numerical stability by specifically considering the spreading of non-Newtonian power-law fluids in a variable-width confined channel-like geometry (a `Hele-Shaw cell') subject to a given mass conservation/balance constraint. We show that this constraint can be implemented by re-expressing it as nonlinear flux boundary conditions on the domain's endpoints. Then, we show numerically that the scheme achieves its full second-order accuracy in space and time. Furthermore, we also highlight through numerical simulations how the proposed scheme respects, with high accuracy, the mass conservation/balance constraint.


A modeling and simulation study of anaerobic digestion in plug-flow reactors

A mathematical model for anaerobic digestion in plug-flow reactors is pr...

Constraint-Aware Neural Networks for Riemann Problems

Neural networks are increasingly used in complex (data-driven) simulatio...

Refined RBF-FD analysis of non-Newtonian natural convection

In this paper we present a refined RBF-FD solution for a non-Newtonian f...

Error Inhibiting Methods for Finite Elements

Finite Difference methods (FD) are one of the oldest and simplest method...

A moving-boundary model of reactive settling in wastewater treatment

Reactive settling is the process of sedimentation of small solid particl...

High accuracy power series method for solving scalar, vector, and inhomogeneous nonlinear Schrödinger equations

We develop a high accuracy power series method for solving partial diffe...

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 20th 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 concentration-dependent diffusivity Barenblatt (1952).111More specifically, Barenblatt (1952) (see also (Ostriker et al, 1992, p. 13)) credits the observation of finite-speed of propagation in a nonlinear diffusion equation to a difficult-to-find 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, all-encompassing 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 traveling-wave and ‘wavefront’ solutions now rests upon solid mathematical foundations Gilding and Kersner (2004); Vázquez (2007), including the case of gradient-dependent 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 finite-speed wave-like 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 flow-wise 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 self-similarity transformation. If the similarity variable can obtained by a scaling (dimensional) analysis, then the solution is termed a self-similar solution of the

first kind (Barenblatt, 1996, Ch. 3). Specifically, the transformation is ( in units222Throughout the chapter, we use SI units for all dimensional quantities. of meters), where is the similarity variable (dimensionless), is the self-similar 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 one-dimensional (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 m2/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 self-similarity 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 .

Figure 3: Spreading via 1D linear diffusion and the approach to the universal intermediate self-similar asymptotics. (a) An arbitrary wavy IC spreads and levels until reaching a flat steady state . (b) A first-kind self-similar transformation (obtained from dimensional analysis) yields a universal profile (highlighted in gold) towards which the solution evolves in the intermediate period after the IC is forgotten (but prior to leveling).

Having illustrated the notion of first-kind self-similarity 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 CO2 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 self-similar 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 non-Newtonian gravity currents, specifically ones for which the denser fluid obeys a power-law rheology. This tractable model of non-Newtonian rheological response is also known as the Oswald–de Weale fluid Bird et al (1987). In unidirectional flow, the power-law 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 power-law 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 power-law fluids under confinement (i.e., in Hele-Shaw cells) using the lubrication approximation. These works have contributed to the use of a modified Darcy law to model the flow of non-Newtonian fluids in porous media using the analogy to flow in Hele-Shaw cells. Aronsson and Janfalk (1992) were perhaps the first to combine a Darcy law for a power-law 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 two-dimensional axisymmetric spreading in media with uniform porosity but variable permeability. All these flows are of interest because exact analytical self-similar 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 self-similar 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 self-similar 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 self-similarity transformation may no longer be possible simply by scaling (dimensional) arguments. Gratton and Minotti (1990)classified a number of such situations, including the so-called ‘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 self-similar regime, even though a self-similar transformation cannot be obtained by scaling arguments alone. The exponents and in the transformation are unknown a priori, hence this situation represents a self-similarity 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 phase-plane 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 ‘pre-wetting film’ ahead of the current’s sharp wavefront ( where ) was required to avoid numerical instabilities. The scheme therein was also first-order accurate in time only. In this chapter, we propose a modern, high-order-accurate implicit numerical method for use in such problems.

Figure 4: A summary of the gravity current flows and domains considered in this work. (a) Flow away from the origin in a completely porous () HS cell of variable width given by (, ). (b) Flow in uniformly porous () passage of variable width given by the same as in (a). (c) Flow in a uniform-width slab (i.e., ) with horizontally heterogeneous porosity and permeability given by and , respectively. The effective permeability of the medium in (a) and (b) is set by the Hele-Shaw analogy via the width: . Figure reproduced and adapted with permission from [Zheng et al, Influence of heterogeneity on second-kind self-similar solutions for viscous gravity currents, J. Fluid Mech., vol. 747, p. 221] © Cambridge University Press 2014.

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 non-Newtonian viscous gravity currents in variable geometries (specifically, Hele-Shaw 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 finite-difference scheme for simulating the spreading of a finite mass of Newtonian fluid in a variable-width Hele-Shaw cell. Owing to its accuracy and stability, this finite-difference scheme has been recently applied by Alhashim and Koch (2018) to study hydraulic fracturing of low-permeability 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 finite-difference 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 finite-difference scheme (up to a specified physical time) against an analytical solution obtained through a self-similar 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 non-Newtonian (power-law) displacing fluid in a variable-width channel. As a benchmark, we use previously derived first-kind self-similar 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 non-Newtonian fluids. We study their spreading in a fixed- or variable-width channel geometry (also known as a “Hele-Shaw 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 Hele-Shaw (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 uniform-thickness 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 power-law 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 quickly-widening 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 non-Newtonian. Specifically, it obeys the power-law rheology. In unidirectional flow, the one unique non-trivial 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 shear-thinning (e.g., blood), and fluids with are termed shear-thickening (e.g., dense particulate suspensions). In the special case , the power-law 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 m3), is a dimensionless exponent, and is an injection pseudo-rate (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 no-slip along the bottom of the cell and zero shear stress at the fluid–fluid interface (Leal, 2007, Ch. 6-C). The resulting velocity field is combined with a depth-averaged continuity equation to derive the nonlinear parabolic PDE for . We propose to summarize all gravity current propagation along horizontal surfaces through a single ‘thin-film’ Oron et al (1997) equation:


According to Engelbrecht (2015, Ch. 5), eq. (1) can be classified as an ‘evolution equation.’ The term in the parentheses on the right-hand side of eq. (1), roughly, represents a fluid flux balanced by the change in height on the left-hand side. The multiplicative factor arises due to (i) geometric variations of the flow passage in the flow-wise direction, (ii) porosity variations in the flow-wise 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 non-Newtonian fluid (as in the third and fifth rows of table 1).333Interestingly, an ‘-Laplacian’ PDE, similar to eq. (1) for a power-law fluid in a HS cell (third row of table 1), arises during fluid–structure interaction between a power-law fluid and an enclosing slender elastic tube Boyko et al (2017). This PDE can also be tackled by the proposed finite-difference scheme.

As stated in §1, several versions of eq. (1) will be explored herein, incorporating geometric variations, porosity variations, non-Newtonian 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 space-time 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


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


Taking a time derivative of eq. (3) and employing eq. (1), we obtain


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


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 variable-width 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 re-written as


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 self-similarity analysis, leading to solutions of the second-kind (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:


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 finite-speed 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 fully-implicit 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 variable-coefficient terms arising in the PDE (1), due to the non-Galilean 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 no-flux boundary condition (5a). To this end, let the IC be given by


where is a ‘release-gate’ 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 ‘point-source’ 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:


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 self-similar solution of Ciriello et al (2016) for a power-law fluid in a uniform-width () HS cell (see also the Appendix). The reasoning behind this particular choice is further expounded upon in §4.1.

Case/Variable [ms] [–] [–] [m]

Newtonian fluid,
fixed-width HS cell: .
(see Huppert and Woods (1995))
0 0
Newtonian fluid,
variable-width HS cell: .
(see Zheng et al (2014))
Power-law fluid: ,
variable-width HS cell: .
(see Di Federico et al (2017); Longo (2017))
Newtonian fluid,
2D porous medium,
variable porosity: ,
variable permeability: .
(see Zheng et al (2014))
Power-law fluid: ,
2D porous medium,
variable porosity: ,
variable permeability: .
(see Ciriello et al (2016))

Table 1: Selected models of the propagation of viscous gravity currents herein simulated by a finite-difference scheme.

3 The Numerical Scheme

The proposed numerical method is a finite-difference scheme using the Crank–Nicolson approach toward implicit time-stepping. 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 time-dependency 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 cell-centers, which are offset by with respect to the equispaced grid points. As a result, there is a node lying a half-grid-spacing 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 second-order accuracy via central differences, by default, using only two cell-centered 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 .

Figure 5: A sample twelve-node equispaced but staggered 1D grid. The grid nodes are staggered by half a grid step from the cell faces. The boundary conditions are implemented at the ‘real’ domain boundaries (here marked by x). The two grid points outside the physical domain (i.e., or and in this example) are used to implement the Neumann BCs, which require computing a derivative at the ‘real’ domain boundaries (i.e., or and in this example).

3.2 The Nonlinear Crank–Nicolson Scheme

Let us denote by the continuous spatial operator acting on on the right-hand side of eq. (1), i.e.,


Since is a second-order spatial operator and, thus, eq. (1) is a diffusion equation, we are inclined to implement a second-order-accurate time-stepping 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 time-discrete version of eq. (1) is


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 flux-conservative central differencing using two cell-face values, while staggering the nonlinear terms:


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 mid-time-step , for both and .

Substituting eqs. (12) into eq. (11) results in a system of nonlinear algebraic equations because is evaluated at mid-time-step 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 fixed-point 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:


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 three-point central finite-difference formula for the second derivative Strikwerda (2004). Therefore, can be discretized using a two-point central-difference approximation on the staggered grid. For example, at any time step:


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


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.

Figure 6: Representative stencil of the proposed scheme. After performing internal iterations, the nonlinear terms are computed at the intermediate stage ‘’ (highlighted in blue) from the known quantities and . The unknown quantity at the next internal iteration, stage ‘’ (highlighted in red), is found by solving the linear system in eq. (18). The process continues until the convergence criterion in eq. (19) is met, yielding the (initially unknown) solution at .

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 power-law non-Newtonian viscous gravity current spreading in a variable-width 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 time-stepping, 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 three-point stencil for . In particular, for interior grid points, we use a central-difference formula, giving rise to the expression (at any time step):


This choice of approximation ensures second-order accuracy at all interior grid nodes. However, at the second () and the penultimate () nodes, the second-order 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) three-point difference approximations:


Finally, substituting the discretization for from eq. (14) into eq. (13), it is possible to re-arrange the scheme into a tridiagonal matrix equation:


for the interior grid points . In eq. (18), the right-hand side and the variable coefficients in brackets on the left-hand 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,


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 finite-difference 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 ():


If the grid function represents the temperature field along a 1D rigid conductor situated on , eq. (20) is then the original second-order (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 ‘fully-implicit’ sense (to further endow numerical stability and accuracy to the scheme Christov and Deng (2002)) as follows:


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


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


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 finite-difference scheme is conservative (i.e., it accurately maintains the imposed time-dependency of the fluid volume set by eq. (3)) and has second-order 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% glycerol-water 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/m3
Consistency index ()
or dynamic viscosity (), 0.62119 Pas
Table 2: Summary of the simulation parameters used in convergence and conservation studies. The fluid was assumed to be a 95% glycerol-water mixture at 20°C. The width exponent and fluid’s rheological index were varied on a case-by-case basis to simulate different physical scenarios.

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 self-similar solution is provided by Ciriello et al (2016). Specifically the solution is for a power-law fluid in uniform HS cell (). The derivation of the self-similar 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 self-similar 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 self-similar solution, we start the simulation with the exact self-similar solution evaluated at some non-zero initial time (

). Then, we let the current propagate up to a final time , with the expectation that the current will remain in the self-similar regime for all . Comparing the final numerical profile with the exact self-similar 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 function-space norms Evans (2010):