Thin-gap averaging of variable-viscosity flows: application to thermoviscous fingering

Several physico-chemical problems offer instances of transport phenomena taking place in confined geometries, such as a Hele-Shaw cell, wherein one of the spatial dimensions is much smaller than the other two. In such cases, a convenient and often-used strategy is to employ a thin-gap averaged model based on a parabolic cross-gap velocity profile to describe the underlying physics. This yields a simplified Darcy-like averaged model, which then serves as an efficient tool for understanding the behavior of such systems (Homsy, 1987). However, one often encounters instances wherein the cross-gap velocity profile may significantly deviate from the Poiseuille parabolic flow. As an example, this is true when the flow is coupled to the transport of a scalar, such as temperature. This scenario occurs in the case of basaltic magma flow through fissures, where cross-gap variations in temperature result in a non-uniform viscosity across the thin gap (Wylie & Lister, 1995). Even in such cases, though, ad hoc thin-gap averaging based on a parabolic velocity profile is often employed (Pearson et al., 1973; Whitehead & Helfrich, 1991; Bercovici, 1994; Helfrich, 1995). Other situations where this issue arises include Hele-Shaw flows of suspensions of colloids, non-Brownian particles, or microswimmers with a number concentration that varies across the thin gap (Lyon & Leal, 1998; Marmet et al., 2017; Vennamneni et al., 2020), thus resulting in a non-uniform effective viscosity and a non-parabolic flow profile.

There are of course special cases where no cross-gap variations occur, for example, when adiabatic or no-flux conditions are imposed on the sidewalls (Pritchard, 2009). However, in general, this is not true, as non-uniformities can arise due to transport across the boundaries, as when the fluid loses heat through the side walls (Helfrich, 1995; Nagatsu et al., 2009), or due to gradient-producing physical effects, such as shear-induced migration in suspensions (Leighton & Acrivos, 1987; Nott & Brady, 1994), wall effects (Marmet et al., 2017), and phoretic motion (Buongiorno, 2005; Choudhary et al., 2020). The inherent inconsistency in using Darcy models in such situations has been recognized previously; however, in the past, one has had to either embrace the complexity of the full three-dimensional equations (Wylie & Lister, 1995) or treat the problem in a linearized limit (Viola et al., 2017). In this work, we address this issue by developing a low-dimensional modelling strategy that takes advantage of the thin-gap geometry while systematically accounting for cross-gap variations in viscosity and the consequent non-parabolic velocity profile.

Given the wide range of possible applications, we first present the averaging method in a general setting in §2 for an arbitrary cross-gap variation of viscosity. Our approach is based on the weighted residual integral boundary layer (WRIBL) method, which was developed to model the inertial dynamics of thin liquid films (Ruyer-Quil & Manneville, 2002), and has been subsequently used in a variety of free-surface and two-phase flow problems (Kalliadasis et al., 2012; Trevelyan & Kalliadasis, 2004; Oron & Heining, 2008; Dietze & Ruyer-Quil, 2013), including cases where the fluid is non-Newtonian (Ruyer-Quil et al., 2012; Pradas et al., 2014). This averaging method has also been used to incorporate inertial corrections into the Darcy model for constant-viscosity Hele-Shaw flows (Ruyer-Quil, 2001). The first step of the WRIBL procedure is to obtain simplified long-wave partial differential equations from the full incompressible Navier-Stokes equations, based on the disparity of cross-gap and in-plane length scales. These equations are then subsequently averaged across the thin gap using appropriate weight functions. While several weighting strategies are possible, the Galerkin weighting method was shown by Ruyer-Quil & Manneville (2002) to be the most efficient and we employ the same. The key step in the method is choosing the weight functions so that an asymptotically consistent, closed system of averaged equations are obtained; it is here that the effect of a non-uniform viscosity is incorporated.

After developing the general averaged flow model, we then apply it to the problem of thermoviscous fingering in a Hele-Shaw geometry (Helfrich, 1995; Wylie & Lister, 1995). This problem serves as an example to elucidate the significance of incorporating cross-gap variations during the averaging procedure. Figure 1 shows a schematic of the system, in which a pressure drop drives the flow of a hot fluid through a narrow slot, whose side-walls are held at a fixed lower temperature. This mimics the upward flow of pressurized hot magma through a narrow fissure whose walls, composed of relatively cold rock, cool the magma as it flows (Wylie et al., 1999). In a simplified model the magma is taken to be a Newtonian fluid with a viscosity that decreases sharply as it cools (Huppert, 2002; Whitehead & Griffiths, 2001)—a property that nonlinearly couples flow and heat transport. For a sufficiently strong temperature-dependence of viscosity, the laterally-uniform steady flow is linearly unstable over a range of flow rates, giving rise to laterally alternating bands of low-viscosity (fast-moving) and high viscosity (slow-moving) magma (Helfrich, 1995; Wylie & Lister, 1995).

