A quadratic viscous fluid law for ice deduced from Steinemann's uni-axial compression and torsion experiments

Abstract

The response of ice to applied stress on ice-sheet flow timescales is commonly described by a non-linear incompressible viscous fluid, for which the deviatoric stress has a quadratic relation in the strain rate with two response coefficient functions depending on two principal strain-rate invariants I₂ and I₃. Commonly, a coaxial (linear) relation between the deviatoric stress and strain rate, with dependence on one strain-rate invariant I₂ in a stress formulation, equivalently dependence on one deviatoric stress invariant in a strain-rate formulation, is adopted. Glen's uni-axial stress experiments determined such a coaxial law for a strain-rate formulation. The criterion for both uni-axial and shear data to determine the same relation is determined. Here, we apply Steinemann's uni-axial stress and torsion data to determine the two stress response coefficients in a quadratic relation with dependence on a single invariant I₂. There is a non-negligible quadratic term for some ranges of I₂; that is, a coaxial relation with dependence on one invariant is not valid. The data does not, however, rule out a coaxial relation with dependence on two invariants.

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, 1992), 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, 1989; Faria and others, 2014). 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, 1960). 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 (2020), a substantial update of Hutter's earlier book ‘Theoretical Glaciology’ (Hutter, 1983).

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 (1979). Glen (1958) (acknowledging Fritz Ursell) presented this general quadratic viscous relation for the strain rate, but adopted the simple coaxial form proposed by Nye (1953) 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, 1952, 1953, 1955, 1958) 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, 1981) gave closer correlations, and avoid the infinite viscosity at zero stress arising with the power law. Morland and Staroszczyk (2020) 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 (1954) 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 Steinemann's (1958) 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 (2020).

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 σ ₀ = $10^{5}{\, \hbox {Pa}}$ and strain-rate unit D ₀ = $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 σ ₀ and D ₀ are those used by Morland and Johnson (1980), Morland (1984) and Smith and Morland (1981) to obtain a normalised dimensionless viscous relation at the melt temperature in the conventional coaxial form.

Now let σ denote the Cauchy stress tensor, then

(1)$$\hat{\boldsymbol \sigma} = {\boldsymbol \sigma} + p{\boldsymbol I},\; \quad p = -{1\over 3}{\rm trace}( {\boldsymbol \sigma}) ,\; \quad {\rm trace}( \hat{\boldsymbol \sigma}) = 0,\; $$

where p is the mean pressure and I is the unit tensor. The principal invariants of $\hat{{{\boldsymbol \sigma }}}$ are defined by

(2)$$J_1 = {\rm trace}( \hat{{{\boldsymbol \sigma}}}) = 0,\; \quad J_2 = {\rm trace}( {\hat{{{\boldsymbol \sigma}}}}^{2}) /2,\; \quad J_3 = {\rm det}( \hat{{{\boldsymbol \sigma}}}) ,\; $$

where for convenience J ₂ 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 (1979), this is the assumption of a thermorheologically simple response (Schwarzl and Staverman, 1952; Morland and Lee, 1960), 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

(3)$${{\boldsymbol D}} = a( T) {\bar{{\boldsymbol D}}},\; \quad I_1 = {\rm trace}( {{\boldsymbol D}} ) = a( T) \bar{I}_1 = 0,\; $$

(4)$$I_2 = {\rm trace}( {{\boldsymbol D}}^{2}) /2 = [ a( T) ] ^{2}\bar{I}_2,\; \quad I_3 = {\rm det}( {{\boldsymbol D}}) = [ a( T) ] ^{3}\bar{I}_3.$$

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 ₁ is the incompressibility condition, and $\bar {I}_2$ is a measure of the normalised strain-rate magnitude squared. Note that this definition of I ₂, used for convenience, has the opposite sign to the usual definition of the second principal invariant.

Smith and Morland (1981) 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 (1969), and that with the widest validity is

(5)$$a( T) = 0.7242\exp( 11.9567\, \bar{T}) + 0.3438\exp( 2.9494\, \bar{T}) ,\; $$

(6)$$T = T_0 + T_{\rm s}\, \bar{T},\; \quad T_0 = 273.15{\, \hbox{K}},\; \quad T_{\rm s} = 20{\, \hbox{K}},\; $$

