of a thixotropic or antithixotropic fluid in a slowly varying channel : the

A general formulation of the governing equations for the slow, steady, two-dimensional ﬂow of a thixotropic or antithixotropic ﬂuid in a channel of slowly varying width is described. These equations are equivalent to the equations of classical lubrication theory for a Newtonian ﬂuid, but incorporate the evolving microstructure of the ﬂuid, described in terms of a scalar structure parameter. We demonstrate how the lubrication equations can be further simpliﬁed in the weakly advective regime in which the advective Deborah number is comparable to the aspect ratio of the ﬂow, and present illustrative analytical and semi-analytical solutions for particular choices of the constitutive and kinetic laws, including a purely viscous Moore–Mewis–Wagner model and a regularised viscoplastic Houška model. The lubrication results also allow the calibration and validation of cross-sectionally averaged, or otherwise reduced, descriptions of thixotropic channel ﬂow which provide a ﬁrst step towards models of thixotropic ﬂow in porous media, and we employ them to explain why such descriptions may be inadequate.


Introduction
Recent years have seen increasing interest in thixotropic flow. This interest stems both from applications, which include the flow of muds, processed foods, polymer solutions and waxy crude oils, and from the challenge that thixotropic fluids present to the modeller. Typically, the macroscopic rheological properties of such a fluid depend on its microscopic structure (for example, a network of flocculated colloidal particles or a tangle of long-chain polymers [1] ) and thixotropy arises because the microstructure gradually breaks down under shear and rebuilds through Brownian motion. The theoretical modeller is faced with two problems: the rheometric problem of describing this build-up and breakdown, along with the corresponding relationship between the structure and the rheology; and the fluid-dynamical problem of describing the resulting flows.
Most attention has been paid to the rheometric problem. In the simplest models of thixotropic fluids, the state of the microstructure is described by a scalar "structure parameter" λ, which evolves according to an advection-kinetic equation. Many such models have been developed over the last fifty years and calibrated against rheometric data [1,2] . However, less research has been carried out on non-rheometric flows, and it is still uncer-In the flow realised, for example, in cone-and-plate rheometers, the shear rate is uniform in both the streamwise and the transverse directions. Thus, even for a thixotropic fluid, the relationship between shear stress and shear rate can be described in terms of ordinary differential equations. Several studies [20][21][22] have extended this to configurations in which the shear rate may vary in the transverse direction but remains uniform in the streamwise direction: this represents a limiting case of the lubrication regime.
In a preliminary study of such a configuration, Coussot et al. [12] modelled the acceleration of a uniform layer of fluid on an inclined plane in terms of a layer-averaged streamwise velocity and a layer-averaged structure parameter. This further reduction of the equations recovers the simplicity of a purely time-dependent problem, at the cost of the ad hoc assumption that the dynamics are well represented by layer-averaged quantities.
Similar ad hoc reductions have been employed to model flows that were evolving both in the streamwise direction and in time: Chanson et al. [23] considered dam-break flow on an inclined plane, while Pritchard and Pearson [24] considered flow in a narrow fracture, taken to be equivalent to Darcy flow in a porous medium. Both studies reduced the governing equations on the assumption that the rheological state of the fluid in a given crosssection could be characterised by a single quantity: [23] employed a "vertically averaged" value of the structure parameter, while [24] employed a "cross-sectionally averaged" value of the fluidity in a version of Bautista et al.'s [25] model.
The study by Livescu et al. [26] , who considered the levelling of a thin film of thixotropic fluid on a horizontal substrate, represents a bridge between lubrication theory and reduced models. They simplified the governing hydrodynamic equations using a lubrication approximation, then integrated them numerically, and proposed a reduced model based on these numerical results. This approach is an advance on that of [23] and [24] , because it does not postulate in advance that the transverse variation of the structure is known. However, the weakness of this approach is that the transverse variation must be obtained by numerical simulations of a non-reduced system, and there is no guarantee that the approximate profiles for λ obtained in this way will be equally applicable to different rheologies or to different problems.
With this in mind, our goal in this paper is to systematically develop the governing equations for lubrication flow of thixotropic and antithixotropic fluids in a slowly varying geometry. Given the uncertainties involved in the rheological characterisation of thixotropic fluids [2,14] , we will develop this lubrication theory as generally as possible, instead of following most previous studies by restricting our discussion to a specific rheology from the start.
One category of behaviour exhibited by structure-parameter models will not be discussed here, although our approach could in principle be extended to include it. For certain choices of the kinetic model that determines the evolution of λ, even in steady uniform flow λ may have multiple equilibrium values for a given shear rate [20,27] . This non-uniqueness in turn causes nonmonotonicity in the equilibrium stress-strain-rate curve and nonuniqueness of the equilibrium flow profiles. When the structure response time is very short, this behaviour may be described by considering a non-unique stress-strain-rate relation and tracking which branch of this relation applies at each point in the flow. If local flow conditions alter so that a solution on a given branch is no longer available, a "viscosity bifurcation" occurs and the structure is assumed to adjust immediately to another branch. (Here we use the term "viscosity bifurcation" in the sense of Hewitt and Balmforth [27] , who incorporated this behaviour in a model of thin-film flow and tracked the surfaces in the flow where viscosity bifurcations occured.) While non-uniqueness is certainly worthy of further study and may be associated with important physical phenomena such as shear banding [28,29] , we do not regard it as the defining feature of thixotropic flow and so will not discuss it here. Moreover, ruling out non-uniqueness allows us to formulate our leading-order solutions in a convenient analytical form. For similar technical reasons we will not consider true yield-stress behaviour, although we will consider the behaviour of a regularised yieldstress model. We note that although in some materials thixotropy and yield stress are intimately linked phenomena, each may occur without the other [2,30] , so this is also not a fundamental restriction on the present analysis.
In Section 2 we present the governing equations for thixotropic and antithixotropic fluids, and a systematic expansion of these equations for lubrication flow. In the course of this derivation we define an advective Deborah number D, and we show that different regimes may be identified in terms of the relative magnitudes of this Deborah number and the small aspect ratio δ ≪ 1 employed in the lubrication expansion. In Section 3 we specialise to the "weakly advective" regime D = O(δ) , and develop semi-analytical solutions for general constitutive laws and structure evolution laws. In Section 4 we present illustrative results for two thixotropic models: the purely viscous Moore-Mewis-Wagner model and a regularised version of the viscoplastic Houška model. In particular, we discuss the flow profiles across the channel, and consider pressure gradients and pressure drops in channels of specified shape. In Section 5 we investigate the behaviour of a reduced Darcy model for channel flow, and show how lubrication theory can be used both to calibrate such models and to assess their validity. Finally, in Section 6 we summarise our results and discuss directions for the further development of our approach.