This problem serves as a useful testing ground for the proposed averaging method because the Dirichlet boundary conditions applied at the cold walls cause the temperature to vary across the thin gap, which in turn produces a non-trivial variation in viscosity. We consider two cases, distinguished by the temperature of the walls which is assumed to be either constant or linearly decreasing along the direction of flow. While the latter scenario has no geological basis, it serves to emphasize the importance of cross-gap viscosity variations, possibly being relevant to industrial processes such as injection moulding and glass manufacturing in which the walls are actively cooled. In § 3, the WRIBL method is used to average the advection-diffusion equation for temperature, which is then coupled to the averaged flow equations of §2 to yield the final coupled WRIBL model.

The Darcy model of Helfrich (1995) is presented in § 4, where we discuss the ad hoc assumptions that must be made to obtain it from the WRIBL model.

In §5, we compare the predictions of the coupled WRIBL model with that of the Darcy model, as well as with calculations of the full three dimensional (3D) equations. We focus on the multiplicity of the uniform steady base state and on its linear stability. The WRIBL model is found to be much more accurate than the Darcy model, with the latter significantly over-predicting the range of multiplicity and instability. However, the key qualitative features of the 3D results are captured by both averaged models. Therefore, while the Darcy model is found to be adequate for qualitative analysis, incorporating cross-gap variations of viscosity during the averaging procedure is important for making quantitative estimates.

In this section, we present the depth-averaging strategy for the flow of a fluid with variable viscosity in a thin-gap geometry. For the sake of generality, we consider an arbitrary cross-gap variation of viscosity. The averaged model can be then used in different physical settings by coupling it to an additional conservation equation for the field variable on which the viscosity depends, as is done for temperature-dependent viscosity in §3.

A Cartesian co-ordinate system is adopted, with $(u,v,w)$ being the components of the velocity vector ($\vec{v}$) in the ($x,y,z$) directions (figure 1). The depth ($d$) of the channel in the wall-bounded $y$ direction is much smaller than its length ($L$) along the primary flow direction ($z$), as well as its width in the lateral direction ($x$). We assume that all field variables, including the velocity and viscosity, vary much more slowly in the $x$ and $z$ directions than in the thin-gap $y$ direction. Exploiting this disparity in length scales, a nonlinear reduced-order model, based on the WRIBL method, is developed for the in-plane, gap-averaged, flow rates $q_{z}(z,x,t)$ and $q_{x}(z,x,t)$.

The flow is governed by the equation of continuity, and the Navier-Stokes equation with variable viscosity, $\mu(y;z,x,t)$, given below in dimensional form:

where, $\rho$ is the fluid density, $\textbf{S}^{*}=\mu(\nabla^{*}\vec{v}^{*}+\nabla^{*}\vec{v}^{*T})$ is the deviatoric part of the Cauchy stress tensor, and $p^{*}$ is the isotropic fluid pressure. The asterisk (^∗) denotes dimensional variables.

At the walls, we apply no-slip boundary conditions:

The following scaling is used to non-dimensionalize the governing equations:

where the subscript $c$ denotes a characteristic scale, and $\epsilon$ is the long-wave parameter defined as $\varepsilon\equiv d/L\ll 1$. Also $\mu_{0}$ is a characteristic value of viscosity for the system, for example an extremal or averaged value, and $U$ is a characteristic velocity. The scaled governing equations and boundary conditions, retaining terms only up to $\mathcal{O}(\varepsilon)$ and ignoring the higher order terms, reduce to

In obtaining the above equations, the Reynolds number $\Rey=\rho dU/\mu_{0}$ is assumed to be small enough for inertial terms to be disregarded. Also, we assume that $\partial_{y}\mu\sim O(1)$ while $\partial_{x}\mu\sim O(\varepsilon)$ and $\partial_{z}\mu\sim O(\varepsilon)$, which requires the field variable (e.g., temperature) on which the viscosity depends to also vary slowly in the $x$ and $z$ directions. Typically, this field variable will be governed by an additional conservation equation which must also admit this separation of cross-gap and in-plane length scales.

