2.1 Rheology

For a broad picture of the full range of complex fluids used in hydraulic fracturing, we refer the reader to the recent review of Barbati et al. (Reference Barbati, Desroches, Robisson and McKinley2016). As mentioned in the introduction, we focus on viscous frac fluids as opposed to slickwater slurries, and we target dispersive phenomena in proppant transport. Thus, shear rheology is more important than viscoelasticity for these aspects, although viscoelasticity is ever present in frac fluids. A wide range of fluids are used, according to operation and company, often with proprietary formulation, e.g. typically aqueous polymer gels (guar, hydroxypropyl guar (HPG), etc.), either linear gels or cross-linked (e.g. with borate). These fluids are quite viscous and often strongly shear thinning. A common oilfield characterization is via a power-law model, with consistency $\hat{\unicode[STIX]{x1D705}}_{0}\approx 0.01{-}5~\text{Pa}~\text{s}^{n}$ and power-law index $n\approx 0.2{-}1$ . Some fluids may have a modest yield stress, $\hat{\unicode[STIX]{x1D70F}}_{Y0}\approx 0{-}15~\text{Pa}$ . For a typical frac fluid based on guar, we might expect $0.35<n<0.6$ and $0.1<\hat{\unicode[STIX]{x1D705}}_{0}<1$ $(\text{Pa}~\text{s}^{n})$ . One also sees frac fluids with shear-thinning exponents in the 0.2–0.3 range, often with larger $\hat{\unicode[STIX]{x1D705}}_{0}$ . Rheological parameters are not independently variable. At high shear rates the effective viscosity tends to plateau, e.g. we would not expect the effective viscosity of viscous frac fluids to fall significantly below $0.01~\text{Pa}~\text{s}$ .

To incorporate the range of shear rheologies, we will assume that the effective viscosity of the pure frac fluid is modelled reasonably well as a Herschel–Bulkley fluid, with effective viscosity $\hat{\unicode[STIX]{x1D702}}_{f,0}(\hat{\dot{\unicode[STIX]{x1D6FE}}})$ , possibly bounded below at high shear. Ovarlez and co-workers have recently developed a general framework for shear-thinning and yield stress fluid suspensions that has been validated at least partially (see Ovarlez et al. Reference Ovarlez, Bertrand and Rodts2006, Reference Ovarlez, Bertrand, Coussot and Chateau2012, Reference Ovarlez, Mahaut, Deboeuf, Lenoir, Hormozi and Chateau.2015; Chateau et al. Reference Chateau, Ovarlez and Luu Trung2008; Mahaut et al. Reference Mahaut, Chateau, Coussot and Ovarlez2008; Coussot et al. Reference Coussot, Tocquer, Lanos and Ovarlez2009; Vu et al. Reference Vu, Ovarlez and Chateau2010; Dagois-Bohy et al. Reference Dagois-Bohy, Hormozi, Guazzelli and Pouliquen2015). The bulk suspension viscosity $\hat{\unicode[STIX]{x1D702}}$ , is decomposed as

where $\hat{\unicode[STIX]{x1D702}}_{f}$ is referred to as the liquid-phase viscosity and $\unicode[STIX]{x1D702}_{r}(\unicode[STIX]{x1D719})$ is the dimensionless relative viscosity, modelled for example by the Krieger–Dougherty law

or a close variant (see Maron & Pierce Reference Maron and Pierce1956; Dougherty Reference Dougherty1959; Krieger Reference Krieger1972; Quemada Reference Quemada, Casas-Vazquez and Lebon1982). Apart from conceptual simplicity, there exist generalizations of such laws to particles of different shapes, e.g. rods/fibres (see Wierenga & Philips Reference Wierenga and Philips1998). Here $\unicode[STIX]{x1D719}_{m}$ denotes the maximal packing fraction. The relative viscosity is further decomposed into $\unicode[STIX]{x1D702}_{r}(\unicode[STIX]{x1D719})=\unicode[STIX]{x1D702}_{p}(\unicode[STIX]{x1D719})+1$ , where $\unicode[STIX]{x1D702}_{p}(\unicode[STIX]{x1D719})$ is called the particle-phase viscosity, satisfying $\unicode[STIX]{x1D702}_{p}(\unicode[STIX]{x1D719})\sim \unicode[STIX]{x1D719}$ as $\unicode[STIX]{x1D719}\rightarrow 0$ and $\unicode[STIX]{x1D702}_{p}(\unicode[STIX]{x1D719})\rightarrow \infty$ as $\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}_{m}$ . The combination $\hat{\unicode[STIX]{x1D702}}_{f}\unicode[STIX]{x1D702}_{p}(\unicode[STIX]{x1D719})$ usually appears as part of the solids-phase stress.

In the liquid phase the particles act to amplify effects of a bulk shear rate imposed on the suspension, i.e. since the particles themselves do not deform. This amplification can be crudely estimated as

which is verified reasonably well by experimental results (Ovarlez et al. Reference Ovarlez, Bertrand and Rodts2006; Mahaut et al. Reference Mahaut, Chateau, Coussot and Ovarlez2008; Dagois-Bohy et al. Reference Dagois-Bohy, Hormozi, Guazzelli and Pouliquen2015) and agrees with the theoretical considerations of Chateau et al. (Reference Chateau, Ovarlez and Luu Trung2008). When the strain rate $\hat{\dot{\unicode[STIX]{x1D6FE}}}$ is imposed on the suspension, $\hat{\dot{\unicode[STIX]{x1D6FE}}}_{loc}$ is representative of that felt by the local fluid and particle, and assuming a Herschel–Bulkley type closure,

This can be interpreted as assuming the liquid within the suspension obeys a Herschel–Bulkley type law, but with $\unicode[STIX]{x1D719}$ -dependent consistency and yield stress defined by

When $\unicode[STIX]{x1D702}_{r}(\unicode[STIX]{x1D719})$ is included, the effective viscosity of the suspension increases with $\unicode[STIX]{x1D719}$ , but note that the fluid viscosity ( $\hat{\unicode[STIX]{x1D702}}_{f}$ above) decreases with $\unicode[STIX]{x1D719}$ .

2.2 Dimensional analysis