Governing equations and boundary conditions
We consider steady, two-dimensional flow of an incompressible thixotropic or antithixotropic fluid at zero Reynolds number. This flow is governed by the continuity equation  More specifically, we consider an ideal thixotropic fluid (in the sense of Larson [31] ) and take the shear stress tensor to be of generalised Newtonian form, for an apparent viscosity ˆ η that depends on both the total shear rate ˙ γ and on the local state of the microstructure represented by λ. The momentum equations (2) thus become and The steady evolution equation for the structure parameter must represent the advection of microstructure along with its build-up Please cite this article as: D. Pritchard August 19, 2016;9:12 ] x yû λ(xŷ)( xŷ) and breakdown; the latter are also taken to depend on both the total shear rate ˙ γ and the local state of the microstructure represented by λ. A general form for such an equation is where for convenience we have defined It will be useful to refer to the equilibrium structure parameter, For simplicity of presentation we consider channel flow between symmetric impermeable walls, and choose our co-ordinate system such that ˆ y = 0 is the centreline of the channel, and the walls are at ˆ Fig. 1 . (Note that, because the corrections associated with a non-straight channel centreline enter at the same order as the other geometrical corrections to lubrication theory, the following analysis also applies to non-symmetric channels.) We assume symmetry of ˆ u and λ, and antisymmetry of ˆ v , about the channel centreline ˆ y = 0 . Working in the lower half-channel where ∂ ˆ u /∂ ˆ y > 0 , we may then impose nopenetration and no-slip boundary conditions on the lower wall and symmetry conditions on the centreline, In general we also require an upstream boundary condition on λ; we will discuss this point in more detail below. Finally, we require a boundary condition related to the pressure gradient that drives the flow. This gradient is not necessarily the same at each streamwise location ˆ x ; rather, it is natural to specify the volume flux, ˆ Q , along the channel, noting that by continuity this must be the same at each cross-section, so At each streamwise location ˆ x , the pressure gradient must be determined so that it is consistent with the flux condition (9) . If a comparison is then required with experiments in which the net pressure drop between two cross-sections, such as the ends of a pipe, has been specified, then this may be achieved by adjusting ˆ Q until the required pressure drop is achieved. (We will return briefly to the question of pressure drops in Section 4.4 .)

Non-dimensionalisation and lubrication scalings
We take the typical width of the channel to be ˆ H , so quantities vary over a typical distance ˆ H in the transverse ˆ y -direction, while in the streamwise ˆ x -direction they vary over a typical distance ˆ H /δ, where the small parameter δ ≪ 1 is the aspect ratio of the flow.
As in classical lubrication theory, our approach will be based on an asymptotic expansion in the limit δ → 0.
We define the dimensionless quantities where we have written the general dimensional viscosity function as ˆ η( ˙ γ , λ) = ˆ μ 0 η(Ŵ, λ) (11) in which ˆ μ 0 is a dimensional viscosity parameter. As usual, these nondimensionalised quantities are implicitly assumed generally to be of order unity. With the scalings (10) , the governing hydrodynamic equations (1) , (4) and (5) become ∂u ∂x ∂ p ∂y = δ 2 2 ∂ ∂y η ∂v ∂y while from (7) Note that only even powers of δ appear in these expressions. The significance of this observation in the present context will become clear in Section 2.4 . In nondimensional form the boundary conditions (8) are simply and η ∂u ∂y while the flux condition (9) becomes To rescale the structure evolution equation, we assume that we may write  August 19, 2016;9:12 ] where the constant ˆ f 0 has dimensions of inverse time, and where the dimensionless structure evolution rate f ( Ŵ, λ) is assumed to be of order unity when Ŵ and λ are of order unity. With the scalings (10) and (18) , the structure evolution equation (6) becomes where the advective Deborah number D is defined by

General expansion scheme
To simplify the equations, we follow the usual procedure and expand all variables in powers of δ, We will pursue the expansions only to first order, although higher orders could readily be developed if necessary. Note that, from equation (14) , the pressure is independent of y up to and including O(δ) .