where T ₀ is the melting point, T _s is a temperature scale and a (T ₀) = 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 (1969) 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 ₀.) 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:

(7)$$\hat{{{\boldsymbol \sigma}}} = \phi_{1}( \bar{I}_{2},\; \bar{I}_{3}) \, {\bar{{\boldsymbol D}}} + \phi_{2}( \bar{I}_{2},\; \bar{I}_{3}) \left[{{\bar{{\boldsymbol D}}}}^{2} - {2\over 3}\bar{I}_{2}{\boldsymbol I}\right],\; $$

(8)$${\bar{{\boldsymbol D}}} = \psi_{1}( J_{2},\; J_{3}) \, \hat{{{\boldsymbol \sigma}}} + \psi_{2}( J_{2},\; J_{3}) \left[{\hat{{{\boldsymbol \sigma}}}}^{2} - {2\over 3}J_{2}{\boldsymbol I}\right].$$

The vanishing of J ₁ is by the definition of the deviator, so the response coefficient functions ϕ ₁, ϕ ₂, ψ ₁, ψ ₂ each depend on only two non-trivial invariants. $\bar {I}_2$ and J ₂ are measures of shear strain-rate and shear stress magnitudes squared, but $\bar {I}_3$ and J ₃ have no physical description. Although expansions (7) and (8) are physically equivalent, there is no explicit algebraic inversion. The coaxial form requires that both ϕ ₂ and ψ ₂ 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 ϕ ₂ = 0, and then only for the strain-rate formulation (8) since the arguments $\bar {I}_{2}$ and $\bar {I}_{3}$ of ϕ ₁ 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 (1980), Morland (1984), equivalent to the shallow ice approximation (SIA) formulated by Hutter (1983). 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 (2019) applied the simplified stress formulations

(9)$$\hat{{{\boldsymbol \sigma}}} = \phi_{c1}( \bar{I}_2) \, \bar{{\boldsymbol D}},\; \quad \hat{{{\boldsymbol \sigma}}} = \phi_{s1}( \bar{I}_2,\; \bar{I}_3) \, \bar{{\boldsymbol D}},\; \quad \hat{{{\boldsymbol \sigma}}} = \phi_{q1}( \bar{I}_{2}) \, \bar{{\boldsymbol D}} + \phi_{q2}( \bar{I}_{2}) \left[{\bar{{\boldsymbol D}}}^{2} - {2\over 3}\bar{I}_{2}{\boldsymbol I}\right],\; $$

to uni-axial stress and simple shear deformation responses to show the correlations with the conventional coaxial response coefficients $\phi _{c1}( \bar {I}_2)$, Steinemann's (1954) 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)₁. We will show the criterion on the uni-axial and shear responses which determines a unique coaxial form (9)₁, and when this is not satisfied determine what can be deduced for the coaxial form (9)₂ 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)₃. Morland (2007) 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 σ ₀ and D ₀ respectively, the stress and deviatoric stress tensors are

(10)$${{\boldsymbol \sigma}} = \left(\matrix{{0} & {0} & {0}\cr {0} & {0} & {0}\cr{0} & {0} & {-\sigma}}\right),\; \quad p = {1\over 3} \sigma ,\; \quad \hat{{{\boldsymbol \sigma}}} = {1\over 3} \left(\matrix{{ \sigma} & {0} & {0}\cr{0} & {\sigma} & {0}\cr{0} & {0} & {-2\sigma }}\right).$$

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

(11)$$\bar{{\boldsymbol D}} = {1\over 2} \left(\matrix{{\dot{\bar{\epsilon}}} & {0} & {0}\cr{0} & {\dot{\bar{\epsilon}}} & {0}\cr{0} & {0} & {-2\dot{\bar{\epsilon}}}}\right),\; \quad {\bar{{\boldsymbol D}}}^{2} = {1\over 4} \left(\matrix{{\dot{\bar{\epsilon}}^{2}} & {0} & {0}\cr{0} & {\dot{\bar{\epsilon}}^{2}} & {0}\cr{0} & {0} & {4\dot{\bar{\epsilon}}^{2}}}\right),\; $$

(12)$$\bar{I}_2 = {3\over 4}\dot{\bar{\epsilon}}^{2} ,\; \quad \bar{I}_3 = -\dot{\bar{\epsilon}}^{3} = -( 4\bar{I}_2/3) ^{{3}/{2}},\; $$