From equation (4b) we see that neglecting variations of $\mu$ in the cross-gap ($y$) direction, as is done in Darcy models, introduces $O(1)$ errors via $\partial_{y}\mu\partial_{y}u$ and $\partial_{y}\mu\partial_{y}w$ terms. Here, we consistently incorporate all leading order terms using the weighted residual integral boundary layer (WRIBL) technique, to obtain the reduced-order nonlinear model up to $O(\varepsilon)$. We follow the procedure described in Dietze & Ruyer-Quil (2013), and begin by decomposing the velocity field into $\mathcal{O}(1)$ and $\mathcal{O}(\varepsilon)$ parts:

The leading-order contributions are required to satisfy (4b) truncated to $\mathcal{O}(1)$, but with additional $\mathcal{O}(\varepsilon)$ forcing terms, $K_{u}$ and $K_{w}$, which ensures that the local gap-averaged flow rates, $q_{x}$ and $q_{z}$, are determined entirely by $\hat{u}$ and $\hat{w}$. We thus have the following definitions:

Here, $K_{u}$ and $K_{w}$ are $y$-independent functions, which are determined by solving the boundary value problems in (6) and applying the integral conditions. Thus, we obtain the leading order velocity profiles $\hat{u}$ and $\hat{w}$, which only depend on $x$, $z$ and $t$ through their dependence on $q_{x}$ and $q_{z}$, respectively. Note that $\hat{u}$ and $\hat{w}$ will be parabolic only if $\mu$ is uniform in the $y$ direction. Conversely, when $\mu$ is a function of $y$ then (6) will capture the consequent non-parabolic variation of velocity across the thin gap.

The specification of the leading order velocity field is completed by using the continuity equation (4a) to determine the leading-order cross-gap velocity $\hat{v}$:

where we have used the no-penetration boundary condition, $\hat{v}|_{y=0}=0$.

Returning to the decomposition (5), we see that the definitions in (6) imply that the $\mathcal{O}(\varepsilon)$ corrections to the velocity field satisfy the following boundary and integral conditions:

The domain equations for the $\mathcal{O}(\varepsilon)$ corrections may be obtained by substituting (5) and (6) into (4), but they are not needed to close the averaged model.

We now proceed to average the governing equations, beginning with the continuity equation (4a), which on integrating across the thin gap and using the no-slip and no-penetration conditions yields the following exact mass balance equation:

Next, we must average the momentum equations in order to obtain relations between the flow rate and the pressure gradient. In order to consistently close the $\mathcal{O}(1)$ viscous term, however, we carry out a weighted average by multiplying (4b) with a weight function $F_{v}$ and then integrate across the gap, to yield (after substituting the velocity decomposition (5)):

The $O(\epsilon)$ corrections, $\tilde{u}$ and $\tilde{w}$, may be entirely eliminated from (10) by using a Galerkin weight function defined as follows (Ruyer-Quil & Manneville, 2002):

Now, on applying integration by parts twice to the right-most integrals of (10) and using the boundary conditions (8) and (11), we obtain the closed WRIBL flow model:

where we have repeated the mass balance equation (9) for completeness.

The cross-gap variation of viscosity, and the consequent non-parabolic local flow profile, enter the WRIBL model (12) through the weight function $F_{v}$, as well as through the leading-order velocity profiles $\hat{u}$ and $\hat{w}$. In fact, if $\mu$ is constant, then $\hat{u}$, $\hat{w}$, and $F_{v}$ become parabolic, and (12) reduces to the Darcy flow model. Note that while $\hat{u}$, $\hat{w}$ and $\hat{v}$ are not seen in (12), they will appear in the convective terms of the conservation equation for the field variable that determines the viscosity (as it does in the energy balance equation in §3).

In this derivation, we have, for simplicity, retained terms only up to $\mathcal{O}(\varepsilon)$. However, the model can be further improved by retaining $O(\varepsilon^{2})$ terms and using a similar weighted residual procedure to include the effects of mechanical inertia (Ruyer-Quil et al., 2012) and in-plane viscous diffusion, as demonstrated by Dietze & Ruyer-Quil (2013, 2015) for interfacial flows. (The leading-order velocity profiles $\hat{u}$ and $\hat{w}$ will then make an appearance in the averaged flow equations via these additional terms.) The resulting $\mathcal{O}(\varepsilon^{2})$ WRIBL model would be an improved version of the inertial Darcy-Brinkmann equation, with corrections due to cross-gap variations in viscosity. We leave such extensions for future work, however, and turn to a specific application of the WRIBL flow model (12).