Figure 1 schematically illustrates a fracture, somewhat distant from the well. Fundamentally, this is a shear flow along a channel of varying width and orientation. We assume three characteristic lengths: width $\hat{D}$ (5–25 mm), length $\hat{L}$ (100–600 m) and height ${\hat{H}}$ (20–100 m). We assume that in-plane variations occur on a length scale $\hat{L}_{t}\lesssim \min \{\hat{L},{\hat{H}}\}$ , where $\hat{L}_{t}\gg \hat{D}$ . On a scale smaller than $\hat{D}$ , one might also include a roughness scale  $\hat{d}_{w}$ .

Figure 1. Schematic of a channel-like fracture. The solid phase, i.e. proppant with volume fraction $\unicode[STIX]{x1D719}_{in}$ , is added at the inlet to the fracturing fluid in a cyclic fashion to produce a network of open channels.

The pumping process is modelled by a representative inflow velocity $\hat{U} _{0}$ ( $0.1{-}3~\text{m}~\text{s}^{-1}$ ), or alternatively a flow rate. Additionally we specify the inflowing proppant solids fraction $\unicode[STIX]{x1D719}_{in}$ , (typically 5 %–40 %). These variables will be time-varying in a complex pumping operation. As the slurry travels along the fracture, some of the frac fluid may leak off into the formation. This leads to higher solid volume fractions in the fracture. Although at times this is a significant effect, for now we focus only on transport along the fracture. Although a range of proppants may be used, common is a 20/40 mesh sand with approximately spherical particles, a representative diameter $\hat{d}_{p}$ in the range 0.42–0.84 mm and density $\hat{\unicode[STIX]{x1D70C}}_{s}=2650~\text{kg}~\text{m}^{-3}$ . Slurry rheology has been discussed in § 2.1. We assume a frac fluid density $\hat{\unicode[STIX]{x1D70C}}_{f}$ ( ${\approx}1000~\text{kg}~\text{m}^{-3}$ ) and a representative viscosity  $\hat{\unicode[STIX]{x1D707}}_{f}$ .

Even with the simplistic description above, a minimum of 11 dimensional parameters arise, depending on three independent dimensions. Thus, we expect eight dimensionless groups (plus $\unicode[STIX]{x1D719}_{in}$ ) to govern these flows. Five of these groups are geometric. The roughness $\unicode[STIX]{x1D6FF}_{w}=\hat{d}_{w}/\hat{D}\ll 1$ we shall neglect for simplicity. On adopting $\hat{D}$ and $\hat{L}_{t}$ as natural length scales for transverse and streamwise directions, we have two groups, $L=\hat{L}/\hat{D}$ and $H={\hat{H}}/\hat{D}$ , that simply indicate the extent of the fracture. More relevant to the transport processes are $\unicode[STIX]{x1D6FF}_{p}=\hat{d}_{p}/\hat{D}$ and $\unicode[STIX]{x1D6FF}_{t}=\hat{D}/\hat{L}_{t}$ i.e. a scaled particle diameter (ratio of particle diameter to fracture width) and the local fracture aspect ratio, respectively. Note that throughout this paper we write all dimensional quantities with a $\hat{\cdot }$ symbol and dimensionless parameters without.

The remaining dimensionless groups include $\unicode[STIX]{x1D719}_{in}$ , the density ratio $s=\hat{\unicode[STIX]{x1D70C}}_{s}/\hat{\unicode[STIX]{x1D70C}}_{f}\approx 2.65$ and two others, namely the densimetric Froude number ( $Fr$ ) and Reynolds number ( $Re$ ):

In characterizing the rheology, we would expect at least two dimensionless groups, e.g.  $n$ and a Bingham number $B_{l}=\hat{\unicode[STIX]{x1D70F}}_{Y0}/[\hat{\unicode[STIX]{x1D705}}_{0}(6\hat{U} _{0}/\hat{D})^{n}]$ . Here, $\hat{\unicode[STIX]{x1D707}}_{f}$ is the effective viscosity of the suspending fluid and the Bingham number $B_{l}$ is based on the nominal laminar wall shear rate, $8\hat{U} _{0}/\hat{D}$ . Further dimensionless groups may arise in characterizing the suspension (e.g.  $\unicode[STIX]{x1D719}_{m}$ ) and in operational parameters.

2.3 Characteristics of fracture flows

In the above two subsections we have indicated the main ranges of slurry rheologies and other dimensional parameters that arise in hydraulic fracturing. In order to model these flows effectively, we need to consider the behaviour at both the bulk scale and the particle scale.

Regarding the bulk flow, pure liquid Reynolds numbers up to a maximum of $(5{-}10)\times 10^{3}$ are possible, but would be reduced via the presence of particles. More typically $Re\sim 10^{2}$ would be common in fracture transport regimes, reducing to zero in cases of screen-out. In considering asymptotic models of flow in long thin geometries, the leading-order inertial terms would be of size $\unicode[STIX]{x1D6FF}_{t}Re$ . Thus, at the outset we deduce that, while non-inertial lubrication-type models have validity in suitably straight fractures, inertia of the bulk flow is important in more complex fractures even at moderate flow rates.

To understand the shear-rate regime with the fluids considered, note that a uniform channel flow of pure frac fluid ( $\unicode[STIX]{x1D719}=0$ ) would have a layered structure consisting of an unyielded plug flow occupying a central fraction $y_{Y}$ of the channel:

( $y_{Y}\in [0,1]$ is readily calculated numerically). Outside of the plug we have yielded shear layers approaching the wall. This simple model can be adapted here by adjusting $B_{l}$ to include the effects of $\unicode[STIX]{x1D719}$ and typical variations in rheology, to give some idea of the range of plug variation. We find that the plug could occupy between 10 % and 85 % of the width. Now we need to consider that the fracture is not uniform. For such fluids, slow geometric variations in the flow direction induce extensional stresses that break the plug. In such situations (e.g. Putz et al. Reference Putz, Frigaard and Martinez2009) we may assume the deviatoric stress remains just above the yield stress with the shear-rate scale determined by geometric variations, i.e. approximately $\hat{\dot{\unicode[STIX]{x1D6FE}}}\gtrapprox \hat{U} _{0}/\hat{L}_{t}$ here. Considering the various effects of extensional stresses, inertial stresses, irregular streamwise variations in geometry and potential transverse settling of dense particles, it is sensible to assume that any predicted plug region is effectively a low-shear pseudo-plug, rather than a rigid plug. Finally, let us note that this type of velocity profile (low-shear-rate pseudo-plug plus sheared wall layers) is found in various models for pressure-driven suspension flows, simply due to particle migration effects and regardless of non-Newtonian effects.