so $\bar {I}_3$ and $\bar {I}_2$ are not independent. Thus

(13)$$U\left(2\left[{\bar{I}_2\over 3}\right]^{1/2}\right) = \left[3\bar{I}_2\right]^{1/2} \phi_{c1}( \bar{I}_2) \ \hbox{or}\ = \left[3\bar{I}_2\right]^{1/2}\phi_{c1}( \bar{I}_2,\; \bar{I}_3) \ \hbox{or}\ = \left[3\bar{I}_2\right]^{1/2}\phi_{q1}( \bar{I}_2) - \bar{I}_{2}\phi_{q2}( \bar{I}_2) $$

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

(14)$$\bar{{\boldsymbol D}} = \left(\matrix{{0} & {0} & {\dot{\bar{\gamma}}}\cr{0} & {0} & {0}\cr{\dot{\bar{\gamma}}} & {0} & {0}}\right),\; \quad {\bar{{\boldsymbol D}}}^{2} = \left(\matrix{{\dot{\bar{\gamma}}^{2}} & {0} & {0}\cr{0} & {0} & {0}\cr{0} & {0} & {\dot{\bar{\gamma}}^{2}}}\right),\; $$

(15)$$\bar{I}_2 = \dot{\bar{\gamma}}^{2},\; \quad \bar{I}_3 = 0,\; $$

so no dependence on $\bar {I}_3$ could be inferred. The corresponding stress is

(16)$$\hat{{{\boldsymbol \sigma}}} = {{\boldsymbol \sigma}} = \left(\matrix{{0} & {0} & {\tau}\cr{0} & {0} & {0}\cr{\tau} & {0} & {0}}\right),\; \quad p = 0.$$

Thus

(17)$$S\left(\bar{I}_{2}^{1/2}\right) = \bar{I}_{2}^{1/2}\phi_{c1}( \bar{I}_2) \ \hbox{or}\ = \bar{I}_{2}^{1/2}\phi_{s1}( \bar{I}_2,\; 0) \ \hbox{or}\ = \bar{I}_{2}^{1/2} \phi_{q1}( \bar{I}_2) ,\; $$

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)₁ and (17)₁ implies

(18)$$C1\, \colon \, \ 3^{{{{1}/{2}}}}S\left(\bar{I}_{2}^{{{{1}/{2}}}}\right) = U\left(2\left[{\bar{I}_2\over 3} \right]^{{{{1}/{2}}}}\right),\; $$

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)₁ to hold. Responses not satisfying (18) imply that the coaxial relation (9)₁ cannot be valid. For the coaxial relation (9)₂, 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)₃, for which (13)₃ and (17)₃ apply. Immediately from (17)₃,

(19)$$\phi_{q1}( \bar{I}_2) = \bar{I}_2^{-{{{1}/{2}}}}S( \bar{I}_2^{{{{1}/{2}}}}) ,\; $$

relating $\phi _{q1}( \bar {I}_2)$ directly to the shear response, independent of the uni-axial response. Then from (13)₃,

(20)$$\phi_{q2}( \bar{I}_2) = \bar{I}_2^{-1}\left[3^{{{{1}/{2}}}}S\left(\bar{I}_2^{{{{1}/{2}}}}\right)- U\left(2\left[{\bar{I}_2\over 3}\right]^{{{{1}/{2}}}}\right)\right],\; $$

which is zero if (18) applies as required. Thus, the quadratic relation (9)₃ would be determined by given uni-axial and shear responses.

Analogous analyses apply to the strain-rate formulations

(21)$$\bar{{\boldsymbol D}} = \psi_{c1}( J_2) \, \hat{{{\boldsymbol \sigma}}},\; \quad \bar{{\boldsymbol D}} = \psi_{s1}( J_2,\; J_3) \, \hat{{{\boldsymbol \sigma}}},\; \quad \bar{{\boldsymbol D}} = \psi_{q1}( J_2) \, \hat{{{\boldsymbol \sigma}}} + \psi_{q2}( J_2) \left[{\hat{{{\boldsymbol \sigma}}}}^{2} - {2\over 3}J_{2}{\bf I}\right],\; $$