Regimes of the Deborah number
From equation (31) , it is evident that the behaviour of the fluid depends on the relative magnitudes of the Deborah number D and the aspect ratio δ. Several regimes are possible.
When the Deborah number is very small, specifically when D = O(δ 2 ) or smaller, the effect of the advection of the structure parameter is of the same order as, or smaller than, the neglected O(δ 2 ) terms in classical lubrication theory, and so in this "very weakly advective" regime the fluid simply behaves like a generalised Newtonian fluid to within the usual accuracy of lubrication theory. When the Deborah number is somewhat larger, but still small, specifically when D = O(δ) , to leading order the structure of the fluid is determined by a local balance between build-up and breakdown, f ( Ŵ, λ) ≈ 0, and the effect of the advection of the structure parameter enters as a correction at O(δ) . We will refer to this as the "weakly advective" regime, and consider it in detail in Section 3 below. When the Deborah number is large, but not very large, specifically when D = O(1 /δ) , the opposite occurs: to leading order the structure parameter is simply advected downstream, and the effects of build-up and breakdown enter at O(δ) . Consequently, the structure parameter evolves over long dimensionless streamwise distances of O(1 /δ) and for a full treatment a multiple-scales analysis is required. We will refer to this as the "strongly advective" regime. When the Deborah number is very large, specifically when D = O(1 /δ 2 ) or larger, the effects of buildup and breakdown are of the same order as, or smaller than, the neglected terms in classical lubrication theory, and so in this "very strongly advective" regime the structure parameter is entirely de- Deborah number in their asymptotics implicitly treat this regime, and so the present work provides a natural complement to these studies.

Estimates of typical Deborah numbers
Although the emphasis of the present study is on developing a general approach to thixotropic lubrication flow rather than on modelling specific fluids or experiments, it is useful to consider briefly the conditions under which the theory may be applicable. The crucial parameter is the advective Deborah number D given by (20) which may be interpreted as D = (response time of structure) · (typical velocity) (streamwise distance) .
We may obtain rough estimates for D based on some recent experimental studies. The most elusive quantity is the typical response time of the fluid microstructure. The figures presented by Coussot et al. [12] and Huynh et al. [21] suggest that the clay suspensions with which they worked had response times of the order of 10 s to 10 0 s, although these authors did not report direct estimates of these times. Dullaert and Mewis [32] reported response times ranging from about 0.1 s to 10 s for a suspension of fumed silica particles. Boek et al. [33] found response times of roughly 1 s for a wormlike micellar solution. Ardakani et al. [34] reported response times of the order of 10 s for toothpaste. Finally, Wachs et al. [15] refer to the structure of waxy crude oils breaking down over a timescale of seconds but requiring hours to build up.
Estimates of speeds and distances are more readily available, though these are highly contingent on the experimental setting. For example, the mud dam-break releases of Chanson et al. [23] travelled distances of the order of 1 m, with flow depths generally in the range 0.01 m to 0.1 m and maximum speeds of the order of 0.1 m s −1 . In the toothpaste simulations of Ardakani et al. [34] , channel diameters varied from 0.015 m to about 0.25 m, while channel lengths were typically about 20 times the maximum diameter. Shear rates as high as U/R = 640 s −1 are quoted, where R is the minimum radius, and these correspond to velocities of around 5 m s −1 , although one might guess that domestic toothpaste dispensing occurs at rather lower speeds. The diameters of oil pipelines are usually of the order of 10 −1 m to 1 m, although some may be rather larger or smaller, while flow speeds are typically a few metres per second; distances travelled can be many kilometres, though these may greatly exceed the distances over which the diameter of the pipeline changes. Finally, as a prototype for flow in a permeable rock, we may consider a fracture of typical width 10 −4 m to 10 −3 m, varying over streamwise distances of around 10 −2 m, and Darcy velocities varying from 10 −5 m s −1 (as considered by Boek et al. [33] ) up to perhaps 10 −3 m s −1 . In summary, we may conclude that all regimes of the Deborah number relative to the aspect ratio described in Section 2.4 may occur in practice in particular applications. In this study, however, we will concentrate only on the weakly advective regime, in which further simplification of the governing equations is possible.

General solutions in the weakly advective regime, D = O(δ)
We now focus on the weakly advective regime in which D = O(δ) and hence the convective derivative in (31) is smaller than the structure evolution term by a factor of O(δ) . Accordingly, we

The problem at O(1) : generalised Newtonian behaviour
The problem at O(1) is given by to be solved subject to the conditions Eqs. (36) -(40) constitute a standard problem for the flow of a generalised Newtonian fluid [35,Section 8.4] , with a constitutive equation defined implicitly by f (Ŵ 0 , λ 0 ) = 0 , in a slowly varying channel. As we might expect in this regime of small Deborah number, at leading order there are no thixotropic effects, and the structure parameter departs from its equilibrium value only at O(δ) .
The details of the O(1) solution are somewhat intricate, and are given in Appendix A : the only point that need concern us here is that these solutions can be computed for general forms of the functions f, η and h , as long as certain inverse functions exist. It is this invertibility condition that rules out non-monotonic stress-strain-rate curves of the kind considered by [20] and [27] . In principle it also rules out yield-stress behaviour, since for a yieldstress fluid the stress-strain-rate relation becomes non-unique below the yield stress. However, our approach can readily accommodate a regularised yield-stress model: we discuss one of these in Section 4.3 .

The problem at O(δ) : thixotropic effects
The problem at O(δ) is given by  August 19, 2016;9:12 ] to be solved subject to the conditions We can rearrange (43) to give λ 1 in terms of ∂ u 1 / ∂ y and a known function of x and y determined by the leading-order solution, namely Substituting this into (42) yields where for convenience we have defined We now have a linear boundary-value problem for the velocity perturbation u 1 ( x, y ) at each value of x , with the pressure gradient perturbation G 1 ( x ) acting as an eigenvalue.
Integrating (47) once and applying the boundary condition at y = 0 , we obtain and thus Substituting this expression into the flux condition (45) and integrating we obtain and having derived this expression we can substitute it back into (51) to obtain u 1 ( x, y ) in terms of previously computed quantities.
Once the variables of the O(1) solution and their streamwise derivatives have been calculated (see Appendix A ), the O(δ) solution can be obtained by evaluating these integrals using numerical quadrature. The results presented in the following sections were obtained by this means using the computer algebra package Maple 18, and the implementation was validated against the explicit solutions presented in Section 4.1 .
We note that in this regime no upstream boundary condition can be imposed upon λ, because at O(1) the structure is completely determined by the local flow conditions, and the changes to the structure at O(δ) reflect only changes in the leading-order quantities. In this sense, the weakly advective regime is a singular perturbation of the full system, although the singularity is not severe. It means, however, that to apply our theory in the vicinity of an upstream boundary, such as an inlet with conditions prescribed upon the structure, additional care and analysis would be required.