The structure of pseudo-plug and sheared wall layer is helpful in evaluating the local behaviour within the suspension, to which we now turn. First, given the relatively high shear rates and nature of the solids-phase, Brownian and colloidal interactions may be safely ignored. Key interactions are likely to be hydrodynamic and possibly collisional. Within the wide operational ranges of dimensional parameters outlined in § 2.2, we now consider combinations of parameters that produce characteristically different effects in the flow at both bulk and particle scales. For simplicity we categorize these as slow, normal and fast, as defined by the first four rows in table 1.

For each flow category we now use the rheological models outlined in § 2.1 together with (2.7) to estimate: (i) the width of the low-shear pseudo-plug; (ii) characteristic low-shear viscosities within the pseudo-plug; (iii) typical particle Stokes numbers within the pseudo-plug; (iv) minimal viscosities, found near the walls; (v) maximal local shear rates, found near the walls; and (vi) maximal particle Stokes numbers within the wall layers. These parameters are tabulated in the lower rows of table 1. This type of analysis is crude in not explicitly including particle migration effects, but does begin to paint a clearer picture of the flows encountered.

In general, the shear rate is relatively constant across the pseudo-plug and then increases continuously but sharply through the shear layers. Shear thinning together with the local strain-rate amplification of (2.3) result in significant centre-to-wall variations in local particle Stokes number ( $St_{p}$ ) given by

which is interpreted as the ratio of the time scale for viscous drag to affect particle momentum and the characteristic time scale of the suspension: ${\sim}1/\hat{\dot{\unicode[STIX]{x1D6FE}}}_{loc}$ . Note that, since $\unicode[STIX]{x1D6FF}_{p}\ll 1$ , transverse variations in both effective viscosity and shear rate (due to the bulk imposed flow) are felt and evaluated locally on the particle length scale.

The Stokes number gives a measure of fluid–solid coupling. As expected, considering the full range of likely process and geometric conditions, the pseudo-plug regions are completely non-inertial in all flow types: $St_{p}\ll 1$ . The main differences are reflected in the size of the central low-shear region $y_{Y}$ . In moving from the pseudo-plug towards the wall, $St_{p}$ increases due to both the increase in $\hat{\dot{\unicode[STIX]{x1D6FE}}}_{loc}$ and the decrease in $\hat{\unicode[STIX]{x1D702}}_{f}$ . For slow fracture flows we find $St_{p}\lesssim 1$ throughout the fracture, whereas for normal flows $St_{p}\geqslant 1$ for approximately half of the sheared layer, approaching ${\sim}10$ at the walls. For fast flows we would have $1\lesssim St_{p}\lesssim 20$ over approximately half the fracture width, for typical 20/40 proppant.

A representative particle Reynolds number is

which closely follows the variations in $St_{p}$ across the fracture as $3\unicode[STIX]{x03C0}/s\sim O(1)$ . While $St_{p}$ and $Re_{p}$ have similar size, the physical interpretation of these quantities is different. Larger $Re_{p}\geqslant 1$ values imply that we enter a nonlinear regime of drag (viscous to inertial), whereas $St_{p}\geqslant 1$ indicates that the particle velocity is not fully relaxed. The latter implies that the fluctuating components of the particle velocity becomes important. Commonly, this is represented via the granular or particle temperature $\hat{\unicode[STIX]{x1D6E9}}\geqslant 0$ , which is the root-mean-square fluctuating velocity of the particle phase. Temperature $\hat{\unicode[STIX]{x1D6E9}}$ is present throughout the flow, but in Stokesian regimes (such as within the pseudo-plug) $\hat{\unicode[STIX]{x1D6E9}}$ becomes very small, serving primarily to regularize diffusive fluxes (see e.g. Nott & Brady Reference Nott and Brady1994). Particle temperature is more commonly used in kinetic theories of particle dynamics and in granular flows.

2.4 Modelling flow in the fracture

The main modelling challenge is to allow for a broad range of both bulk-scale flow effects (arising from geometry and process variations) and particle-scale effects as we traverse the fracture width. We adopt a continuum approach in which variables are interpreted as being volume-averaged over a suitably chosen local averaging volume, but not time-averaged. (Note that for typical $\unicode[STIX]{x1D6FF}_{p}\approx (2{-}4)\times 10^{-2}$ , we have 25–50 particle diameters across the fracture, so that we are close to the limit of validity of a continuum approach.) Operational choices, cyclical pumping (e.g. CFT), particle migration, the possibility of screen-out, etc. all combine to imply that models for transport along the fracture should consider a full range of $\unicode[STIX]{x1D719}$ , from relatively dilute to concentrated suspensions.