where the principal deviatoric stress invariants J ₂ and J ₃ are defined in (2). For the quadratic form (21)₃, the shear relation determines

(22)$$\psi_{q1}( J_2) = J_2^{-{{{1}/{2}}}}S^{-1}( J_2^{{{{1}/{2}}}}) ,\; $$

then the uni-axial relation determines

(23)$$\psi_{q2}( J_2) = J_2^{-1}\left[3^{{{{1}/{2}}}}S^{-1}\left(J_2^{{{{1}/{2}}}}\right)-{3\over 2}U^{-1}\left(( 3J_2) ^{{{{1}/{2}}}}\right)\right],\; $$

where S ⁻¹ and U ⁻¹ are the inverse functions of S and U; that is, x = U(X), X = U ⁻¹(x), y = S(Y), Y = S ⁻¹(y). The criterion for a coaxial form (21)₁ is then

(24)$$2S^{-1}( J^{{{{1}/{2}}}}) = 3^{{{{1}/{2}}}}U^{-1}\left(3^{{{{1}/{2}}}}J^{{{{1}/{2}}}}\right)\quad \hbox{for}\ J^{{{{1}/{2}}}} \geq 0,\; $$

with the limit as J ^1/2 → 0 also necessary for the form (21)₂. 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 Steinemann's (1958) experimental data for uni-axial compression and torsion of a hollow cylinder with the quadratic law (9)₃, removing the subscript q for convenience. The former is analysed above, and leads to relation (13)₃ for ϕ ₁ and ϕ ₂. The simple shear analysis must now be replaced by a torsion analysis to obtain a second relation for ϕ ₁ and ϕ ₂. Note that there is no analogous correlation with the strain-rate formulation to determine ψ ₁ and ψ ₂.

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 (1958) in his laboratory measurements; namely

(25)$$H = 3{\, \hbox{cm}},\; \quad R_{\rm i} = 1.5{\, \hbox{cm}},\; \quad R_{\rm e} = 4{\, \hbox{cm}}.$$

Fig. 1. (a) Hollow ice cylinder of height H and internal and external radii R _i and R _e, respectively, subjected to the action of torque M applied on the upper horizontal surface z = H. Θ is an initial (at t = 0) azimuth angle, and θ is a current (at t > 0) azimuth angle of an ice particle. (b) Vertical cross section through an ice sample showing its physical dimensions.

It is assumed that the rotation is linear in elevation z, so the deformation is defined by

(26)$$r = R,\; \quad \theta = {\it \Theta} + z\kappa( t) /H,\; \quad z = Z,\; $$

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 ₀ is

(27)$$\dot{\bar{\kappa}}( t) = \dot{\kappa}( t) /D_{0}.$$

The velocity is

(28)$$v_{\rm r} = 0,\; \quad v_\theta = rz\dot{\kappa}( t) /H,\; \quad v_z = 0,\; $$

with one non-zero velocity gradient

(29)$$\partial v_{\theta}/\partial z = r\dot{\kappa}( t) /H = D_{0}r\dot{\bar{\kappa}}( t) /H,\; $$

so that $r\dot {\bar {\kappa }}( t) /H$ is the dimensionless velocity gradient with unit D ₀. Thus, the dimensionless strain rate and invariants are given by

(30)$$\bar{{\boldsymbol D}} = {r\over 2H} \left(\matrix{{0} & {0} & {0}\cr{0} & {0} & {\dot{\bar{\kappa}}}\cr{0} & {\dot{\bar{\kappa}}} & {0}}\right),\; \quad \bar{{\boldsymbol D}}^{2} = {r^{2}\over 4H^{2}} \left(\matrix{{0} & {0} & {0}\cr{0} & {\dot{\bar{\kappa}}^{2}} & {0}\cr{0} & {0} & {\dot{\bar{\kappa}}^{2}}}\right),\; $$

(31)$$\bar{I} = {r^{2}\over 4H^{2}}\dot{\bar{\kappa}}^{2},\; \quad \left[\bar{{\boldsymbol D}}^{2} - {2\over 3}\bar{I}{\boldsymbol I}\right] = {r^{2}\over 12H^{2}} \left(\matrix{{-2\dot{\bar{\kappa}}^{2}} & {0} & {0}\cr{0} & {\dot{\bar{\kappa}}^{2}} & {0}\cr{0} & {0} & {\dot{\bar{\kappa}}^{2}}}\right),\; $$

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)₃ shows immediately that σ _θθ = σ _zz, so the stress components and invariant are