Illustrative results for specific rheological models
To illustrate the general approach described above, and to gain insight into the generic behaviour of thixotropic and antithixotropic lubrication flow, we now present some solutions for specific models. First we present results in a specific case of the Moore-Mewis-Wagner model ( Section 4.1 ) in which the O(1) and O(δ) solutions may be written explicitly; the rheological model used in this case suffers from pathological behaviour at zero shear rate, but the results are nevertheless informative and provide a useful benchmark for the more general solution method. Then we present equivalent results for the full Moore-Mewis-Wagner model ( Section 4.2 ) and a regularised version of the Houška model ( Section 4.3 ). Finally , we examine the pressure gradients and pressure drops along channels of specified shape ( Section 4.4 ).

The Moore-Mewis-Wagner model with d = 0
The first specific choice of constitutive law and evolution equation which we will use to illustrate our approach is the so-called Moore-Mewis-Wagner (MMW) model. The MMW model has recently been investigated by McArdle et al. [36] in the context of unsteady rectilinear flow, where it exhibits a variety of thixotropic and antithixotropic behaviours depending on the values of the exponents a, b, c and d .
The constitutive law is a version of that introduced by Moore where ˆ η 0 is a dimensional viscosity parameter.
In conjunction with the constitutive law (53) we will employ the structure evolution model presented by Mewis and Wagner [2] , which contains many previous models as special cases (see their Table 3). In this model, the structure evolution rate is given by where ˆ k 1 and ˆ k 2 are dimensional constants which control the breakdown and build-up rates respectively, while a, b, c and d are non-negative dimensionless exponents. In the limiting case d = 0 , the build-up rate is independent of λ, which may therefore increase unboundedly and may take any non-negative value; this is the case, for example, in the model used by Coussot et al. [12] and Liu and Zhu [22] . In the more general case d > 0, the value of λ is restricted to lie between 0 and 1, as in most thixotropic models [2] .
For the MMW model we take ˆ μ 0 = ˆ η 0 and η = λ in (11) . In ad- In equilibrium, the nondimensionalised form of the structure equation (55) yields and we conclude that λ eq ( Ŵ) is uniquely defined if b > 0. We implicitly take κ = O(1) ; if κ ≫ 1 or κ ≪ 1, as when build-up times are much longer than breakdown times [15] , then (56) predicts that λ may become asymptotically large or small and, depending on the constitutive law employed, some further rescaling may be required [36] .
With the particular parameter choice d = 0 , the equilibrium rheology η = λ eq , where λ eq is defined by (56) , is simply that of a power-law fluid, At O(1) we recover the classical solution for channel flow of a power-law fluid [39] , where for brevity we have introduced the notation At O(δ) the perturbation to the streamwise pressure gradient is The perturbations to the streamwise velocity and the structure parameter are and The first-order velocity perturbation u 1 given by (63) is singular at y = 0 when (2 n − c) /n < 0 , i.e. when n < c /2. When (2 n − c) /n ≤ −1 , i.e. when n ≤ c /3, this singularity is non-integrable, giving an unbounded perturbation to the flux. These singularities reflect the breakdown of the power-law model at low shear rates [40] .
Despite this singularity in λ 0 , for the range of parameters chosen in Fig. 2 the velocity perturbation u 1 remains finite (see Fig. 2 (b)). For thixotropic fluids, u 1 is negative close to the walls and positive in the centre of the channel, whereas for antithixotropic fluids the converse is true. As noted above, the structure parameter perturbation λ 1 (see Fig. 2 (d)) unfortunately inherits the singular behaviour of λ 0 , except in the most strongly antithixotropic case plotted (for which n = 1 . 5 and 2 n − c = 2 , so all powers of | y | in (64) are non-negative). Apart from this singular behaviour, the crucial feature is that for thixotropic fluids the structure parameter perturbation is negative, λ 1 < 0, whereas for antithixotropic fluids it is positive, λ 1 > 0. We note also that for thixotropic fluids the pressure gradient perturbation is negative, G 1 < 0, whereas for antithixotropic fluids it is positive, G 1 > 0.
To explain the form of these perturbations physically, it is sufficient to consider an expanding channel, h ′ > 0, since all perturbation quantities are proportional to the local value of h ′ . We may also discuss only a thixotropic fluid, n < 1, since the explanation for an antithixotropic fluid is simply the converse of that for a thixotropic fluid.
In an expanding channel, the shear rate is higher upstream and lower downstream, so the microstructure of the fluid tends to be more broken down upstream. Since the microstructure is advected with the fluid, this broken-down structure is carried downstream by the flow. The result of this is that, at any location, the thixotropic fluid is less structured than the corresponding shearthinning generalised Newtonian fluid would be, so its apparent viscosity is lower. (In terms of our asymptotic expansion, this corresponds to the condition λ 1 < 0, apparent in Fig. 2 (d).) We also note from Fig. 2 (d) that the reduction in viscosity is more pronounced (i.e. | λ 1 | is largest) near the centre of the channel where the rate of downstream advection is highest, and the reduction in viscosity is least pronounced near the walls where advection is lowest.
The velocity perturbation u 1 must reflect both the changes in the viscosity due to thixotropy and the requirement that the net flux is unchanged. In particular, unless u 1 is identically zero, the flux condition requires that u 1 should be positive in some regions and negative in others. Near the centre of the channel, the larger reduction in the viscosity due to thixotropy makes the fluid easier to shear, so it is in this region that the fluid moves faster ( u 1 > 0). This faster flow near the centre of the channel must be compensated for by slower flow near the channel walls ( u 1 < 0). This is indeed the pattern that can be observed in Fig. 2 (b).
Finally, because u 1 is negative near the channel walls, the velocity gradients and the corresponding viscous shear stresses at the walls are reduced. Since these viscous stresses must be balanced by the driving force exerted by the streamwise pressure gradient, the effect of thixotropy is to reduce the magnitude of this pressure gradient. In terms of our asymptotic expansion, this corresponds to G 1 < 0.
This qualitative argument relies only on whether the fluid is thixotropic or antithixotropic; it does not rely on the fine details of the choice of parameters, although these do have a quantitative effect. Fig. 3  so that n = 0 . 6 in each case and thus the leading-order solutions are identical. Although there are some differences because different values of these exponents make the structure parameter more or less sensitive to changes in the shear rate, the overall pattern is unchanged.