The description emerging is of a Stokesian pseudo-plug layer bounded by shear layers in which inertial and unsteady fluid–particle coupling effects are progressively important towards the walls. Variables specific to solid or liquid phases are phase-averaged, defined using the characteristic function of the phase, the local instantaneous variable (e.g. solids-phase velocity) and a suitable smoothing or weighting function ${\mathcal{G}}$ (see e.g. Drew Reference Drew1983, [). Our dimensional variables will be denoted with a $\hat{\cdot }$ throughout. The suspension mass and momentum balances are

where $\hat{\boldsymbol{u}}$ is the volume-averaged suspension velocity, $\hat{\boldsymbol{b}}$ and $\hat{\unicode[STIX]{x1D72E}}$ are the volume-averaged suspension body force and stress tensor, respectively (see Batchelor Reference Batchelor1970). The suspension stress $\hat{\unicode[STIX]{x1D72E}}$ is often decomposed into individual contributions from both fluid and solid phases:

The material derivative in (2.11) is

These equations are obtained by phase-averaging the individual equations for solid and liquid phases, then summing to eliminate the inter-phase terms. Thus, we also need to consider conservation equations for one phase, taken here as the solids phase:

Here $\hat{\boldsymbol{u}}_{p}$ is the phase-averaged particle velocity and $\unicode[STIX]{x1D719}$ is the solids volume fraction. On the left-hand side of (2.15), $\text{D}_{p}/\text{D}\hat{t}$ denotes the material derivative using $\hat{\boldsymbol{u}}_{p}$ for the advective term. The first term on the right-hand side is the particle stress, i.e. particle contribution to the suspension stress. The second term is the volume-averaged traction on the particle surfaces and the last term is the average body force.

2.5 Developing the suspension balance framework

Diffusive and dispersive effects combine with the averaged forces acting on the solids phase to distribute the particles. The solid-phase mass conservation equation (2.14) is typically manipulated to give a transport equation for evolution of $\unicode[STIX]{x1D719}$ . Note that the phase averaging adopted preserves the instantaneous dynamics of each configuration, so that (2.14) has no diffusive flux. At least two general approaches have been taken to model particle-phase diffusion: (i) the diffusive flux approach of Leighton & Acrivos (Reference Leighton and Acrivos1987), modified by Phillips et al. (Reference Phillips, Armstrong, Brown, Graham and Abbott1992); and (ii) the SBM of Nott & Brady (Reference Nott and Brady1994).

Following the SBM approach here, the relative velocity ( $\hat{\boldsymbol{u}}_{r}=\hat{\boldsymbol{u}}_{s}-\hat{\boldsymbol{u}}_{f}$ ) is substituted into (2.14) to give

where $\hat{M}$ denotes the particle mobility and $\hat{\boldsymbol{f}}_{D}$ is the phase-averaged particle drag force. The usual approach now is to consider (2.15) in the limit of small inertia, whereby $\hat{\boldsymbol{f}}_{D}$ may be substituted into (2.16), leading to the divergence of the remaining terms on the right-hand side of (2.15).

More recently Lhuillier (Reference Lhuillier2009) and Nott et al. (Reference Nott, Guazzelli and Pouliquen2011) have shown that the volume-averaged traction (second term on the right-hand side of (2.15)) may be expressed as the sum of a particle-averaged force and the negative divergence of the particle stress. In this sense, the particle stress does not contribute to the particle momentum equation. The only terms remaining on the right-hand side of (2.15) are the external body force and the particle-averaged force, which consists of the hydrodynamic forces, contact traction forces and interparticle forces. This questions the foundation of the SBM model in attributing the shear-induced particle migration to the divergence of particle-phase stress. However, Nott et al. (Reference Nott, Guazzelli and Pouliquen2011) showed that the particle-averaged force can itself be expressed in terms of an interphase drag force and the divergence of a stress tensor, which includes the effects of hydrodynamic, contact and interparticle forces. Consequently, the form of particle-phase momentum equation used in the SBM approach is correct, but using conventional closures for the interphase drag force and the stress tensor is currently under debate (see Lhuillier Reference Lhuillier2009; Nott et al. Reference Nott, Guazzelli and Pouliquen2011; Brady Reference Brady2015). Thus, we continue below to follow the classical SBM approach in developing our model, but incorporate the local viscosity of the shear-thinning yield stress suspension as outlined in § 2.1, and other necessary extensions.

One advantage of the SBM approach is that the particle temperature is included, which is helpful for modelling non-Stokesian particle effects present in fracture flows, as we do below, and in interpreting the results. We therefore replace (2.15) with the solid momentum equation:

The phase-averaged forces acting on the particles contribute to two terms: $\hat{\boldsymbol{f}}_{B}$ representing the net solids-phase body force and $\hat{\boldsymbol{f}}_{P}$ representing the hydrodynamic forces. The phase-averaged net solids-phase body force is

In this paper we simplify by assuming that the phase-averaged net particle force is given primarily by the drag force, $\hat{\boldsymbol{f}}_{P}\approx \hat{\boldsymbol{f}}_{D}$ , neglecting all other hydrodynamic particle forces. This in turn is modelled via a hindered settling closure:

We adopt the same framework as Ovarlez and co-workers. In Ovarlez et al. (Reference Ovarlez, Bertrand, Coussot and Chateau2012) the authors study particle settling in yield stress fluids, perpendicular to the main direction of shear. They advocate using a Newtonian hindering function $h(\unicode[STIX]{x1D719})$ ,

to modify the single-particle settling speed. They identify two limits (a plastic regime and a viscous regime) in each of which they use a drag law incorporating $\hat{\unicode[STIX]{x1D702}}_{f}(\hat{\dot{\unicode[STIX]{x1D6FE}}},\unicode[STIX]{x1D719})$ , as defined in (2.4). A quite similar usage of $\hat{\unicode[STIX]{x1D702}}_{f}$ for fitting a drag coefficient is described in Tabuteau, Coussot & de Bruyn (Reference Tabuteau, Coussot and de Bruyn2007). The drag coefficient $C_{D}(Re_{p,l})$ is extended to cover inertial regimes, for which there are a number of similar closure laws. There is an advantage to having a drag coefficient closure law that is easily invertible, and consequently we adopt

where $A_{CD}=24/1.4^{0.375}$ . Note that the range of $Re_{p,l}>500$ is unlikely to be attained. The Reynolds number is based on $\hat{\unicode[STIX]{x1D702}}_{f}(\hat{\dot{\unicode[STIX]{x1D6FE}}},\unicode[STIX]{x1D719})$ and $|\hat{\boldsymbol{u}}_{r}|$ , i.e.

Inverting (2.19), with the drag coefficient (2.21), leads straightforwardly to the relation $\hat{\boldsymbol{u}}_{r}=-\hat{M}(\unicode[STIX]{x1D719},\hat{\dot{\unicode[STIX]{x1D6FE}}},|\hat{\boldsymbol{f}}_{D}|)\hat{\boldsymbol{f}}_{D}$ , required for (2.16). We specify a dimensionless version of the mobility $\hat{M}$ later.

Returning to the SBM derivation, it is assumed that the left-hand side of (2.17) is relatively small (which would hold, for example, in conditions where the flow is steady, rectilinear and fully developed). In this case, we may write

and the solids mass conservation equation becomes

The particle stress tensor is usually modelled by the expression

(see e.g. Fang et al. Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002). The second term in (2.25) is the particle shear stress term. The bulk rate-of-strain tensor for the suspension is $\hat{\dot{\unicode[STIX]{x1D738}}}=\hat{\dot{\unicode[STIX]{x1D6FE}}}_{ij}$ . Note that the bulk suspension strain rate is

The first term in (2.25) is a product of the particle pressure $\hat{\unicode[STIX]{x1D6F1}}$ and the tensor $\unicode[STIX]{x1D655}$ , through which we account for normal stress differences. It is argued that a reasonable approximation to the tensor $\unicode[STIX]{x1D655}$ is

where the directions $x_{1}$ and $x_{2}$ are in the plane of shear (here this would be locally aligned with the mean flow along the fracture). According to the scaling, $\unicode[STIX]{x1D706}_{1}+\unicode[STIX]{x1D706}_{2}+\unicode[STIX]{x1D706}_{3}=3$ , and from Brady & Morris (Reference Brady and Morris1997), a common choice for shear flow is $\unicode[STIX]{x1D706}_{1}=\unicode[STIX]{x1D706}_{2}=2\unicode[STIX]{x1D706}_{3}=6/5$ . Other choices can be made if there is specific knowledge of the normal stress differences. The particle pressure is modelled using the particle temperature in place of the strain rate, considering that $\hat{\dot{\unicode[STIX]{x1D6FE}}}^{2}\propto \hat{d}_{p}^{-2}\hat{\unicode[STIX]{x1D6E9}}$ in a homogeneous suspension. The form selected is

An evolution equation for $\hat{\unicode[STIX]{x1D6E9}}$ is derived in Nott & Brady (Reference Nott and Brady1994) as a simplified form of the mechanical energy balance for the particle phase:

There are slight differences due to our preference to work with the particle diameter rather than radius, and a factor of $1/3$ absent in Nott & Brady (Reference Nott and Brady1994), compared to the usual definition of particle temperature. The heat-capacity term $C(\unicode[STIX]{x1D719})$ is not specified, and not used below. The functions $\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D719})$ and $\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D719})$ are specified semi-empirically by considering different limits of the model. The field $\hat{\boldsymbol{q}}$ represents the phase-averaged fluctuating component of the solids-phase dissipation, i.e.

