High-shear rate rheometry of micro-nanofibrillated cellulose (CMF/CNF) suspensions using rotational rheometer

Suspensions of cellulose micro- and nanofibrils are widely used in coatings, fibre spinning, 3D printing and as rheology modifiers where they are frequently exposed to shear rates > 104 s−1, often within small confinements. High-shear rate rheological characterisation for these systems is therefore vital. Rheological data at high-shear rates are normally obtained using capillary and microfluidic rheometers, which are found in relative scarcity within research facilities compared to rotational rheometers. Also, secondary flows and wall depletion prevalent at such high-shear rates often go unnoticed or unquantified, rendering the measurement data unreliable. Reliable high shear rate rheometry using rotational rheometers is therefore desirable. Suspension of TEMPO-oxidised CMF/CNF was tested for its high-shear rate rheological properties using parallel plate geometry at measurement gaps 150–40 µm and concentric cylinder at 1 mm gap. The errors from gap setting, radial dependence of shear stress and wall depletion were quantified and accounted for. Viscosity data from 0.1 to 30,000 s−1 shear rates was constructed using both geometries in agreement. Possibilities of secondary flows, radial migration of fluid and viscous heating were ruled out. Steady shear flow data of CMF/CNF suspension from 0.1 to 30,000 s−1 obtained using rotational rheometer