The Moore-Mewis-Wagner model with d > 0
A natural question, which arises from the pathological centreline behaviour of the explicit solutions with d = 0 described in Section 4.1 , is how the solutions will change if d > 0 so that the structure parameter λ is constrained to lie between 0 and 1. This Restricting the value of the structure parameter by setting d > 0 makes only minor differences to the O(1) velocity profiles ( Figs. 4 (a) and (b)). The effect is felt only near the centreline, where it makes the signatures of shear-thinning and shearthickening behaviour less pronounced: the velocity profile for a thixotropic fluid is less plug-like than the corresponding shearthinning power-law profile ( Fig. 4 (a)), while that for an antithixotropic fluid is less angular than the corresponding shearthickening power-law profile ( Fig. 4 (b)). In contrast, the O(1)  ( Fig. 4 (e)), while the converse holds for antithixotropic fluids ( Fig. 4 (f)). Fi-nally, although the magnitudes of the structure-parameter perturbations change dramatically when d > 0 ( Figs. 4 (g) and (h)), the signs of these perturbations are unaltered. As might be expected, the perturbations to the structure parameter remain largest in the centre of the channel, where advection is highest. Although u 1 remains finite at the centreline, there is a weak singularity in λ 1 at the centreline, which is discussed in Appendix B .

The regularised Houška model
The rheological model introduced by Houška [41] has become a popular and relatively tractable description of waxy crude oils [14,15,18,43] . It is a thixotropic yield-stress model which comprises a Herschel-Bulkley constitutive law [11] with a structuredependent yield stress τ y and consistency parameter η H , coupled to a structure evolution equation which is a special case of the Mewis-Wagner model (54) . The constitutive law may be written as where the yield stress ˆ τ y (λ) and the viscosity ˆ η H (λ) are linear functions of λ, given by ˆ τ y (λ) = ˆ τ y0 + λˆ τ y1 and ˆ η H (λ) = ˆ η H0 + λˆ η H1 .
The structure evolution rate is given by setting b = 1 , c = 0 and d = 1 in (54) , to obtain The advantages of the Houška model, which it inherits from the Bingham and Herschel-Bulkley models, are its relative simplicity and the ease with which it may be fitted to experimental data. However, it also inherits some disadvantages of these models: the constitutive relation is non-differentiable at the yield stress ˆ τ = ˆ τ y , and the equilibrium stress-strain-rate relation is not well-defined at ˙ γ = 0 . This makes it impossible to apply our semianalytical approach because not all the required inverse functions exist. To circumvent the mathematical difficulties of introducing a genuine yield stress, we therefore regularise [42] the constitutive law (65) to where ˆ k is a parameter which is chosen to be large relative to the reciprocal of typical shear rates in the fluid. We note that the regularisation of yield-stress models may be surprisingly problematic [10] , and that particular caution would be required if this regularised model were used to track plugs or pseudo-plugs in the flow. For simplicity, we follow Wachs et al. [15] and consider only the case n = 1 , so for any constant value of λ the fluid has a Bingham rheology [35,Section 2.3] .
Non-dimensionalising (68) using the scale ˆ μ 0 = ˆ η H0 , and nondimensionalising (67) in the same way as (54) , we obtain Note that in the limit k → ∞ yield-stress behaviour is recovered; we will consider finite but numerically large values of k so that the regularised model approximates the unregularised Houška model (65) . The equilibrium structure parameter is given by (72) Fig. 5 illustrates the effect of the regularisation on the equilibrium apparent viscosity ( Fig. 5 (a)) and the equilibrium shear stress ( Fig. 5 (b)). Note that the deviation from the pure yield-stress behaviour of the unregularised model is very small for k = 10 0 , and almost imperceptible for k = 10 0 0 . Fig. 6 illustrates the leading-order solutions and perturbations for the streamwise velocity and the structure parameter for the regularised Houška model (compare with Figs. 4 (a), (c), (e) and (g) for the thixotropic MMW model). All dimensionless parameters have been set equal to unity so that all the physical effects (viscosity, yield stress and the dependences of both on λ) come into play.
The leading-order velocity profile ( Fig. 6 (a)) is practically independent of the regularisation parameter k , and very similar to the corresponding profile for pipe flow [43] . There is a pseudoplug region in the centre of the channel, which becomes more exactly plug-like as k increases. The corresponding leading-order result for the structure parameter ( Fig. 6 (c)) demonstrates, as one might expect, that the fluid is most broken down near the walls and most built up in the centre. The effect of the regularisation on the pseudo-plug is much more evident here than in the velocity profile; it is clear that as k → ∞ the fluid becomes completely structured, λ = 1 , throughout the plug.
The velocity perturbation u 1 ( Fig. 6 (b)) is, as for the MMW model, rather small ( Fig. 4 (e)). The pattern is the same as in the other thixotropic cases we have seen, with u 1 < 0 close to the walls and u 1 > 0 near the centreline; throughout the pseudo-plug region, the velocity perturbation is effectively constant. Finally, the structure parameter perturbation λ 1 ( Fig. 6 (d)) is more sensitive to the regularisation than the velocity perturbation. It follows a similar pattern to that for the MMW model ( Fig. 4 (g)), with the exception that λ 1 ≈ 0 in the pseudo-plug region, so the fluid remains fully structured here. The physical explanation for the shape of these perturbations is the same as that for the MMW fluid presented in Section 4.1 , with the exception that it is now the  regions at the edges of the pseudo-plug, rather than on the centreline, where the viscosity is most reduced and the shear rate is thus most readily increased.

