A comprehensive constitutive law for waxy crude oil: a thixotropic yield stress ﬂ uid † Society of

Guided by a series of discriminating rheometric tests, we develop a new constitutive model that can quantitatively predict the key rheological features of waxy crude oils. We ﬁ rst develop a series of model crude oils, which are characterized by a complex thixotropic and yielding behavior that strongly depends on the shear history of the sample. We then outline the development of an appropriate preparation protocol for carrying out rheological measurements, to ensure consistent and reproducible initial conditions. We use RheoPIV measurements of the local kinematics within the ﬂ uid under imposed deformations in order to validate the selection of a particular protocol. Velocimetric measurements are also used to document the presence of material instabilities within the model crude oil under conditions of imposed steady shearing. These instabilities are a result of the underlying non-monotonic steady ﬂ ow curve of the material. Three distinct deformation histories are then used to probe the material's constitutive response. These deformations are steady shear, transient response to startup of steady shear with di ﬀ erent aging times, and large amplitude oscillatory shear (LAOS). The material response to these three di ﬀ erent ﬂ ows is used to motivate the development of an appropriate constitutive model. This model (termed the IKH model) is based on a framework adopted from plasticity theory and implements an additive strain decomposition into characteristic reversible (elastic) and irreversible (plastic) contributions, coupled with the physical processes of isotropic and kinematic hardening. Comparisons of experimental to simulated response for all three ﬂ ows show good quantitative agreement, validating the chosen approach for developing constitutive models for this class of materials.


Introduction
The term "yield stress uid" is used to describe so rheologically-complex materials that behave like a solid at low levels of imposed stress, yet ow when subjected to stresses that exceed a critical value. 1,2 These types of materials are ubiquitous in nature and in industrial applications, so predicting their rheology with a constitutive model is an important challenge. The term yield stress uid is a rather broad classication for materials which exhibit this general behavior of transitioning from solid to liquid, and these materials can exhibit additional complexities in their rheological behavior. Møller et al. 3 proposed using thixotropy as a distinguishing behavior that separates "ideal" yield stress uids from "non-ideal" yield stress uids. Thixotropy describes the ability of a material to exhibit a continuous and reversible time-dependent change (usually a decrease) in viscosity at a particular constant imposed shear rate. [4][5][6] These time-dependent changes in viscosity are due to a gradual change in the microstructure of the material resulting from the application of shear; this mechanism is commonly referred to as shear rejuvenation. Thixotropic materials also exhibit aging, which describes the ability of the material to reacquire its initial structure under the absence of shear (hence the reversibility of the initial time-dependent change). The aging phenomenon is commonly the result of thermally activated Brownian motion causing a sufficient rearrangement of the material's microconstituents. When this rearrangement occurs, the attractive forces between these microconstituents can cause a reformation of the sample spanning network in the material structure. 5,6 The classication used by Møller et al. is a simplied view of yield stress uids. For example, Divoux et al. 7,8 showed that even Carbopol, widely considered to be an "ideal" yield stress uid, 3 can exhibit some non-idealities (e.g. shear banding) that are associated with thixotropic yield stress uids. Given the ubiquity of thixotropic behavior in yield stress uids, it is important to have an appropriate constitutive framework for modeling the rheology of these uids. The development of such a framework is a critical issue which remains to be addressed. 6 There is a signicant body of literature which has strived towards this goal. [9][10][11][12][13][14][15][16][17][18][19] Many of the models in the literature adopt a scalar structure parameter (frequently, but not always, denoted as l) to describe some appropriate measure of the internal material structure. This material structure can build up or break down over time due to the application of shear, and the evolution equation for the structure parameter therefore typically contains an aging term, and a shear rejuvenation term. The structure parameter is then related to macroscopic rheological parameters in the uid, such as the yield stress, shear viscosity, or elastic modulus.
Other models in the literature take a different approach by modeling the underlying microstructural physics in these materials. Rather than prescribing a constitutive behavior based on bulk rheological measurements of the material, the constitutive relations are formally derived by considering interactions between the material's microstructural elements. Some of the more prominent examples of such models include the So Glassy Rheology (SGR) model, 13,20 the Shear Transformation Zone (STZ) model 14,21 and the Kinetic Elastoplastic (KEP) model. 18,22 As a testament to the capability of these models, several important experimental observations are predicted including shear banding/nonlocal effects, 23,24 power law dependency of linear viscoelastic moduli on frequency, and stress overshoots. 13 Together with their generality and the added intuition developed from the microstructural approach, these models are quite appealing. However considerable mathematical manipulation is oen required in order to express them in terms of macroscopic variables such as stress or strain, and they have not oen been used to quantitatively describe experimental data sets for real elastoplastic materials. On the other hand, models such as those of Quemada, 15 Mujumdar et al. 17 as well as many others, are explicitly written in terms of stress and strain. This form facilitates the tting and use of such models in engineering applications.
The particular subclass of thixotropic yield stress uids studied here are model waxy crude oils. Wax is a commonly encountered component of crude oils that can be present as a solid precipitate at low temperatures. These solid precipitates have adverse effects on the rheology of the crude oil during production scenarios, 25,26 so engineers and scientists must be aware of the thixotropic nature of waxy crude oils. 4,27 Detailed studies of the rheology of these uids emphasized the importance of sample preparation protocol (or beneciation 28 ) in order to ensure repeatable rheological data. [29][30][31] However, one of the issues presented when measuring the rheology on these uids is the variability in sample composition across different crudes (asphaltene content, wax content 32 ). For these reasons, we utilize a model waxy crude oil in the present study, the composition of which can be tightly controlled by formulating the uid in the laboratory.
The thixotropic behavior of these uids also makes reproducible rheological data difficult to obtain. Comprehensive, quantitative comparisons with model predictions are therefore scarce. Many of the models employed in the literature are simplied inelastic models, e.g. the frequently used Houska model 9,33 does not account for material elasticity, and the Pedersen model 28,34 assumes a generalized Newtonian uid with temperature dependency. A more complete and experimentallysubstantiated constitutive law would be useful for developing ow assurance strategies. 35 The goal of ow assurance is to ensure continuous ow of production uids (crude oils) under adverse scenarios, such as the presence of precipitates or non-Newtonian behavior. 26 As an illustration, we note that the pressure drop DP required to restart a gelled pipeline typically scales with the material yield stress s y . In the simplest case, a force balance on a gelled plug of wax gives DPpR 2 $ s y pDL, where D is the diameter of the pipe and L is the length. Designing a pipeline network with a xed pumping capacity DP therefore requires a priori knowledge of how the yield stress s y varies in the gelled crude oil (especially when aging processes and shear history cause large variations of s y with time). This necessitates a constitutive law that can be both t to rheological data, and has predictive capabilities.
A more comprehensive approach to the rheology of elastoviscoplastic materials, both in terms of experiments and modeling, is therefore required. Coussot et al. 10 and Mujumdar et al. 17 propose several different rheological experiments, among them startup of steady shear and large amplitude oscillatory shear (LAOS) as methods to validate and discriminate between constitutive models. There have also been many recent experimental studies which have combined bulk rheological measurements with velocimetry. These are important experimental efforts, because thixotropic yield stress uids oen exhibit shear banding and other types of heterogeneities in simple viscometric ows. This behavior is commonplace and is observed in a number of specic materials, including suspensions of star polymers, 36 Laponite gels, 37,38 emulsions, 24 Carbopol micro gels 7 and other colloidal gels. 39 Shear banding or inhomogenous ow can signicantly affect interpretation of experimental data, so the use of velocimetric techniques has become an important part of the experimentalist's toolbox.
Hence, combined with a detailed preparation protocol for our model uid (which is validated using direct measurements of ow kinematics), we will use a suite of well-dened viscometric ows to evaluate a constitutive model. The model we develop will be veried across a wide range of data. Once the model parameters are xed, the constitutive equations are solved for an additional independent test protocol in order to explore the predictive capabilities of the model. This gives further experimental validation to future numerical studies which would utilize the model to predict pipeline restart. 40,41 We endeavor to expand the approaches for describing thixotropic materials by introducing ideas from plasticity that have not been widely used in prior studies of such phenomena. The term "yield stress uid" generally describes a yielding material with rate-dependent plastic ow above a critical stress. This behavior is generic, and is even exhibited by metals at elevated temperatures. 42,43 The strengthening/soening processes that characterize thixotropy are similar to hardening mechanisms that are used in plasticity, which also utilize an internal microstructural variable of the form of the scalar parameter l. 42 These similarities suggest that there are underlying links between the models employed to describe thixotropic yield stress uids, and the models used in plasticity. 44 Early work hinted at the utility of employing concepts from plasticity to model yield stress uids, 45 but only recently have these mechanisms been used to successfully predict the behavior of a (non-thixotropic) yield stress uid. 46 In the present contribution, we show for the rst time how a constitutive model can quantitatively predict the rheology of waxy crude oil (a specic thixotropic yield stress uid) across a comprehensive set of experimental data. The data includes steady and transient tests, as well as linear and nonlinear rheological measurements. This is accomplished by utilizing detailed experiments on a model crude oil to motivate the development of the constitutive equations. We take a modeling approach that allows us to rapidly express our elastoviscoplastic equations of state in terms of macroscopic variables (stress and strain), which are the variables that are of interest to the practicing engineer or scientist who deals with ow of a real crude oil. While this approach may not give as detailed an account of the microstructural yielding and thixotropic mechanisms as other approaches, 13,14,18 we will still take care to provide a general discussion of the microstructural interpretation of these plasticity mechanisms. We will also illustrate the links between this model, and other models that have been employed to predict thixotropic behavior. The constitutive framework can describe the material response to a number of different rheometric test protocols, and also has a general 3D continuum mechanics formulation. It will therefore ultimately be useful for numerical simulations of complex ows of waxy crude oils that can then be used to guide ow assurance strategies.

Model uids & rheometry
The class of thixotropic yield stress uids used here has been characterized in detail in a previous study. 47 They are created by combining a heavy mineral oil (Sigma Aldrich 330760) with a paraffin wax (Sigma Aldrich 327212). The components are combined together at a high temperature ($60 C) and stirred continuously to fully dissolve the wax in the oil. The model uid behaves as a Newtonian liquid in this state, and can be loaded into the rheometer and cooled in situ.
Two different rheometers are used to probe the response of this uid to imposed deformations. A controlled strain rheometer (ARES LS, TA Inst.) is used for large amplitude oscillatory shear (LAOStrain) measurements. The ARES rheometer is equipped with a lower Peltier plate for temperature control of the sample. The upper xture is a cone with a diameter of 50 mm and an angle of 2 . The upper and lower surfaces are also roughened (root mean squared roughness R q x 30 mm) in order to suppress wall slip.
A stress controlled rheometer (AR-G2, TA Inst.) is also used. This rheometer is utilized in two different congurations. The rst of these is the standard conguration, which incorporates a lower Peltier plate for temperature control, and an upper cone (diameter 60 mm and angle of 2 ). This conguration also utilizes roughened upper and lower surfaces, with R q x 30 mm. The second conguration is a custom-designed Rheo-PIV conguration, which allows in situ measurements of the ow eld within the uid to be extracted. This apparatus is described in detail in Section 2.2.
In Fig. 1, we illustrate one of the characteristic thermorheological features of these model waxy crude oilsthe wax appearance temperature, T wa . Fig. 1 presents measurements of the viscosity of the heavy oil, and a 5% and 10% wax in oil mixture over a range of temperatures (obtained using the standard AR-G2 conguration). In Fig. 1, a shear rate of _ g ¼ 50 s À1 is imposed on the uid, while lowering the temperature at a rate of 1 C min À1 from 55 C to 25 C. As seen in Fig. 1, the mineral oil exhibits an Arrhenius-like exponential dependence of viscosity on temperature. Both the 10% and 5% wax in oil systems exhibit a similar Arrhenius-like variation of viscosity on temperature for large values of T. However, they rapidly diverge from this dependency at a temperature known as the wax appearance temperature, T wa . This temperature corresponds to the point at which wax crystals rst begin to precipitate in a sufficient amount to affect the viscosity of the uid. The morphology of these crystals is typically sheet-like or needle like, 30,48 and the high aspect ratio of the crystallites results in mechanical percolation and formation of a gel at low concentrations of precipitated wax. 49