which is modelled via a Fourier-type ‘conductive heat flux’ law,

Although phenomenological in many respects the above model has proven successful in predicting solids-phase distributions in pressure-driven shear flows and many others.

2.6 Unsteady inertial effects

The model outlined above is the standard SBM approach, extended to include non-Newtonian frac fluids (§ 2.1) and steady inertial effects via the nonlinear drag regime. In fracturing flows, within the shear layers approaching the walls, unsteady inertial effects also arise. We now look in more detail at the response of the suspension in the different parts of the flow, with the aim of developing a rational framework that covers these effects.

Following Cassar, Nicolas & Pouliquen (Reference Cassar, Nicolas and Pouliquen2005) and Andreotti, Forterre & Pouliquen (Reference Andreotti, Forterre and Pouliquen2013) we consider a particle that accelerates within the suspension, under the action of a pressure $\hat{P}$ and drag force $\hat{F}_{d}$ . The particle responds on at least three different microscopic time scales related to rearrangement of the suspension. Either $\hat{P}$ is balanced by pure acceleration (leading to $\hat{t}_{1}\sim \hat{\unicode[STIX]{x1D70C}}_{s}\hat{\unicode[STIX]{x1D6E9}}^{1/2}\hat{d}_{p}/\hat{P}$ ), or $\hat{P}$ is balanced by $\hat{F}_{d}$ . The latter leads to $\hat{t}_{2}~({\sim}\hat{\unicode[STIX]{x1D702}}_{f}/\hat{P})$ if the drag is principally viscous, or to $\hat{t}_{3}~({\sim}\hat{d}_{p}\sqrt{\hat{\unicode[STIX]{x1D70C}}_{f}C_{D}/\hat{P}})$ if the drag is principally inertial. When the suspension is itself under shear, there is also a macroscopic time scale arising from the strain rate. We have seen that the bulk strain rate ( $\hat{\dot{\unicode[STIX]{x1D6FE}}}$ ) is felt locally as $\hat{\dot{\unicode[STIX]{x1D6FE}}}_{loc}$ , which is still macroscopic on the particle scale.

Comparison of the macroscopic time scale $(\hat{\dot{\unicode[STIX]{x1D6FE}}}_{loc}^{-1})$ , with each of $\hat{t}_{1}$ – $\hat{t}_{3}$ leads to three dimensionless groups, the relative size of which determines the dominant characteristic behaviour of the suspension, i.e. at different positions across the fracture and in different types of bulk flow. For example, when $\hat{t}_{1}\gtrsim \hat{t}_{2}$ and $\hat{t}_{1}\gtrsim \hat{t}_{3}$ , the rate-limiting microscopic time scale is $\hat{t}_{1}$ and suspension behaviour is described primarily by $I=\hat{\dot{\unicode[STIX]{x1D6FE}}}\hat{t}_{1}$ (case I). In this regime particle inertia dominates. If $\hat{t}_{2}\gtrsim \hat{t}_{1}$ and $\hat{t}_{2}\gtrsim \hat{t}_{3}$ , then $\hat{t}_{2}$ is the rate-limiting microscopic time scale and suspension behaviour is described primarily by $J=\hat{\dot{\unicode[STIX]{x1D6FE}}}\hat{t}_{2}$ (case II: viscous drag regime). Finally, when $\hat{t}_{3}$ is rate-limiting, we expect $J_{I}=\hat{\dot{\unicode[STIX]{x1D6FE}}}\hat{t}_{3}$ to describe the suspension behaviour (case III: inertial drag regime).

These three microscopic regimes are plotted schematically in figure 2. Note that the drag coefficient $C_{D}$ is used on the vertical axis in order to capture both viscous and inertial balances. The fluctuating component of velocity has characteristic size $\hat{\unicode[STIX]{x1D6E9}}^{1/2}$ , which is used to define a (fluctuating) particle Reynolds number:

Using $Re_{p,\unicode[STIX]{x1D6E9}}$ in the usual Stokesian drag law allows us to evaluate $\hat{t}_{1}/\hat{t}_{2}=\sqrt{sRe_{p,\unicode[STIX]{x1D6E9}}/24}\propto \sqrt{St_{p,\unicode[STIX]{x1D6E9}}}$ . One can similarly evaluate $\hat{t}_{1}/\hat{t}_{3}\propto s^{1/2}Re_{p,\unicode[STIX]{x1D6E9}}^{5/16}$ , over the inertial range likely in this process (and eventually $\hat{t}_{1}/\hat{t}_{3}\propto s^{1/2}$ as $Re_{p,\unicode[STIX]{x1D6E9}}\rightarrow \infty$ ).

In each microscopic regime it is possible to derive order-of-magnitude estimates for the key terms in the solids-phase conservation laws, as was done by Jenkins & McTigue (Reference Jenkins, McTigue, Joseph and Schaeffer1990) for the granular regime; see table 2. In table 2, ${\hat{h}}=\hat{d}_{p}[\unicode[STIX]{x1D702}_{r}(\unicode[STIX]{x1D719})/(1-\unicode[STIX]{x1D719})]^{-1/2}$ denotes a representative particle separation, estimated consistently with our earlier treatment of the local strain rate. Particle conductivity and dissipation terms are described later with the transport equation for $\hat{\unicode[STIX]{x1D6E9}}$ (see appendix A).