Pressure gradients and pressure drops
In many practical situations, the details of the flow within a channel are not readily measured, and are in any case of relatively little practical interest: the key physical quantity is the pressure gradient required to drive a given flux along a section of a channel, or the pressure drop required to drive a given flux along a channel of known geometry and length. The present approach allows these quantities to be calculated up to O(δ) . An important point is that both the leading-order pressure gradient G 0 and the first-order perturbation G 1 depend in a complicated and generally nonlinear manner on the local channel width h , while G 1 is also proportional to h ′ . Fig. 7 (a) illustrates the dependence of the pressure gradients on h for a thixotropic case of the MMW model ( a = 1 . 4 , d = 1 ; compare with the solid lines in Figs. 4 (a), (c), (e) and (g)). Unsurprisingly, as the channel becomes wider the pressure gradient required to drive a given flux decreases, because lower velocities and lower shear rates are required. As h becomes small, so the fluid is highly sheared, the upper limit on λ becomes irrelevant and the behaviour of both G 0 and | G 1 | is asymptotic to that of the case when d = 0 ( Eqs. (60) and (62) ). As h becomes large, on the other hand, the fluid is less sheared and more completely structured, λ → 1; in this limit, G 0 asymptotically approaches its value for a Newtonian fluid, G 0 = 12 /h 3 , while | G 1 | decays much faster than the model with d = 0 would predict. (For these parameter values, the model with In contrast to the case where d = 0 , the ratio G 1 / G 0 remains of roughly the same magnitude as h varies ( Fig. 7 (b)), although it does vary by a factor of about two over the range of h considered.
For a channel with a given width h ( x ), it is straightforward to construct the corresponding contributions to the pressure gradient at O(1) and O(δ) . Fig. 8   and (e)), the channel width varies sinusoidally, as in the standard "wavy-wall" test case [9,10] ; in the right column ( Figs. 8 (b), (d) and (f)), it has a more "saw-tooth" form, with more rapid decrease of h towards the constriction and more gradual increase of h afterward ( Figs. 8 (a) and (b), respectively). In each case, the leading-order pressure gradient G 0 is positive ( Figs. 8 (c) and (d)), and it is highest in the constriction, where the resistance to flow is greatest; the wide range of variation of G 0 reflects the strongly nonlinear relation between G 0 and h ( Fig. 7 (a)). The pressure gradient perturbation G 1 changes sign when h ′ changes sign ( Figs. 8 (e) and (f)), so that just upstream of the constriction (a contracting channel) the pressure gradient is increased, whereas just downstream of the constriction (an expanding channel) the pressure gradient is decreased. The asymmetrical profile of the "saw-tooth" channel is reflected in the asymmetry of G 1 ( Fig. 8 (f)): the region where G 1 > 0 is smaller than that where G 1 < 0, but the maximum of G 1 is greater in magnitude than the minimum.
It is also interesting to consider the total pressure drop p along the channel. We may write where L is the length of the channel ( L = 2 π for the cases plotted in Fig. 8 ). For the sinusoidally varying channel ( Figs. 8 (a), (c) and (e)), p 0 ≈ 104.51, whereas for the "saw-tooth" channel ( Figs. 8 (b), (d) and (f)), p 0 ≈ 85.846. The lower pressure drop in the "saw-tooth" case reflects the fact that h ( x ) is more rapidly varying around its minimum, so the section that is most strongly constricted is slightly shorter: this can be seen in the difference between the widths of the peaks in G 0 in Figs. 8 (c) and (d). Meanwhile, the net pressure drop associated with the perturbation, p 1 , is zero in each case. This last fact is not an artefact of the channel shape, but a consequence of the fact that G 1 is linear in h ′ . Writing In other words, the thixotropic contribution to the pressure drop along the channel depends only on the difference between the values of the function ˜ G (h ) at the ends of the channel, regardless of the behaviour of h ( x ) between these end points. Since the channels plotted in Fig. 8 satisfy h (L ) = h (0) , this contribution is identically zero. More generally, it is straightforward, having constructed a plot such as Fig. 7 (a), to evaluate ˜ G (h ) by quadrature and thus to calculate the thixotropic contribution to the pressure drop along any given channel without calculating any more specific details of the flow. This contrasts strikingly with the task of evaluating p 0 , which remains a hard problem requiring detailed knowledge of the geometry of the channel. It does, however, mean that the inverse problem of determining the flux ˆ Q that will be driven through a channel by a given pressure drop ˆ p is not significantly harder for weakly advective thixotropic flow than it is for the flow of a generalised Newtonian fluid.

