1. Introduction
Ice-sheet flow plays a significant role in climate change, and flow solutions require a constitutive law for the stress dependence on ice deformation and strain rate, and on temperature. To a very good approximation ice is incompressible, and the pressure is a workless constraint. The response to an applied stress first shows a primary creep, over hours or less, in which the final deformation is very small, so on the large timescales, years, of ice-sheet flow, this is neglected. The subsequent stationary creep is described by a viscous law for a fluid of grade 1 (Truesdell and Noll, Reference Truesdell and Noll1992), for which we will use the common description viscous fluid law. This is followed by a tertiary accelerating creep depending on the deformation history, describing the evolution of induced anisotropy (Budd and Jacka, Reference Budd and Jacka1989; Faria and others, Reference Faria, Weikusat and Azuma2014). Ice-sheet modelling has commonly adopted the viscous fluid law for the complete response, neglecting tertiary creep. Non-simple fluid or solid laws to describe the tertiary creep will still incorporate the viscous fluid response as an initial condition. Here, we focus on that viscous fluid response. We also assume the ice is thermorheologically simple, for which the strain rate at a given stress can be normalised by a temperature-dependent rate factor (Morland and Lee, Reference Morland and Lee1960). Detailed discussion of the theoretical mechanics underlying the various constitutive theories for ice deformation and flow, and their role in the flow of glaciers and ice sheets, are presented in the recently published book ‘Fundamental Glaciology’ by Hutter (Reference Hutter2020), a substantial update of Hutter's earlier book ‘Theoretical Glaciology’ (Hutter, Reference Hutter1983).
A viscous fluid law, necessarily isotropic by material frame indifference, has the general Rivlin–Ericksen quadratic representation, with two response coefficient functions depending on two principal invariants of the strain rate, in a stress formulation, or on two principal invariants of the deviatoric stress in a strain-rate formulation, as discussed by Morland (Reference Morland1979). Glen (Reference Glen1958) (acknowledging Fritz Ursell) presented this general quadratic viscous relation for the strain rate, but adopted the simple coaxial form proposed by Nye (Reference Nye1953) with dependence on one, the second principal, deviatoric stress invariant, which is a measure of the shear stress magnitude. His pioneering uni-axial stress experiments at Cambridge (Glen, Reference Glen1952, Reference Glen1953, Reference Glen1955, Reference Glen1958) provided data to determine the one response coefficient function dependence on the one, the second principal, deviatoric stress invariant, for which he assumed a power law. This coaxial law, ‘Glen's law’, has been used extensively. Subsequent correlations with the same data assuming a polynomial law (Smith and Morland, Reference Smith and Morland1981) gave closer correlations, and avoid the infinite viscosity at zero stress arising with the power law. Morland and Staroszczyk (Reference Morland and Staroszczyk2020) determined an accurate fractional power expansion for the equivalent coaxial stress formulation.
Simple shear response data would also determine a coaxial law with one response coefficient function depending on one invariant. For both uni-axial and shear responses to determine the same law requires a relation between these responses: the criterion for a unique coaxial law with one response coefficient function depending on one invariant. This criterion is determined for both strain-rate and stress formulations. Glen noted that Steinemann (Reference Steinemann1954) claimed that his compression and shear data were not consistent with a common coaxial form with dependence on only one invariant, but suggested retaining the coaxial form with additional dependence on the third principal deviatoric stress invariant. Here, the coaxial strain-rate formulation and equivalent stress formulation with dependence on one invariant are analysed to determine required criteria on the uni-axial and shear responses. Note although that these two responses are not sufficient to determine the response coefficient dependence on two invariants. The first principal strain-rate and deviatoric stress invariants are zero by incompressibility.
Given the lack of simple shear data, we now extend the analysis and criteria to Reference SteinemannSteinemann's (Reference Steinemann1958) uni-axial compression and torsion data. The torsion of the hollow cylinder is dominated by shear, but this is non-uniform through the cylinder, depending on height above a fixed base and on radial distance from the axis. The applied torque then involves an integral of a response coefficient function between the inner and outer radii, but an explicit correlation is deduced for the stress formulation – that needed in the momentum balance. It shows that the conventional coaxial relation with one response coefficient function depending on one invariant does not apply at all stress levels. Again, the two responses cannot determine a coaxial response coefficient function with dependence on two invariants, but a necessary criterion is derived. Steinemann's data do not rule out this form. We correlate the two responses with a quadratic relation involving two response coefficient functions, both depending on only one, the second principal strain-rate invariant, which shows that a non-trivial quadratic term arises at both high and low stress level ranges; that is, the coaxial form with dependence on one invariant is not valid at all stress levels, including shear stress levels arising in ice-sheet flows which is where the viscous law is required. Note that the ice-sheet stress levels are much smaller than the experimental data stress levels, so this response must be inferred by extrapolating the experiments correlated relation to these lower stress levels. This background of Steinemann's largely ignored work is presented in the Preface of the book by Hutter (Reference Hutter2020).
2. The simple viscous fluid relation
We adopt normalised dimensionless deviatoric stress $\hat {{\boldmath {\sigma }}}$ and strain rate D, the symmetric part of the velocity gradient tensor, based on a stress unit σ 0 = $10^{5}{\, \hbox {Pa}}$ and strain-rate unit D 0 = $1{\, \hbox {a}}^{-1}$ = ${3.17}\times 10^{-8}{\, \hbox {s}}^{-1}$, where ‘a’ denotes year, which are the deviatoric stress and strain-rate magnitudes expected, at melt point, in ice-sheet flow. The units σ 0 and D 0 are those used by Morland and Johnson (Reference Morland and Johnson1980), Morland (Reference Morland1984) and Smith and Morland (Reference Smith and Morland1981) to obtain a normalised dimensionless viscous relation at the melt temperature in the conventional coaxial form.
Now let σ denote the Cauchy stress tensor, then
where p is the mean pressure and I is the unit tensor. The principal invariants of $\hat{{{\boldsymbol \sigma }}}$ are defined by
where for convenience J 2 has the opposite sign to the usual definition. Ice response to stress exhibits a strong dependence on temperature T, which is assumed to be described by applying a rate factor a(T) to the strain rate, where a(T) is a rapidly increasing function of T; that is, the actual strain rate at a given stress and temperature increases rapidly with temperature. As noted by Morland (Reference Morland1979), this is the assumption of a thermorheologically simple response (Schwarzl and Staverman, Reference Schwarzl and Staverman1952; Morland and Lee, Reference Morland and Lee1960), in which the same processes occur, but on a timescale factored by a(T). This is satisfied by introducing a temperature-normalised strain rate ${\bar {{\boldsymbol D}}}$ with principal invariants $\bar {I}_1$, $\bar {I}_2$, $\bar {I}_3$ defined by
That is, for a given ${\bar {{\boldsymbol D}}}$ the actual strain rate D is magnified by the temperature-dependent rate factor a(T). The vanishing of I 1 is the incompressibility condition, and $\bar {I}_2$ is a measure of the normalised strain-rate magnitude squared. Note that this definition of I 2, used for convenience, has the opposite sign to the usual definition of the second principal invariant.
Smith and Morland (Reference Smith and Morland1981) constructed exponential representations for the rate factor a(T) over different temperature ranges from the constant uni-axial stress data obtained by Mellor and Testa (Reference Mellor and Testa1969), and that with the widest validity is
where T 0 is the melting point, T s is a temperature scale and a (T 0) = 1.068 is approximately unity, normalising the factor at the melt point. (We appreciate the advice given to us by Dr David Cole who has examined many datasets for temperature variation and judges the Mellor and Testa (Reference Mellor and Testa1969) data to show a consistent variation with temperature, although with higher strain rates. The above a(T) variation is therefore consistent, here normalised at the melt temperature T 0.) At $2{\, \hbox {K}}$ below melting a (271.15) = 0.4751, less than half that at the melt point, and at $30{\, \hbox {K}}$ below the melt point, a temperature magnitude found in cold ice sheets, a (243.15) = 0.0041, implying much smaller strain rates than those near melting.
Now the most general thermorheologically simple frame-indifferent incompressible viscous law of grade 1 is a relation between $\hat {{{\boldsymbol \sigma }}}$ and ${\bar {{\boldsymbol D}}}$ which can be expressed in two alternative, but physically equivalent, forms of the Rivlin–Ericksen representation between tensors with zero trace:
The vanishing of J 1 is by the definition of the deviator, so the response coefficient functions ϕ 1, ϕ 2, ψ 1, ψ 2 each depend on only two non-trivial invariants. $\bar {I}_2$ and J 2 are measures of shear strain-rate and shear stress magnitudes squared, but $\bar {I}_3$ and J 3 have no physical description. Although expansions (7) and (8) are physically equivalent, there is no explicit algebraic inversion. The coaxial form requires that both ϕ 2 and ψ 2 are identically zero. Furthermore, note that the strain rate ${{\boldsymbol D}} = a( T) {\bar {{\boldsymbol D}}}$ is not recovered by simply applying the factor a(T) to $\hat {{{\boldsymbol \sigma }}}$, except in the coaxial form with ϕ 2 = 0, and then only for the strain-rate formulation (8) since the arguments $\bar {I}_{2}$ and $\bar {I}_{3}$ of ϕ 1 depend on ${\bar {{\boldsymbol D}}}$.
These expansions for a simple viscous fluid are necessarily isotropic in all reference configurations, and cannot describe induced anisotropy associated with the fabric developed as the ice elements deform and crystal glide planes are re-oriented. However, this viscous response describes the initial isotropic response as ice is first formed at an ice-sheet surface, and is a crucial part of the constitutive behaviour as the ice deforms and induced anisotropy develops. It is the stress formulation (7) which is the constitutive law required for substitution in the momentum and energy balances of a general ice-sheet flow. The equivalent strain-rate formulation (8) is that required in the reduced model, Morland and Johnson (Reference Morland and Johnson1980), Morland (Reference Morland1984), equivalent to the shallow ice approximation (SIA) formulated by Hutter (Reference Hutter1983). The reduced model, SIA, requires longitudinal gradients to be very small, so is applicable only when the bed undulations are very small, which is not a realistic situation. However, the reduced model, SIA, is useful for determining solutions of idealised flows against which solutions from the more complex numerical schemes needed for the full equations can be compared to validate numerical models for polar ice-sheet flows.
3. Uni-axial stress and simple shear deformation relations
The general relations (7) and (8) involve two response coefficient functions depending on two invariants, which could only be determined by experiments yielding four independent data relations, and those must cover a wide invariants domain arising in ice-sheet flow. Such data do not exist.
Morland and Staroszczyk (Reference Morland and Staroszczyk2019) applied the simplified stress formulations
to uni-axial stress and simple shear deformation responses to show the correlations with the conventional coaxial response coefficients $\phi _{c1}( \bar {I}_2)$, Reference SteinemannSteinemann's (Reference Steinemann1954) coaxial conjecture $\phi _{s1}( \bar {I}_2,\; \bar {I}_3)$, and the quadratic pair $\phi _{q1}( \bar {I}_{2})$ and $\phi _{q2}( \bar {I}_{2})$. Each response determines one response coefficient function of one argument, so will in general determine two distinct $\phi _{c1}( \bar {I}_2)$ in the conventional coaxial form (9)1. We will show the criterion on the uni-axial and shear responses which determines a unique coaxial form (9)1, and when this is not satisfied determine what can be deduced for the coaxial form (9)2 with dependence on two invariants. Then we show that the two responses determine the two coefficients $\phi _{q1}( \bar {I}_{2})$ and $\phi _{q2}( \bar {I}_{2})$ of one argument in the quadratic relation (9)3. Morland (Reference Morland2007) analysed the reduced model, SIA scaling with all forms to show that a leading order (in surface slope or dimensionless viscosity magnitude) simplification still follows.
Given a dimensionless uni-axial experimental response $\sigma = U( \dot {\bar {\epsilon }})$, where σ is the axial compressive stress and $\dot {\bar {\epsilon }}$ is the axial compressive strain rate, in units σ 0 and D 0 respectively, the stress and deviatoric stress tensors are
These give rise to the axial vertical compressive strain rate $\dot {\epsilon }$ and equal horizontal strain rates ${{1\over 2}}\dot {\epsilon }$ for which the strain-rate tensor, strain-rate squared tensor and invariants are
so $\bar {I}_3$ and $\bar {I}_2$ are not independent. Thus
for the three respective relations (9).
Similarly, given a simple shear response $\tau = S( \dot {\bar {\gamma }})$, where τ is the shear stress and $\dot {\bar {\gamma }}$ is the simple shear strain rate, the strain-rate tensor and invariants are
so no dependence on $\bar {I}_3$ could be inferred. The corresponding stress is
Thus
for the respective relations (9), with no dependence on $\phi _{q2}( \bar {I}_{2})$ in the shear relation $\tau = S( \dot {\bar {\gamma }})$.
Hence, eliminating $\phi _{c1}( \bar {I}_2)$ between (13)1 and (17)1 implies
which is the criterion on the uni-axial response $U( \dot {\bar {\epsilon }})$ and shear response $S( \dot {\bar {\gamma }})$ necessary for the first coaxial relation (9)1 to hold. Responses not satisfying (18) imply that the coaxial relation (9)1 cannot be valid. For the coaxial relation (9)2, only the line $\bar {I}_3 = 0$ arises, so there are no relations over a finite domain $( \bar {I}_2,\; \bar {I}_3)$ with $\bar {I}_3 \neq 0$.
Given the two responses $\sigma = U( \dot {\bar {\epsilon }})$ and $\tau = S( \dot {\bar {\gamma }})$ not satisfying (18), we will now focus on the quadratic relation (9)3, for which (13)3 and (17)3 apply. Immediately from (17)3,
relating $\phi _{q1}( \bar {I}_2)$ directly to the shear response, independent of the uni-axial response. Then from (13)3,
which is zero if (18) applies as required. Thus, the quadratic relation (9)3 would be determined by given uni-axial and shear responses.
Analogous analyses apply to the strain-rate formulations
where the principal deviatoric stress invariants J 2 and J 3 are defined in (2). For the quadratic form (21)3, the shear relation determines
then the uni-axial relation determines
where S −1 and U −1 are the inverse functions of S and U; that is, x = U(X), X = U −1(x), y = S(Y), Y = S −1(y). The criterion for a coaxial form (21)1 is then
with the limit as J 1/2 → 0 also necessary for the form (21)2. Criteria (18) and (24) are identical relations between the responses $S( \bar {I}^{{{{1}/{2}}}})$ and $U( 2[ {\bar {I}_2}/{3}] ^{{{{1}/{2}}}})$.
Given a lack of uni-axial and shear data for the same ice specimen, we now correlate Reference SteinemannSteinemann's (Reference Steinemann1958) experimental data for uni-axial compression and torsion of a hollow cylinder with the quadratic law (9)3, removing the subscript q for convenience. The former is analysed above, and leads to relation (13)3 for ϕ 1 and ϕ 2. The simple shear analysis must now be replaced by a torsion analysis to obtain a second relation for ϕ 1 and ϕ 2. Note that there is no analogous correlation with the strain-rate formulation to determine ψ 1 and ψ 2.
4. Torsion experiment
Let (r, θ, z) be dimensional cylindrical polar coordinates, with (R, Θ, Z) the initial particle positions, in a vertical hollow cylinder of height H, with internal and external radii R i and R e. The torsion experiment determines the shear strain rate on the upper surface of a hollow cylinder due to an applied torque on the upper surface, with the bottom surface fixed, and the upper surface held at a fixed height, as illustrated in Figure 1. Figure 1a defines the notations and coordinates adopted in the analysis, and Figure 1b shows the physical dimensions of the ice samples used by Steinemann (Reference Steinemann1958) in his laboratory measurements; namely
It is assumed that the rotation is linear in elevation z, so the deformation is defined by
where t denotes time, with no rotation on the base Z = 0 and a dimensionless twist (rotation angle) at the upper surface κ(t). The surface twist rate is $\dot {\kappa }( t)$, with dimension reciprocal of the time dimension, and the corresponding dimensionless surface twist rate $\dot {\bar {\kappa }}( t)$ with unit D 0 is
The velocity is
with one non-zero velocity gradient
so that $r\dot {\bar {\kappa }}( t) /H$ is the dimensionless velocity gradient with unit D 0. Thus, the dimensionless strain rate and invariants are given by
where the subscript 2 is now omitted from the single invariant $\bar {I}_2$. The experiments relate the stress, through the torque, to strain rate at a sequence of strain rates $\dot {\bar {\kappa }}$, and hence to $\bar {I}$; there is no explicit time dependence.
The non-zero stress components are the applied shear stress σ zθ and constraints σ rr, σ θθ and σ zz. Applying the viscous law (9)3 shows immediately that σ θθ = σ zz, so the stress components and invariant are
The diagonal components of the viscous law (9)3 give only one independent relation
so σ zz and σ rr will be determined by the radial and vertical momentum (equilibrium) balances and boundary conditions. Finally, apply the stress formulation (9)3, omitting the subscript q, to obtain the shear relation
which is independent of ϕ 2, analogous to relation (19). Note that σ zθ is independent of z due to the assumption of a rotation linear in elevation.
The applied torque on the surface z = H at a given $\dot {\kappa }$ is
where σ zθ is expressed as a function of $\bar {I}$ by (34)1, and $\bar {I}_{\rm i}$ and $\bar {I}_{\rm e}$ are the strain-rate invariant limits corresponding to R i and R e. Note that M given by (35) has dimension force times length as required, that is, stress times length cubed, but has dependence on $\dot {\kappa }$ through the integral limits $\bar {I}_{\rm i}$ and $\bar {I}_{\rm e}$, as well as the proportionality to $\dot {\kappa }^{-3}$. Note also that M is independent of z: each surface $z = \hbox {constant}$ feels the same torque. For the later correlations relation (35) is expressed in dimensionless form
where
which is an integral equation to determine $\phi {_1}( \bar {I})$ from the given torque $\bar {M}( \dot {\bar {\kappa }})$. This has no algebraic solution, unlike the explicit relation (17)3 from the simple shear test. The correlation is completed by determining $\phi _{2}( \bar {I})$ from the uni-axial compression relation (13)3, thus
where ϕ 1(I) is determined by the torsion response (34).
We assume that the uni-axial and torsion viscosities at zero strain rate are finite and strictly positive, so
in the limit as $\dot {\bar {\epsilon }} \rightarrow 0$ and $\dot {\bar {\kappa }} \rightarrow 0$ respectively. It then follows from (36) and (37), with the definition (31)1, that
which, given m 1 by data correlation and the dimensions (25), determines
Then (38) shows that as $\bar {I} \rightarrow 0$, $\phi _2 \sim \bar {I}^{-{{{1}/{2}}}}$, and hence a sensible bounded representation of the quadratic response coefficient is
where, by (39)1, the latter term is bounded in the limit of zero $\bar {I}$, determining the limit relation
Here, (42) is an algebraic relation for ${\it \Phi}_{2}( \bar {I})$, depending on the correlated uni-axial response $\sigma = U( \dot {\bar {\epsilon }})$.
Note that the strain-rate formulation (21)3 does not allow the extraction of σ zθ and construction of M. Furthermore, an alternative torsion configuration with zero constraining stresses, but then needing vertical, radial and azimuthal strain rates, yields more direct relations, although still not allowing direct correlation. Essentially, a strain-rate formulation cannot be correlated with a torsion integral.
5. The quadratic viscous fluid law
We will now apply relations (36) and (42) to determine the response coefficient functions $\phi _{1}( \bar {I})$ and $\phi _{2}( \bar {I})$ in the quadratic viscous fluid stress formulation (9)3 by correlation with data points $\bar {M}( \dot {\bar {\kappa }})$ and $U( \dot {\epsilon })$. Steinemann's experiments give six torsion data points and 16 uni-axial stress data points, shown in both physical and dimensionless variables in Tables 1 and 2, respectively. Note that the torsion experiment applied the torque on a fixed upper plate consistent with our formulation. The experiments were carried out at a temperature $T = -1.9^{\circ }\hbox {C}$, for which, by (5) and (6),
The conversion of strain rates in units s−1 to the required unit ${\, \hbox {a}}^{-1}$ has a factor 3.15 × 107, but a further factor 1/a(T) = 2.04 is required to convert to the normalised strain rate D defined in (3)1; that is a total factor 6.4 × 107.
M and $\dot {\kappa }/H$ are measured torques and twist angle rates per unit height in physical units, $\bar {M}$ and $\dot {\bar {\kappa }}$ are corresponding dimensionless quantities and ${{\bar {D}}}_{z\theta }( R_{\rm e})$ and ${{\bar {I}}}_2( R_{\rm e})$ are dimensionless shear strain rates and second principal invariants at the cylinder external radius R e.
σ z and $\dot {\epsilon }_z$ are measured axial stress and axial strain rates in physical units, σ and $\dot {\bar {\epsilon }}$ are corresponding dimensionless quantities and ${{\bar {I}}}_2$ are dimensionless second principal invariants.
The data in the tables show that $U( \dot {\bar {\epsilon }})$ and $\bar {M}( \dot {\bar {\kappa }})$ have similar shapes: zero at zero strain rate, then increasing with increasing strain rate with decreasing gradient. In turn, analogous to the inversion of Glen's uni-axial law (Morland and Staroszczyk, Reference Morland and Staroszczyk2019) when $\phi _{2}( \bar {I})$ is zero, it is supposed that $\phi _{1}( \bar {I})$ is positive with strictly positive ϕ 1(0), and decreases with increasing strain rate with decreasing negative gradient. The data correlation supports this property, although here $\phi _{2}( \bar {I})$ is non-zero. Thus
However, it is $\bar {I}^{1/2}$ which defines the strain-rate magnitude, so the correlation with data is improved by expressing ϕ 1 in terms of $\bar {I}^{{{{1}/{2}}}}$. Each term of the following adopted expansions satisfy the conditions (45):
where ϕ 1(0) is related to m 1 by (41). $U( \dot {\bar {\epsilon }})$, $\bar {M}( \dot {\bar {\kappa }})$ and $\phi _1( \bar {I})$ have respectively 3K, 3L and 3N free parameters to be determined by the correlations.
Note that inserting the representation (48) for $\phi _1( \bar {I})$ in the integral (36) allows integration by parts for each term to determine an explicit algebraic representation for $\bar {M}( \dot {\bar {\kappa }})$, which can be correlated with the data correlation (47). The analytic integration requires three successive integrations by parts with a corresponding lengthy result. In view of (18) and (37), adopt
then the integration yields
Various correlation strategies and expansions have been explored, with the most accurate, adopted, strategy as follows. Representations of the functions $U( \dot {\bar {\epsilon }})$ and $\bar {M}( \dot {\bar {\kappa }})$ are determined by least squares correlations with the 16 and 6 data points respectively, shown in Tables 1 and 2. The determined continuous $\bar {M}( \dot {\bar {\kappa }})$ is then used to generate a set of 25 points regularly spread over the curve $\bar {M}( \dot {\bar {\kappa }})$, which are then correlated by least squares with the integral representation (36) to determine the representation (48) for $\phi _1( \bar {I})$. This requires numerical integration over the range $\bar {I}_{\rm i}$ – $\bar {I}_{\rm e}$ at each trial set of coefficients in (48). Then the response coefficient function ϕ 2 is determined from the algebraic relation (42) involving the known functions $U( \dot {\bar {\epsilon }})$ and $\phi _1( \bar {I})$. Alternatively, the correlation for $\phi _1( \bar {I})$ can use the analytic representation (50) for $\bar {M}( \dot {\bar {\kappa }})$. Here the former method was adopted, and then the resulting $\phi _1( \bar {I})$ was substituted in (50) to determine the consequent $\bar {M}( \dot {\bar {\kappa }})$, which, if the correlated $\phi _1( \bar {I})$ was accurate, will match the data correlated $\bar {M}( \dot {\bar {\kappa }})$.
Figure 2 shows the continuous function $U( \dot {\bar {\epsilon }})$ obtained by correlating (46) with Steinemann's 16 uni-axial compression data points shown by squares in the figure. The chosen range $0 \leq \dot {\bar {\epsilon }} \leq 200$ extends a little beyond the final data point $\dot {\bar {\epsilon }} = 164$ (see Table 2). Good agreement between the theoretical curve and measured data has been achieved by using only two terms in expansion (46), with the coefficients given by
Figure 3 shows the continuous function $\bar {M}( \dot {\bar {\kappa }})$ over the range $0 \leq \dot {\bar {\kappa }} \leq 800$, extended a little beyond the last data point $\dot {\bar {\kappa }} = 750$ (see Table 1), represented as the solid line in the plots, calculated by correlating (47) with 6 torsion data points shown by squares in the figure. Again, good accuracy has been achieved by using only two terms in expansion (47), with the coefficients
The continuous function $\bar {M}( \dot {\bar {\kappa }})$ defined by (47) with the coefficients (52) allows the choice of a closer spaced set of point for correlation purposes, and such a set of 25 correlation points are shown as solid circles in Figure 3. These 25 points are used to determine by correlation expansion (48) for the response coefficient function $\phi _1( \bar {I})$ which arises in the integral in relation (36) for $\bar {M}( \dot {\bar {\kappa }})$. The resulting expansion coefficients are
with ϕ 1(0) given by m 1 in (48)2, and Φ 2(0) by (43). These coefficients of the function $\phi _1( \bar {I})$ can now be substituted in the analytic relation (50) to determine a consequent $\bar {M}( \dot {\bar {\kappa }})$ shown as the dash-dot line in the figure. The very close match of this $\bar {M}( \dot {\bar {\kappa }})$ with the continuous $\bar {M}( \dot {\bar {\kappa }})$ (the solid line) obtained by the data correlation (47) is seen in Figure 3, which demonstrates the accuracy of the correlated $\phi _1( \bar {I})$.
Given the response coefficient function $\phi _1( \bar {I})$ established by the above correlation procedure, the second, bounded, response coefficient function ${\it \Phi}_2( \bar {I})$ can be calculated from the algebraic relation (42), with its limit value Φ 2(0) = 2.536 defined by (43). This, in turn, determines the response coefficient function $\phi_2( \bar {I}) = \bar {I}^{-1/2}\, {\it \Phi}_2( \bar {I})$, unbounded at $\bar {I} = 0$. The functions ϕ 1 and ϕ 2, or alternatively Φ 2, determine the deviatoric stress $\hat {{{\boldsymbol \sigma }}}$ in terms of the strain rate D by the law (9)3, which can be rewritten in the form
Although the response function ϕ 1 is defined analytically by the series (48) with the coefficients and ϕ 1(0) given by (53), the function ϕ 2 must be calculated from the algebraic equation (42), so it cannot be expressed by a simple analytic formula.
The three response coefficient functions $\phi _1( \bar {I})$, ${\it \Phi}_2( \bar {I})$ and $\phi _2( \bar {I})$ are illustrated in Figure 4 over the strain-rate invariant ranges $0 \leq \bar {I}^{1/2} \leq 200$ (Fig. 4a) and with more details of the low strain-rate range response in $0 \leq \bar {I}^{1/2} \leq 5$ (Fig. 4b). Recall the uni-axial and torsion correlation ranges are $\bar {I}^{1/2} = 140$ and 500, respectively.
Note that inequalities (45) applied to the ϕ 2 relation (42) do not impose any sign to ϕ 2 or to its derivative. Specifically, ϕ 2 is not necessarily positive nor monotonic over the full range $\bar {I} \geq 0$, as seen in Figure 4a. Figure 5 shows again the variation of the uni-axial dimensionless stress σ with the axial strain rate $\dot {\bar {\epsilon }}$, over the ranges of $0 \leq \dot {\bar {\epsilon }} \leq 200$ (Fig. 5a), $0 \leq \dot {\bar {\epsilon }} \leq 30$ (Fig. 5b), $0 \leq \dot {\bar {\epsilon }} \leq 1$ (Fig. 5c) and $0 \leq \dot {\bar {\epsilon }} \leq 0.01$ (Fig. 5d), with the data points shown as squares. Here, in addition, are shown the uni-axial stress contributions σ (1) and σ (2) of the linear and quadratic terms in the law (9)3, namely
Note that σ (2) is negative over a low strain rate range. In magnitude it makes a significant contribution to the total stress at the higher stress data level $\sigma _z = 1.55 \times 10^{6} {\, \hbox {Pa}}$ (the highest experimental stress data point, see Table 2), corresponding to a deviatoric stress $2\sigma _z /3 = 1.03 \times 10^{6} {\, \hbox {Pa}}$, where σ (2)/σ (1) is ~0.228.
However, the maximum deviatoric stress in a large ice sheet, estimated from a steady plane flow based on the SIA, implies a deviatoric stress magnitude $10^{4}{\, \hbox {Pa}}$, corresponding to a dimensionless uni-axial stress σ = 0.15, which is much smaller than the lowest experimental stress data point σ = 1.86 (see Table 2). (This comparison of shear stress levels in the data range with those in an ice-sheet flow was prompted by Professor K. Hutter, who had queried the stress levels in which σ (2) is relevant.) It is the response at these lower stress levels arising in ice-sheet flows that are actually required, implying lower stress experiments over longer times are needed. Here, that behaviour can only be estimated from the response coefficient functions extrapolated from the higher stress data correlation. The details of this lower stress behaviour can be seen in Figure 5d. The zero stress limit behaviour for uni-axial stress is given by the vertical component (z component) of (55), which, with the calculated limits from the correlations, is
from which the quadratic to linear term ratio is
which shows a Φ 2(0) contribution $12\percnt$ of the ϕ 1(0) contribution, not very significant, but not negligible. Figure 6a shows the ratio s r = |σ (2)/σ (1)| as a function of the axial stress σ, over a wide dimensionless axial stress range $0 -16.1$, with the upper limit value being the stress at the strain rate $\dot {\bar {\epsilon }} = 200$ (see Fig. 5a). This ratio, which measures the significance of the quadratic term, varies dramatically over the displayed range, but indicates significance at both high and low stress levels. In contrast, Figure 6b shows the ratio over a low stress range, which includes an ice-sheet stress magnitude, over which it is uniform. Its value is ~− 0.12, see the limit estimation (57); that is, the quadratic term is ~$12\percnt$ of the linear term in the ice-sheet flow range, not very significant, but not negligible.
To complete the calculations from the torsion experiments, Figure 7 shows the distributions of dimensionless shear stresses τ = σ zθ on the internal (r = R i) and external (r = R e) surfaces of the hollow cylinder, and the distribution of the corresponding shear stress on the mid-surface (R i + R e)/2. These stresses were not measured in the torsion test performed by Steinemann, in which only the resultant torques M for applied twist angle rates $\dot {\bar {\kappa }}$ were recorded.
Finally, consider the coaxial form (9)2, which, analogous to (41), requires ϕ s1(0, 0) = 0.4110 m 1, whereas, with the limit relation (39)1, (13)2 shows ϕ s1(0, 0) = 2u 1/3, which imposes a necessary zero strain-rate limit relating u 1 and m 1 for the coaxial form (9)2 to be possible, namely
In practice, the correlated limit derivatives u 1 and m 1 given in (51) and (52) are not expected to be accurate, and a sensible test of (58) is the condition
Here, this ratio is 0.1412, so Steinemann's conjecture is not ruled out by his uni-axial stress and torsion experiments. However, these experimental configurations do not provide response data over a 2-D (I 2, I 3) domain which is necessary to determine a response coefficient depending on two invariants.
6. Conclusions
We have examined three different particular stress formulations of the general quadratic viscous fluid relation: coaxial with dependence on one strain-rate invariant, coaxial with dependence on two strain-rate invariants, and quadratic with dependence on one strain-rate invariant. Correlations with uni-axial stress and simple shear relations are analysed, and the criteria between these responses necessary for the coaxial relations are determined. Lacking simple shear data, an alternative analysis is made of Steinemann's torsion experiment configuration, which combined with his uni-axial stress data determines the analogous criteria for coaxial relations. Uni-axial and shear experiments do not cover a finite domain of independent $\bar {I}_2$ and $\bar {I}_3$, which is required to determine response coefficient functions depending on two invariants. However, the uni-axial and torsion experiments can be correlated with the quadratic relation with both response coefficients depending on only one invariant, $\bar {I}_2$. This correlation shows that the quadratic term is significant over the large stress and strain-rate range, implying that the conventional coaxial relation with the one response coefficient depending only on $\bar {I}_2$ would not be valid in that range. Furthermore, examining the low stress behaviour extrapolated from Steinemann's data implies that the quadratic term is not negligible relative to the linear term at the lower shear stress levels arising in ice-sheet flows. This does not reject a coaxial relation with dependence on two invariants. Experimental configurations which cover a finite domain of independent $\bar {I}_2$ and $\bar {I}_3$ are required to determine a general coaxial relation.
Acknowledgements
Professor L. W. Morland is grateful for the award of a Leverhulme Trust Emeritus Fellowship (reference number 61410) supporting his collaboration with Dr R. Staroszczyk, and we appreciate Professor K. Hutter reminding us of Reference SteinemannSteinemann's (Reference Steinemann1958) experimental research, and Dr J. Glen loaning us his Steinemann papers. Dr R. Staroszczyk is grateful to the Dean of the School of Mathematics of the UEA for the invitation to spend 1 month in Norwich in August 2019, when the work on this research was started.