Figure 2. Variation of the local suspension regimes with slow, normal and fast fracture flows. The arrows indicate the variation across the fracture from the pseudo-plug in the centre to more inertial regimes at the fracture walls (styled after figure 7.13 in Andreotti et al. (Reference Andreotti, Forterre and Pouliquen2013)).

It remains to interpret the different fracture flows in terms of the localized suspension characteristics described above. To do this we balance $\hat{\dot{\unicode[STIX]{x1D6FE}}}_{loc}\sim \hat{\unicode[STIX]{x1D6E9}}^{1/2}/d_{p}$ , to enable us to relate $St_{p,\unicode[STIX]{x1D6E9}}$ to $St_{p}$ . Note that $\hat{\dot{\unicode[STIX]{x1D6FE}}}_{loc}$ is imposed by the macroscopic flow, whereas $\hat{\unicode[STIX]{x1D6E9}}^{1/2}/d_{p}$ is a measure of the local response of the fluctuating particle velocity to this strain rate. Figure 2 shows schematically where the (process-scale) fracture flows ‘live’ in terms of the $(I,J,J_{I})$ cartoon that indicates microscopic regimes. The arrows indicate the variation across the fracture from the pseudo-plug in the centre to more inertial and unsteady regimes at the fracture walls.

This schematic indicates to us how non-Stokesian effects enter as we move across the fracture width, and table 2 suggests the appropriate scaling for extending the different terms in the SBM. Although we take this approach below, it is not the only way to model and include inertial/unsteady effects in the SBM. An alternative approach is to introduce an $O(Re_{p,\unicode[STIX]{x1D6E9}})$ ‘correction’ into the energy equation of the SBM (and into the other closures), drawing from kinetic theory. Closure models in these regimes have been developed by Koch and co-workers (e.g. Sangani et al. Reference Sangani, Mo, Tsao and Koch1996; Wylie, Koch & Ladd Reference Wylie, Koch and Ladd2003), targeted primarily at dense inertial gas suspensions, but also applicable to liquid suspensions with sufficiently large density ratio, $s$ . An application of this style of model to slurry flows in fractures is developed by Eskin & Miller (Reference Eskin and Miller2008). The gas kinetic theory closures are not strictly applicable for the $s$ typical of fracturing flows, but this body of work does indicate the order of $Re_{p,\unicode[STIX]{x1D6E9}}$ correction. We outline this approach and compare with the direct scaling approach in appendix B.

In the Stokesian regime (case II) these scales are well known and those in table 2 are consistent with those of the SBM in Nott & Brady (Reference Nott and Brady1994). Comparing the scalings for case I with those for case II gives us the extension to the SBM closures as we transition from Stokesian to the granular regime. Comparing the scalings for case III with those for case II gives us the extension to the SBM closures as we transition from Stokesian to the inertial hydrodynamic regime. The results of making these comparisons are illustrated in table 3. We observe that (dimensionally) the first column has a prefactor $sRe_{p,\unicode[STIX]{x1D6E9}}$ compared to the Stokesian regime (case II), and the third column has a prefactor $Re_{p,\unicode[STIX]{x1D6E9}}$ . The case III scalings are based on a fully inertial regime (in which the drag coefficient becomes Reynolds-invariant), but in this application we are still in the weakly inertial regime at the particle scale. Typically $s\approx 2.65$ (sand–water), or sometimes $s\approx 3{-}4$ (ceramic proppants), so that the relevant extensions to the Stokesian closures may be regarded as being of order  $Re_{p,\unicode[STIX]{x1D6E9}}$ .

Although we have used table 2 directly, the same order of magnitudes result on using the micro–macro time scale balances outlined above, which has some current popularity. A recent example of this is from Trulsson et al. (Reference Trulsson, Andreotti and Claudin2012), who use a 2D computational approach to model the transition between case I and case II regimes. They show that the $sRe_{p,\unicode[STIX]{x1D6E9}}$ prefactor is multiplied by a constant ( ${\approx}0.6$ for the ranges studied). This suggests that the scaling approach followed here is valid and that we may also expect $O(1)$ constant multipliers to the prefactors listed in table 3. Determination of these multipliers, however, requires extensive further experimentation and/or simulation, which is not our intent here. Instead, we assume constant unit multipliers for simplicity, which leads to the following extensions to the SBM closures:

3.1 Scaled governing equations

The dimensionless suspension mass and momentum equations are:

The terms identified as $CT$ and $SSG$ indicate the next order of terms that enter the equations. Here $CT$ denotes curvature terms, i.e. due to the streamwise variations in the local coordinate system. Similarly, $SSG$ denote the next-order terms in the shear stress gradients.

We now consider the solids-phase momentum balance. The aim is to use the momentum balance to estimate the relative size of the drag force components and hence the relative velocity, which is needed in order to evaluate the solids-phase transport via the mass conservation equation, i.e. in (2.16). The $x$ - and $z$ -momentum equations are

where we have neglected terms of $O(\unicode[STIX]{x1D6FF}_{t})$ related to curvature and terms of $O(\unicode[STIX]{x1D6FF}_{t}^{2})$ related to shear stress gradients. The term $\boldsymbol{f}_{\!D^{\prime }}$ denotes all forces other than drag and buoyancy that act on the particles. Assuming that $\unicode[STIX]{x1D6FF}_{t}Re\lesssim O(1)$ , we see that the leading-order contributions to the drag force are