Introduction
Cellulose is the most abundant polymer in nature, and has been used historically for its mechanical properties (wood), as a fuel and as a raw material for paper making (Silva et al. 2012). Cellulose nanofibrils (CNFs) and microfibrils (CMFs) are the flexible elongated structural constituents of cellulose which are of widespread interest across various sectors of the industry owing to their mechanical properties, biodegradability and low density (Jonoobi et al. 2015).
The production of CNFs involves disintegration of cellulose fibres into their constituents microfibrils and consecutively, nanofibrils. Various mechanical processes such as high-pressure homogenisation, processing through Microfluidiser and grinding have been employed to break down the hydrogen bonds and cell wall structure present within cellulose. These methods involve application of high-shear forces to swollen cellulose fibres under small confinements (Abdul Khalil et al. 2014). High intensity ultrasonic treatment is also used for isolating CNF (Masruchin et al. 2015).
Although the energy input for CNF production using a given method chiefly depends upon the botanical source of cellulose as well as the pulping process, the abovementioned methods in general consume vast amounts of energy, making the industrial scale production of CNFs often prohibitively expensive. Various pre-treatment approaches have therefore been developed over the years to reduce the energy consumption. TEMPO mediated oxidation is one of these pre-treatment processes in CNF processing that has become popular in literature because of the ease of preparation and the quality of CNFs, although the cost of reagents has limited its industrial applications. Other pre-treatment routes and different processing methods have been extensively compared (Sacui et al. 2014;Onyianta et al. 2017) and reviewed (Klemm et al. 2011;Lavoine et al. 2012;Nechyporchuk et al. 2016).
High aspect ratio, high surface area and the hygroscopic nature of CNFs and CMFs (due to the presence of accessible hydrogen bonds) contribute towards the ability to form high-viscosity aqueous suspensions at relatively low solid contents of 1-5% (w/w) (Nechyporchuk et al. 2016). The nanofibrils absorb water within their nanoporous structure (Dimic-Misic et al. 2013a), which increases the effective solid volume fraction and form a sample spanning network that rheologically manifests as a viscoelastic solid. The deformation/break-down of this network during steady shear flow causes viscosity to reduce, giving rise to rheological properties such as shear-thinning (Pääkko et al. 2007;Lasseuguette et al. 2008;Saarikoski et al. 2012;Martoïa et al. 2015;Kumar et al. 2017b) and thixotropy (Iotti et al. 2011;Shafiei-Sabet et al. 2016). At small applied stresses, the network is considered to be elastically deformed, i.e. the strain is recovered upon removal of stress. The network supposedly undergoes irreversible deformation at stresses r [ r y , referred to as the yield stress. Yield stress behaviour has been reported for both CNF and CMF suspensions (Karppinen et al. 2011;Dimic-Misic et al. 2013b;Kumar et al. 2016c;Taheri and Samyn 2016).
Such structural and rheological properties make CNF a material of choice in several applications such as fibre spinning (Lundahl et al. 2016(Lundahl et al. , 2017, coating (Gruneberger et al. 2014;Kumar et al. 2016a;Oh et al. 2017) as well as 3D printing (Leppiniemi et al. 2017). All these applications involve processing of CNF suspension at shear rates between 10 3 and 10 6 s -1 within small confinements (Kumar et al. 2016c;Lundahl et al. 2016;Leppiniemi et al. 2017). CNFs and CMFs are used as rheology modifiers in several products which may get exposed to shear rates greater than 10 3 s -1 during the operations such as pumping and filling containers (Barnes et al. 1993). Moreover, the processing of CNF using the Microfluidiser involves high-pressure pumping of cellulose suspension through flow channels with diameters generally less than 200 lm (Saarinen et al. 2009;Naderi et al. 2014aNaderi et al. , b, 2015aTaheri and Samyn 2016;Onyianta et al. 2017).
An accurate knowledge of the rheological behaviour of CNF/CMF suspensions at high shear rates is therefore necessary to study such processes of industrial relevance. For example, CNF hydrogels are used for spinning filaments where a high shear rate in the spinneret promotes nanofibril orientation which improves the mechanical properties of the filaments. However, the flow regime shifts from laminar to turbulent flows (defined by erratic motion and swirls) at even higher shear rates, thereby lowering the overall orientation and leading to weaker mechanical properties (Lundahl et al. 2017). Designing such industrial processes is therefore only possible with reliable high shear rate rheology data. Nevertheless, in-depth rheological characterisation of CNF/CMF suspensions at high shear rates with microstructural interpretations has not been extensively reported in the literature. Most rheological data available for such systems are limited in the order of 10 3 s -1 of shear rate.
One of the major reasons for lack of such data is the scarcity of 'conventional' high-shear rate rheometers (e.g. capillary extrusion rheometers and microfluidic rheometers) within the research facilities accessible to the wider scientific community. Most of the high shear rate rheological data reported for industrial systems in general and CNF/CMF suspensions in particular have been obtained by so-called conventional high-shear rate rheometers. Such instruments require specialist training for operation and are far less available in laboratories compared to more common rheometers such as controlled-stress and controlled-strain rotational rheometers. A method to achieve high shear rate measurement data using rotational rheometers could therefore prove to be a very useful tool.
Since shear rate is the ratio of flow velocity to the measurement gap, high shear rates can be achieved by (i) increasing the flow velocity within the geometry and/or (ii) decreasing the flow gap of the geometry. Nonetheless, each of these approaches have their limitations. At high flow velocity, the fluid inertia often dominates the flow kinematics. Whereas slip layer (the solvent rich and particle-depleted layer near the wall) thickness as well as slip velocity become more relevant at narrow flow geometry gaps. Ignoring the quantitative contribution of such factors may lead to implausible interpretations of high shear rheometry data. For example, Samyn and Taheri (2016)   studied the flow behaviour of fibrillated cellulose suspensions at up to 2000 s -1 shear rate using concentric cylinder geometry. They observed apparent shear thickening at shear rates near 1000 s -1 for many of their samples which was attributed to the dynamic interactions between the cellulose fibrils and other organic nanoparticles present in the system. However, a more plausible explanation for such behaviour is the presence of Taylor vortices, i.e. an axisymmetric secondary flow induced by inertia which cause the measured torque to increase in a couette flow, thereby exhibiting an apparent shear thickening. The onset criteria for Taylor vortices is defined by Taylor number Ta: Here q is the fluid density, X is the angular velocity, R o and R i are the radii of outer and inner concentric cylinders respectively and g is the viscosity at the shear rate corresponding to X. The systems with Ta [ 3400 are considerably prone to the Taylor vortices (Macosko 1994) even though the conservative limit is Ta [ 1700 (Ewoldt et al. 2015). Based on the experimental conditions reported by , we calculated Ta for the measurement made at 1000 s -1 to be 3752, i.e. their measurement at that shear rate was highly prone to the presence of Taylor vortices. Similar comments can be made on the data reported by .
There have been some reports on high shear rate rheometry of CMF/CNF suspensions using slot die rheometers (Kumar et al. 2016a(Kumar et al. , b, d, 2017b and vane rotor (Nazari et al. 2016). Many of these reports correctly identified the possibility of heterogeneous flow profiles across the measurement gap owing to wall depletion phenomenon (i.e. formation of solvent rich and particle depleted layer over the boundary) and its effect on the data. The presence of wall depletion cannot be avoided in many processes involving colloidal dispersions. In coating systems, formation of a low viscosity layer depleted in particles at the boundary affects lubricating flows at high stresses by causing unusually high velocities (Triantafillopoulos et al. 2001). Wall depletion may also affect pumping through pipes. Estimation of slip layer thickness and slip velocity becomes necessary as a result.
It is possible (although often difficult) to eliminate/ minimise the heterogeneous flow in particulate systems within the context of rheological measurements using geometry with roughened walls. In the absence of such measures, the data must be resolved into two separate components, i.e. bulk flow and slip (wall depletion) flow (Barnes 1995). However, the contribution of slip on the viscosity data of CNF/CMF suspensions at high shear rates has never been reported to be quantified.
Using capillary rheometer, Iotti et al. (2011) (Iotti et al. 2011) reported the viscosity of CMF suspensions between 1.9 9 10 5 and 3.5 9 10 5 s -1 shear rate range and Ming et al. (2016) (Ming et al. 2016) reported the viscosity of CMF containing clay suspensions between 3 9 10 5 and 1.8 9 10 6 s -1 shear rate range. Both studies observed shear thickening within abovementioned shear rate ranges. However, no satisfactory microstructural explanation was provided for this effect. These effects are more likely to be a result of turbulent flow at such high shear rates rather than any other phenomenon (Barnes 1989;Macosko 1994;van Vliet 2013).
The present work reports high-shear rate rheometry of TEMPO-oxidised CMF/CNF suspension using parallel plate geometry at several measurement gaps, with some as low as 40 lm. Many different corrections to the data were made owing to the characteristics associated with the measurement technique (e.g. non-parallelism and non-flatness of plates and the gapzeroing procedure), the nature of fluids (e.g. viscosity being a function of shear rate/stress) and discontinuities within the shear flow profile (e.g. slip, also known as wall depletion).
The systematic error in rheometer gap settings has been determined by measuring the viscosity of polydimethylsiloxane oil at different measurement gaps. Since shear rate in parallel plate geometry varies with the distance from the rim, the apparent shear rate values calculated from the angular velocity are applicable only to Newtonian fluids. These values were corrected in the present work for a non-Newtonian sample such as the CMF/CNF suspension. Wall depletion effects become even more dominant at low measurement gaps, leading to an apparent reduction in viscosity with decreasing gap at the same shear rate. Such trend was identified and corrected by extrapolating the observed shear rate to an 'infinite' measurement gap. Following all the corrections, a steady shear flow curve with measurement data at shear rates up to [ 30,000 s -1 was generated. As measurements performed using parallel plates at low gaps have relatively higher ratio for sample surface to sample bulk, heat exchange is more efficient and therefore viscous heating at high shear rates is less likely to be present (Pipe et al. 2008).
The abovementioned corrections have been explained in detail in the Theory section.

Determination of gap error
Rotational rheometry requires the gap between both upper and lower components of the geometry to be 'zeroed' prior to sample loading, i.e. the point of their contact is established as the 'zero gap' which can then be used as a reference for setting non-zero gaps. The point of solid-solid contact between both plates is determined by an abrupt increase in the normal force. However, such change in normal force can even be observed before the actual contact owing to the viscosity of air resisting the squeezing gap between two plates resulting into an underestimation of set gap. Such an error can also stem from either or both the plates being non-flat and/or them being non-parallel to each other which is almost unlikely to be noticed visually.
Such errors typically have a magnitude of the order of 10 lm. Their significance can often be ignored when working at gaps of 1000 lm or greater but must be taken into account at gaps lower than 200 lm. An underestimation of measurement gap leads to overestimation of apparent shear rate and lower apparent viscosity at gaps typically used in high shear rate measurements. It is therefore essential to determine the gap error in the present study in which gaps much lower than 200 lm are to be used. The method to determine the gap error has been reported earlier (Kramer et al. 1987;Davies and Stokes 2008;Pipe et al. 2008) which involves measuring the viscosity of an apparently Newtonian fluid such as silicone oil (polydimethylsiloxane).
The shear stress r and the apparent shear rate _ c app are calculated using the following relationships in parallel plate rheometry: Here M is the torque and R is the radius of the upper plate. d is the commanded measurement gap. As discussed, the true shear rate _ c true (calculated at d þ e gap instead where e is the unknown gap error) is different to the apparent shear rate: Based on Eqs.
(2), (3) and (4), apparent (g app ) and true (g true ) viscosities can be defined as: Equations (3), (4), (5) and (6) lead to the following expression: Equation (7) can be rearranged to that of a straight line: Hence when plotting d/g app as a function of d, the true viscosity and the gap error can be determined using the slop 1/g true and the intercept e/g true of the fitted straight line respectively. Following the estimation of e, Eq.
(3) can be converted to Eq. (4) by multiplying the former with d/d ? e.
Non-Newtonian nature of the sample Equation (3) indicates that the shear rate (and as a consequence shear stress) in a parallel plate geometry is not constant across the geometry at a given angular velocity but instead varies linearly from 0 at the centre to the maximum at the rim. Therefore Eq.
(2) is valid only for Newtonian fluids which exhibit linear relationship between shear rate and shear stress. The shear stress versus shear rate data for non-Newtonian fluids obtained using parallel plate must be corrected for the heterogeneity of shear stress in the parallel plate geometry by the expression below, which is widely known as Weissenberg-Rabinowitsch correction (Kramer et al. 1987;Davies and Stokes 2008;Pipe et al. 2008;De Souza Mendes et al. 2014): Here r corrected is the corrected shear stress. Equation (9) can be seen as a generalised expression of Eq.
(2) as the term d ln M/d ln X is unity for Newtonian fluids and therefore Eq. (9) reduces to Eq. (2). For power-law fluids, the term can be taken as the exponent of the power-law and for any other type of non-Newtonian fluids, the term varies with angular velocity and needs to be calculated as the first derivative of ln M versus ln X data. The significance of this often-ignored correction can be realised by an example in which a corrected shear stress r Corrected for a power-law fluid with the exponent value of 0.4 is 15% lower compared to the recorded shear stress. Such difference can be as significant as 25% where d ln M/d ln X is zero. It is very common for shear thinning CMF and CNF suspensions to have powerlaw exponent values around 0.4-0.6. Despite this, majority of the published reports on the rheology of CNF or CMF suspensions which use parallel plates do not perform the above correction according to the best of our knowledge (Iotti et al. 2011;Missoum et al. 2012;Dimic-Misic et al. 2013a, b;Rezayati Charani et al. 2013;Kumar et al. 2016bKumar et al. , d, 2017aNazari et al. 2016).

Quantifying wall-depletion effects
A depletion in the concentration of dispersed particles near the geometry wall compared to that in the bulk is a commonly encountered phenomenon in colloidal dispersions (Mewis and Wagner 2012). As a consequence, the shear gradient across the measurement gap is discontinuous (higher in the solvent-rich outer surface and lower in the bulk). The slip (more appropriately wall-depletion for particulate systems) effect is expected to become more dominant at lower gaps, leading to apparent lowering of viscosity with decreasing gap. The contribution of wall-depletion effects within apparent viscosity at different measurement gaps can be considered as a systematic error and therefore can be calculated by the following expression referred to as the Mooney method (Davies and Stokes 2008): Equation (10) suggests plotting the measured shear rate _ c meas against the reciprocal of corrected measurement gap (at a constant shear stress for all gaps) will allow one to calculate the shear rate in bulk sample _ c bulk and the slip velocity t s if a straight line is obtained, i.e. the slip velocity remaining constant at all gaps for a given shear stress. This could also be seen as extrapolating the measured shear rate to 'infinite' gap where the contribution of wall depletion effects tends to diminish.
The slip velocity can subsequently be used to calculate the thickness of the depleted layer h dep using the following relationship (Soltani and Yilmazer 1998;Bécu et al. 2005): Here g medium is the viscosity of the continuous medium found in the slip layer.

Turbulent flow
Majority of flows are laminar at moderate shear rates, i.e. they contain virtual layers of fluid in relative motion with respect to one another (sliding over) at a velocity gradient in an orderly fashion. However, the inertia of the molecular constituents of the fluid may dominate over the laminar nature of flow at sufficiently high shear rates. The laminar structure of the flow is no longer sustained and instead is replaced by chaotic motion characterised by swirls which are more energy intensive compared to the laminar flow. This generally results in overestimation of viscosity.
The onset of such turbulent flows depends on the relative importance of inertial forces to viscous forces, often quantified as the Reynolds number: For rheometry using parallel plates, Re value of 100 is generally regarded as the limit beyond which the turbulent flow can exist (Davies and Stokes 2008). It is necessary to consider the possibility of turbulent flow in the present study as high shear rates will be approached and presence of such flow deteriorates the data quality which leads to misdirected conclusions.

Radial migration of fluid by centrifugal force
During high angular velocities encountered in high shear rate measurements using parallel plates, the fluid experiences considerable centrifugal force due to conservation of momentum. This force grows with the angular velocity as well as the measurement gap and it tends to eject the fluid away from the measurement geometry. Nevertheless, ejection of fluid generates more surface area and hence is resisted by the surface tension of the fluid. For Newtonian fluids, the critical shear rate for the ejection of the fluid by centrifugal forces is defined by following expression: Here C is the surface tension of the sample. This equation is only valid for Newtonian fluids as viscoelastic stresses present in non-Newtonian fluids tend to prevent the fluid migration well above the critical shear rate as defined by Eq. (13) (Davies and Stokes 2008;Pipe et al. 2008).

Materials
For CMF/CNF production, never-dried hardwood bleached sulphite cellulose pulp of commercial grade was used for the pre-treatment. Sodium hydroxide (NaOH), ethanol, hydrogen chloride (HCl), sodium chloride (NaCl) and 2,2,6,6-tetramethylpiperidine-1oxyl (TEMPO) were all of analytical grade and were used as received from Sigma Aldrich (UK). Ultrapure water (Purelab Option-Q, Class 1 water system, ELGA) was used throughout the experiments. 1000 cSt polydimethylsiloxane (silicone oil) from Sigma Aldrich (UK) was used as received for gap error determination.
CMF/CNF production TEMPO-oxidation of the cellulose was carried out following the general procedure described by Saito and Isogai (2004) (Saito and Isogai 2004). The cellulose pulp was dispersed in ultrapure water at 1% (w/w) concentration and transferred to a glass reactor. Appropriate amounts of TEMPO and NaBr were added and stirred until completely dissolved. The reaction was started by adding NaClO (5000 lmol g -1 ) in dropwise manner. The reaction pH was maintained at 10 ± 0.2 by adding 0.5 M NaOH. After 1 h, the reaction was quenched with 20 ml of ethanol. 1% (w/w) aqueous dispersion was prepared from the insoluble cellulose obtained after oxidation and passed once through the Microfluidiser (Microfluidics Inc.) at 172 MPa pressure using auxiliary chamber with 200 lm diameter and interaction chamber with 100 lm diameter. The total surface charge of the TEMPO-oxidised cellulose were determined by conductometric titration according to the method used by Saito and Isogai (2004) to be 1063 lmol of carboxyl groups per gram of cellulose. The degree of oxidation DO was determined to be 0.18, using the equation: where C is the concentration of base (mol/L) and V 1 and V 2 are the amount of base (L), w (g) is the weight of the oven-dried sample. The values of 162 and 36 corresponds to the molecular weight (g/mol) of anydroglucose repeating unit as well as the difference between the molecular weight of an anhydroglucose unit and that of the sodium salt of a glucuronic acid group, respectively (Habibi et al. 2006).

Optical microscopy
The sample was diluted using ultrapure water and homogenised using a sonication probe (GEX 130, ultrasonic processor, 130 W, Cole-Parmer, UK) for 180 s at 60% intensity, forming a 0.02% (w/w) solution. This mild ultrasonication treatment allowed dispersion of the gel-like material in water at low concentration and was not believed to have significantly altered the fibrillar structure. The material was stained using small amount of methylene blue solution in order to increase the contrast between the background and cellulose material. The stained sample was observed using an optical microscope (Leica DMRX) in transmission mode.

Scanning electron microscopy
For SEM observation, the sonicated sample was further diluted forming the final solution of 0.001% (w/w). The solution was dropped on a fresh cleaved mica disc (muscovite, 9.9 mm diameter and 0.22-0.27 mm thickness, Agar Scientific, UK) that was attached on a SEM aluminium stub. Subsequently, the sample was dried in a vacuum oven (Gallenkamp) at 35°C and 700 mbar in vacuum overnight. The dried sample was gold coated for 90 s using a sputter coater (EMITECH K550X, Quorumtech, UK) in order to provide adequate conductivity prior to image capture using field emission scanning electron microscope (FE-SEM, S-4800 Hitachi, Japan). 3 kV and 8.5 mm were used as the acceleration voltage and observation distance respectively in order to prevent damage on samples while being observed.

Rheometry
In order to determine and account for the gap error occurring due to non-parallelism and non-flatness of plates and the gap-zeroing procedure, the viscosity of 1000 cSt polydimethylsiloxane was tested at different measurement gaps (10-600 lm) using the 60 mm parallel plate, intended to be used for testing the TEMPO-oxidised CMF/CNF suspension sample on AR-G2 controlled stress rheometer by TA Instruments. The shear rate range was 1-10 s -1 with total of 7 data points logarithmically equidistant to one another. Following each measurement, the geometry was removed, cleaned with acetone and reattached. Following reattachment, the gap was 'zeroed' at 5 N normal force. For reference, the sample was also tested using a 60 mm 1°cone geometry.
1% (w/w) aqueous suspension of TEMPO oxidised CMF/CNF was tested using the 60 mm parallel plate with solvent trap at seven different commanded gaps ranging from 150 to 40 lm, where shear stress was applied logarithmically over 350-1 Pa with 10 data points per decade. The steady state was considered to have been achieved when the average measurements from three consecutive periods of 20 s were within 5% range of one another. The maximum time allowed for measurement per point was 240 s. For comparison, the sample was also tested using serrated concentric cylinders geometry with 1 mm gap (12 mm bob radius, 60 mm bob length and 13 mm cup radius), in which 100-0.1 s -1 of shear rate was applied logarithmically with 4 data points per decade. The criteria for steady state were the same as those with the experiments with 60 mm parallel plate.
All rheological tests were performed in duplicates to ensure reproducibility and at 20°C temperature.

Results and discussion
Optical microscopy Figure 1 shows the image of processed sample that has been stained using methylene blue, being observed under the optical microscope. It is clear that CMFs are present, with diameter of several tens of microns. Blue clouds are obvious around the CMFs, which could be attributed to the presence of fibrils of small diameter including partially delaminated and completely detached nanofibrils. However, these structures are too small to be identified individually using optical microscope.

Scanning electron microscopy
More detailed microstructure of the sample can be seen from the SEM micrographs shown in Fig. 2a, b. It is clearly manifested that the surface of fibre fragment is not smooth and fibril delamination occurred around it, which is believed as the fibrillation from outer layer of cellulose cell wall structure. This is more pronounced at the two ends of the fragment. These fibrils are tightly attached to one another like a 'mat', indicating the incomplete delamination/fibrillation process, which are commonly seen in samples of lower number of passes using a Microfluidiser (early stage of fibrillation and delamination). At higher number of passes, fibre fragments are expected to reduce their size and eventually disappear while more individual fibrils get released.
It is worth noting that individual fibrils can indeed be seen in the image of higher magnification Fig. 2c, although the amount is very small at this stage of the process. Due to strong hydrogen bonding, individual fibrils tend to attach to one other forming a complex, highly entangled, web-like network structure which results in the 'clouds' seen in optical microscope. It is clear that the fibrils are of approximately 20 nm in diameter and the length appears in the order of several hundred nm. However, it is not possible to determine the accurate length due to the great aspect ratio of nanofibrils which are overlapping each other significantly. From both optical and scanning electron microscopic images, the largest dimension of the large fragments in the sample appear to be between 40 and 50 lm.
It is evident that the processed sample contains CMFs and CNFs, although the amount of latter is insignificant when comparing with that of the former. A separate work has shown that as the number of passes through the Microfluidiser increase, the amount of CMFs is reduced till it vanishes completely at 5 passes (Onyianta et al. 2017). Therefore, from the observation of the optical and electron microscopic images, a better description of such material would be CMF/CNF.

Determination of gap error
The viscosity of 1000 cSt polydimethylsiloxane was independent of applied shear rate at each measurement gap (data not shown). Newtonian model was fitted to the data to obtain viscosity values, which appeared to decrease with decreasing gap. This observation was consistent with earlier findings (Davies and Stokes 2008;Pipe et al. 2008;Ewoldt et al. 2015). As suggested in the Theory section, this decrease was believed to be a consequence of the systematic gap error. The data was plotted according to Eq. (8) as shown in Fig. 3.
By fitting the data in Fig. 3 with a straight line, the gap error was calculated to be 19.1 lm and the viscosity (which is theoretically independent of the gap) was calculated to be 1.12 Pa s. The viscosity measured using 60 mm 1°cone was 1.07 Pa s.

Rheometry of TEMPO oxidised CNF suspension
No ejection of the sample at high shear stresses was observed. It was not possible to achieve commanded measurement gaps lower than 40 lm owing to the normal stress transducer getting saturated frequently while lowering the gap.
The steady shear stress sweep data for 1% (w/w) TEMPO-oxidised CMF/CNF suspension as measured have been plotted in Fig. 4 for all the gaps. It is apparent from Fig. 4 that the measurement data are in more agreement at higher shear rates/stresses compared to low shear rates/stresses whilst being shear thinning owing to deformation/break-down of the network in response to shear. Also apparent is the lack of agreement between the data measured using concentric cylinders and that measured using parallel plates. The difference however, decreases with increasing shear stress/shear rate and is particularly noticeable in the shear stress versus shear rate and viscosity vs shear stress graphs compared to the viscosity vs shear rate graph.  The reason for providing shear stress versus shear rate and viscosity versus shear stress plots in addition to the viscosity versus shear rate plot is because the latter graph often conceals the data artefacts arising due to slip (or wall depletion), shear banding and shear fracture. For example, the shear stress can either be seen to decrease, remain nearly constant, increase moderately or increase strongly with increasing shear rate. All these possibilities except for the last are seen in the viscosity versus shear rate plot as shearthinning, whereas only the last one will be expressed as shear-thickening owing to viscosity being the ratio of shear stress and shear rate (Barnes 1995). The suspensions of CMF and CNF, among vast majority of other commercially relevant fluids are predominantly shear thinning over a broad range of shear rates. Hence observing a viscosity versus shear rate plot which appears to shear thin (while possibly concealing other artefacts due to abovementioned reasons) is not seen with due scepticism.
Nevertheless, such lack of scepticism has staggering consequences for the scientific validity of the data. Since tangential viscosity at any given shear rate or shear stress is obtained by the slope of the tangent on shear stress versus shear rate curve, shear stress decreasing with increasing shear rate will imply negative viscosity. This would either violate the first law of thermodynamics (by assuming that more work of shear is done with lower energy input) or necessitates the existence of a high-shear rate and a low-shear rate region (Shaw 2012). Another implication is that the deformation in one direction is caused by the forces acting in the opposite direction (Malkin and Isayev 2017). Even the shear stress being constant as a Fig. 4 a Viscosity versus shear rate, b shear stress versus shear rate and c viscosity versus shear stress for 1% (w/w) TEMPOoxidised CMF/CNF suspension at various commanded measurement gaps prior to any correction. Also included are the data obtained using concentric cylinders function of shear rate (i.e. a shear stress plateau) implies faster flow without additional resistance being met by the fluid, essentially an excess flow without energy costs. Despite this, the scientific literature on the rheological properties of suspensions/dispersions in general and those of CMF and CNF in particular is replete with viscosity data which, if replotted as shear stress versus shear rate would reveal that the shear stress either decreases with shear rate or remains constant (Pääkko et al. 2007;Lasseuguette 2008;Karppinen et al. 2011Karppinen et al. , 2012Missoum et al. 2012;Dimic-Misic et al. 2013b;Rezayati Charani et al. 2013;Fukuzumi et al. 2014;Naderi and Lindström 2014;Naderi et al. 2014bNaderi et al. , 2015a. Such unnoticed errors do not only mislead towards unsuitable microstructural interpretations being propagated across the scientific community but also false empirical relations. For example, Naderi et al. (2014a, b) (Naderi et al. 2014b) reported viscosity as a function of shear rate for CNF suspensions at several concentrations, which revealed the abovementioned artefact upon replotting. Without realising this, the data was used to derive a power-law relationship between viscosity at 0.1 s -1 shear rate and concentration. All shear rate/ shear stress dependent rheological data in the present study has therefore been plotted in all three forms.
In Fig. 4, there is a clear trend of decreasing viscosity with decreasing gap which is similar to the trend observed in the viscosity data of silicone oil. Prior to any corrections, shear rates [ 10 5 s -1 have been attained. However, the Weissenberg-Rabinowitsch correction must be applied according to Eq. (9) for the radial dependence of shear stress and shear rate. This correction was applied both manually and using TRIOS software by TA Instruments. Identical results were obtained using both the methods which have been plotted in Fig. 5.
While the relative trends observed in the data following Weissenberg-Rabinowitsch correction (i.e. agreement between the curves at different shear rates and at different gaps) are largely the same as those found in the data prior to the correction, the entire family of curves of the data from parallel plates has shifted lower owing to the correction at the same shear rate values.
As shown in the Theory section, the shear stress can reduce up to 25% from its apparent value following the correction. For the present data, the shear stress values at highest shear rates have dropped between 5 and 7% while those at the lowest shear rates have dropped by more than 20%. This is expected as shear thinning is more profound at lower shear rates compared to that at higher shear rates as the viscosity approaches plateau, and the magnitude of correction in Eq. (9) is based on the first derivative of logarithmic torque and logarithmic angular velocity plot. The magnitude of change remained nearly independent of commanded gap.
The data from parallel plates measurement in Fig. 5 were further corrected for the gap error which was calculated from Fig. 1 (19.1 lm). In order to correct for the gap error, the shear rate values were multiplied by a correction factor, i.e. d/d ? e. Again, the results were identical when the correction was performed within TRIOS. The results have been plotted in Fig. 6.
Since the gap correction lowers the shear rate values while keeping the shear stress values the same, the maximum shear rate achieved is now under 10 5 s -1 following the correction. As the gap error is the same in magnitude for all commanded gaps, its impact in percentage was the lowest for the highest gap (11% reduction in shear rate) and the highest for the lowest gap (32% reduction in shear rate). For a given gap, the impact has to remain constant across the entire range of shear rates.
It appears as if the data at shear rates [ 10 3 s -1 at various measurement gaps has slightly better agreement compared to that prior to the gap correction, while the data at lower shear rates largely remain in disagreement. In addition, the gap correction has shifted the family of curves obtained using parallel plate geometry slightly towards the curve obtained using the concentric cylinders.
The impact of all successive transformations on the data can be seen in the representative example given in Fig. 7. For the sake of completeness, these two transformations were also performed on 'as measured data' in the opposite order, i.e. gap correction followed by Weissenberg-Rabinowitsch correction. The results obtained were the same as obtained from performing the Weissenberg-Rabinowitsch correction followed by gap correction (data not shown).
From Fig. 6, the gradual reduction in the apparent viscosity with decreasing gap indicates growing importance of wall depletion at low gaps. The dimensions of suspended particles become more significant compared to the gap as the gap decreases, thereby rendering the wall depletion increasingly dominant on the measurement data (Barnes 1995;Saarinen et al. 2009). As seen from the optical and scanning electron microscopy, the largest fragments in the sample are between 40 and 50 lm long, being of the same order as the smallest measurement gap 59.1 lm. For an accurate determination of viscosity, the data therefore needs to be resolved into components of particle-rich bulk and particle-depleted wall.
The data in Fig. 6 were analysed using Eq. (10) for slip velocity and bulk shear rate. Since the wall depletion effects get weaker with increasing gap, the observed shear rate is extrapolated to the 'infinite' gap, i.e. the reciprocal of the gap being zero in order to determine the bulk shear rate as the y-intercept. Representative curves have been shown in Fig. 8 on a semi-log plot.
It was possible to determine bulk shear rates and slip velocities using Eq. (10) from these data for shear stresses from 325 to 38.4 Pa. For shear stresses \ 38.4 Pa, the bulk shear rate values were negative and were therefore regarded as unphysical. With the bulk shear rates obtained at corresponding shear stresses, a slip corrected flow curve has been generated ranging from 0.1 to 30,000 s -1 shear rate.
Following due consideration of the errors ranging from gap setting to wall depletion, the data presented in Fig. 9 is highly shear thinning from low to high shear rates, where it appears to be approaching a plateau value. The agreement between the data obtained from the concentric cylinders and the parallel plates seems excellent. A notable exception to this is a data point at 38.4 Pa shear stress which was obtained using parallel plates.
It has been reported that at measurement gaps around 5 times that of the biggest particle size present Fig. 5 a Viscosity versus shear rate, b shear stress versus shear rate and c viscosity versus shear stress for 1% (w/w) TEMPOoxidised CMF/CNF suspension at various commanded measurement gaps following the Weissenberg-Rabinowitsch correction. Also included are the data obtained using concentric cylinders in a fluid sample, 'jamming' effects occur when the sample behaves like a solid under elastic deformation under stress rather than a fluid. Such state originates primarily from geometrical confinement of particles at lower shear rates and resolves at higher shear rates (Davies and Stokes 2008). The inconsistency of data point at 38.4 Pa shear stress and lower seems to be stemming from these 'jamming' effects mainly due to the largest particles present in the sample with length between 40 and 50 lm being sheared within a measurement gap of 59.1 lm. Even though the measurement point is not consistent with the remaining data, such measurements can be useful in specific instances where gap-dependent rheological properties of nanoparticle suspensions are to be investigated specifically at small measurement gaps owing to the nature of the process which is to be simulated using the said rheological experiment.
The data points at highest shear rates seem to approach a plateau, most possibly indicating so called infinite-shear rate viscosity commonly observed in suspensions and solutions at very high shear rates when the particles move adjacent to one another in the laminae (Barnes 1989). However, high shear rate flows are also prone to turbulence which is known to increase the apparent shear stress and mimic the infinite-shear rate viscosity plateau. The likelihood of turbulence must therefore be determined. The Reynolds number (Re) calculated for the highest gap at the highest shear rate was found to be 86, lower than the critical value 100. Hence the flow is likely to be laminar at such shear rates.
It is also remarkable that the sample did not eject from the parallel plate geometry at such high shear rates despite the presence of centrifugal force. At the highest measurement gap (commanded-150 lm and Fig. 6 a Viscosity versus shear rate, b shear stress versus shear rate and c viscosity versus shear stress for 1% (w/w) TEMPOoxidised CMF/CNF suspension at various measurement gaps following the Weissenberg-Rabinowitsch correction and the gap correction. Also included are the data obtained using concentric cylinders corrected-169.1 lm) and at the assumed density and surface tension of 1000 kg m -3 and 0.072 N m -1 (similar to those of the continuous medium water) respectively, the critical shear rate according to Eq. (13) was calculated to be 9972 s -1 . Much higher shear rates were nevertheless achieved due to viscoelastic stresses present within the fluid preventing the migration of the fluid to the edge of the plate even beyond the critical shear rate.
Although a power-law relationship between shear stress and shear rate can effectively describe the profound shear thinning observed in the data, powerlaw would fail at higher shear rates when the viscosity appears to approach the so-called 'infinite-shear rate viscosity'. Sisko model has been considered as an Fig. 7 a Viscosity versus shear rate, b shear stress versus shear rate and c viscosity versus shear stress 1% (w/w) TEMPOoxidised CMF/CNF suspension at 40 lm commanded measurement gap as measured, following the Weissenberg-Rabinowitsch correction and Weissenberg-Rabinowitsch ? Gap correction. Also included are the data obtained using concentric cylinders Fig. 8 Observed shear rate plotted against the reciprocal of the measurement gap (corrected) at several representative shear stresses for 1% (w/w) TEMPO-oxidised CMF/CNF suspension. The curved black lines indicate linear fits on a semi-log plot appropriate empirical expression for such data (Sisko 1958;Barnes et al. 1993): Equation (15) is referred to as Sisko equation, where g is the shear viscosity, g ? is asymptotical viscosity at very high shear rates, n is the power-law index and k S (despite of not having the dimensions of viscosity) is the viscosity additional to g ? at the shear rate _ c.= 1 s -1 . This equation seems to represent the data quite well when the data point at 38.4 Pa is disregarded. The fit parameters are g ? = 9.12 9 10 -3 Pa s, n = 0.176 and k S = 14.5 Pa s n . The slip velocity values calculated using Eq. (10) can be seen plotted against the corresponding shear stress values in Fig. 10. It can be represented using a Fig. 9 a Viscosity versus shear rate, b shear stress versus shear rate and c viscosity versus shear stress for 1% (w/w) TEMPOoxidised CMF/CNF suspension following correction for slip effects. Also included are the data obtained using concentric cylinders. The black line in (a) is the fit to Eq. (15) Fig. 10 The slip velocity plotted against shear stress. The black line represents power-law fit power-law fit with 1.35 exponent. It is evident how the slip velocity increases faster than shear stress and therefore how vital it becomes to take wall depletion effects into consideration at high shear rates.
Similarly, the thickness of the depleted layer calculated using Eq. (11) is also seen to increase with shear stress as shown in Fig. 11. With increased layer thickness as well as increased velocity at higher shear stresses/rates, wall depletion governs the flow kinematics at high shear rates.
Particularly at the highest shear stress, the depletion layer thickness is 10.2 lm, i.e. 17% of the smallest measurement gap 59.1 lm. At this point, large fragments present in the sample dominate the entire bulk flow layer with their largest dimension ranging between 40 and 50 lm.

Conclusions
The present work reports high shear-rate rheometry for cellulose micro/nanofibril suspension using a rotational rheometer at very small confinements. 1% (w/ w) suspension of TEMPO-oxidised CNF was tested using parallel plate geometry at commanded gaps ranging between 150 and 40 lm. Setting these gaps is prone to error owing to the plates being non-parallel and/or non-flat as well as the 'squeezing' viscosity of air affecting the zero-gap setting. The gap error was determined by measuring the viscosity of a Newtonian silicone oil at different commanded gaps. Also accounted for were the factors such as radial dependence of shear stress in parallel plates and the presence of particle-depleted layer close to the geometry walls rendering the shear gradient discontinuous across the measurement gap.
The possibility of turbulent flows and radial migration of fluid affecting measurement data was also ruled out. Viscous heating at high shear rates is less likely due to small measurement gaps enabling efficient temperature control. Consequently, a steady shear flow curve ranging from 0.1 to 30,000 s -1 of shear rate was constructed in which excellent agreement was seen between the measurements performed by two different type of geometries. The viscosity was seen approaching a so called 'infinite shear rate' viscosity plateau.
Accurate determination of viscosity and the ability to predict the onset of turbulence at high shear rates using rotational rheometers has positive consequences for operations such as coating, fibre spinning, pumping and filling containers. On the other hand, use of CNFs as a rheology modifier in several products leads to its suspensions being pumped and filled in containers which may expose it to shear rates around 10 4 s -1 . A precise knowledge of shear stresses expressed as pressures at high shear rates in such processes is therefore of commercial interest.