We now take up the problem of thermoviscous fingering which arises in the Hele-Shaw flow of a fluid with a viscosity that decreases sharply with increasing temperature. Figure 1 presents a schematic of the Hele-Shaw flow, in which hot fluid, entering the channel at a mean temperature $T_{H}$, is driven primarily in the $z$ direction by a fixed pressure drop $P_{in}-P_{out}$. The temperature of the colder walls is assumed either to be constant, with a value $T_{C}$, or linearly decreasing along the flow direction, as $T_{H}-(T_{H}-T_{C})z/L$. The former case has been studied before in the context of magma flows (Wylie & Lister, 1995; Helfrich, 1995), while the latter case, which is more relevant to industrial processes such as injection moulding, is considered here for the first time. In both cases, the fluid will lose heat to the colder walls and become cooler and more viscous as it flows. Importantly, the temperature and therefore the viscosity will vary across the gap and result in a non-parabolic cross-gap velocity profile (Wylie & Lister, 1995). It is this leading-order feature of the problem that is neglected by the Darcy model used in previous work (Helfrich, 1995), and which we seek to include using the WRIBL averaging scheme.

This problem serves as an ideal test case because the dependence of viscosity on temperature is not merely incidental, but rather plays a key physical role in thermoviscous fingering. This channelling instability arises because advective heat supply, which competes with heat loss to the cold side-walls, depends on the viscosity and thus on the local temperature. A positive perturbation to the local temperature inside the fissure decreases the local viscosity, and therefore increases the velocity. This, in turn, decreases the residence time of the fluid parcel in the conduit, thus reducing heat loss to the cold walls and further increasing the local temperature. It is this positive feedback that underlies the fingering instability, which ultimately produces channels of low-viscosity, fast-moving magma separated by regions of high viscosity slow-moving magma. On the surface, at the outlet of the fissure, this instability manifests as a transition in the spatial form of the emerging flow: from a uniform sheet of magma to well-separated spouts.

To model this coupled flow and heat transport problem, we derive an averaged transport equation for temperature, which is then coupled to the averaged flow model. We use the WRIBL method for this purpose, and the derivation parallels that of the flow model (§2), except that the situation is now simpler as we are dealing with a single scalar field and a constant diffusion coefficient. We begin with the energy equation in dimensional form:

where, $\kappa$ is the thermal diffusivity.

In case 1, we consider the temperature of both side walls $T_{w}^{*}$ to be constant,

and assume that the fluid enters with a mean temperature $T_{H}$ and a parabolic cross-gap profile,

Wylie & Lister (1995) showed that if the fluid enters with a uniform temperature, then the temperature profile rapidly evolves, over a short entrance length, to a near-parabolic profile. This initial rapid evolution may be seen as a relaxation of the system from its initial condition to the low-dimensional slow-manifold on which the dynamics may be described by reduced-order models (Roberts, 2015; Roberts & Bunder, 2017). In general, this initial relaxation cannot be captured by averaged models, although error-minimizing strategies have been proposed (Balakotaiah & Ratnakar, 2010; Roberts, 1992, 2001). Here, because our goal is to compare between the WRIBL and Darcy models, we avoid the singular initial condition of a uniform temperature and instead use the parabolic initial condition (15).

In case 2, we enforce a linearly decreasing temperature profile along the walls,

and consider the temperature of the fluid at the inlet to be the same as the local wall-temperature, which implies that

In contrast to case 1, the fluid in case 2 starts out in thermal equilibrium with the walls, and so the issue of a rapid initial relaxation does not arise. This contradistinction is why we find it instructive to consider both cases.

To scale these equations, we use (3) along with the nondimensional temperature,

and on dropping terms of order higher than $\varepsilon$ we obtain the following scaled heat transport equation:

Here, we have chosen the characteristic velocity to be $U=\kappa/d$, which amounts to setting the Péclet number $Pe=Ud/\kappa$ to unity, in contrast to $\Rey$ which was taken to be $\mathcal{O}(\varepsilon)$ or smaller in §2. Thus, we retain thermal inertia while neglecting momentum inertia, a choice that is motivated by the fact that in magma flow through basaltic fissures, $Pe$ is typically much larger than $\Rey$ (Wylie & Lister, 1995; Wylie et al., 1999). Moreover, as discussed above, convective heat transport plays a key role in the positive feedback between flow rate and temperature perturbations which underlies the fingering instability.