Although we ignore $\boldsymbol{f}_{\!D^{\prime }}$ in our reduced model later, we include here some discussion, as it is helpful for future modelling. Other than viscous drag and gravitational/buoyancy forces, a number of other forces act on particles immersed in a fluid (see e.g. Maxey & Riley Reference Maxey and Riley1983). These include the Archimedes force, added mass, Basset force, components of lift force and various smaller terms (e.g. Faxen corrections). The added mass and Basset forces are related to fluid and particle accelerations and to the time evolution of the relative motion. Although closure expressions for these forces can be found in the literature, we ignore them as there is no clear directional bias to these contributions. Similarly, we ignore the Magnus component of the lift force, which is related to systematic rotation of solid particles. The leading-order Archimedes force is likely to be transverse to the main flow direction, driven by streamline curvature (see below). There has been considerable work on the lift forces experienced by particles in shear flows at moderate $Re$ (e.g. Hogg Reference Hogg1994; Asmolov Reference Asmolov1999; Asmolov, Lebedeva & Osiptsov Reference Asmolov, Lebedeva and Osiptsov2009; Asmolov & Osiptsov Reference Asmolov and Osiptsov2009; Lebedeva & Asmolov Reference Lebedeva and Asmolov2011), with the last specifically focused on fracturing. These studies do require small particle Reynolds numbers and are aimed at dilute suspensions. Although this limits the accuracy of applying these closures directly, they are valuable in assessing the magnitude of the lift force. In particular, when the principal components of the velocity are in the plane of the fracture, the main lift force components acts transversely, in the $y$ -direction. In summary, the main contributions to $\boldsymbol{f}_{\!D^{\prime }}$ will be in the $y$ -direction, which means that the main contribution to $(f_{D,x},f_{D,z})$ comes from the term $\unicode[STIX]{x1D6FF}_{p}^{2}\unicode[STIX]{x1D719}Re/Fr^{2}$ , aligned with gravity.

The solid-phase $y$ -momentum balance is important in that the particle forces here contribute to the diffusive particle fluxes:

with the same meaning as before, regarding the terms identified as $CT$ and $SSG$ . The terms marked $IT$ in (3.20) are inertial terms of the specified order. At lower order (moved to the right-hand side) is a centrifugal term that may become important in cases where the fracture has streamwise curvature: here the scaled radius of curvature is $r_{\unicode[STIX]{x1D705}}$ and $u_{s,str}$ indicates the solids-phase speed in the local $(x,z)$ -plane, along the fracture. As commented above, the two leading-order terms in $f_{D^{\prime },y}$ that have a consistent directional bias are a transverse Archimedes force and lift force. The transverse Archimedes force is similar to the centrifugal term but involves the fluid density and speed in the local $(x,z)$ -plane. Combining this with the centrifugal term in (3.20) leads to

with relative error of the size of the relative velocity. We denote the scaled lift forces by $f_{L,y}$ and note simply that prescribing $f_{L,y}$ fully is not possible from the current literature. The expressions deriving from Hogg (Reference Hogg1994) and Asmolov (Reference Asmolov1999) and similar relate to Newtonian fluids and dilute regimes. Recent studies by Lashgari et al. (Reference Lashgari, Picanoa, Breugemc and Brandt2014, Reference Lashgari, Picanoa, Breugemc and Brandt2016) show that lift affects the particle distribution only at smaller solid volume fractions (see also Segre & Silberberg (Reference Segre and Silberberg1961)): for $\unicode[STIX]{x1D719}>0.2$ , the particle distribution is governed increasingly by particle–particle interactions and particle-phase stress.

In summary, from (3.20) we find that

In substituting for the relative velocity in (2.16), we will take the divergence of the relative velocity. Owing to the scaling, the $y$ -derivative in the divergence operator is multiplied by $1/\unicode[STIX]{x1D6FF}_{t}$ . Thus, in neglecting small terms in $f_{D,y}$ , to be consistent with the $(x,z)$ -directions we should neglect only the terms of $O(\unicode[STIX]{x1D6FF}_{p}^{2}\unicode[STIX]{x1D6FF}_{t})$ and $O(\unicode[STIX]{x1D6FF}_{p}^{2}\unicode[STIX]{x1D6FF}_{t}^{2}Re)$ . We cannot say whether or not the lift and Archimedes terms are consistently smaller/larger than the stress gradient term. In all likelihood there will be parameter regimes where these forces are significant, e.g. tortuous fractures at high speeds for $f_{A,y}$ and dilute regimes for $f_{L,y}$ . We continue by assuming that the three contributions to $f_{D,y}$ are comparable in (3.22). Note, however, that we have neglected terms that are strictly $O(\unicode[STIX]{x1D6FF}_{p}^{2})$ in estimating $f_{D,x}$ and $f_{D,z}$ , retaining only the term $\unicode[STIX]{x1D6FF}_{p}^{2}\unicode[STIX]{x1D719}Re/Fr^{2}$ . Consistent with this approximation is that the magnitude of the drag force vector does not involve $f_{D,y}$ at leading order, i.e. at leading order we have $|\boldsymbol{f}_{\!D}|=\unicode[STIX]{x1D6FF}_{p}^{2}\unicode[STIX]{x1D719}Re/Fr^{2}$ and consequently

The leading-order relative velocities are now inserted into the solids mass conservation law:

Finally, the granular temperature equation scales as follows:

where neglected terms are lower-order terms in the heat flux (HFT) and viscous dissipation functional (VDT). The convective terms are small and to leading order the balance is simply

As in Nott & Brady (Reference Nott and Brady1994) the heat flux term is largely inactive unless we have large gradients in $\unicode[STIX]{x1D6E9}$ .

3.2 Viable modelling frameworks

The $(x,z)$ -domain has size $L\times H$ . We regard $\unicode[STIX]{x1D6FF}_{t}$ as characteristic of streamwise geometrical variations. The principal equations are (3.13)–(3.16), (3.24) and (3.26), with all terms $CT$ and $SSG$ neglected. Evidently, if we consider all $O(\unicode[STIX]{x1D6FF}_{t})$ terms, we invite suffocating geometric complexity, quite apart from the additional normal stresses and shear stress contributions. Neglecting terms strictly of $O(\unicode[STIX]{x1D6FF}_{t})$ , we see that the Hele-Shaw type of model is recovered as $\unicode[STIX]{x1D6FF}_{t}\rightarrow 0$ at fixed $Re$ . If $Re$ is $O(1)$ , then anyway the Hele-Shaw type of averaging is valid, but in practice many fracturing flows have large $Re$ . Insofar as the bulk suspension flow is concerned, the terms of $O(\unicode[STIX]{x1D6FF}_{t}Re)$ are probably the most important to include at high flow rates. These terms change the form of the momentum balance. Equally, at higher flow rates we would expect larger gradients in shear rates and thus inertial effects to appear at the particle scale in sheared layers nearer the walls. In summary, it appears that two types of model are feasible within this framework, according partly to the fracture geometry, fluid rheology and flow rates.

4.1 Computation of $F(\bar{\unicode[STIX]{x1D719}}_{0},\bar{U}_{0})$