Rheo-PIV conguration
A Rheo-PIV apparatus was used to obtain in situ measurements of the velocity/displacement eld within the model uids as they are sheared in the rheometer. This Rheo-PIV apparatus is based on a previous design, 47,50 and a schematic diagram is shown in Fig. 2(a). A laser light sheet (10 mW, 635 nm wavelength), created by a plano-concave lens, is directed vertically downward into the sample at a location r l . Displacements of illuminated tracer particles are then tracked to extract the deformation eld of suspended particles across the gap between the upper plate and lower cone (with angle f). The uid is seeded with reective TiO 2 particles of average size 3 mm, at a volume fraction of 2 Â 10 À5 (this volume fraction is low enough not to affect the rheology of the uid being studied). Cross correlation PIV soware (Digiow, Dalziel Research Partners, Cambridge UK) was used to process sequences of images and obtain the 2-dimensional velocity eld within the uid from a sequence of particle displacements. The 2-dimensional eld is averaged along the direction of ow, providing a velocity prole of the form v(y). In this case, y is the position between the lower cone and upper plate, with y ¼ 0 at the lower cone, and y ¼ H h r l tan f at the upper plate. This conguration differs from the standard conguration of the AR-G2 by using a lower cone and an upper plate. A stepped lower Peltier plate (TA Instruments part #531052.901) is utilized that allows for interchangeable lower geometries to be used. All measurements made using the Rheo-PIV conguration use either a 4 or 2 lower machined aluminum cone. The cone is black anodized in order to suppress reection of the laser sheet from this surface. The upper and lower geometries are smoother than the geometries used in the standard conguration (R q $ 0.1 mm)a polished upper geometry is necessary to provide a clear optical path for the reected light from the tracer particles.
The RheoPIV xture shown in Fig. 2 has one important alteration from previous versionsthe upper rotating plate, which is made from a transparent acrylic disc of diameter 50 mm, has a polished 45 bevel machined into the upper corner. This allows the laser-illuminated plane in the sample to be imaged through the at beveled surface, rather than through the curved uid meniscus at the edge of the geometry. The uid meniscus may change in shape over the course of an experiment due to sample shrinkage (which waxy crude oils are known to exhibit 51 ), or edge instabilities. 52 The at beveled edge avoids local image distortion which arises from refraction of light across the irregularly shaped air-uid meniscus. The index of refraction of the acrylic and model oils (n ¼ 1.49 and n ¼ 1.47 respectively) is very close. This refractive index matching, combined with the bevel, avoids optical aberrations and systematic errors that could arise from the lens/camera axis not being positioned orthogonally to an interface between two media with different values of refractive index n 1 s n 2 .
The apparatus shown in Fig. 2(a) allows for the location of the laser light sheet r l and the angle of the reective mirror q m to be adjusted and optimized. For the experiments presented in this work, the value of r l was kept constant at 20 mm, while q m was held constant at 45 (creating a vertical laser light sheet). It is possible to decrease q m so that the plane of illuminated particles is orthogonal to the camera's imaging axis. However, this was not necessary because the camera/lens assembly had a high enough focal depth (>2 mm), and all of the illuminated particles across the gap were sufficiently in focus. The laser light sheet is also very thin (approximately 50 mm), resulting in radial variations in particle location having a negligible effect on the velocity prole. A velocity calibration prole of a heavy mineral oil undergoing steady shear of _ g ¼ 0.5 s À1 in a 50 mm diameter, f ¼ 4 cone-plate geometry is shown in Fig. 2(b). The prole is averaged from 100 frames of video taken over a course of 3 seconds. It shows excellent agreement with the linear, theoretical prole predicted for a Newtonian uid with no wall slip at either xture surface.
2.2.1 Spatially averaged vs. local stresses and shear rates. Throughout this work, we carefully distinguish between spatially averaged measurements of stress and shear rate, and the local values of the stress and shear rate within the uid. The model waxy crude oils here are known to exhibit wall slip and other shear localization phenomena. 47 When these phenomena are present in cone-plate geometries, the stress and shear rate may no longer be uniform across the radial extent of the geometry. 53, 54 We will therefore dene a spatially averaged stresŝ s and shear rateĉ g as follows:  Protocol used to prepare the model waxy crude oil to its "slurry" state (black curve) and the "strong gel" state (blue curve). The rate of cooling is the same in both cases, however the slurry is formed under the application of a high shear rate (50 s À1 ).
where T is the torque that the rheometer imposes on the sample, U is the angular rotation rate of the geometry, R is the radial size of the geometry, and f is the cone angle. We also dene an average strain accumulated asĝh ð t2 t1 UðtÞ f dt. Under homogenous shearing with no wall slip, these average quantities reduce to the local true stress s, the local true shear rate _ g and the local true strain g. Due to the presence of wall slip and other shear localization events in the uids, we will reserve the symbols s and _ g to refer to the stress and shear rate under homogenous deformations only.

Preparation methods of the model uidslurry vs. strong gel
In order to carry out reproducible bulk rheological measurements on waxy crude oil, an appropriate sample preparation protocol must be selected. Waxy crude oils are extremely sensitive to variations in their thermal and shear history as they are cooled to below their wax appearance temperature. 55 Here we justify the selection of one particular preparation protocol which minimizes the impacts of wall slip and material instabilities. This provides reproducible rheometric data, which can subsequently be used for the development of a constitutive model.
We envision two distinct ow assurance scenarios under which waxy crude oils are transported through pipelines. For the rst scenario, consider a length of pipeline containing waxy crude at a temperature T > T wa , under quiescent (zero ow-rate) conditions. If the temperature of this uid slowly and uniformly drops to below T wa , then a gel network forms in the material due to the presence of solid wax precipitates. The cooling of this uid occurs under no ow conditions with zero shear rate. The crystallite aggregates therefore grow to large sizes and the percolated gel structure formed in the crude oil will be stiff (high viscosity and modulus). Restarting the ow in the pipeline aer cooling will therefore require large applied pressure drops, with the possibility that the uid undergoes an adhesive failure (or wall slip) at the pipe wall. 51 The second scenario involves cooling the waxy crude oil under non-quiescent (owing) conditions. For these conditions, the gel structure formed in the oil will be soer (lower modulus and viscosity), due to the disruptive effect of shearing the aggregates of wax crystals as they form. Smaller pressure drops will be required in order to maintain the ow of the uid through the pipeline, and wall slip will be less likely.
We refer to the state of the crude oil in the rst scenario as a "strong gel" and the state of the material in the second scenario as a "slurry". The strong gel and slurry states of the oil are reproduced in a laboratory scale setup by cooling the model wax-oil mixture from a single phase mixture at T > T wa in the rheometer either quiescently (under zero imposed shear rate) or under a high shear rate ðĉ g ¼ 50 s À1 Þ. These two sample preparation protocols are illustrated in Fig. 3. In previous work 47 we focused on the behavior of the wax-oil system in its strong gel state, which is dominated by wall slip and other interfacial failure processes. Here, we restrict our focus to the rheology of the material in its slurry state. In the slurry state, bulk rheological behavior is easier to interpret, due to the absence of these failure processes. Fig. 4 illustrates the critical differences between the response of a 5% wax-oil mixture in the two different states to a step in the shear rate toĉ g ¼ 1:2 s À1 . The response of the model uid in the slurry state is shown in Fig. 4(a). The model uid is brought to this state by cooling the sample from 55 C to 27 C at a cooling rate of 1 C min À1 , under an applied shear rate of c g ¼ 50 s À1 . Once the uid reaches the nal temperature of 27 C, the shear rate ofĉ g ¼ 50 s À1 is held for an additional 3 minutes, so that any thermal transients in the system die out. Aer these 3 minutes, the shear rate is immediately reduced tô c g ¼ 1:2 s À1 . This sequence of steps is administered using the AR-G2 in its RheoPIV conguration with a 2 lower cone (discussed in Section 2.2). Both the spatially averaged values ofŝ andĉ g , and the local ow eld within the sample are measured.
The response of the model uid in its strong gel state is shown in part (b) of Fig. 4. The uid reaches this state by undergoing the same protocol for the slurry state, however for the cooling step and the subsequent 3 minute holding time, the shear rate is held atĉ g ¼ 0 s À1 . This results in an essentially unperturbed gel network before the startup of steady shear of The RheoPIV conguration enables measurements of the average stressŝ in the uid, as well as measurements of the local velocity eld during the 300 second application of the shear rateĉ g ¼ 1:2 s À1 . In the slurry state, the spatially average stressŝ shows a constant measured value of 0.4 Pa. The strong gel state shows an initial peak in the stress, followed by a slow decay in the measured value ofŝ, occurring over the course of several minutes. The velocity eld within the uid in the two states (illustrated through the use of spatiotemporal diagrams for the rst 10 seconds, as well as the average velocity proles shown in (ii)) is markedly different for the two cases. For the strong gel state, larger stresses are needed to break the gel structure in the bulk of the uid. The stiff material therefore ows through the mechanism of wall slip at the lower surface, and remains adhered to the upper surface. The velocity prole in this state is plug-like, in contrast to what is seen for the material in the slurry state. In the slurry state, a velocity prole closer to linear is observed, with a non-zero shear rate in the bulk.
We quantify the long term evolution of the velocity proles in these two states by utilizing a non-dimensional metric F, dened in a previous study. 47 This dimensionless scalar measures deviations from homogenous linear shear and varies from 0, for a perfectly linear velocity prole with vðyÞ ¼ĉ g y, to 1, for a plug-like prole with velocity v(y) ¼ V p . F is evaluated using the following equation: A single value of F can be determined for each frame of video, so measured values of F corresponding to each frame of video are plotted for the rst 10 seconds of the experiment. Subsequently, values of F which have been averaged over 10 seconds of video are plotted at 30 second intervals (and the time axis on the gure is expanded in this region). The material in the slurry state shows constant values of F # 0.25 over the entire 300 second duration of the test. Measurements of F ¼ 0 are difficult to obtain for wax-oil uids below T wa , due to some residual wall slip, and due to the turbidity of the uid. This Illustrating the difference between the wax-oil system in its "slurry state" (a) and in its "strong gel" state (b). For each of the states, a step in (spatially averaged) shear rate toĉ g ¼ 1:2 s À1 is imposed at t ¼ 0, andŝ is measured over the course of 300 seconds (iv). In (i), we show spatiotemporal diagrams of the velocity field within the fluid for the first 10 seconds. Plots of the parameter F for the entire 5 minutes are given in (iii), while in (ii) we give the average velocity field within the fluid over the 300 seconds. Error bars are equivalent to 2 standard deviations of the velocity measurement at each height. turbidity introduces errors in the velocity measurements that always positively contribute towards F. However this behavior can still be compared to the behavior of the material in the strong gel state. In the strong gel state, F x 0.7 initially, then subsequently decreases. This process has been observed and documented previously in these uids 47 and is the result of erosion 37 whereby the wax crystal structure breaks down into smaller aggregates and the gel weakens.
The measurements ofŝ and F in Fig. 4 show that the strong gel response to an imposed shear rate is dominated by localized failure events at the wall. This is followed by long transients in the evolution of the heterogenous microstructure. These circumstances are different from the uniform distribution of stress and shear rate typically assumed for a uid being deformed in a cone and plate rheometer. In order to avoid the impact that these complicated yielding scenarios have on interpreting bulk rheological data, we choose to utilize the slurry preparation protocol for all experiments that follow. As is seen in Fig. 4(a), the 5% wax-oil system in its slurry state exhibits signicantly less wall slip. It also does not exhibit a prolonged decay in the stress signal due to uctuating stick-slip events and erosion of the microstructure. Some deviations of the velocity prole from the linear form for values of y/H < 0.4 are observedthese are due to the turbidity of the material decreasing the precision of the velocity measurements deeper into the uid (and consequently the error bars are much larger for these values of y/H < 0.4). Bulk rheological data for material prepared in this slurry state is more congruent with a continuum mechanics constitutive modeling framework (especially when roughened surfaces are used to eliminate residual wall slip). Fig. 5 graphically illustrates the difference between the "strong gel" and "slurry" states. The gure contains two bright-eld microscopy images of the wax crystal microstructure for a 5% wax-oil mixture at 25 C in its slurry state (a), and in the strong gel state (b). The wax crystals appear as dark particles in a light background. These crystals are discotic, but appear needlelike when viewed edge-on in the two-dimensional plane. The wax network formed in the slurry state is looser with more vacant (light gray) regions where no precipitates are present. The strong gel state consists of fewer vacant regions with a more closely packed network of wax precipitates. The wax crystals in the strong gel are non-Brownian in character, and form a sample-spanning jammed structure. However in the slurry the discrete crystals (with typical sizes 10-50 mm) can be aligned by the shearing ow. This results in strong shear thinning in the apparent viscosity observed in the steady-state ow curve (see Fig. 6).

Rheology of wax-oil system in slurry state
We propose three different rheometric tests as canonical ows to probe the material response of the wax-oil mixture, and provide a set of data for tting to a constitutive model. In what follows, we describe these three different test protocols, and discuss implications of the experimental results for constitutive modeling.

Steady state ow curve
For thixotropic systems, steady state ow curves (that is, plots of viscosity vs. shear rate, or shear stress vs. shear rate) are difficult to obtain. 5 The measured viscosity will depend on the duration of the experiment, due to the thixotropic transients exhibited by the material. However, ow curves are commonly used in tting of constitutive models to thixotropic yielding materials. 31,56 de Souza Mendes and Thompson 57 argue that the steady state ow curve describes the equilibrium locus of the dynamic thixotropic system. For constitutive models which incorporate an evolving internal structural parameter (commonly denoted generically as l(t) 5 ), this implies that the steady state ow curve can be t to the constitutive model predictions for the special case when _ l ¼ 0. This signicantly simplies the tting procedure.
In Fig. 6 we plot the measured ow curves for the heavy mineral oil, and the 5% and 10% wax-oil systems in their slurry state (T ¼ 27 C). These measured ow curves are obtained by using the AR-G2 in its standard conguration. The measurements are carried out by shearing at a given globally averaged rateĉ g, and are determined to have reached a steady state by  using a "steady state sensing" option on the rheometer. Using this option we measure the temporally averaged torque T over successive 30 second periods, and determine that the measurement has reached a steady state when three successive sampling periods are within 5% of the same value. Each measurement point usually requires 5 minutes to reach a steady state, and the measurement of the entire ow curve requires 2-3 hours.
Several key features can be identied from the ow curves shown in Fig. 6. First, both the 5% and 10% model uids exhibit Newtonian-like behavior at high shear rates, with the average stressŝ being linearly proportional to average shear rateĉ g . Second, the 10% system has a larger steady state viscosity than the 5% system, due to the higher volume fraction of solid precipitates. Third, the ow curves of the model crude oils exhibit a non-monotonicity, i.e. there is a region at lower applied shear rates ðĉ g # 0:1 s À1 Þ where the average stress decreases as the shear rate is increased.
Non-monotonic ow curves (NMFC's) point towards the presence of a material instability. 58-61 Some shear banding uids, for example wormlike micellar solutions, exhibit an underlying non-monotonic ow curve. 62 This non-monotonic ow curve cannot be measured, because at imposed shear rates in the decreasing region of the ow curve, the material bifurcates into two 63 or more 64 coexisting regions with different local shear rates. We refer to this as a "standard shear banding scenario". 61 There are several examples in the literature of measured NMFC's. Dijksman et al. 65 measured NMFC's in granular media in a number of different geometries, while Michel et al. 66 measured non-monotonic ow curves in microemulsions of oil in water which form transient networks. In the latter, the authors observe "wavy deformations of the sample surface" and "progressive cavitation-like proliferation of bubbles in the bulk". These manifestations are indicative of an unstable ow within the material, even though the imposed deformation is steady. This differs from the standard shear banding scenario envisioned for wormlike micellar solutions, in which a steady (banded) velocity eld is established in the material. The underlying non-monotonic constitutive curve is then manifested as a stress plateau with an approximately constant value of the average stressŝ. 61 For our model crude oil, the data shown in Fig. 6 indicates that a NMFC ofŝ vs.ĉ g is measurable. We therefore emphasize that the results shown in Fig. 6 correspond to temporally converged values of average stressŝ in the material at an imposed global shear rateĉ g. However, as pointed out by Michel et al., 66 substantial temporal and spatial uctuations in microstructure and local stress must persist.
3.2.1 Rheo-PIV in the non-monotonic region of the ow curve. The AR-G2 in its Rheo-PIV conguration was used to provide experimental evidence for a material instability associated with the NMFC. The 5% wax-oil system in its slurry state at a temperature of T ¼ 27 C was utilized for this (higher wax content systems are too turbid to provide discriminating velocimetric data). The particular protocol that is used is as follows: rst; preshearing the material atĉ g ¼ 50 s À1 for 5 minutes. Then, a sequence of three steps in shear rate is imposed. The rst step is toĉ g ¼ 0:02 s À1 (at a time denoted as t ¼ 0 s), followed by a step toĉ g ¼ 1:2 s À1 at t ¼ 1432 s, and nally a step back toĉ g ¼ 0:02 s À1 at t ¼ 3232 s. The measured stress and velocity proles for this test are shown in Fig. 7.
From Fig. 7(a) it is apparent that the rheometer imposes an apparent shear rate in the two cases which differs by a factor of 60, however the steady state value of the average stressŝ is roughly equal (within AE12%) for both values ofĉ g . This is indicative of the extreme shear thinning exhibited by the model waxy crude oil, and is consistent with the non-monotonicity in the measured ow curve of Fig. 6. Transient responses are also exhibited in the response of the material to the steps in shear rate. At t ¼ 1432 s, when the shear rate is stepped tô c g ¼ 1:2 s À1 , an instantaneous increase in the stress occurs, followed by a subsequent decrease with the stress reaching a minimum at approximately t ¼ 1452 s. Then, there is a gradual increase (from t ¼ 1452-2800 s) in the stress towards a steady value ofŝ ¼ 0.71 AE 0.01 Pa. At t ¼ 3232 s, when the shear rate is suddenly dropped back toĉ g ¼ 0:02 s À1 , the stress instantaneously decreases, followed by a subsequent increase in the stress reaching a maximum (and steady value) at t ¼ 3280 s.
These long transient responses are consistent with a thixotropic material undergoing cycles of structuring and destructuring. 67 In particular, from t ¼ 1432 s to t ¼ 1452 s, the material evolves from a structured state (due to the low value of the shear rate imposed for t # 1432 s) to a destructured state. The high shear rate aer the jump results in a sudden increase in the stress, which subsequently breaks the wax structure. The structure breaks down completely and the stress reaches a minimum at t ¼ 1452 s. From t ¼ 3232 s to t ¼ 3280 s when the shear rate is decreased, the material starts from a destructured state due to the high value of the shear rate imposed previously. A smaller stress is required to drive the ow, but this stress progressively increases due to a build up in the structure at the lower imposed shear rates. From t ¼ 1452 s to t ¼ 2800 s we also observe the additional effect of a long-time structuring which results in an increase in the measured or apparent viscosity within the uidĥðtÞhŝðtÞ=ĉ g . This presence of multiple timescales for restructuring of the uid is to be expected. Fig. 5(a) indicates that for the material in the slurry state, there is a distribution of characteristic sizes of the wax precipitates and aggregates. The characteristic rotary diffusion time for rearrangement of the discotic microstructure is a strong function of the size of these crystals and aggregates. 5 The evolution of the bulk uid rheology should therefore reect this distribution of timescales.
The measurements ofŝ in Fig. 7(a) provide clear evidence that the model waxy crude oil is highly thixotropic. The larger temporal uctuations in the average stressŝ(t) at lower shear rates is also notable. However, Fig. 7(a) provides no indication whether the ow eld within the material is heterogenous or shear-banded. Rheo-PIV evidence supporting a material instability is given in Fig. 7(b) and (c). In these two sub-gures we contrast the measured velocity eld within the uid as two different values ofĉ g are imposed. At the lower imposed shear rate ofĉ g ¼ 0:02 s À1 , we observe stochastic uctuations in the average shear rate across the gap, in the slip velocity at the upper plate, and in the general form of the velocity prole. From t ¼ 0-100 s, the velocity prole indicates a region of higher shear rate near the center of the gap. This is reminiscent of the 3-banded velocity proles which have been observed in wormlike micellar solutions under imposed steady shear in cone-plate geometries. 50,64,68 However, this structure in the velocity prole is ephemeral. The velocity proles at later times (e.g. t ¼ 1150 s and t ¼ 1300 s plotted in Fig. 7(b)(ii)) show a more uniform shear rate across the gap with larger slip velocities at the upper and lower plates. The velocity prole plotted for t ¼ 1300 s is nearly plug-like (F ¼ 0.73) with scaled slip velocities of (1 À v(y)/V w ) ¼ 0.48 at the upper surface and (v(y)/V w ) ¼ 0.28 at the lower surface. We quantify these uctuations by considering a Reynolds decomposition of the velocity eld v(y, t) such that: where v(y) is the average part of the velocity dened as follows: and v 0 (y, t) is the uctuating part of the velocity. We dene a temporally and spatially averaged dimensionless shear rate measured across the gap as follows: where v(H) and v(0) are the temporally-averaged velocities of the upper and lower walls. To quantify the effect that the uctuating part of the velocity v 0 ( y) has on the shear rate, we dene a standard deviation of the shear rate as follows: For the velocity data in Fig. 7(b), the average shear rate measured across the gap is h _ gi AE 2S _ g , where the AE represents the 95% condence bound in this average over time, and is indicative of the extent of the uctuations in the local shear rate _ g over time. For Fig. 7(b), we compute h _ gi ¼ 0.59 and S _ g ¼ 0.13. In Fig. 7(c), we show measured proles for the case when c g ¼ 1:2 s À1 . At this higher shear rate, images are acquired at a higher frame rate due to the increased velocity of the tracer particles. Due to limited digital storage space, sampling of the velocity eld is limited to several smaller intervals in time. The velocity eld is therefore sampled for 5 different 10 second intervals from t ¼ 1670 s to t ¼ 2040 s, which coincides with the region of increasing measured stress. The velocity is also sampled for 30 seconds at t ¼ 3030 s, where the spatially averaged stressŝ has approached and remains at a constant value. In these regions, the average shear rate across the gap uctuates considerably less, with h _ gi ¼ 0.66 and S _ g ¼ 0.068. The 95% condence bounds for the velocity proles have been reduced by a factor of 2 compared to the lower shear rate. The velocity eld also appears linear for all times, with some wall slip still occurring at both the upper ((1 À v(H)/V w ) ¼ 0.14) and lower (( v(0)/V w ) ¼ 0.21) surfaces. The uctuations in shear rate and slip velocity that are observed in Fig. 7(b) are therefore markedly reduced when compared to Fig. 7(b). Fig. 7(b) (ii) and (c)(ii) also include envelopes which represent the magnitude of the uctuations in the velocity v 0 ( y). These envelopes are dened as v( y) AE 2S v ( y), where S v ( y) is the standard deviation of the velocity as a function of position across the gap, and is dened as follows: We nd that S v ð yÞj^c g ¼1:2 \S v ð yÞj^c g ¼0:02 for all values of y, indicating that there are more uctuations in the velocity when the applied shear rate isĉ g ¼ 0:02 s À1 . The difference in the magnitude of these uctuations is highest near the upper surface (y ¼ H), where uctuations are 4.8 times larger for the case whenĉ g ¼ 0:02 s À1 . Incidentally, at y ¼ H the velocity measurements are most accurate due to diminished effects of sample turbidity. The RheoPIV data and associated statistical analysis presented in Fig. 7 provides strong evidence of spatiotemporal uctuations in the velocity eld associated with a material instability in the wax-oil mixture. This unstable ow is most pronounced at low shear rates (which lie in the decreasing region of the NMFC).

Stress overshootsstartup of steady shear.
The second rheometric test that we utilize to characterize thixotropy is the measurement of overshoots in the average stressŝ under startup of steady shear. These tests probe the effect of the "aging time" that elapses between sample preparation, and imposing the step in shear rate. The initial sample preparation is typically carried out by shearing (or rejuvenating) the material at a high and constant deformation rate (well within the stable region of the ow curves shown in Fig. 6) for an extended period of time. Then, during the period of aging, the sample is le unperturbed in the rheometer so that the material microstructure starts to build up. When the step in shear rate is subsequently imposed, the gelled and structured material will go through a yielding transition which depends on the initial structure, which in turn depends on the aging time. These tests have been used previously to quantify the effect that internal microstructure has on the rheology of a thixotropic system, and to quantify the characteristic time scales which may be associated with the formation of this structure. 69 The tests presented here are carried out on the 10% model waxy crude oil prepared to its slurry state at a temperature of 27 C. The AR-G2 in its standard conguration is used, which utilizes roughened cone-plate surfaces to eliminate wall slip. We adopt the nomenclature of Fielding et al. 70 and denote the aging/waiting time as t w . Aer the sample is prepared, an initial waiting time of t w ¼ 2 s is applied, and then a step in shear rate toĉ g ¼ 2 s À1 is imposed for 10 minutes. Following this step, the material is le unperturbed for a waiting time t w ¼ 5 s, and an applied steady shear rate ofĉ g ¼ 2 s À1 is imposed for 10 minutes. In total, seven startup tests of steady shear are imposed successively, with the waiting time between these steps being incremented in a logarithmic fashion up to a maximum value of t w ¼ 200 s. For each value of t w , measurements of the stressŝ are obtained for the 10 minute step in shear rate. For consistency, the measurements of stress and strain rate in this section will be indicated as averages (ŝ andĉ g ). However, the average strain rateĉ g here is high enough, and wall slip is suppressed by the roughened geometry, so deviations of the average values from the local shear rate _ g and local stress s will be small.
In Fig. 8 we plot the measured response of the average stresŝ s as a function of the imposed strainĝ, for seven different values of t w . An overshoot is seen at applied deformations of approximatelyĝ ¼ĉ g Dtx0:3 dimensionless strain units. The magnitude of this overshoot increases with t w , however it appears to saturate for t w $ 20 s. The overshoot is clearly the result of structuring within the uid that occurs over the course of the aging time t w . There is also a gradual monotonic increase in the long time steady state value ofŝ which occurs over much longer time scales, and over multiple steps in the shear rate. This slower increase inŝ is a manifestation of the same slow aging process that is observed in Fig. 7(a) even at a (stable) applied shear rate of 1.2 s À1 it takes a long time (on the order of tens of minutes) for the value ofŝ(t) (and henceĥ(t)) to stabilize. This slow structuring occurs on time constants that are much larger than the waiting time of t w ¼ 20 s which is approximately required for the peak stress to saturate. This is additional experimental evidence of a distribution of processes and time scales associated with the internal material structuring.
Our constitutive model will primarily focus on predicting the transient yield peak in the stressŝ(t) at intermediate strainsĝ x 0.3. This is frequently of greater relevance to ow assurance than the very long time response of the material stress. These initial yield peaks observed in Fig. 8 would be responsible for relatively large additional pressure drops being required to restart a pipeline where ow has ceased for a short period of time. Several other constitutive models have had success in predicting these types of overshoots (e.g. see the SGR model, 13 or the thixo-elastic Maxwell model, 15 among any others). In the present work, we develop a model capable of predicting these overshoots by incorporating mechanisms from plasticity theory. The model will quantitatively predict these stress overshoots, as well as several other key experimental observations. Although lacking the more detailed microstructural physics interpretation of some other models, 13,14,18 the model introduced here will be written in terms of a set of evolution equations in macroscopic variables, rather than in integral form. This differential form is more suitable for implementation in common Eulerian ow solvers. Furthermore, the mechanisms utilized in this model are grounded in a robust continuum mechanics framework, so the model can also be generalized to tensorial form in a straightforward manner (see the table in the ESI †).
To quantify the magnitude of the yield peak, which will in turn assist in the tting of the constitutive model, we dene our stress overshoot Ds 0 as the difference between the peak stress during each test and the measured stress 10 seconds (i.e. 20 strain units) aer the step in shear rate is imposed. At this level of imposed strain it is clear that the initial transients associated with the thixotropic yield peak have died out. The measured value of Ds 0 is plotted as a function of t w in Fig. 9. In this gure, the magnitude of the yield peak saturates with Ds 0 approaching a value close to 1 Pa, or approximately 40% of the nal ow stress. Some further increases at longer times (between 100-200 s) are also evident, and these smaller increases are evidence of the same long time restructuring behavior that was discussed in Section 3.2.1.
3.2.3 Large amplitude oscillatory shear. The third rheometric test utilized here to probe the time-varying nonlinear constitutive response of the material is large amplitude oscillatory shear (LAOS). LAOStrain, is generally carried out by imposing an oscillatory deformation on the material of the following form: For the measurements presented here an oscillatory average strain will be imposed on the material,ĝ ¼ g 0 sin ut. Even though we suppress wall slip, we will utilize the averaged quantitiesŝ,ĝ andĉ g to represent the measured rheological data. This is to account for the possibility of small uctuations in the local variables s, g and c g. Fig. 8 Stress overshoots for different waiting times t w in the 10% model waxy crude oil system. Fluid is prepared to its slurry state at 27 C. The applied shear rate isĉ g ¼ 2 s À1 . Two independent test parameters can be modied in LAOS; the oscillation frequency u and the strain amplitude g 0 . When the strain amplitude is small enough, the stress in the material s(t) is linear in strain and can be decomposed into in phase and out of phase components as follows: where the storage and loss moduli G 0 and G 00 are both functions of the frequency u. When strains are large, most materials exhibit a nonlinear response with a periodic stress that can be expressed in terms of a Fourier series as follows: where G 0 n (u, g 0 ) and G 00 n (u, g 0 ) are higher harmonic moduli that now depend on both the frequency of oscillation and the strain amplitude. This approach, using Fourier transform rheology, 71 is the basis for many techniques which analyze the response of a material to LAOS deformations. 72 We do not analyze our LAOS data using the Fourier transform representation, because this decomposition can only be applied to periodic stress waveforms. Because of the thixotropic nature of the model crude oil, the material response to a sinusoidal shearing deformation will contain prolonged stress transients due to the aging and shear rejuvenation of the uid that can occur over multiple periods of oscillation. Our goal is to quantify these transients and then develop a constitutive model that predicts them. Any measures based on Fourier transform rheology cannot be applied to these transients. Instead, an instantaneous approach such as that introduced by Rogers 73 would be required.
The LAOStrain measurements presented here will be used primarily as a tting tool in order to guide the development of a suitable elasto-viscoplastic model, and then to compare the model response to the response of the real uid. One benet of using LAOS tests to t constitutive models is that several key aspects of the rheological behavior of the uid are revealed under LAOS forcing. At low values of g 0 , the linear viscoelastic behavior of the uid is probed. Larger values of g 0 can then probe the plastic yielding transition of the test material. It is also possible to observe transients in the LAOS response which result from the thixotropic behavior of these uids. These transients are represented clearly in the Lissajous curves as decays towards a periodic attractor. Hence, if a constitutive model can predict the response of a material to LAOS, then it is likely to predict other time-and strain-varying responses that occur over a wide range of deformations.
Several researchers have used cyclic stress-strain loading curves, or Lissajous-Bowditch curves, as signatures that are used to t material responses to a constitutive model. 17,46,74 We represent the response of our model crude oil to LAOS through the use of these Lissajous-Bowditch curves (with strain as the abscissa and stress as the ordinate). A series of these curves can be plotted for a number of different strain amplitudes g 0 and frequencies u. We carry out a series of LAOS tests on the 10% model crude oil uid in its slurry state at a temperature of 27 C.
These measurements are done on the ARES rheometer, with a roughened cone-plate geometry conguration. A series of measurements at progressively larger strain amplitudes g 0 and at a single xed frequency of u ¼ 1 rad s À1 is shown in Fig. 10. Between each measurement, a waiting time of t w ¼ 100 s is applied so that the structure associated with the transient yield peaks shown in Fig. 8 has enough time to build up. For each measurement at a given strain amplitude, four complete oscillation periods are applied.
In Fig. 10(b) we represent the Lissajous curves as 3D trajectories with the measured stressŝ plotted against the sinusoidally varying strainĝ, and the strain rateĉ g which varies in quadrature. We are primarily interested in the stress-strain projection of these 3D curves, as shown in Fig. 10(a), because this reveals the elastoviscoplastic character of the material response. The curves in Fig. 10(a) and (b) are only shown for their stable limit cycle or "alternance state". 75 This alternance state is reached aer multiple cycles of oscillation. In Fig. 10(c), we illustrate the transients which occur before reaching the alternance state by plotting these Lissajous curves on independent (rescaled) axes for each strain amplitude g 0 . The magnitude of the maximum stress s m and strain at each strain amplitude is shown by each gure. In Fig. 10(c), overshoots and prolonged transients are clearly evident in the response of the material.
From Fig. 10(c) we can see that at low values of g 0 the material response is that of a linear viscoelastic solid, with G 0 > G 00 . This can be distinguished from the elliptical shape of the loading curve, which implies a linear response. As g 0 is progressively increased, the maximum stress s m also increases. However s m then saturates at moderate strain amplitudes, a behavior typically associated with a yield-like response. For moderate strain amplitudes between 5% and 50% the characteristic shape of the cyclic loading curves changes and indicates yielding. At these values of g 0 the material undergoes a sequence of processes which involve elastic loading, followed by saturation of the stress and subsequent viscoplastic ow, before elastically unloading when the direction of straining is reversed. This sequence of events occurs within each cycle (i.e. it occurs on an intra-cycle basis), and is responsible for the transition from elliptical to rhomboidal shaped loading curves. 76,77 There is additional valuable information contained in the material's inter-cycle response, i.e. in the prolonged transients exhibited, especially for intermediate strain amplitudes of 20% # g 0 # 100%. These transients are manifested in the form of an initial stress overshoot which occurs as the material is rst strained. This stress overshoot is evidence of the same transient yield peak that was documented in Section 3.2.2 for the startup of steady experiments. Several cycles are required (spanning a time period on the order of 10 seconds) for the stress to settle into its steady alternance state.
Several key features of Fig. 10 can be used to guide the development of our constitutive model (the detailed formulation of which will be given in Section 4). First, there is no observable inter-cycle change in the linear viscoelastic behavior between each cycle of loading. This implies that G 0 and G 00 remain constant with time for the material. If both G 0 and G 00 are to be non-constant functions of some microstructure parameter l, which in turn depends on t w , then deformations in the linear regime cannot cause a change in the microstructural parameter l. This will be accounted for by specifying a plastic strain as being responsible for changing the structure parameter l, and this plastic strain is only allowed to accumulate when the material yields at larger stresses and strains. Thus, a central component of our modeling approach will be to additively decompose the total strain into plastic and viscoelastic contributions, with plastic strain only accumulating when the shear stress in the sample exceeds some critical value. As a mechanical analog representation, this model can be visualized as a viscoelastic element in series with a nonlinear plastic yielding element. The yielding element will then utilize two important concepts from plasticity theoryisotropic hardening and kinematic hardeningin order to capture both the nonlinear intra-cycle behavior of the material, as well as the transient inter-cycle behavior (i.e. stress overshoots) which is exhibited under LAOS.

Formulation of the isotropic-kinematic hardening model (IKH model)
In this section we outline the development of a constitutive model that quantitatively predicts the salient features of the model crude oil's rheological response. As noted above, a central component of our approach is to additively decompose the total shear strain in our material into linear viscoelastic (g ve ) and plastic (g p ) components, such that Irreversible plastic strain accumulates in the material when the applied stress is above a certain critical value, so for low stresses g ¼ g ve and the response of the material is linear and viscoelastic. This decomposition differs slightly from classical descriptions of viscoplasticity, where strain is additively decomposed into plastic and purely elastic components. Here we rename our elastic strain a viscoelastic strain in order to account for some linear viscoelastic effects that are present in our material below the yielding point. This linear viscoelastic response can specied further by decomposing the viscoelastic strain g ve into elastic and viscous components, such that g ve ¼ g v + g e . This would correspond to specifying a Maxwell-like linear viscoelastic behavior below the yield stress, with g e ¼ s/G and _ g v ¼ s/h. Alternatively, a Kelvin-Voigt viscoelastic behavior can be specied. In fact any LVE-type behavior can be specied within this general constitutive framework, depending on the characteristics of the so gel and the desired delity of the model prediction. We show three canonical forms in Fig. 11 in the form of corresponding mechanical analog elements. Because the ows of interest are nonlinear, we will predominantly focus on prescribing the correct behavior of the nonlinear yielding element that is characterized by the plastic strain g p . As a result, we will not discuss tting of linear viscoelastic behavior in detail. Interested readers can consult any of the classical textbooks on linear viscoelasticity 78-80 and replace the simple linear viscoelastic elements shown in Fig. 11 with the elements of their choice (and this will not signicantly affect the large strain nonlinear behavior).
For most of the calculations below, the simplest canonical form of LVE behavior is considered, i.e. a linear Maxwell-like behavior as shown in Fig. 11(a) with elastic modulus G and damping h. Incidentally, this MIKH form of the model is compatible with the classical description of yielding materials with g ¼ g e + g p , since the plastic strain can be redened to include the linear dampling element with viscosity h. In addition to this, two assumptions are made about the nature of the linear viscoelastic behavior. First, we observed in Section 3.2.3 that under small amplitude oscillatory shear, G 0 and G 00 remain constant in timei.e. imposed deformations in the linear regime do not affect material structure. Therefore, only the plastic accumulated strain g p can drive a thixotropic restructuring of the uid. The second, stronger assumption is based on measurements of G 0 and G 00 as a function of t w . These measurements show that both G 0 and G 00 are in fact at most weak functions of waiting time t w , with G 0 only varying by 20% and G 00 varying by even less as t w is changed (see ESI †). For this data, the LVE parameters G and h are therefore assumed to be constant, and do not vary signicantly with changes in material structure. This assumption has been discussed by Quemada,15 and is utilized here to reduce the number of constitutive parameters associated with this model. However, the IKH model will still quantitatively capture the three nonlinear ows discussed in Section 3.
4.1.1 Specifying the plastic strain g p . Having discussed the linear viscoelastic behavior, the form of the plastic shear strain g p must now be specied. A common rst order assumption 4,6,81 is to prescribe a Bingham-like behavior, in which plastic ow only occurs above a critical stress s y : where n p ¼ s/|s| is the direction of the applied stressthis is the codirectionality hypothesis, 42 which ensures that plastic deformations always occur in the same direction as the applied stress. The parameter m p is the plastic viscosity.
The simple Bingham equation given in eqn (14) captures some of the basic features of the ow curve in Fig. 6, however it lacks the ability to predict the transient behavior and nonlinear stress overshoots that were documented in Section 3. To predict these transients, the ow rule in eqn (14) must be modied to account for the constitutive processes of kinematic hardening and isotropic hardening/soening. 42 We have used the concept of kinematic hardening in previous work to model the response of simple (non-thixotropic) yield stress uids, 46 Kinematic hardening accounts for a movement of the center of the yield surface in stress space. Isotropic hardening, on the other hand, accounts for an expansion or contraction of the yield surface, and is frequently used to model the strengthening of polycrystalline metals under imposed deformations. 43 These hardening mechanisms are abundantly used in the plasticity literature, and full three-dimensional versions of models with combined isotropic and kinematic hardening have been developed 82,83 that are both frame invariant and thermodynamically consistent. For simplicity, we only describe here a one-dimensional version of the model. To capture these two additional microstructural mechanisms, the Bingham equation given in eqn (14) is modied as follows: where the total stress s has been replaced by the effective stress, s eff , which is dened as follows: here s back is the back stress in the material, which effectively tracks the location of the center of the yield surface. The yield surface refers to a surface dened in the three-dimensional "stress space" formed by the principal stress axes. Many different denitions of yield surfaces exist, e.g. von Mises, Tresca. However, we deal here with a one-dimensional form of the IKH model, so our stress space consists of a single axis (s), and our yield surface is thus dened by three points. The rst is the back stress, s back , which is the center of the yield surface, and the other two are the yield limits in the positive and negative directions of stress, s back + s y (l) and s back À s y (l), respectively. In this model, plastic ow only occurs when the applied stress s lies outside of the yield surface dened by these three points. The presence of both the back stress and yield stress in the material can be more thoroughly illustrated in Fig. 11 by decomposing the stress s acting on the nonlinear yielding element into individual stress branches. In addition to the presence of the back stress in eqn (15), the yield stress s y is now a function of a single internal structure parameter l. Furthermore, the direction of plastic strain n p is now codirectional with the effective stress, so that: With the formulation in eqn (15)- (17), the center of the yield surface can change through the variation of the back stress s back with the imposed plastic strain (kinematic hardening) and the yield surface itself can contract/expand through the variation of s y with microstructure (isotropic hardening). The incorporation of two yield parameters (center of yield surface s back and size of yield surface s y ) bears some similarity to the Houska model used by Sestak et al. 84 and Cawkwell and Charles, 9 in which the yield stress is decomposed into two components (a constant term and a varying component). This model is commonly used to predict the rheology of crude oils, however as an inelastic model, it cannot capture the full behavior of waxy crude oils, especially in large amplitude oscillatory ows.
The variation of the back stress s back and the yield stress s y for a given deformation must be specied. The back stress is related to an internal structural variable A which was utilized in previous work to model ideal (non-thixotropic) yield stress uids: 46 The linear relationship between s back and A results from taking the derivative with respect to A of a quadratic form of a defect energy in the material that accumulates with the plastic strain g p . 46,83,85 This is analogous to an additional elastic free energy for the defects, so C can be thought of as the back stress modulus, and has units of stress. 85 The internal variable A can be physically interpreted as a dimensionless strain-like variable, or a back strain. The defect energy accumulates because of a rearrangement of the material microstructure, which in our case consists of the interlocking wax crystal network. For threedimensional versions of constitutive models that include kinematic hardening, the variable A necessarily takes on a tensorial form. 83 The back strain A follows an evolution equation that is a more general form of the Armstrong-Frederick equation 86 discussed by Jiang and Kurath 87 that can be used to capture the Bauschinger effect, 88 The Bauschinger effect refers to a decrease in the yield strength of a material upon reversal of the direction of stress. 42 In the above equation, when f (A) ¼ qA (with q a dimensionless material constant) the classical Armstrong-Frederick equation is recovered, which is typically used to describe cyclic loading in metals, but was more recently used to describe the LAOS behavior of so Carbopol microgels. 46 To more accurately capture the intra-cycle behavior under LAOS model waxy crude oil over a wide range of imposed strain amplitudes g 0 , a power law function of f (A) is specied: where both m and q are material constants. In the tting which will follow in Section 4.2, eqn (19) is always evaluated with an initial condition of A ¼ 0. This is the approach taken in previous workit presupposes that the material is initially in a "virgin state 82 with no directional preference for yielding (i.e. a back stress of zero). It also follows from the assumption that the deformation is beginning from a state where the stress and shear rate in the material is zero. Under such conditions, the back stress and back strain typically relax back to their zero values. When the plastic strain rate _ g p is nonzero, the back strain A will initially increase in the direction of the plastic strain g p . This results in an increase in the back stress, or a "hardening" of the material in the direction of ow. However, the back strain A eventually saturates due to the second term in eqn (19). The saturation value of A is AE1/q, with the back stress therefore saturating at AEC/q (with the sign depending on the direction of the plastic straining _ g p ). The value C/q therefore corresponds to the center of the yield surface under steady owing conditions ( _ g p s 0) and for the case of a material with s y ¼ 0, and purely elastic behaviour below yielding, then C/q would also correspond to the steady state yield stress of the material as the ow curve is ramped to zero.
Finally, the variation of the additional contribution to the yield stress s y (l) is specied. Here we select the simplest linear relationship between the yield stress s y and the dimensionless microstructural variable l where the parameter k 3 is a material constant with units of stress, and l is an additional evolving internal variable that characterizes the material microstructure. A differential equation which determines how l evolves over time is then specied. The following canonical form, which is frequently encountered in thixotropic material models, will be used here 5,17,89 The change in l is therefore determined by the competition between a microstructural build up term (rst term in eqn (22)) and a breakdown term (second term of eqn (22)). The buildup term is assumed to be the result of Brownian motion or other internal aging processes causing a rearrangement in the local wax microstructure which increases the stiffness of the gel. The breakdown term results from the irreversible plastic shearing of this microstructure. The two new material constants in eqn (22) are the rate constant k 1 (with units s À1 ) and the coefficient k 2 which is dimensionless. The constant k 1 sets the scale for the waiting time t w $ 1/k 1 required for a build up of the wax microstructure, while the constant k 2 determines the rate at which the wax microstructure breaks down for a given shear rate. An important distinction of the evolution equation proposed in eqn (22) compared to other proposed models is that only the plastic strain rate, and not the total strain rate, is responsible for destroying the material structure. This reects the fact that only large and irreversible strains should cause a signicant change in the material microstructure. Furthermore, small strains which result in primarily elastic behavior should not perturb the material microstructure, and should not have an effect on the microstructural parameter l or on the back strain A. This distinction is frequently made in the plasticity literature, 42 but is less common in thixotropic rheological models.
The evolution equation for l, also requires specifying an initial condition. The initial value of l will depend on how long the material has been aging (l can build up under the absence of any type of deformation). Typically an initial condition of l ¼ 0 is used when a deformation begins immediately aer a long history of preshearing at a high shear rate. If a deformation begins immediately aer a long aging time, then an initial condition of l ¼ 1 is used.
Eqn (22) allows for both isotropic hardening and soening. The isotropic hardening only results from a build up in structure in the material which occurs even when the shear rate is zero. This contrasts to some of the hardening equations that are used in the plasticity literature, where shear can cause an increase in the yield stress of the material. 82,85 For our model crude oil, shear only results in a nonzero soening term (i.e. the second term in eqn (22)) which causes the yield stress s y (l) to decrease.
4.1.1.1 Physical interpretation of the internal parameters. In most thixotropic models, a single variable l is frequently used to characterize the material structure. 5 The model introduced here bears more similarities to the approach taken in plasticity, 82 due to the introduction of two variables, l and A. Here we discuss the signicance of these two variables with respect to the material microstructure.
Mujumdar et al. 17 argued that the scalar parameter l quan-ties the ratio of the number of links in a transient network at a particular level of structure, to the number of links in the network when it is fully structured. For the wax microstructure, we envision l as being representative of the total number of links per unit volume between the wax particles that form aggregates. Assuming that the force required to break such links is always the same, then the yield stress s y can be approximated as a linear function of l. Because the variable l quanties the number of links, it will necessarily be a positive scalar parameter, even in a 3-dimensional version of the model. Furthermore, this scalar parameter can only be used to quantify an isotropic resistance to deformations (i.e. a resistance which is the same in all directions of loading).
However, an additional parameter is necessarily required to characterize anisotropy in the local material microstructure and the resulting non-isotropic resistance, i.e. a stiffening of the material along the direction of imposed deformation. Nonisotropic resistance to deformation is most important when considering oscillatory loading, where the direction of deformation changes (e.g. LAOS), or more complex 3-dimensional deformations. Capturing this non-isotropic resistance is the purpose of the back strain A and the resulting back stress s back , which are allowed to take on both negative and positive values for the one-dimensional case of this model (unlike l and s y which are always positive scalar quantities). Furthermore, a fully three-dimensional version of this model would generalize A to a tensor valued quantity (A), and the back stress s back to a tensor valued quantity (T back ). 46 Thus, in contrast to the scalar variable l which only measures the total number of links between wax particles, this tensor valued internal variable A additionally characterizes the orientation of the wax microstructure. This aspect of the microstructure will be affected by the orientation of the wax sheets observed in Fig. 5 (a property which cannot be described by a simple scalar), which will in turn play a role in determining how these particles resist deformation along a certain direction. The inclusion of both of these internal variables therefore provides a more complete picture of the material microstructure.

Summary of the IKH model equations.
Here we summarize the constitutive model in a concise form, and provide a list of the material constants required for tting. The decomposition of shear strain into linear viscoelastic strain and nonlinear plastic strain is as follows: The linear viscoelastic strain is specied through the introduction of an appropriate LVE constitutive element (Maxwell, Kelvin-Voigt, or simple elastic solid as shown in Fig. 11), while the rate of change of the plastic strain is given by the following: i f js À CAj\k 3 l s À CA js À CAj js À CAj À k 3 l m p ifjs À CAj $ k 3 l ; (24) where A and l are the internal variables that vary through the following differential equations: With the Maxwell model specied as the LVE behavior shown in Fig. 11(a), this model has 9 tting parameters, which are: G, h, m p , k 1 , k 2 , k 3 , C, q, m. This is one additional parameter compared to the Houska model, which is commonly used to describe the rheology of crude oils. 27 For the reader interested in more general tensorial forms of this model, the additive strain decomposition for the case of the EIKH model shown in Fig. 11(c) can be generalized to a multiplicative Kroner Decomposition 90 of the deformation gradient F such that where F e is the elastic part of the deformation gradient, and F p is the plastic part of the deformation gradient. Anand et al., 82 Henann and Anand 83 and Ames et al. 85 all give examples of frame invariant forms of constitutive models with both isotropic and kinematic hardening (although their evolution equations for the internal variables differ).

Quantitative predictions
The IKH model can be used to predict the model crude oil rheological response in the three canonical ows given in Section 3. We use data sets for the 10% model waxy crude system in its slurry state at a temperature of 27 C. First the steady state ow curve given in Fig. 6 is considered. Under steady owing conditions, the IKH model reduces to the limiting case where _ l ¼ 0 and _ A ¼ 0. The relationship between the magnitude of the plastic ow rate | _ g p | and the magnitude of the stress |s| in eqn (15) simplies to the following: In the case where the KIKH or EIKH variants of the model are used, _ g p is the only contribution towards the strain rate. Eqn (28) therefore gives a functional form of the IKH model which can be tted to steady state ow curves such as that in Fig. 6. In Fig. 12 we show such a t, and it correctly predicts the important features of the ow curve. It predicts a Newtonian limit at high shear rates, as well as locally nonmonotonic decreasing region of the ow curve at the intermediate shear rates. As the shear rate approaches zero, the EIKH and KIKH variants predict a constant stress |s| ¼ C/q + k 3 . The MIKH variant however will predict a high viscosity Newtonian region at lower shear rates with s x h _ g (which our model crude oil does not appear to exhibit). These types of high viscosity, Newtonian regions at very low imposed shear rates oen appear for thixotropic yield stress uids, 1 but are difficult to measure precisely. Tuning the LVE element in the IKH model allows for exibility in the model in capturing this type of behavior.
The values of the model parameters used in Fig. 12 are C/q ¼ 0.85 Pa, m p ¼ 0.42 Pa s, k 3 ¼ 1.5 Pa, k 1 /k 2 ¼ 0.033 s À1 and h ¼ 500 Pa s. The tting to the steady state ow curve only depends on the ratios k 1 /k 2 and C/q, and not the specic values of the individual constants. The tting also does not depend on the modulus G or the soening exponent m (or the viscosity h in the case of the MIKH and KIKH variants). It is therefore necessary to t the model to additional deformation histories in order to determine the values of these parameters.
Although the simple equation given in eqn (28) predicts a non-monotonic dependency of stress s on shear rate _ g, it does so assuming the case of homogenous shear, i.e. a spatially and temporally uniform shear rate within the material. Such a nonmonotonic response can also lead to shear banding. 91 In Section 3.2.1 we indicated experimentally that for shear rates where we observe a decrease in the average stressŝ as the average shear rateĉ g is increased, unstable ow occurs within the material and the local shear rate is not homogenous across the gap. Therefore, eqn (28) should only be viewed as a rst order approximation of what the local stress is within the material as a function of the globally applied shear rate. A more complete prediction for the stress in Fig. 12 would account for these heterogeneities and material instabilities in the ow. In Section 4.3 we briey discuss possible avenues for predicting this complex ow behavior using the IKH model.
The stress overshoots that were measured in Section 3.2.2 can also be quantitatively predicted by the IKH model using the same values of the model parameters used in Fig. 12. Values of  Fig. 11 is used). Initial conditions for the internal variables are A(t ¼ 0) ¼ 0 and l(t ¼ 0) ¼ l(t w ). Fig. 13(a) shows that the IKH model does indeed capture the key features of the response of the 10% wax oil system to startup of steady shear. Stress overshoots are predicted, and these increase as a function of waiting time t w due to the magnitude of the microstructure parameter l increasing during the waiting time as a result of rheological aging. The value of the stress overshoot also saturates beyond waiting times of t w $ 1/k 1 ¼ 10 s, in accordance with our experimental observations. Fig. 13(b) gives a quantitative comparison between the predicted value of the stress overshoot Ds 0 h s m À s ss with its measured value. At a shear rate of 2 s À1 , which is far from the unyielded state, the steady state stress is approximately s ss x C/q + m p _ g. The maximum in the stress scales as s m $ (C/q + k 3 + m p _ g), so the stress overshoot will scale as Ds 0 $ k 3 . The IKH model correctly predicts the overshoot saturating at approximately 1 Pa, and it predicts the critical waiting time t w $ 1/k 1 required for this overshoot to saturate.
One aspect of the material response that the IKH model does not predict is the long term increase in the steady measured stressŝ(t) which is observed in Fig. 8. This slow transient aging inŝ(t) can be attributed to the presence of additional longer time scales for restructuring in the wax. The IKH model can be readily modied to account for this behavior by introducing a third internal variable, and making the high shear rate plastic viscosity m p a function of this variable (similar dependencies are also adopted for some isotropic hardening models 83 ). An evolution equation for this third variable that is similar to eqn (22) could be specied. However the associated parameters determining time scales (e.g. the aging rate k 1 ) would have to be modied in order to reect this slow restructuring. This would result in the introduction of several new material parameters, making the model signicantly more complex. For the sake of simplicity and compactness, we choose to focus on predicting the relative stress overshoot Ds 0 (which is more important for ow assurance restart procedures), and therefore avoid introducing a third internal variable. The introduction of a third internal variable also requires a subtle distinction in the microstructural interpretation for the l variable used here. This can be justied by considering rearrangement of the individual wax crystals vs. rearrangement of ocs of crystals. Due to the difference in size, the ocs rearrange much slower than the individual crystals, and have different impacts on the uid rheology (e.g. the slow oc rearrangement will impact long term changes in steady state viscosity whereas the local rapid rotary diffusion of the wax crystals will affect initial yielding transients).
We also assess the predictive capabilities of the IKH model with respect to the measured LAOS response of the waxy crude oil presented in Fig. 10. For the LAOS response, the values of the model parameters used for the ts in Fig. 12 and 13 are kept constant, other than k 3 and C/q which are both set to 0.7 Pa (while keeping C constant at 70 Pa). With a value of k 3 ¼ 1.5 Pa, the general features under LAOS are still captured (see supplemental information), however a lower value of k 3 more accurately predicts the absolute value of s m under LAOS. In Fig. 14(a), the t of the MIKH variant to the LAOS response is shown for the 10% wax oil system in its slurry state at 27 C. The simulated LAOS response of the IKH model always uses the initial conditions of A ¼ 0 and l ¼ 1. For visual clarity Fig. 14 only illustrates the limit cycle, or alternance state of these LAOS tests.
There is good agreement between the experimental data and the IKH model for the limit cycles observed under LAOS. Both material and model exhibit a transition from linear viscoelastic behavior (corresponding to elliptical Lissajous curves) to nonlinear behavior at strain amplitudes of approximately g 0 ¼ 0.2%. This very small linear range of strains implies that at moderate strain amplitudes (1-10%), the IKH response is dominated by the behavior of the nonlinear yielding element with strain g p .
In Fig. 14(b) we compare the values of the maximum stress s m determined from the experimental data and from the IKH model. The maximum stress s m is dened as the maximum stress over all cycles of oscillation, so it includes the initial Fig. 13 Prediction of stress overshoots under startup of steady shear for the IKH model. In (a), the imposed shear rate isĉ g ¼ 2 s À1 . In (b), we plot the predicted stress overshoot (red line) from the IKH model simulated in (a) vs. the experimental data from Fig. 9 (circles). The model parameters are the same for (a) and (b), and are as follows: transients (overshoots) that are not included in the alternance plots of Fig. 14(a) (but are shown below in Fig. 15). Fig. 14(b) includes the predicted values of s m for a Maxwell LVE element (with h ¼ 500 Pa s and G ¼ 250 Pa), and a Bingham model (with m p ¼ 0.42 Pa s and s y ¼ 1.4 Pa). The MIKH model reduces to the LVE model for small strains, and to the Bingham model for large strains, while also maintaining the ability to predict the behavior at intermediate strains. Fig. 14(c) decomposes the "sequence of physical processes" 77 that the IKH model undergoes as the stress and strain vary under LAOS at moderate strain amplitudes. Starting from the lower le point of the curve, the material undergoes an initial (primarily) elastic loading, characterized by a slope/modulus of G ¼ 250 Pa. When the stress has increased by a value of 2k 3 , the material yields in the forward direction. As soon as yielding starts, there is an onset of kinematic hardening (relative to a perfect yield plateau), leading to an increase in the stress with ds dg xC initially. The recovery term in the evolution equation for A (eqn (19)) causes this initial slope to decrease as the imposed deformation continues to increase. A rhomboidal shape in the curve arises, and this is indicative of the presence of the Bauschinger effect. 42 For low shear rates (i.e. m p ug 0 ( (k 3 + C/q)), the maximum stress s m never exceeds the value of k 3 + C/q, corresponding to the limiting value of static yield stress in Fig. 12 at low shear rates. For larger strain amplitudes (and hence larger shear rates), the plastic viscosity contribution to the stress m p _ g becomes more prominent, resulting in Lissajous curves with a diminished rhomboidal shape and increasing convexity.
The IKH model is also capable of predicting the initial transient response of the material as it undergoes large amplitude deformations. Before each test is carried out at a given amplitude g 0 , a waiting time t w ¼ 100 s is applied, resulting in the material being fully structured before the oscillatory deformation is imposed. For this reason, the IKH simulation begins with the initial condition l ¼ 1, and over the course of several oscillations l will decrease as a result of a nonzero plastic strain rate _ g p , resulting in a decrease in the size of the Lissajous orbits. Fig. 15 shows two representative transient Lissajous curves of the material undergoing an initial three cycles of deformation at two different amplitudes, g 0 ¼ 0.1% and g 0 ¼ 20%. For the small amplitude case in Fig. 15(a), the material requires one half of a cycle of oscillation before the behavior settles into a steady state. Furthermore, no overshoots are observed in the stress when the oscillatory deformation is initialized. This behavior is consistent with a standard linear viscoelastic  response, where no thixotropy is present. The IKH model predicts this behavior, validating our argument that only the plastic contribution to the strain g p should appear in the evolution equation of l in eqn (22). Thixotropic effects only become apparent at larger strain amplitudes when g p is no longer zero. These thixotropic effects are apparent in the LAOS measurements at a higher strain amplitude in Fig. 15(b). An initial stress overshoot is observed, and multiple cycles of oscillation are required for the material to settle into an alternance state. The IKH model quantitatively predicts the magnitude of this overshoot, and also predicts the slow contraction of the Lissajous orbits towards the ultimate limit cycle. At a given frequency, the ratio of parameters k 2 /k 1 ¼ 30 s controls the rate of the contraction; as k 2 /k 1 increases the orbits will contract more rapidly.
The LAOS tests in Fig. 14 only show the dependency of the material response on the imposed strain amplitude g 0 . The oscillation frequency u was also varied at a moderate strain amplitude (g 0 ¼ 2%) to further illustrate the nature of the material's response in the nonlinear regime of deformations. The specic strain amplitude of 2% was chosen because it lies in the intermediate regime of Fig. 14(b) between linear viscoelastic (LVE) behavior and Bingham-plastic behavior. In each of these limiting regimes, the effect of varying u is well understood. In the LVE regime, varying u enables the appropriate viscoelastic mechanical model for the material properties below yield to be probed. In the Bingham-plastic regime, the dependency on u is dominated by the plastic viscosity m p , since the instantaneous shear rates are high leading to increasing circular trajectories (as shown in Fig. 14(a)).
However, in the intermediate regime where strains are between 1 and 20%, the material exhibits a rate-independent behavior. Fig. 16(a) illustrates this by showing the cyclic Lissajous curves of the rst 3 cycles of oscillation for the material undergoing LAOS at 4 frequencies spanning an order of magnitude (u ¼ 0.5 rad s À1 to u ¼ 5 rad s À1 ). At this moderate strain amplitude, the material response is clearly nonlinear, with a small amount of thixotropic contraction occurring in the orbits over multiple cycles. However, there is little dependency of the material response on the frequency u. This is indicative of a deformation rate-independent response at moderate strain amplitudes. The predicted model response shown in Fig. 16(b) captures this type of behavior. This is because of the rateindependent nature of the kinematic hardening behavior prescribed in eqn (19). This equation can be rewritten in terms of increments in A and g p as follows: The back strain A therefore evolves independently of the rate of plastic strain _ g p . For the strain amplitude shown in Fig. 16, s m x 1 Pa, so g ve $ (s m /G + s m /(hu)) x 0.006 when u ¼ 1 rad s À1 . A large proportion of the total imposed strain, g 0 ¼ 0.02, is therefore taken up by the plastic strain g p ¼ g À g ve . As a result, the behavior at this strain amplitude is primarily rate independent. The other internal variable, l, evolves in a ratedependent fashion set by the ration k 2 /k 1 , but on a time scale longer than the period of a single orbit 2p/u. Hence, the overall response of the material at this range of frequencies is characterized by rate-independent kinematic hardeninga behavior which is most appropriately described by equations of the Armstrong-Frederick type 86 used in eqn (29).
The IKH model predicts all of the important features of the material response to the three canonical rheological ows outlined in Section 3. To end this subsection, we will discuss one additional ow, and compare it to the predicted response of the IKH model.
In Fig. 6, the particular experimental protocol that was used to measure the ow curve of the material controlled the applied global (or spatially averaged) deformation rateĉ g on the material. An alternative way to control the ow in the material would be to ramp the applied average stressŝ on the material at a very slow rate, and measure the corresponding average shear rate at subsequent points in time. This "stress ramp" test was conducted with the 10% wax-oil system in its slurry state at 27 C using the AR-G2 in its standard conguration. The material is initially brought to its fully structured state by applying a wait time of t w ¼ 300s, and then a stress ramp from 0.1 Pa to 5 Pa is imposed over the course of 2 hours (this results in a constant value of dŝ/dt ¼ 0.00136 Pa s À1 ). Once the stress reaches the maximum value of 5 Pa, the ramp is immediately reversed and brought down to 0.1 Pa at a rate of dŝ/dt ¼ À0.00136 Pa. The average shear rate within the uid (as measured by the rate of rotation of the conical xture) is sampled at regular intervals over the course of the increasing and decreasing ramps in stress.
In Fig. 17, we show the experimental results of this stress test, as well as the predicted response of the IKH model to such a stress ramp using the same model parameters that were determined from Fig. 12 and 13. For comparison with the average shear rate measurements, the ow curve obtained in Fig. 6 is also included on this gure. The non-monotonicity in the IKH model results in hysteresis being captured in the stress ramp experiments. Specically, a larger imposed stress (asymptotic value of s s x 2 Pa) is required to start the ow from the structured state, and as the stress is decreased, the same shear rates can be sustained at lower values of the applied stress (asymptotic value of s d x 1.1 Pa). This essentially corresponds to different measurements for the static and dynamic yield stress within the material, 92 and the nonzero value of the microstructural parameter l is responsible for this difference. While s approaches the limit of C/q + k 3 ¼ 2.35 Pa in the limit _ g / 0 for the rate-controlled experiment, the asymptotic values of s determined from the stress ramp experiments are both different from this value. The difference between the static and dynamic yield stress is responsible for the "avalanche effect" that is frequently observed in thixotropic yield stress uids 93the supplementary information shows how the IKH model can predict this avalanche effect.
The microstructural physics captured by the IKH model is thus extremely versatile and is able to capture the complex rheology of a structured and aging material such as an elastoviscoplastic crude oil. It can capture the salient features of the response exhibited by the uid to a number of different steady and time-varying rheological ows. The numerical values of the model parameters were consistent among these different ows, so the measurements in Fig. 17 constitute a true test of the predictive nature of the model.

Further discussion
To conclude our discussion of the IKH model, we expand on a few important points related to the current form of the model. We will also comment on how this model is related to other constitutive models that are commonly encountered in the literature.
4.3.1 Fitting procedure. We outline a general tting procedure used to determine values of the model parameters in the IKH model. This procedure can be summarized as follows: We rst t the parameters m p , k 3 and the ratios C/q and k 1 /k 2 to the steady state ow curve in Fig. 6. Next frequency sweeps at low strains below the yield point can be used to t the linear viscoelastic behavior. If the material is Maxwell-like, then G, h and the corresponding viscoelastic relaxation time h/G can be determined. Subsequently, LAOS can be used to determine appropriate values of the back stress modulus C, and the constants q and m which control the shape of the Lissajous curves. The rate constants k 1 and k 2 which control the rate of thixotropic contraction towards a nal periodic alternance state can also be determined from the LAOS measurements. Finally, stress overshoot measurements can be used to verify that the value of the aging parameter 1/k 1 agrees with the critical waiting time t w required for the saturation of the yield overshoot or yield peak. However, the tting of the parameters to the steady state ow curve of Fig. 6 may proceed differently for different classes of thixotropic yield stress uids. Specically, the nature of the unstable ow that was observed in Section 3.2.1 may differ among different classes of so solids and gel-like materials. For example, an alternative shear banding scenario may be observed where the non-monotonic region of the ow curve is inaccessible to time-averaged measurement.
4.3.2 Accounting for shear heterogeneities. The steady state ow curve of the IKH model show in Fig. 12 represents a homogenous shear ow, however, at the lowest globally applied shear rates, PIV measurements indicate that the ow of the model crude oil is not homogenous. The nonmonotonicity of the IKH model means that it can predict these types of shear heterogeneities. One way to account for these would be to initially specify a spatially heterogenous value of l(y, t ¼ 0), and then allow for the material elements to evolve under application of a constant global shear rate. The IKH model predicts transient shear banding under such a scenarioregions with initially high values of l will not deform, but regions with low values of l will exhibit high shear rates. This shear banding would be extremely sensitive to initial spatial variations of l(y). A heterogenous distribution of the l parameter is a reasonable assumption, given the observations in this work, as well as previous observations that have shown thixotropic yield stress uids exhibiting an "erosion" behavior that is oen spatially heterogenous 38 Shear banding is also predicted by the IKH model under deformations in geometries with sufficiently large stress gradients. We illustrate this in Fig. 18 by simulating startup of shearing ow for the IKH model in a Taylor-Couette cell under a steady applied torque T . For the simulations in Fig. 18, the inner radius of the cell is set to R in ¼ 23.7 mm and the outer radius is R out ¼ 25 mm, corresponding to a typical experimental cup and bob. The dimensionless radial position is dened as r h (r À R in )/(R out À R in ), and thus varies from 0 to 1. The torque is set so as to impose a stress which varies (in a 1/r 2 fashion) from s ¼ 2.05 Pa at r ¼ 0 to s ¼ 1.85 Pa at r ¼ 1. The IKH simulation presented in Fig. 18 assumes that the uid is initially in a fully structured state with spatially uniform values of l ¼ 1 and back strain A ¼ 0. The model parameters used The data in Fig. 18 is plotted against dimensionless time t/t c , where t c ¼ k 2 /(k 1 q). The timescale t c arises from equating the natural timescales for the kinematic hardening and isotropic soening/hardening processes.
At time t/t c ¼ 0, the step increase in the torque is imposed and the material starts to deform away from its fully structured state (corresponding to l ¼ 1). There is initially a nonzero shear rate across the gap, which decreases up until t/t c ¼ 1. During this initial deformation stage, the material deforms elastically, begins to yield, and then kinematically hardens, resulting in an increase in A and a strengthening of the material in the ow direction. Consequently, the shear rate across the gap decreases during this time period.
For t/t c > 1, the uid begins to shear band. For the region of the uid located closest to the inner wall, i.e. for values of r < 0.4, the isotropic soening process (accounted for by the term Àk 2 | _ g p |l in the evolution equation for l) begins to dominate over the kinematic hardening process. The value of l decreases in this region, and a high shear rate band forms in the material closest to the inner wall. For times t/t c [ 1 the material near the inner wall has fully yielded, and at r ¼ 0 the velocity reaches its maximum value, v ¼ v m ¼ 1.52 mm s À1 . The material closer to the outer wall experiences lower shear rates and continues to harden and age, until it exhibits a negligibly small shear rate at long times. These steady shear banded proles (with one strongly sheared band and one stationary, unsheared band) have been observed in numerous thixotropic yield stress uids; Fig. 17 Rheological measurements for an imposed stress ramp shown in (a), and measured (points) and simulated (lines) shear rates in (b). The different static and dynamic yield stresses are also annotated. The resulting flow curves of the model crude oil (red data points) and the IKH model (solid and dashed lines) are given in (c). The data is for the same model system at the same temperature shown in Fig. 12-16, with the same model parameters as in Fig. 12 and 13. for example they are carefully documented in the work by Martin and Hu. 94 4.3.3 Genealogy of constitutive models. The general IKH framework we have described in this work can be simplied to several other constitutive models under various limits. It can be reduced to the Bingham generalized Newtonian uid model by setting k 3 ¼ 0, and then taking the limit of G / N and C / N, while keeping the ratio C/q constant and equal to the yield stress s y . This results in the back stress s back immediately responding to a given stress level s y , by either saturating at AEs if |s| < s y or saturating at AEs y if |s| > s y . An alternative way to reduce the IKH model to a Bingham model would be to set the back stress to zero s back ¼ 0 and set l to a constant.
The IKH model can also be readily simplied to the kinematic hardening model used by Dimitriou et al. 46 to describe non-thixotropic yield stress uids. This is accomplished by setting k 3 ¼ 0. While the KH model was shown in that work to predict the correct LAOS behavior for a "simple yield stress uid", 91 it is not able to predict stress overshoots or prolonged transients under multiple cycles of large amplitude deformation.
The limit of the inelastic Houska model studied by Cawkwell and Charles 9 and others 33,41 can be realized by taking the limit of C / N while keeping the ratio C/q constant and equal to their "xed" component of the yield stress. The dependence of s y on l however remains with an evolution equation for l resulting in the second component of the yield stress not varying. The Houska model also does not account for any elastic behavior below yield, corresponding to the limit G / N.
Finally, the IKH model can be reduced to the "toy" l-model employed by Coussot et al. 95 This limit can be achieved by setting both yield parameters s y and s back equal to zero, taking the limit of G / N, and making the viscosity m p a function of l. Coussot employs a different evolution equation for l, which results in l growing unbounded in time when the shear rate is zero.
The IKH model is more broadly related to the Mujumdar 17 and de Souza Mendes 19 models. These models exhibit key differences from the IKH model. The Mujumdar model does not explicitly decompose strain into elastic and plastic components as a result the authors necessarily introduce a critical strain g c into the model which is used to determine whether or not the structure of the material is breaking. The Mujumdar model does utilize the same form of the evolution equation for l that is employed here.
The de Souza Mendes family of models 19,57,96 use a different form for the evolution equation of l, specically the stress s enters this equation and is assumed to be the driving force for disrupting the material microstructure. de Souza Mendes and Thompson take this approach in order to avoid disruption of the material microstructure in any elastic, reversible regime of deformation. 57 The use of _ g p in our evolution equation in l achieves this same goal. de Souza Mendes also writes the linear viscoelastic moduli G and h as functions of l, resulting in additional material parameters that must be determined by regression to data. 19 This may be required for a more complete form of the IKH model, however for simplicity and compactness these dependencies are neglected in the IKH model. Fig. 18 Shear banding predicted by the IKH model in a Taylor-Couette concentric cylinder geometry. A steady imposed torque T is applied, resulting in a steady imposed, heterogenous stress in the material. The material is allowed to evolve from its initial configuration of being fully structured, i.e. l ¼ 1. In (a), a space-time profile of the velocity field is shown. In (b) and (c) several velocity profiles are given for different instants in time. Velocity is plotted in a non-dimensional form using the maximum velocity for long times at near the inner wall, v m ¼ 1.52 mm s À1 .

Conclusions
The primary goal of this work has been to develop a suitable elasto-viscoplastic constitutive model that can predict the rheology of a thixotropic yield stress uid. The uid studied here is a waxy crude oil that is relatively easy to formulate and prepare in a repeatable manner, but still captures the rheological complexity of a real crude oil.
To develop this model, Rheo-PIV measurements were rst used to distinguish between the material behavior under different sample preparation protocols. We distinguished between the material in a "slurry" state vs. a "strong gel" state. In the former, the material exhibits a weaker gel network and is less likely to ow through the mechanism of interfacial wall slip. In the latter, the material can undergo irreversible breakdown of an initially strong gel network. Furthermore, the material in its "strong gel" state is more likely to exhibit localized fracture and failure at the uid/solid interface.
Second, we outlined three canonical ow scenarios to probe the rheology of the model uid in its slurry state. These are measurements of the steady (spatially averaged) ow curve of the material (ŝ vs.ĉ g ), startup of steady shear ow following different waiting times t w , and large amplitude oscillatory shear strain (LAOStrain). Furthermore, we used Rheo-PIV measurements to elucidate the nature of ow instabilities that occur in the material at low shear rates. These coincide with a nonmonotonic region in the material's ow curve. The material instabilities consist of spatially and temporally uctuating values of local shear rate and shear stress within the material, with a constant globally applied shear rateĉ g and average shear stressŝ.
Third, we developed an elastoviscoplastic constitutive model (referred to for simplicity as the IKH model) which quantitatively captures the material response to these three ow scenarios. The key features of this model are the additive decomposition of strain into two components (g ve and g p ) and the use of the isotropic and kinematic hardening mechanisms adapted from plasticity theory. An added benet of this approach is that the model is extendable to a frame invariant, 3dimensional tensorial form that can be utilized in simulations of more complex ow scenarios. The IKH model was tted to the 10% wax in oil model system in its "slurry state" at 27 C. Under these isothermal conditions it was capable of capturing the rheological aging behavior exhibited by the uid, and the subsequent effects of this aging process on the material response to deformations. Values of the model parameters were generally consistent across the three different rheological ows. The IKH model was also shown to be capable of predicting transient shear banding behavior under deformations in geometries with sufficiently large gradients in the stress.
There are several avenues for a continuation of this work on the IKH modelsome of which we have already highlighted. These include accounting for additional slower restructuring timescales which affect the high shear rate plastic viscosity. In the present study we have taken this viscosity to be a constant, m p . The IKH model could also be improved by introducing nonlocal terms into the constitutive model, to better describe the nature of the ow instabilities that are observed at low globally-applied shear rates. 18,97 Other than the results presented in Section 4.3, all of the predictions of the IKH model were obtained by assuming homogenous ow. The work by Moorcro et al. 98,99 suggested that instances of inhomogenous ow may be more ubiquitous than previously thought (in particular where stress overshoots are present). Thus, understanding how the IKH framework might be able to predict these instabilities in more detail would be a worthwhile endeavor. Finally, the model could be modied to capture the continuous, temperature-dependent transition that the material exhibits as it is cooled from above T wa (where it behaves as a Newtonian liquid) to below T wa where it begins to exhibit thixotropic, elastoviscoplastic behavior. The ultimate goal of such work would be to develop a temperature-dependent IKH model, which can then be implemented into simulations of nonisothermal ow in pipelines for ow assurance applications. This would substantially build on the temperature-dependent models that are available today, which typically assume a generalized Newtonian uid and do not account for elastic or yielding behavior. 28,34 Practical applications for such a model would be widespread, and would greatly assist in optimizing operating conditions for pipelines that transport waxy crude oil.