Equation (18) is accompanied by the following non-dimensional inlet and boundary conditions, which differ for the two cases.
Case 1:

We now employ the WRIBL technique to obtain an averaged heat transfer equation. The derivation closely follows that in §2, and so we only present the key steps here. First, the temperature field is split into the wall temperature $T_{w}$ and an internal variation which in turn is decomposed into $\mathcal{O}(1)$ and $\mathcal{O}(\varepsilon)$ components:

The leading order component of the internal variation is defined as:

where $K_{\theta}$ is completely determined in terms of the thin-gap averaged temperature $\theta(z,x,t)$ and the wall temperature $T_{w}(z)$. Thus, the leading order cross-gap temperature profile $\hat{T}(y;\theta)$ arises from a competition between thermal conduction and convection and takes on a parabolic form, parameterized by $\theta$ and $T_{w}$. The $\mathcal{O}(\varepsilon)$ correction in (20) must satisfy the following conditions:

Next, we substitute the temperature decomposition (20) along with the velocity decomposition (5) from §2 into (18) and neglect all terms of $O(\varepsilon\tilde{u})\sim O(\varepsilon^{2})$ and smaller to obtain

In order to derive a closed averaged model, we take the Galerkin weighted residual of (23) using a weight function $F_{\theta}$, defined as

On using integration by parts, along with the boundary and integral conditions in (21), (22) and (24), we exactly close the diffusive term on the right-hand-side of (23) and thereby obtain

It remains to specify the functional dependence of the viscosity on temperature. A simple relationship that has been used in previous work (Helfrich, 1995) is $\mu^{*}(T^{*})=\mu_{0}\exp[\gamma^{*}(T_{H}-T^{*})]$, where the characteristic viscosity $\mu_{0}=\mu^{*}(T_{H})$. In non-dimensional terms, we have

Now, replacing $T$ with its leading order approximation, we obtain the following relationship which couples the averaged flow  (12) and heat transport (25) equations:

Using (27), we solve (6) for $\hat{u}$ and $\hat{w}$, (7) for $\hat{v}$, as well as (11) for the weight function $F_{v}$. We then evaluate the integrals appearing in the averaged equations of flow (12) and heat transport (25) to obtain the coupled WRIBL model in terms of the gap-average temperature and the fluid flow rates:

The inlet condition for the cross-gap averaged temperature is

Equations (28) constitute the coupled WRIBL model. It is distinguished from previously used Darcy models, such as the one described in the next section, by the nonlinear functions $F_{i}$ and $G_{i}$ through which the cross-gap variations in temperature and viscosity enter the problem.

We now briefly describe the Darcy-like model used by Helfrich (1995), which unlike our WRIBL model neglects cross-gap variations in viscosity. The cross-gap averaged flow rates are described by Darcy’s law,

To better appreciate the differences between this Darcy model and our WRIBL model, we shall attempt to reduce the WRIBL model (28) to the Darcy model (29). Firstly, one must neglect the $O(1)$ variation of temperature and therefore of viscosity which leads to $F_{i}=1/12$ and reduces (28a) to (29a). The disregard of cross-gap temperature variations amounts to assuming that $\hat{T}=\theta$ and choosing the Galerkin weight function $F_{\theta}=1$, both of which result in $G_{1}=G_{2}=1$. The convective terms of (28b) then simplify to the convective terms on the left-hand-side of (29b). However, we also obtain $G_{3}=0$ which implies an absence of heat transfer to the side walls; the loss of this key physical ingredient is a direct outcome of assuming the temperature to be uniform across the thin gap. In Helfrich (1995), the Darcy model is salvaged by introducing a heuristic heat loss term, $-\pi^{2}(\theta-T_{w})$, which then yields (29b). This sink term is the conductive heat loss that would occur if the cross-gap temperature profile is that of the leading eigenfunction of the heat diffusion equation, that is, $(\theta-T_{w})(\pi/2){\rm sin}(\pi y)$ [note that $T_{w}=0$ in Helfrich (1995)]

In summary, to obtain the Darcy-like model (29), two asymptotically unjustified simplifications must be made: (a) neglect the leading order, cross-gap variation of temperature, (b) artificially introduce a heat sink term to model the conductive heat loss to the wall. The WRIBL strategy, presented in §2 and §3, allows us to retain all leading order effects and yet derive an averaged model whose computational simplicity is comparable to that of the Darcy-like model.