Our leading-order model is based on the conservation law (4.13), which is governed by $F(\bar{\unicode[STIX]{x1D719}}_{0},\bar{U}_{0})$ . In order to construct $F$ , let us assume that $\bar{\unicode[STIX]{x1D719}}_{0}$ and $\bar{U}_{0}$ are given. If we integrate (4.7) once and impose a symmetry condition on $\unicode[STIX]{x1D6F1}_{0}$ at $y=0$ , we see that $\unicode[STIX]{x1D6F1}_{0}$ is constant across the fracture. Thus, we have three variables to find, $U_{0}(y)$ , $\unicode[STIX]{x1D719}_{0}(y)$ $\unicode[STIX]{x1D6E9}_{0}(y)$ , as well as two constants, the pressure gradient $|G|$ and the leading-order solids pressure  $\unicode[STIX]{x1D6F1}_{0}$ .

To solve our nonlinear system of three variables and two constants, we first note that (4.8) relates $\unicode[STIX]{x1D6E9}_{0}$ algebraically to $\unicode[STIX]{x1D719}_{0}$ and $|U_{0}^{\prime }|$ (denoting the $y$ -derivative of $U_{0}(y)$ ). It may be verified that this expression has a unique solution for given $|U_{0}^{\prime }|$ and $\unicode[STIX]{x1D719}_{0}\in [0,\unicode[STIX]{x1D719}_{m})$ . Thus, below, we assume that $\unicode[STIX]{x1D6E9}_{0}$ is a known function of $|U_{0}^{\prime }|$ and $\unicode[STIX]{x1D719}_{0}$ . Since we have $\unicode[STIX]{x1D6F1}_{0}=\text{const.}$ , we set this value from values at $y=0$ , i.e. we use

This uniquely defines $\unicode[STIX]{x1D6F1}_{0}$ in terms of the centreline value of $\unicode[STIX]{x1D719}_{0}$ , i.e.  $\unicode[STIX]{x1D719}_{0}(0)$ . On evaluating $\unicode[STIX]{x1D6F1}_{0}$ in this way, we may then calculate $\unicode[STIX]{x1D719}_{0}(y)$ across the fracture, with the distribution depending only on $\unicode[STIX]{x1D719}_{0}(0)$ and $|U_{0}^{\prime }|$ .

Next, for any $(\unicode[STIX]{x1D719}_{0},\unicode[STIX]{x1D6E9}_{0})$ with $\unicode[STIX]{x1D719}_{0}\in [\!0,\unicode[STIX]{x1D719}_{m}\! [$ , the function $\unicode[STIX]{x1D702}_{f}[1+\unicode[STIX]{x1D702}_{s}]|U_{0}^{\prime }|$ is strictly monotone with respect to the velocity gradient $|U_{0}^{\prime }|$ . This is effectively the flow curve (constitutive law) of the suspension under steady shear. As discussed above, $(\unicode[STIX]{x1D719}_{0},\unicode[STIX]{x1D6E9}_{0})$ are themselves determined by $|U_{0}^{\prime }|$ and $\unicode[STIX]{x1D719}_{0}(0)$

Thus, for given $|G|$ and $\unicode[STIX]{x1D719}_{0}(0)$ , (4.5) uniquely determines the velocity gradient at each $y$ . We now iterate with respect to $\unicode[STIX]{x1D719}_{0}(0)$ at fixed $|G|$ in order to satisfy (4.10). Finally, on integrating out from the wall, $U_{0}(y)$ is determined from $|U_{0}^{\prime }|(y)$ , and a further integration yields the flow rate. The applied pressure gradient $|G|$ is adjusted in order to satisfy the flow rate constraint (4.9).

The various nonlinear equations to be solved are either algebraically solved, or involve monotone continuous functions, for which various robust root-finding procedures can be used. We regularize the relative viscosity to remove the singular behaviour of $\unicode[STIX]{x1D702}_{r}(\unicode[STIX]{x1D719}_{0})$ as $\unicode[STIX]{x1D719}_{0}\rightarrow \unicode[STIX]{x1D719}_{m}^{-}$ . The relative viscosity is regularized as

where $\unicode[STIX]{x1D716}$ is the regularization parameter. The results of this paper are invariant for $\unicode[STIX]{x1D716}\leqslant 10^{-5}$ . The equations are then solved iteratively in the nested fashion indicated in the above description, i.e.  $|G|\rightarrow U_{0}(y)\rightarrow \unicode[STIX]{x1D719}_{0}(0)\rightarrow \unicode[STIX]{x1D719}_{0}(y)\rightarrow \unicode[STIX]{x1D6E9}_{0}(y)$ . Finally, having computed $U_{0}(y)$ and $\unicode[STIX]{x1D719}_{0}(y)$ , we subtract $\bar{U}_{0}$ and $\bar{\unicode[STIX]{x1D719}}_{0}$ , respectively, and integrate (4.14) to give $F(\bar{\unicode[STIX]{x1D719}}_{0},\bar{U}_{0})$ .

Figure 3. Flux function for a range of flow parameters: $B=1,2,4,6,8,10$ , $Re=1,10,100,200,400$ and $n=0.3,0.4,0.5,0.6,0.7,1.0$ . The solid red line corresponds to a representative set of dimensionless parameters: $Re=10$ , $B=1$ , $n=0.5$ , $Fr=1$ , $s=2.65$ , $\unicode[STIX]{x1D6FF}_{p}=0.05$ and $\unicode[STIX]{x1D6FF}_{t}=0.0001$ . The flow profiles are associated with the case of $\bar{\unicode[STIX]{x1D719}}_{0}=0.4$ .

4.2 Proppant propagation

To study proppant propagation behaviour, we solve (4.13) numerically for a given $F(\bar{\unicode[STIX]{x1D719}}_{0})$ under different operational scenarios. This is a first-order quasi-linear hyperbolic equation, for which many numerical methods exist. Here we have used a second-order extension of the staggered Lax–Friedrichs scheme introduced by Nessyahu & Tadmor (Reference Nessyahu and Tadmor1990) with Superbee limiter (see Leveque Reference Leveque2002). The implementation has been tested with a range of test problems, omitted here for brevity. In fact, for the numerical integration, we use an interpolated $F(\bar{\unicode[STIX]{x1D719}}_{0})$ evaluated by solving for $F$ at a number of discrete values. Here we use spacing $\unicode[STIX]{x1D6FF}\bar{\unicode[STIX]{x1D719}}_{0}\approx 0.01$ .