(32)$${{\boldsymbol \sigma}} = \left(\matrix{{\sigma_{rr}} & {0} & {0}\cr{0} & {\sigma_{zz}} & {\sigma_{z\theta}}\cr{0} & {\sigma_{z\theta}} & {\sigma_{zz}}}\right),\; \quad \hat{{{\boldsymbol \sigma}}} = {1\over 3} \left(\matrix{{2\sigma_{rr}-2\sigma_{zz}} & {0} & {0}\cr{0} & {\sigma_{zz}-\sigma_{rr}} & {3\sigma_{z\theta}}\cr{0} & {3\sigma_{z\theta}} & {\sigma_{zz}-\sigma_{rr}}}\right).\; $$

The diagonal components of the viscous law (9)₃ give only one independent relation

(33)$$\sigma_{zz}-\sigma_{rr} = \phi_{2}( \bar{I}) ( r \dot{\bar{\kappa}}) ^{2}/[ 4H^{2}] ,\; $$

so σ _zz and σ _rr will be determined by the radial and vertical momentum (equilibrium) balances and boundary conditions. Finally, apply the stress formulation (9)₃, omitting the subscript q, to obtain the shear relation

(34)$$\sigma_{z\theta} = S( \dot{\bar{\kappa}}) = \phi_{1}( \bar{I}) \bar{I}^{{{{1}/{2}}}},\; \quad \bar{I}^{{{{1}/{2}}}} = { r\over 2H}\dot{\bar{\kappa}},\; $$

which is independent of ϕ ₂, 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

(35)$$M( \dot{\kappa}) = 2\pi \sigma_{0} {\islant10% \int} _{ R_{{\rm i}}}^{ \hskip1pt R_{{\rm e}}} \sigma_{z\theta}( r,\; t) r^{2}\, {\rm d}r = {8\pi \sigma_{0}H^{3}\over \dot{\kappa}^{3}} {\islant10% \int}_{\bar{I}_{{\rm i}}}^{\hskip1pt \bar{I}_{{\rm e}}} \sigma_{z\theta}\bar{I}^{{{{1}/{2}}}}\, {\rm d}\bar{I} = {8\pi \sigma_{0}H^{3}\over \dot{\kappa}^{3}} {\islant10% \int} _{\bar{I}_{{\rm i}}}^{\hskip1pt \bar{I}_{{\rm e}}} \bar{I}\phi_{1}( \bar{I}) \, {\rm d}\bar{I},\; $$

where σ _zθ is expressed as a function of $\bar {I}$ by (34)₁, 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

(36)$$\bar{M}( \dot{\bar{\kappa}}) = {M( \dot{\kappa}) \over \sigma_{0}H^{3}} = {8\pi\over \dot{\bar{\kappa}}^{3}} {\islant10% \int}_{\bar{I}_{{\rm i}}}^{\hskip1.5pt \bar{I}_{{\rm e}}} \bar{I}\phi_{1}( \bar{I}) \, {\rm d}\bar{I},\; $$

where

(37)$$\bar{I}_{{\rm i}} = \left[{R_{i}\dot{\bar{\kappa}}\over 2H}\right]^{2},\; \quad \bar{I}_{{\rm e}} = \left[{R_{e}\dot{\bar{\kappa}}\over 2H}\right]^{2},\; $$

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)₃ from the simple shear test. The correlation is completed by determining $\phi _{2}( \bar {I})$ from the uni-axial compression relation (13)₃, thus

(38)$$\bar{I}\phi_{2}( \bar{I}) = \left[3\bar{I}\right]^{{{{1}/{2}}}}\phi_{1}( \bar{I}) - U\left(2\left[{\bar{I}\over 3} \right]^{{{{1}/{2}}}}\right),\; $$

where ϕ ₁(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

(39)$$U( \dot{\bar{\epsilon}}) \sim u_1\dot{\bar{\epsilon}},\; \quad u_1 > 0,\; \quad \bar{M}( \dot{\bar{\kappa}}) \sim m_{1}\dot{\bar{\kappa}},\; \quad m_1 > 0,\; $$

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)₁, that