Comparison of lubrication theory with a reduced model
As discussed in Section 1 , several previous studies have employed formally or informally reduced versions of the lubrication equations, typically attempting to describe the flow in terms of variables averaged in the direction transverse to the flow [23,24] , or evaluated at a free surface [26] . For practical applications such models are significantly simpler than those that resolve the transverse variation, and they may also provide a basis for generalisations to more complicated scenarios such as flow in porous media. However, it is not apparent a priori whether such models will necessarily agree, or can be calibrated to agree, with the non-reduced lubrication theory.
For channel flow, the key output of a model is a relationship between the flux ˆ Q , the applied pressure gradient ˆ G , and the channel width ˆ h ; using the nondimensionalisation in Section 2.2 , this corresponds to the relationship between the pressure gradient terms G 0 and G 1 and the channel width described in terms of h and h ′ . The prototype for such relationships is the Newtonian relationship familiar from studies of Hele-Shaw flow [44] , The minimal requirement that we can make of a reduced model of thixotropic flow is that in the weakly advective regime it captures both the dependence of G 0 on h and the dependence of G 1 on h and h ′ to reasonable accuracy; in other words, that it captures the non-Newtonian nature of the leading-order flow and the first-order correction to this caused by thixotropy. This is not a sufficient condition for validity, but if a reduced model fails this test then its claim to capture the behaviour of the flow is at best questionable. It will also be helpful to compare the quantities that appear in reduced models with the cross-channel averages of the quantities that appear in lubrication theory. We will denote the latter by an overbar, so that in particular (76)