We evaluate the relative accuracy of the WRIBL model, and the importance of cross-gap temperature and viscosity variations, by comparing the predictions of the WRIBL model with that of the Darcy model of Helfrich (1995). In addition, we solve the the full 3D model, given by (4) and (18)-(19), and consider its results as the ‘truth’ against which both average models are tested. This 3D model was first studied by Wylie & Lister (1995) for the case of constant wall temperature $T_{w}=0$.

We focus on predictions of the uniform base flow state, its multiplicity, and its linear stability. The ordinary differential equations of the two averaged models are integrated using the adaptive-stepping NDSolve subroutine in Mathematica^® v.10.3.. The partial differential equations of the full 3D model are solved numerically by discretizing the thin gap direction using the Chebyshev pseudo-spectral method and using the implicit integration scheme described in Wylie & Lister (1995) for evolving the temperature along the flow direction.

We now proceed to compare the averaged models, first considering a constant cold temperature at the walls (case 1) and next adopting a linearly decreasing wall temperature (case 2). In both cases, we begin with the prediction of the steady and uniform base flow and then move on to its linear stability characteristics.

A thin-gap averaged nonlinear model is presented, using the weighted residual integral boundary layer (WRIBL) method, for flow of a fluid with non-uniform viscosity. The model is applied to the problem of thermoviscous fingering in a Hele-Shaw geometry, wherein the dependence of viscosity on temperature leads to the spontaneous formation of alternate bands of low-viscosity (fast flowing) hot fluid and high viscosity (slow flowing) cold fluid. The cross-gap variation in temperature, which is crucial to the physics of the problem, results in a strongly non-parabolic flow profile across the thin gap. The WRIBL model accounts for these cross-gap variations, which were neglected by the ad hoc Darcy model used in previous work. Comparing the predictions of the two averaged models with fully three-dimensional calculations we find that, although the results are in qualitative agreement, the WRIBL model performs much better than the Darcy model. Therefore, when quantitative accuracy is important, e.g., for comparison with experiments, the WRIBL model ought to be used. Furthermore, the WRIBL model will be crucial when studying issues that depend critically on the shape of the cross-gap velocity profile, such as laminar dispersion of solutes (Datta & Ghosal, 2009) and inertial focusing of particles (Carlo et al., 2007).

The broader implications of this work are twofold. First, it is shown that the averaging of the variable-viscosity transport equation (WRIBL model) is not the same as introducing variable viscosity into an averaged equation (Darcy model). Second, the ideas behind this work go beyond magma flow dynamics in fissures to other thin-gap problems, such as the Hele-Shaw flow of solutions or suspensions, where the viscosity varies due to gradients in solute or particle concentration. In such cases, the energy conservation equation for temperature would have to be replaced by an appropriate transport equation for the solute/particle concentration, which may then be coupled to the WRIBL averaged flow model by specifying the dependence of viscosity on concentration. Moreover, the WRIBL procedure presented here can also be used to develop reduced-order evolution equations for thin film interfacial flows, in which the viscosity varies across the depth of the thin film, e.g., the flow of a thin film of a particle suspension, in which the particles preferentially accumulate at the free-surface due to shear-induced migration (Timberlake & Morris, 2005). Indeed, the WRIBL method was developed for thin film problems and easily accounts for a deforming free surface.

Finally, it is important to note that the primary contribution of this work is the WRIBL averaging of the momentum equations with a variable viscosity, not the averaging of the thermal transport equation. The latter was done as part of the illustrative example of thermoviscous fingering, in order to compare the WRIBL model with the Darcy-like model of Helfrich (1995), wherein both momentum and energy equations were averaged. In general, however, one need not average the transport equation for the scalar that determines the viscosity. Indeed, one could conceive of an iterative/time-stepping solution scheme in which the full transport equation is solved along with the WRIBL averaged flow model. As the transport equation is linear, such an approach would still be far more efficient, computationally, than solving the full set of nonlinear momentum and transport equations. This “partial-averaging” strategy would also be more accurate than a fully averaged model, particularly in regimes where convective transport is important.

D.S.P. and R.N. acknowledge financial support from NASA (NNX17AL27G) and NSF (2025117). R.N. is grateful to the Institute of Advanced Study-Durham University, UK for a fellowship and to Prof. E.W. Llewellin for introducing him to the physics of magma flows. J.R.P. is grateful for funding from the IIT Bombay IRCC Seed Grant. D.S.P. acknowledges support from the IIT Kanpur Initiation Grant, and DST-SERB grant SRG/2020/00242.