(40)$$m_{1} = {\pi\phi_{1}( 0) ( R_{{\rm e}}^{4} - R_{{\rm i}}^{4}) \over 4H^{4}},\; $$

which, given m ₁ by data correlation and the dimensions (25), determines

(41)$$\phi_{1}( 0) = {4H^{4}m_{1}\over \pi( R_{{\rm e}}^{4} - R_{{\rm i}}^{4}) } = 0.4110\, m_{1}.$$

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

(42)$${\it \Phi}_{2}( \bar{I}) = \bar{I}^{{{{1}/{2}}}}\phi{_2}( \bar{I}) = 3^{{{{1}/{2}}}}\phi_{1}( \bar{I}) - U\left(2\left[{\bar{I}\over 3}\right]^{{{{1}/{2}}}}\right)\bar{I}^{-{{{1}/{2}}}},\; $$

where, by (39)₁, the latter term is bounded in the limit of zero $\bar {I}$, determining the limit relation

(43)$${\it \Phi}_{2}( 0) = 3^{{{{1}/{2}}}}\left[\phi_{1}( 0) -{2\over 3}u_1\right].$$

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)₃ 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)₃ 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),

(44)$$T = -1.9^{\circ}\hbox{C},\; \quad a( T) = 0.49.$$

The conversion of strain rates in units s⁻¹ to the required unit ${\, \hbox {a}}^{-1}$ has a factor 3.15 × 10⁷, but a further factor 1/a(T) = 2.04 is required to convert to the normalised strain rate D defined in (3)₁; that is a total factor 6.4 × 10⁷.

Table 1. Torsion test data

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.

Table 2. Uni-axial compression test data