A reduced Darcy model
The simplest way to represent thixotropic effects on pressuredriven channel flow is the reduced ("coarse") Darcy model employed by Pritchard and Pearson [24] ; the approach employed for free-surface flow by Chanson et al. [23] is conceptually rather similar. In this model, the flow is described in terms of representative values at each cross-section for the structure parameter, ( ˆ x ) , and for the shear rate, ˆ ( ˆ x ) . The representative structure parameter then evolves downstream according to where ˆ K is an effective permeability which we assume takes the form ˆ K = βˆ h 2 for some dimensionless coefficient β; for a Newtonian fluid, β = 1 / 12 . Finally, we assume on dimensional grounds that the representative velocity and the representative shear rate can be related to the flux and to the channel width by for some dimensionless coefficients ω and σ . Nondimensionalising this Darcy model using the scales defined in Section 2.2 , we obtain where We now consider this model in the weakly advective regime, in which we can compare the results at each order in δ directly with those from the lubrication expansion. We define D = δD * as before and write At O(δ) , we have where f , f and η represent derivatives evaluated at ( , 0 ). As Figs. 9 (a) and (b) indicate, it is easy to calibrate the Darcy model so that the leading-order pressure gradient G 0 is predicted fairly well and the "typical" structure parameter 0 in the Darcy model agrees fairly closely with the channel-averaged structure parameter λ 0 from lubrication theory. However, as Figs. 9 (c) and (d) indicate, the Darcy model is not capable of predicting accurately how the pressure gradient perturbation G 1 varies with h , and neither can the "typical" structure parameter perturbation 1 in the Darcy model be made to agree well with the channelaveraged quantity λ 1 from lubrication theory. Although the value of ω can be used to alter the overall magnitude of both G 1 and 1 , the Darcy model substantially overestimates how strongly both of these quantities vary with h . Fig. 10 helps to explain why the Darcy model struggles to capture the behaviour of the perturbations. The leading-order velocity profiles for different values of h are nearly similar, so that when hu 0 is plotted against y / h they almost collapse onto the same curve ( Fig. 10 (a)). In contrast, the leading-order structure parameter does not collapse in this way: for large values of h, λ 0 remains close to 1, while for smaller values of h it varies significantly across the channel ( Fig. 10 (c) The gradients in λ 0 become important at O(δ) because they contribute to the Lagrangian derivative of λ. The fact that the gradients are highest near the centreline results in a perturbation λ 1 which is much more strongly localised for small h , as well as having a more pronounced peak ( Fig. 10 (d)). The localisation of the peak counteracts its increasing magnitude, so that λ 1 varies rather weakly with h , decreasing in magnitude from λ 1 ≈−0 . 452 for h = 0 . 25 to λ 1 ≈−0 . 047 for h = 4 ; in contrast, the Darcy model seems to track the peak magnitude of λ 1 rather than the average value, so 1 decreases in magnitude from 1 ≈−1 . 41 for h = 0 . 25 to 1 ≈−0 . 0056 for h = 4 . The consequence of this playoff between the size and the localisation of the perturbation to λ is that the velocity perturbation remains roughly in proportion to the leading-order velocity ( Fig. 10 (b)), rather than varying much more strongly with h as the Darcy model would suggest ( Fig. 9 (c)). The implication of these results seems to be that, because the Darcy model fails to resolve the variation of λ across the channel, it is not capable of representing accurately the advection of structure, and in particular of capturing the tendency of changes to the structure parameter to localise within the channel. Consequently, it may be seriously inaccurate as a predictor of the behaviour of thixotropic channel flows, and we may expect the inaccuracy found in this simplest problem to persist in more complicated scenarios. This in turn suggests that predictions made on the basis of the Darcy model, such as those of Pritchard and Pearson [24] , should be validated against a model that resolves the transverse variation. For some specific rheological models and flow problems, an approach of the kind employed by Livescu et al. [26] , which prescribes the transverse variation of λ, may be practicable, but it is not clear how to formulate such an approach in general.

Summary and conclusions
We have systematically developed a lubrication theory for the slow, steady, two-dimensional flow of a general thixotropic or antithixotropic fluid in a slowly varying channel. The role of microstructural change in such a flow depends on the magnitude of the advective Deborah number D relative to the aspect ratio δ of the channel; in different regimes of the advective Deborah number, different formulations of the theory are required.
In the present study we have concentrated on the weakly advective regime D = O(δ) , in which the fluid behaves as a generalised Newtonian fluid at leading order and thixotropic effects enter as a perturbation at O(δ) . We have presented illustrative results for the Moore-Mewis-Wagner (MMW) model of purely viscous thixotropic and antithixotropic behaviour, which contains many previous models as special cases, and for a regularised version of the Houška model of thixotropic yield-stress behaviour. In a particular case of the MMW model, explicit solutions are available up to O(δ) ; in other cases we have employed a semi-analytical approach in which the solutions are expressed in terms of inverses of known functions and integrals of these. The qualitative pattern of our results is rather similar for all the fluids considered: in an expanding channel, the effect of thixotropy is to decrease the fluid velocity and shear rate near the channel walls and to increase the velocity near the centre of the channel. The net effect of this is to decrease the pressure gradient required to drive a given flux through this section of the channel. Another conspicuous feature of the solutions is that the perturbations to the structure parameter tend to be strongly localised, either near the channel centreline or near the edges of a "pseudo-plug" region. As well as the transverse variation of the variables, we have considered the net pressure drop required to drive a given flow through a channel of known shape, demonstrating that in order to obtain the thixotropic correction to this quantity, unlike the leading-order term, it is sufficient to know the width only at the upstream and downstream ends of the channel.
The lubrication solutions that we have presented can be used as a means to validate the predictions of more strongly reduced models, such as the Darcy model employed by Pritchard and Pearson [24] . If a reduced model cannot capture the first-order corrections due to thixotropy in the weakly advective regime then its ability to represent thixotropy more generally becomes questionable. We have demonstrated that the disagreement between a Darcy model and lubrication theory can be surprisingly large, because the Darcy model is unable to capture the large transverse gradients of the structure parameter and thus to represent the strongly localised changes to the structure parameter induced by thixotropy. It appears that a more sophisticated approach to model reduction may be required in order to develop tractable models of thixotropic flow in confined geometries.
Although we have treated lubrication flow in one simple flow configuration, there are many other configurations to which an equivalent approach could be applied. Classic lubrication problems involving flow driven by a moving boundary include slider and journal bearings [45,Chapter 5] and squeeze flow; the latter has been used as a non-rheometric device for gauging the rheology of slurries [46,47] . Peristaltic pumping [17] is a similar boundarydriven flow of some practical importance. Alternatively, lubrication flows may be driven by a time-dependent pressure gradient, as in pulsatile flow [16] or in start-up problems [15,18,20] . Thixotropic free-surface flow [23,26] would be an interesting and challenging extension to this theory, involving the coupling of the evolving microstructure and the evolving free surface. With an eye to the subtleties of yield-stress behaviour [9,10] , a systematic investigation of how thixotropy affects the lubrication flow of a non-regularised yield stress fluid would also be of interest. Finally, with a view to petrochemical applications [14] , it would be valuable to extend the theory to incorporate additional physical effects such as temperature variation, fluid compressibility, or the transport of suspended particles.
In conclusion, a systematically developed lubrication theory offers a promising approach to a wide range of thixotropic flow problems. The results presented here, although restricted to one class of lubrication flows, suggest that this approach has the potential to provide dynamical insight as well as to guide the construction of models for practical application. advective regime The boundary-value problem defined by (36) -(40) can be solved numerically to give the variation of u 0 and λ 0 with y once h ( x ) is specified. However, to compute the O(δ) solution, we require the derivatives of the leading-order variables with respect to x . We could estimate these derivatives by differencing in x , but such an approach is clumsy and computationally intensive. It is preferable to write the O(1) solution in a semi-analytical form which allows derivatives to be evaluated in the same way as the variables themselves. We give them here in some detail to allow easy replication.

(A.5)
To obtain G 0 ( x ), we first employ Weissenberg's "trick" and integrate the flux condition by parts [48] , applying the boundary conditions at y = −h/ 2 and y = 0 , to obtain where the wall shear rate q w is defined by q w (x ) = q x, − h 2 , so τ (q w (x )) = We now obtain the asymptotic behaviour of the functions A and B defined by (48) and (49)  We can now obtain expressions for the perturbation quantities. From (50) we have (B.11) The sign of 2 α − a determines which term dominates as y → 0. For the parameters used in Fig. 4 we have 2 α − a = −0 . 2 , and so the second term on the right hand side dominates (although only narrowly, making the dependence rather hard to resolve numerically). Both exponents are, however, positive, so in this case ∂ u 1 / ∂ y remains finite as y → 0. The velocity perturbation u 1 , which is obtained by integration, is also finite.
We can now assemble an expression for λ 1 . Substituting the expressions above into (46) yields (B.12) Since α > 0, the dominant term is the first one, and we conclude that λ 1 ∼ y 2 α−a ; for the thixotropic case plotted in Fig. 4 we have 2 α − a = −0 . 2 , and so there is an integrable singularity in λ 1 at the centreline.
Finally, we comment on the case when 2 α − a > 0 . In this case, the dominant term in ∂ u 1 / ∂ y is O(y ) , so the second term in λ 1 is O(y α ) . The dominant term in λ 1 is therefore either O(y α ) or O(y 2 α−a ) ; since both exponents are positive, we conclude that λ 1 → 0 at the centreline.