σ _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, 2019) when $\phi _{2}( \bar {I})$ is zero, it is supposed that $\phi _{1}( \bar {I})$ is positive with strictly positive ϕ ₁(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

(45)$$\eqalign{ & U( \dot{\bar{\epsilon}}) \geq 0,\; \quad U^{\prime}( \dot{\bar{\epsilon}}) > 0,\; \quad U^{\prime\prime}( \dot{\bar{\epsilon}}) < 0,\; \cr & M( \dot{\bar{\kappa}}) \geq 0,\; \quad M^{\prime}( \dot{\bar{\kappa}}) > 0,\; \quad M^{\prime\prime}( \dot{\bar{\kappa}}) < 0,\; \cr & \phi_{1}( \bar{I}) > 0,\; \quad \phi_{1}^{\prime} ( \bar{I}) < 0,\; \quad \phi_{1}^{\prime \prime}( \bar{I}) > 0.}$$

However, it is $\bar {I}^{1/2}$ which defines the strain-rate magnitude, so the correlation with data is improved by expressing ϕ ₁ in terms of $\bar {I}^{{{{1}/{2}}}}$. Each term of the following adopted expansions satisfy the conditions (45):

(46)$$U( \dot{\bar{\epsilon}}) = \sum_{k = 1}^{K} a_k^{2}\left[b_k^{-2c_k^{2}}-\left(b_k^{2} + \dot{\bar{\epsilon}}\right)^{-c_k^{2}}\right],\; \quad u_1 = \sum_{k = 1}^{K} a_k^{2} c_k^{2} b_k^{-2( c_k^{2} + 1) }> 0,\; $$

(47)$$\bar{M}( \dot{\bar{\kappa}}) = \sum_{l = 1}^{L} d_l^{2}\left[e_l^{-2f_l^{2}}-\left(e_l^{2} + \dot{\bar{\kappa}}\right)^{-f_l^{2}}\right],\; \quad m_1 = \sum_{l = 1}^{L} d_l^{2} f_l^{2} e_l^{-2( f_l^{2} + 1) }> 0,\; $$

(48)$$\phi_1( \bar{I}) = \phi_1( 0) + \sum_{n = 1}^{N} g_n^{2}\left[\left(h_n^{2} + \bar{I}^{\, 1/2}\right)^{-k_n^{2}}-h_n^{-2k_n^{2}}\right],\; \quad \phi_{1}( 0) = 0.4110\, m_{1} > 0,\; $$

where ϕ ₁(0) is related to m ₁ 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

(49)$$\theta_{\rm i} = \bar{I}_{{\rm i}}^{{{{1}/{2}}}} = {R_{{\rm i}}\dot{\bar{\kappa}}\over 2H},\; \quad \theta_{\rm e} = \bar{I}_{{\rm e}}^{{{{1}/{2}}}} = {R_{{\rm e}}\dot{\bar{\kappa}}\over 2H},\; $$

then the integration yields

(50)$$\eqalign{& \bar{M}( \dot{\bar{\kappa}}) = {4\pi\over \dot{\bar{\kappa}}^{3}}\phi_{1}( 0) \left[\theta^{4} \right]_{\theta_{{\rm i}}}^{\theta_{{\rm e}}} \cr & \quad + {16\pi\over \dot{\bar{\kappa}}^{3}}\sum_{n = 1}^{N} g_{n}^{2}\left[-{h_{n}^{-2k_{n}^{2}}\, \theta^{4}\over 4} + {{( h_{n}^{2} + \theta) }^{1-k_{n}^{2}}\, \theta^{3}\over ( 1-k_{n}^{2}) } - {3\, {( h_{n}^{2} + \theta) }^{2-k_{n}^{2}}\, \theta^{2}\over ( 1-k_{n}^{2}) ( 2-k_{n}^{2}) }\right. \cr & \quad \left. + {6\, {( h_{n}^{2} + \theta) }^{3-k_{n}^{2}}\, \theta\over ( 1-k_{n}^{2}) ( 2-k_{n}^{2}) ( 3-k_{n}^{2}) } - {6\, {( h_{n}^{2} + \theta) }^{4-k_{n}^{2}}\over ( 1-k_{n}^{2}) ( 2-k_{n}^{2}) ( 3-k_{n}^{2}) ( 4-k_{n}^{2}) } \right]_{\theta_{\rm i}}^{\theta_{\rm e}} .}$$

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 ϕ ₂ 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

(51)$$\matrix{K = 2\, \colon \, & a_1 = 0.7609,\; & b_1 = 0.5350,\; & c_1 = 1.1640,\; \cr & a_2 = 7.5523,\; & b_2 = 2.7181,\; & c_2 = 0.3107,\; \cr & u_1 = 15.546. }$$

Fig. 2. Function $U( \dot {\bar {\epsilon }})$ obtained by correlating (46) with 16 uni-axial compression data points (squares).

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

(52)$$\matrix{L & = 2\, \colon \, & d_1 = 224.80,\; & e_1 = 0.3993,\; & f_1 = 0.0095,\; \cr & & d_2 = 520.31,\; & e_2 = 214.76,\; & f_2 = 77.869,\; \cr & & m_1 = 28.778. }$$

Fig. 3. Dimensionless torque function $\bar {M}( \dot {\bar {\kappa }})$ (solid line) obtained by correlating (47) with six torsion data points (squares). Solid circles show 25 correlation points used for calculating the response function $\phi _1( \bar {I})$ defined by (48) from the integral representation (36) for $\bar {M}$. The dotted line illustrates the torque $\bar {M}$ determined by the analytic formula (50).

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

(53)$$\matrix{N = 3\, \colon \, & g_1 = 1.8768,\; & h_1 = 1.2917,\; & k_1 = 1.7177,\; \cr & g_2 = 1.9507,\; & h_2 = 1.0402,\; & k_2 = 0.9309,\; \cr & g_3 = 0.7792,\; & h_3 = 0.5819,\; & k_3 = 1.5235,\; \cr & \phi_1( 0) = 11.828,\; & {\it \Phi}_2( 0) = 2.536,\; }$$

with ϕ ₁(0) given by m ₁ in (48)₂, and Φ ₂(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 Φ ₂(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 ϕ ₁ and ϕ ₂, or alternatively Φ ₂, determine the deviatoric stress $\hat {{{\boldsymbol \sigma }}}$ in terms of the strain rate D by the law (9)₃, which can be rewritten in the form

(54)$$\hat{{{\boldsymbol \sigma}}} = \phi_{1}( \bar{I}) \, \bar{{\boldsymbol D}} + \phi_{2}( \bar{I}) \left[{\bar{{\boldsymbol D}}}^{2} - {2\over 3}\bar{I}{\boldsymbol I}\right].$$

Although the response function ϕ ₁ is defined analytically by the series (48) with the coefficients and ϕ ₁(0) given by (53), the function ϕ ₂ 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.

Fig. 4. Scaled response functions $\phi _1( \bar {I})$, ${\it \Phi}_2( \bar {I})$ and $\phi _2( \bar {I})$ for $0 \leq \bar {I}^{1/2} \leq 200$ (a) and response functions $\phi _1( \bar {I})$, ${\it \Phi}_2( \bar {I})$ and $\phi _2( \bar {I})$ for $0 \leq \bar {I}^{1/2} \leq 5$ (b).

Note that inequalities (45) applied to the ϕ ₂ relation (42) do not impose any sign to ϕ ₂ or to its derivative. Specifically, ϕ ₂ 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 σ ⁽¹⁾ and σ ⁽²⁾ of the linear and quadratic terms in the law (9)₃, namely

(55)$$\sigma^{( 1) } = ( 3/2) \phi_{1}( \bar{I}) \, \dot{\bar{\epsilon}},\; \quad \sigma^{( 2) } = -\left(3^{{{{1}/{2}}}} /2 \right){\it \Phi}_{2} ( \bar{I}) \, \dot{\bar{\epsilon}}.$$

Fig. 5. Comparison of dimensionless normal stresses σ given by the proposed quadratic flow law (9)₃ with data measured in the uni-axial compression test (squares), for axial strain rates $0 \leq \dot {\bar {\epsilon }} \leq 200$ (a), $0 \leq \dot {\bar {\epsilon }} \leq 30$ (b), $0 \leq \dot {\bar {\epsilon }} \leq 1$ (c) and $0 \leq \dot {\bar {\epsilon }} \leq 0.01$ (d). σ ⁽¹⁾ and σ ⁽²⁾ represent the stresses given by the linear and quadratic terms in the flow law (9)₃.

Note that σ ⁽²⁾ 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 σ ⁽²⁾/σ ⁽¹⁾ 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 σ ⁽²⁾ 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

(56)$$\sigma \sim \dot{\bar{\epsilon}} \left[( 3/2) \phi _{1}( 0) - \left(3^{{{{1}/{2}}}} /2 \right){\it \Phi}_{2}( 0) \right] = \dot{\bar{\epsilon}} \left[17.741 - 2.195\right] = 15.546\, \dot{\bar{\epsilon}},\; $$

from which the quadratic to linear term ratio is

(57)$$\sigma^{( 2) }/ \sigma^{( 1) } = - \left[3^{-{{{1}/{2}}}}\right]{\it \Phi}_{2}( 0) /\phi _{1}( 0) = -0.1237,\; $$

which shows a Φ ₂(0) contribution $12\percnt$ of the ϕ ₁(0) contribution, not very significant, but not negligible. Figure 6a shows the ratio s _r = |σ ⁽²⁾/σ ⁽¹⁾| 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.

Fig. 6. Absolute stress ratio s _r = |σ ⁽²⁾/σ ⁽¹⁾| over a wide range of dimensionless axial stress 0 ≤ σ ≤ 16.1 (a) and over a small range 0 ≤ σ ≤ 0.2 (b).

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.

Fig. 7. Dimensionless shear stresses $\tau = \sigma _{z\theta }( \dot {\bar {\kappa }})$ in torsion given by the proposed flow law.

Finally, consider the coaxial form (9)₂, which, analogous to (41), requires ϕ _s1(0, 0) = 0.4110 m ₁, whereas, with the limit relation (39)₁, (13)₂ shows ϕ _s1(0, 0) = 2u ₁/3, which imposes a necessary zero strain-rate limit relating u ₁ and m ₁ for the coaxial form (9)₂ to be possible, namely

(58)$$u_1 = {3\over 2}\, \phi_{{\rm s}1}( 0,\; 0) = {6H^{4}m_{1}\over \pi( R_{{\rm e}}^{4} - R_{{\rm i}}^{4}) } = 0.6165\, m_{1}.$$

In practice, the correlated limit derivatives u ₁ and m ₁ given in (51) and (52) are not expected to be accurate, and a sensible test of (58) is the condition

(59)$${0.6165\, m_1-u_1\over u_1} \ll 1.$$

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 ₂, I ₃) 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.