Creep and fracture of warm columnar freshwater ice

This work addresses the time-dependent response of 3 m×6 m floating edge-cracked rectangular plates of columnar freshwater S2 ice by conducting load control (LC) mode I fracture tests in the Aalto Ice Tank of Aalto University. The thickness of the ice plates was about 0.4 m and the temperature at the top surface about −0.3 ^∘C. The loading was applied in the direction normal to the columnar grains and consisted of creep/cyclic-recovery sequences followed by a monotonic ramp to fracture. The LC test results were compared with previous monotonically loaded displacement control (DC) experiments of the same ice, and the effect of creep and cyclic sequences on the fracture properties were discussed. To characterize the nonlinear displacement–load relation, Schapery's constitutive model of nonlinear thermodynamics was applied to analyze the experimental data. A numerical optimization procedure using Nelder–Mead's (N-M) method was implemented to evaluate the model functions by matching the displacement record generated by the model and measured by the experiment. The accuracy of the constitutive model is checked and validated against the experimental response at the crack mouth. Under the testing conditions, the creep phases were dominated by a steady phase, and the ice response was overall elastic–viscoplastic; no significant viscoelasticity or major recovery was detected. In addition, there was no clear effect of the creep loading on the fracture properties at crack growth initiation: the failure load and crack opening displacements.

1 Introduction

Understanding the deformation and fracture processes of columnar freshwater ice is important in many engineering problems. For example, freshwater ice sheets fracture when in contact with ships, river ice fractures during interaction with bridge piers, and thermal cracks form in lakes and reservoirs. Deformation and fracture processes of freshwater ice are highly dependent on temperature, strain rate, sample size, grain type, and grain size. Qualitatively, high temperature and low strain rate lead to viscous behavior and ductile fracture; low temperature and high strain rate lead to elastic behavior and brittle fracture (Gharamti et al., 2021). However, quantitatively these relations are not well known.

As the response of freshwater ice is time-dependent, a general constitutive model should incorporate elastic (immediate and recoverable), viscoelastic (or delayed elastic, time-dependent, and recoverable), and viscoplastic (time-dependent and unrecoverable) components (Jellinek and Brill, 1956; Sinha, 1978). The importance of each component depends on the problem studied. For example, thermal deformations of ice in dams can have a timescale of a few days and creep behavior dominates. In ice–structure interaction problems, the timescale of interest is often seconds and hours, so all three components of deformation need to be modeled.

This paper reports results from laboratory experiments which were conducted to study the time-dependent response and fracture of columnar freshwater ice. The work is directly relevant to a number of practical problems (Ashton, 1986) but has also general relevance in ice research by studying the coupled creep and fracture in a quasi-brittle material. Unless just short timescales are involved, where only elastic response is relevant, the creep deformations must be modeled to obtain the true fracture behavior. In materials with time-dependent properties, the fracture and creep responses are coexistent.

The time-dependent behavior of freshwater ice has been addressed with great attention, and several constitutive models were developed (Michel, 1978; Sinha, 1978; Le Gac and Duval, 1980; Ashby and Duval, 1985; Sunder and Wu, 1989; Mellor and Cole, 1983; Cole, 1990; Duval et al., 1991; Sunder and Wu, 1990; Abdel-Tawab and Rodin, 1997; Santaoja, 1990). Constitutive laws can be phenomenological or micromechanical. Micromechanical modeling in ice faces challenges because the characterization of the microscopic mechanisms of ice deformation is still inadequate (Abdel-Tawab and Rodin, 1997).

Phenomenological laws are classified into two groups. The first group is empirically based relations (Sinha, 1978; Schapery, 1969). Their equations relate macroscopic variables: stress/load, strain/displacement, and time. They do not contain state variables that describe the internal state of the material and are valid only for constant stress/load. The functions in these models can be easily calibrated to simulate the experiments. The second group of phenomenological models starts from physically based models involving internal state variables (dislocation density, internal stresses reflecting hardening, etc.); they develop differential equations for the evolution of these variables with time and quantify the dependence of these variables on stress, temperature, and strain (Le Gac and Duval, 1980; Sunder and Wu, 1989, 1990; Abdel-Tawab and Rodin, 1997). These models provide insights into the microscopic mechanisms taking place, and the state variables describe the deformation resistance offered by changes in the microstructure of the material. However, they require a proper identification of the deformation mechanisms.

The effect of time-dependent loading on the strength of freshwater ice has been examined in the literature. Subjecting freshwater ice to cyclic loading apparently leads to a significant increase in the tensile, compressive, flexural strength, and fracture toughness of that ice (Murdza et al., 2020; Iliescu et al., 2017; Iliescu and Schulson, 2002; Jorgen and Picu, 1998; Rist et al., 1996; Cole, 1990). On the other hand, no detailed investigation of the effect of creep and cyclic loading on the fracture properties of freshwater ice has been conducted in the past.

Laboratory experiments were conducted to measure the time-dependent response and fracture behavior of 3 m×6 m floating edge-cracked rectangular plates of columnar freshwater S2 ice, loaded in the direction normal to the columnar grains. The ice studied was warm, and the temperature at the top surface of the samples was about −0.3 ^∘C. Compared to earlier studies with freshwater ice, the samples were large (3 m×6 m) and very warm. A program of five load control (LC) mode I fracture tests was completed in the test basin (40 m square and 2.8 m deep) at Aalto University. Creep/cyclic-recovery sequences were applied below the failure loads, followed by monotonic ramps leading to complete fracture of the specimen. The LC results were compared with the fracture results of monotonically loaded displacement control (DC) tests of the same ice (Gharamti et al., 2021), and the effects of the creep and cyclic sequences on the fracture properties were analyzed.

The constitutive modeling used in this paper was presented by Schapery (1969) and applied to polymers. Schapery's model belongs to the first phenomenological group and originates from the theory of nonlinear thermodynamics. This study presents the first attempt to use Schapery's model for freshwater ice. The choice of this model for freshwater ice is motivated by the fact that the model was successfully applied to saline ice (Schapery, 1997; Adamson and Dempsey, 1998; LeClair et al., 1999, 1996) with encouraging results. The model accurately described the deformation response during load/unload applications over varying load profiles.

The experiments in this study aim to assess the time-dependent nature of warm columnar freshwater S2 ice. In particular, the study aims to examine (1) the extent to which the elastic, viscoelastic, and viscoplastic components contribute to the ice deformation as defined through the crack mouth opening displacement; (2) the effects of the testing conditions on the creep stages (primary/transient and steady state/secondary) present in the ice; (3) the effects that creep and cyclic sequences have on the fracture properties – i.e., failure load and crack growth initiation displacements; and (4) the ability of Schapery's nonlinear constitutive model to predict the experimental response.

The rest of the paper is structured as follows. In Sect. 2, a description of the experimental setup, testing conditions, and the applied loading profile is presented. Section 3 introduces Schapery's model, which is used to analyze the experiments. In Sect. 4, the experimental and model results are summarized and analyzed. Section 5 concludes the paper.

Figure 1 Specimen geometry. Edge-cracked rectangular plate of length L, width H, and crack length A₀.

Figure 2 (a) Temperature profile. Each data point represents the average of measurements taken at the same depth of different ice cores throughout the 1-month duration of the test program. (b) Grain size distribution. Each data point is measured from one thin section.

2 Creep-recovery fracture experiments

2.1 Experimental details

The ice specimens tested were 3 m×6 m rectangular plates, cut from a 40 m×40 m parent sheet, with a thickness of 340–410 mm, and instrumented as shown in Fig. 1. The experiments were conducted at an ambient temperature of −2 ^∘C. The ice was columnar freshwater S2 ice having a mean grain size of 6.5 mm (Fig. 2b). The temperature at the top surface was about −0.3 ^∘C, as shown in Fig. 2a. An edge crack of length A₀ (A₀≈0.7 L) was cut and tip-sharpened in each ice specimen. The response of the ice was monitored by using a number of surface-mounted linear variable differential transducers (LVDTs). LVDTs were placed at five different locations along the crack to measure directly the crack opening displacements. Figure 1 labels these positions as CMOD, COD, NCOD1, NCOD2, and NCOD3 for the crack mouth, intermediate crack, 10 cm behind the initially sharpened tip, 10 cm ahead of the tip, and 20 cm ahead of it, respectively. A hydraulically operated device was inserted in the mouth of the pre-crack to load the specimen horizontally, in the direction normal to the columnar grains, with a contact loading length of 150 mm, denoted by D in Fig. 1. The tests were load controlled by a computer-operated closed-loop system that also recorded the displacement measurements. Creep/cyclic-recovery sequences were applied below the failure loads, followed by monotonic ramps leading to complete fracture of the specimen. The loading rate used is similar to that used in earlier sea ice studies (LeClair et al., 1999; Adamson and Dempsey, 1998) and thus allows for comparison of these two materials. The global behavior of the crack propagation was straight through the gauges. Detailed description of the experimental setup, ice growth, microstructure, and fractographic analysis is provided in Gharamti et al. (2021).

2.2 Creep-recovery and monotonic loading profile

In two tests, ice specimens were subjected to creep-recovery loading followed by a monotonic fracture ramp. The creep-recovery sequences consisted of four constant load applications, separated by zero-load-recovery periods. Each sequence was composed of alternating load/hold and release/recovery periods. Creep phases were applied at load levels of 0.4, 0.8, 1.2, and 0.4 kN, as given by the loading signal in Fig. 3a. The loads were chosen low enough to avoid crack propagation and failure of the specimen. Each load–hold–unload was applied in the form of a trapezoidal wave function to avoid instantaneous load jump and drop; the loading was applied in approximately 10 s and released in approximately 10 s. The slopes of the wave on loading and load release were 0.04, 0.08, and 0.12 kN s⁻¹ for the 0.4, 0.8, and 1.2 kN load levels, respectively.

Figure 3 Loading consisting of (a) creep-recovery and (b) cyclic sequences followed by a monotonic fracture ramp. The number above each segment indicates the duration in seconds.

Once at the desired hold level, the load was kept constant for a predetermined time interval. The load intervals were multiples of the hold interval for the 0.4 kN load level, Δt₁=126 s. For the 0.8 and 1.2 kN load levels, the time interval was doubled and quadrupled: 2Δt₁=252 s and 4Δt₁=504 s, respectively. The four zero-load-recovery periods, separating the creep load periods, were also a function of Δt₁. Three recovery periods were held at the zero-load level for 5Δt₁=630 s, but the last recovery period was maintained for a longer interval of 10Δt₁=1260 s.

Immediately following the creep and recovery loading sequences, the specimen was loaded monotonically to failure on a load-controlled linear ramp. The ramp up to the peak load and unloading were each applied over an interval of Δt₁.

2.3 Cyclic-recovery and monotonic loading profile

In three tests, ice specimens were loaded with cyclic-recovery sequences followed by a fracture ramp, as shown in Fig. 3b. The cyclic-recovery loading consisted of three sequences, with each being composed of four fluctuating loads, at the levels of 0.4, 0.8, and 1.2 kN. Each cyclic sequence continued for a constant time interval Δt₂=480 s. The slopes of the wave on the loading and load release were $\mathrm{1}/\mathrm{150}$, $\mathrm{1}/\mathrm{75}$, and $\mathrm{1}/\mathrm{50}$ kN s⁻¹ for the 0.4, 0.8, and 1.2 kN load levels, respectively. The 0.4, 0.8, and 1.2 kN cyclic load periods were followed by zero-load-recovery periods of 1.25Δt₂=600 s, 1.25Δt₂=600 s, and 2.5Δt₂=1200 s, respectively.

At the completion of the cyclic-recovery loading sequences, the specimen was loaded to failure by a monotonic linear ramp. The ramp up to the peak load and unloading were each applied over an interval of 0.25Δt₂=120 s.

3 Nonlinear time-dependent modeling of S2 columnar freshwater ice

The model applied in this section to characterize the nonlinear viscoelastic/viscoplastic response of S2 columnar freshwater ice was presented by Schapery; it was used to model the time-dependent mechanical response of polymers in the nonlinear range under uniaxial stress-strain histories (Schapery, 1969). Schapery's stress-strain constitutive equations are derived from nonlinear thermodynamic principles and are very similar to the Boltzmann superposition integral form of linear theory (Flügge, 1975). Schapery's model represents the material as a system of an arbitrarily large number of nonlinear springs and dashpots.

The equations in this section are presented in terms of load and displacement instead of the original stress-strain relations. The notations of the original equations in Schapery (1969) are modified to bring out similarity between all the equations in the paper.

When the applied loads are low enough, the material response is linear. For an arbitrary load input, P=P(t) applied at t=0, Boltzmann's law approximates the load by a sum of a series of constant load inputs and describes the linear viscoelastic displacement response of the material using the hereditary integral in a single integral constitutive equation. The Boltzmann superposition principle states that the sum of the displacement outputs resulting from each load step is the same as the displacement output resulting from the whole load input. If the number of steps tends to infinity, the total displacement is given as

$$\begin{array}{}\text{(1)}& \mathit{\delta }\left(t\right)=C_{\mathrm{0}}P+\underset{\mathrm{0}}{\overset{t}{\int }}\mathrm{\Delta }C(t-\mathit{\tau }){\displaystyle \frac{\mathrm{d}P}{\mathrm{d}\mathit{\tau }}}\mathrm{d}\mathit{\tau },\end{array}$$

where C₀ is the initial, time-independent compliance component and ΔC(t) is the transient, time-dependent component of compliance.

Turning now to nonlinear viscoelastic response, Schapery developed a simple single-integral constitutive equation from nonlinear thermodynamic theory, with either stresses or strains entering as independent variables (Schapery, 1969). Using load as the independent variable, the displacement response under isothermal and uniaxial loading takes the following form:

$$\begin{array}{}\text{(2)}& \mathit{\delta }\left(t\right)=g_{\mathrm{0}}{C}_{\mathrm{0}}P+g_{\mathrm{1}}\underset{\mathrm{0}}{\overset{t}{\int }}\mathrm{\Delta }C(\mathit{\psi }-{\mathit{\psi }}^{\prime }){\displaystyle \frac{\mathrm{d}\left(g_{\mathrm{2}}P\right)}{\mathrm{d}\mathit{\tau }}}\mathrm{d}\mathit{\tau },\end{array}$$

where C₀ and ΔC are the previously defined components of the Boltzmann principle; ψ and ψ^′ are the so-called reduced times defined by

$$\begin{array}{}\text{(3)}& \mathit{\psi }=\underset{\mathrm{0}}{\overset{t}{\int }}{\displaystyle \frac{\mathrm{d}{t}^{\prime }}{a_P\left[P\right(t^{\prime }\left)\right]}}\text{ and }{\mathit{\psi }}^{\prime }=\mathit{\psi }\left(\mathit{\tau }\right)=\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}{\displaystyle \frac{\mathrm{d}{t}^{\prime }}{a_P\left[P\right(t^{\prime }\left)\right]}};\end{array}$$

and g₀, g₁, g₂, and a_P are nonlinear functions of the load. Each of these functions represents a different nonlinear influence on the compliance: g₀ models the elastic response, g₁ models the transient response, g₂ models the loading rate, and a_P is a timescale shift factor. These load-dependent properties have a thermodynamic origin. Changes in g₀,g₁, and g₂ reflect third- and higher-order stress dependence of the Gibbs free energy, and changes in a_P are due to similar dependence of both entropy production and the free energy. These functions can also be interpreted as modulus and viscosity factors in a mechanical model representation. In the linear viscoelastic case, $g_{\mathrm{0}}=g_{\mathrm{1}}=g_{\mathrm{2}}=a_P=\mathrm{1}$, and Schapery's constitutive Eq. (2) reduces to Boltzmann's Eq. (1).

Equation (2) contains one time-dependent compliance property, from linear viscoelasticity theory, ΔC, and four nonlinear load-dependent functions, $g_{\mathrm{0}},g_{\mathrm{1}},g_{\mathrm{2}}$, and a_P, which reflect the deviation from the linear viscoelastic response, that need to be evaluated. Schapery's model uses experimental data to evaluate the material property functions in Eq. (2). Lou and Schapery outlined a combined graphical and numerical procedure to evaluate these functions (Lou and Schapery, 1971). In their work, a data-reduction method was applied to evaluate the properties from the creep and recovery data. Papanicolaou et al. proposed a method capable of analytically evaluating the material functions using only limiting values of the creep-recovery test (Papanicolaou et al., 1999). Numerical methods are also employed and are the most commonly used techniques; they are based on fitting the experimental data to the constitutive equation (LeClair et al., 1999). In the current study, a numerical-experimental procedure is adopted. An optimization procedure is applied using the Nelder–Mead (N-M) method (Nelder and Mead, 1965) to back-calculate the values that achieve the best fit between the model and the experimental data. To avoid multiple fitting treatments of data and account for the mutual dependence of the functions, the properties were determined from the full data. This avoided errors that may result from separating the data into parts and estimating the functions independently from different parts.

Schapery later updated his formulation (Schapery, 1997). He added a viscoplastic term to account for the viscoplastic response of the material and stated that the total compliance can be represented as the summation of elastic, viscoelastic, and viscoplastic components. Adamson and Dempsey applied Schapery's updated constitutive equation to model the crack mouth opening displacement of saline ice in an experimental setup similar to the current study (Adamson and Dempsey, 1998). The theory represents the displacement at the crack mouth (δ_CMOD) as the sum of elastic, viscoelastic, and viscoplastic components:

$$\begin{array}{}\text{(4)}& {\mathit{\delta }}_{\text{CMOD}}={\mathit{\delta }}_{\text{CMOD}}^{\mathrm{e}}+{\mathit{\delta }}_{\text{CMOD}}^{\mathrm{ve}}+{\mathit{\delta }}_{\text{CMOD}}^{\mathrm{vp}},\end{array}$$

where

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathit{\delta }}_{\text{CMOD}}^{\mathrm{e}}=g_{\mathrm{0}}{C}_{\mathrm{e}}P,\text{(6)}& {\displaystyle }{\mathit{\delta }}_{\text{CMOD}}^{\mathrm{ve}}=g_{\mathrm{1}}\underset{\mathrm{0}}{\overset{t}{\int }}{C}_{\mathrm{ve}}(\mathit{\psi }-{\mathit{\psi }}^{\prime }){\displaystyle \frac{\mathrm{d}\left(g_{\mathrm{2}}P\right)}{\mathrm{d}\mathit{\tau }}}\mathrm{d}\mathit{\tau },\text{(7)}& {\displaystyle }{\mathit{\delta }}_{\text{CMOD}}^{\mathrm{vp}}=C_{\mathrm{vp}}\underset{\mathrm{0}}{\overset{t}{\int }}{g}_{\mathrm{3}}P\mathrm{d}\mathit{\tau }.\end{array}$$

In the above equations, ψ and ψ^′ are defined in Eq. (3). g₀, g₁, g₂, g₃, and a_P are nonlinear load functions to be determined. The coefficients C_e, C_ve, and C_vp are the elastic, viscoelastic, and viscoplastic compliances, respectively. Schapery's equation has been developed for uniaxial loading. The response of the test specimen is dominated by the normal stresses at the direction normal to the x axis, ahead of the crack (Fig. 1). This stress state can be approximated as uniaxial in the same way as in beam bending; the stress is uniaxial tension at the crack tip and then changes linearly. Thus, Schapery's equations are used to analyze the experimental data. Few assumptions are applied at this point and are based on the choices made in Adamson and Dempsey (1998). For ice, the elastic displacement is linear with load; this immediately leads to g₀=1. Schapery stated that $g_{\mathrm{1}}=a_P=\mathrm{1}$ if the instantaneous jump and drop in the displacement are equal (Schapery, 1969). Examination of the current data shows that this condition is not valid, and the functions need to be evaluated. Accordingly, the following approximations are employed:

$$\begin{array}{}\text{(8)}& g_{\mathrm{1}}\propto P^a;g_{\mathrm{2}}\propto P^{b-\mathrm{1}};g_{\mathrm{3}}\propto P^{c-\mathrm{1}};a_P\propto P^d.\end{array}$$

From Eq. (3),

$$\begin{array}{}\text{(9)}& \mathit{\psi }-{\mathit{\psi }}^{\prime }=\underset{\mathit{\tau }}{\overset{t}{\int }}{\displaystyle \frac{\mathrm{d}{t}^{\prime }}{a_P\left[P\right(t^{\prime }\left)\right]}}.\end{array}$$

The viscoelastic compliance is assumed to follow a power law in time with a fractional exponent n. This gives

$$\begin{array}{}\text{(10)}& C_{\mathrm{ve}}\left(\mathit{\beta }\right)\approx \mathit{\kappa }{\mathit{\beta }}^n.\end{array}$$

Incorporating each of these conditions, the total displacement is expressed as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\mathit{\delta }}_{\text{CMOD}}& =C_{\mathrm{e}}P+\mathit{\kappa }{P}^a\underset{\mathrm{0}}{\overset{t}{\int }}{\left(\underset{\mathit{\tau }}{\overset{t}{\int }}{\displaystyle \frac{\mathrm{d}{t}^{\prime }}{\left[P\right(t^{\prime }){]}^d}}\right)}^n\\ & \cdot {\displaystyle \frac{\mathrm{d}\left[P\right(\mathit{\tau }){]}^b}{\mathrm{d}\mathit{\tau }}}\mathrm{d}\mathit{\tau }+C_{\mathrm{vp}}\underset{\mathrm{0}}{\overset{t}{\int }}{P}^c\mathrm{d}\mathit{\tau },\end{array}\end{array}$$

where δ_CMOD, P, and t are in meters, newtons, and seconds, respectively. It follows from Eq. (11) that two unknown parameters (C_e and C_vp), one unknown constant (κ), and five unknown exponents (a, b, c, d, and n) need to be determined. As previously mentioned, the problem is optimized through the N-M technique, by minimizing the objective function ℱ given by the difference between the model and data, as shown in Eq. (12). The components of the total displacement were computed and optimized using MATLAB. A positive constraint was applied to the model variables. Initial guesses of the exponents on the load and time functions were assumed based on previous work on saline ice. The optimized values were then obtained by comparing the model response and the experimentally measured response over the full length of the test up to crack growth initiation.

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathcal{F}& =\arg {\min }_{C_{\mathrm{e}},C_{\mathrm{vp}},a,b,\mathrm{\dots }}\\ & \cdot \sum _{i=\mathrm{1}}^N\Vert M_i(C_{\mathrm{e}},C_{\mathrm{vp}},\mathit{\kappa },a,b,c,d,n)-D_i{\Vert }_{\mathrm{2}},\end{array}\end{array}$$

where M_i and D_i refer to the CMOD values given by the model (11) and the experimental data, respectively. $\Vert .{\Vert }_{\mathrm{2}}$ is the Euclidean norm of a vector. N is the number of data points ($\approx \mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$ points). This problem is typically called a least-squares problem when using the Euclidean norm. It is a convex problem because ℱ is a convex function and the feasible set is convex. Thus, the optimization algorithm will converge to the global optimal solution.

As mentioned earlier, Schapery's model originated from the thermodynamic theory. The model is not physically based, and its parameters are not linked to the microstructural properties of the ice (dislocation density, grain size, etc.). In addition, the analysis does not account for the formation of a fracture process zone in the vicinity of the crack tip. Schapery's formulation models the experimental response until crack growth initiation and does not account for crack propagation.

4 Experimental and modeling results

This section presents the results measured and computed for the LC tests. The current results are compared with the fracture results of monotonically loaded DC tests of the same ice and same specimen size (3 m×6 m) (Gharamti et al., 2021). The main aim is to elucidate the effect of creep and cyclic sequences on the fracture properties.

Table 1 Measured experimental data and computed results for the LC tests.

4.1 Effect of the creep and cyclic sequences on the fracture properties

Table 1 shows the measured and computed parameters for the LC experiments. P_max is the measured peak load, which is also the failure load. t_f represents the time to failure, computed from the fracture ramp. CMOD is measured at crack growth initiation. $\dot{\mathrm{CMOD}}$ indicates the displacement rate at the crack mouth and is obtained by dividing CMOD by the failure time. Similarly, NCOD1 (see Fig. 1) represents the displacement at crack growth initiation near the initially sharpened crack tip. $\dot{\mathrm{NCOD}\mathrm{1}}$ indicates the displacement rate in the vicinity of the tip and is obtained by dividing NCOD1 by the failure time.

Figure 4 Experimental results for the (a) peak load P_max, (b) crack mouth opening displacement CMOD, and (c) near-crack-tip opening displacement NCOD1 at crack growth initiation, as a function of time to failure t_f for the monotonically loaded DC tests (Gharamti et al., 2021) and the creep/cyclic and monotonically loaded LC tests.

Figure 4 gives the results of the peak load P_max, crack mouth opening displacement CMOD, and near-crack-tip opening displacement NCOD1 as a function of the loading time for the DC tests (Gharamti et al., 2021) and the current LC tests. In these subplots, first-order power-law fits were applied to the data of the DC tests. The LC values lie above, below, and along the DC fit. No clear effect of creep and cyclic loading on the fracture properties was detected.

Figure 5 Measured load versus CMOD for the (a) DC tests (Gharamti et al., 2021), (b) LC tests, and (c) LC tests up to the peak load.

Figure 5a and b show the experimental load versus the crack opening displacement at the crack mouth for the DC and the LC tests, respectively. Figure 5c displays a magnified view of the fracture ramp of the LC tests. Comparing the failure loads of the DC and LC tests indicates that the failure loads, of tests with comparable loading rates, were similar. Therefore, in these experiments, the creep and cyclic sequences had no influence on the failure load.

Table 1 presents several elastic moduli for each test. The elastic moduli were calculated from the load–CMOD record following Sect. 4 of Gharamti et al. (2021). For the creep tests (RP15 and RP16), this procedure is repeated for the four creep cycles, resulting in E₁, E₂, E₃, and E₄, and for the fracture ramp, resulting in E_f. Similarly for the cyclic tests (RP17, RP18, and RP19), the moduli calculation was done for the last cycle of each cyclic sequence, giving steady-state moduli E₁, E₂, and E₃, and for the fracture ramp, resulting in E_f. Some of the values are missing, caused by the fact that the initial portion of the associated load–CMOD curve was very noisy. The values of the elastic moduli calculation for the creep/cyclic sequences and fracture ramps were similarly linear upon load application, as shown by the loading slope in Fig. 5c and Fig. 6a and b. This linearity justifies the choice of g₀=1 in the elastic CMOD component in Eq. (5).

Figure 6 Load versus CMOD over the (a) creep-recovery cycles for RP15 and the (b) cyclic-recovery sequences for RP17. (c) Schematic illustration of the hysteresis load-displacement diagram. The whole of the hysteresis loop area is the energy loss per cycle. The dashed area is the part of that total that is due to the viscoelastic mechanism and the rest is due to viscous processes.

Table 1 in Gharamti et al. (2021) presents the elastic modulus (E_CMOD) calculated at the crack mouth for the DC tests; E_CMOD is similar to E_f in Table 1 here; both values lie within the same range. Therefore, the creep and cyclic sequences preceding the fracture ramp did not affect the load–CMOD prepeak behavior. However, the sequences affected the postpeak response as can be distinguished from Fig. 5b, which displays a longer decay behavior than Fig. 5a. The gradual decay of the load portrays the time dependency in the behavior of freshwater ice.

Figure 7 Experimental results for RP16. (a) Load at the crack mouth; see Fig. 1. (b) Displacement–time records. (c) Load–displacement record. (d) Typical response of a Maxwell model, consisting of a nonlinear spring and nonlinear dashpot, to a constant load step.

Figure 8 Experimental results for RP17. (a) Load at the crack mouth; see Fig. 1. (b) Displacement–time records. (c) Load–displacement record.

4.2 Ice response under the testing conditions

Figure 7 shows the experimental results for RP16: the applied load and the crack opening displacements at the crack mouth (CMOD), halfway of the crack (COD), and 10 cm behind the tip (NCOD1) (see Fig. 1). Similarly, Fig. 8 shows the experimental response for RP17. The time-dependent nature of the ice response is evident. A complete load–CMOD curve was obtained during loading and unloading for each test of Table 1, indicating stable crack growth.

It is clear from Figs. 7b and 8b that the CMOD, COD, and NCOD1 displacements were composed mainly of elastic and viscoplastic components. No significant viscoelasticity was detected in the displacement–time records for all the tests. The primary (transient) creep stage was almost absent or instantaneous. The load sequences were characterized by a non-decreasing displacement rate at all levels. The displacement–time slope was linear and constant, indicating that the secondary/steady-state creep regime dominated during each load application. Although the recovery time was longer than the loading time, ≥1.25Δt₁ (creep test, Fig. 3a and Sect. 2.2) and ≥1.25Δt₂ (cyclic test, Fig. 3b and Sect. 2.3), the recovery (unload) phases consisted mainly of an elastic recovery (instantaneous drop) and unrecovered viscoplastic displacement. The behavior as observed resembles the response of a Maxwell model composed of a series combination of a nonlinear spring and nonlinear dashpot (Fig. 7c). There is no delayed elastic recovery, but there is the elastic response and a permanent deformation.

Figure 6a and b support the same analysis. Unlike the viscoelastic response (Fig. 6c), which displays no residual displacement in the loading and unloading hysteresis diagram, the current load–CMOD plots showed large permanent displacement after each loading cycle. This concludes that the response of columnar freshwater S2 ice in these tests was overall elastic–viscoplastic.

4.3 Nonlinear modeling analysis

The nonlinear theory, outlined in Sect. 3, was used to analyze the experiments. Modeling the viscoelastic term (second term of Eq. 11) proved to be very challenging. Instead, a simplified version was modeled by setting $a_{\mathrm{p}}=g_{\mathrm{2}}=\mathrm{1}$. The results of the initial optimization trials confirmed the previous analysis; the viscoelastic component ${\mathit{\delta }}_{\text{CMOD}}^{\mathrm{ve}}$ had no effect on the final fit between the data and the model. The optimization algorithm fine-tuned κ (Eq. 11) to a very small number (10⁻¹⁸), indicating that the best model–data fit is attained when the viscoelastic term goes to zero.

Table 2 Optimization results computed using Schapery's model.

The final optimization runs were carried out by considering the elastic and viscoplastic components (first and last terms of Eq. 11) only. This resulted in two parameters, C_e and C_vp, and one exponent, c, that need to be optimized. The optimization converged results are given in Table 2: C_e, C_vp, and c. For all the tests, the percent reduction of the objective function exceeded 95 %, and about 110 iterations were needed to reach convergence. A value of c=1 for the viscoplastic load function provided the best fit between the model and the experiment at all load levels over the total experimental time up to the peak load. The final compliance values of the elastic and viscoplastic components were in the ranges 1.8–$\mathrm{3.8}\times {\mathrm{10}}^{-\mathrm{8}}$ mN⁻¹ and 0.2–$\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ mN⁻¹ s⁻¹. respectively.

Figure 9 Experimental and model results for RP16. (a) Load at the crack mouth; see Fig. 1. (b) CMOD–time records.

Figure 10 Experimental and model results for RP17. (a) Load at the crack mouth; see Fig. 1. (b) CMOD–time records.

Figures 9 and 10 give the model results, obtained using Eqs. (4)–(10), and the experimental results for experiments RP16 and RP17, respectively. Figure 9a and b show the measured load and the load applied to the model and the measured CMOD–time record compared to the response of the model, respectively, for RP16. Figure 10 shows similar plots for experiment RP17. Test RP17 showed an excellent model–experiment fit for the three cyclic-recovery sequences over the load and unload periods. The model succeeded in following the increasing and decreasing load levels and the corresponding recovery phases. The experimental response for the creep-recovery test RP16 appeared to generally conform to the model results, but the model overestimated the recovery displacement in the first two cycles. It is unclear to the authors why the model did a better job in fitting the cyclic-recovery than the creep-recovery sequences. This probably hints at some mechanisms that took place in the creep-recovery tests and are not accounted for by Schapery's model. Schapery's model has been tested for creep-recovery sequences of saline ice with an increasing load profile (Schapery, 1997; Adamson and Dempsey, 1998; LeClair et al., 1999, 1996). This is the first application of the model with a load profile of increasing and decreasing load levels (Fig. 3).

Considering the fracture ramp, Schapery's nonlinear equation succeeded in modeling the monotonic displacement response up to crack growth initiation perfectly well for all the tests. As previously mentioned, the model does not account for crack propagation, so modeling was applied until the peak load. The model was also successful in predicting the critical crack opening displacement values at the failure load. Thus, the model gives a very close prediction of the experimental data over the whole loading profile up to the failure load. The other tests displayed the same experiment–model agreement.

Figure 11 Contribution of each individual model component to the total CMOD displacement for (a) RP16 and (b) RP17.

In this study, Schapery's constitutive model is tested for the first time for freshwater ice. The match between the model and the measured data, especially for the cyclic-recovery tests, provides a firm support of the ability of Schapery's constitutive model to describe the time-dependent response of columnar freshwater S2 ice up to crack growth initiation. Figure 11a and b show the contribution of each individual model component, elastic and viscoplastic, to the total CMOD displacement, for RP16 and RP17, respectively. As mentioned earlier, the elastic and viscoplastic components account for the total deformation. For RP16, the viscoplastic component dominated over the elastic component. For RP17, the elastic and viscoplastic components contributed equally to the total displacement.

The applicability of the proposed model and the fitted parameters is limited to the studied ice type, geometry, specimen size, ice temperature, and the current testing conditions. Variation in the operating conditions will change the dominant deformation mechanisms and the ice behavior; and accordingly, new model parameters are needed to adapt to the new response.

5 Discussion

Interestingly, the ice behavior in the current study differs from previous experimental creep and cyclic work on freshwater ice. A large delayed elastic or recoverable component has been previously observed. Several researchers performed creep experiments on granular freshwater ice at lower temperatures (Mellor and Cole, 1981, 1982, 1983; Cole, 1990; Duval et al., 1991) and reported considerable recovery. Duval conducted torsion creep tests on glacier ice at a similar testing temperature of −1.5 ^∘C (Duval, 1978). When unloaded, the ice exhibited creep recovery. According to his analysis, during loading, the internal stresses opposing the dislocation motion increase; upon unloading, the movement of dislocations produced the reversible deformation and is caused by the relaxation of internal stresses. Sinha (1978, 1979) studied columnar-grained freshwater ice and concluded that the high-temperature creep is associated with grain boundary sliding. Cole developed a physically based constitutive model in terms of dislocation mechanics (Cole, 1995) and quantified two mechanisms of anelasticity: dislocation and grain boundary relaxations. He demonstrated that the increased temperature sensitivity of the creep properties of ice within a few degrees of the melting point is due to a thermally induced increase in the dislocation density (Cole, 2020). The question then arises as to why warm columnar freshwater ice tested here showed no significant delayed elastic effect and why the microstructural changes were mainly irreversible upon unloading.

The measured ice response is a novel result for any type of ice. It is important to emphasize that in comparison with earlier freshwater ice studies, the tested samples are very warm and large. Viscoelasticity normally happens due to the elastically accommodated grain boundary sliding. Upon loading, internal stresses build up at local stress concentrations in the grain boundary geometry (triple points and grain boundary ledges). Assuming there is no microcracking, the growing stress impedes further grain boundary sliding and causes sliding in the reverse direction, giving rise to the recoverable component after unloading. However, in the present case, the measurements showed that the grain boundary sliding produced permanent deformation. Several reasons can be discussed, related to the ice temperature, microstructure, and nonlinear mechanisms in the process zone.

Concerning the effect of temperature: the warmer the temperature, the more liquid on the grain boundary. The high homologous test temperature (top ice surface $\approx -\mathrm{0.3}$ ^∘C) causes liquidity on the gain boundaries (Dash et al., 2006). The intergranular melt phase on the grain boundary renders the ice as a two-phase polycrystal and significantly influences the creep and recovery response. In fact, the grain boundary sliding then consists theoretically of two processes: (1) the sliding of grains over one another and (2) the squeezing in/out of the liquid between adjacent grains (Muto and Sakai, 1998). The shear behavior of the liquid film is a function of its properties (thickness and amount). The presence of this liquid at the triple points and the boundary acted as a resisting obstacle for the grains to shear and deform back to their original form, causing the viscoplastic deformation.

The microstructure (grain size, crystalline texture) could be another contributing factor. Sinha (1979) developed a nonlinear viscoelastic model, incorporating the grain size effect, to describe the high-temperature creep of polycrystalline materials. He concluded that delayed elastic strain exhibits an inverse proportionality with grain size. This suggests that the grain size (3–10 mm, Fig. 2b) of the ice samples is coarse enough not to produce any measurable viscoelastic deformation under the testing conditions. It is also probable that for this grain size, there were not enough local concentration points to arrest the grain boundary sliding and drive the recoverable and reverse sliding. In addition, Gasdaska (1994) discussed that regularly ordered and packed microstructures limit the amount of sliding and rearrangement and lead to less anelastic strain. The ice growth in the Aalto Ice Tank was very controlled and resulted in homogeneous ice sheet.

Knauss presented a thorough review of the time-dependent fracture models available to date (Knauss, 2015). The essence of the models is based on modeling the behavior in a finite cohesive/process zone which is attached to the traction-free crack tip. The one-parameter fracture mechanics encompassed by the apparent fracture toughness is not applicable (Dempsey et al., 2018). It is believed that the mechanisms taking place in the process zone play an influencing role in the current tests. The nonlinearity in the fracture zone relieved the internal stresses that would ordinarily accommodate the grain boundary sliding and generate some viscoelastic deformation upon unloading. Thus, any microstructural damage that occurred during loading manifested as permanent deformation at the end of the test.

It is noteworthy that the earlier studies used test sizes which are smaller than the plate size used here. It was shown in the DC fracture tests (Gharamti et al., 2021) that scale had an effect at the tested loading rates. It is probable that the specimen size influenced the time-dependent deformation of freshwater ice. The current tests suggest that for the large sample size and the kind of ice studied (very warm freshwater ice) under the loading applied, viscoelasticity is not an important deformation component. The experimental results support this prediction, but more tests are needed to make more general conclusions.

All the abovementioned factors might have contributed to the measured elastic–viscoplastic response. However, the question as to which factor influenced mostly the behavior is an important research question that requires more testing programs. Testing the effect of each factor separately requires a set of experiments that considers this factor while keeping all the other conditions fixed.

6 Conclusions

In the present work, five 3 m×6 m warm freshwater S2 ice specimens were tested under creep/cyclic-recovery sequences followed by a monotonic ramp. The temperature at the top surface was about −0.3 ^∘C. The tests were load controlled and led to complete fracture of the specimen. The purpose of this study was to examine the time-dependent behavior of freshwater ice using a joint experimental-modeling approach.

In the experimental part, the tests aimed to (1) measure and examine the time-dependent response of columnar freshwater S2 ice through the applied creep/cyclic-recovery sequences and (2) investigate the effect of creep and cyclic sequences on the fracture parameters/behavior through the fracture monotonic ramp. The current tests were compared with other monotonically loaded tests of the same ice. The results showed that the creep and cyclic sequences had no clear effect on the failure load and the crack opening displacements at crack growth initiation. The ice response at the testing conditions was overall elastic–viscoplastic. The loading phases displayed an instantaneous transformation from the primary (transient) stage to the steady-state regime, which resulted in permanent (unrecoverable) displacement. The conducted experiments provided a novel observation for the time-dependent behavior of freshwater ice. Though the delayed elastic component has been reported as a major creep component in freshwater ice, no significant viscoelasticity was detected in this study. Several factors were discussed as possibly contributing to the observed behavior: the very warm columnar freshwater ice, liquidity on the grain boundary, large sample size, coarse grain size, and nonlinear mechanisms in the fracture zone. Testing the effect of each factor on the ice response requires a different set of experiments that varies this factor only while keeping the other conditions fixed.

In the modeling part, Schapery's nonlinear constitutive model was applied for the displacement response at the crack mouth. The elastic–viscoplastic formulation succeeded in predicting the experimental response of columnar freshwater S2 ice over the applied loading profile up to crack growth initiation. The model parameters were obtained via an optimization procedure using the N-M method by comparing the model and experimental CMOD values.

The proposed model parameters are valid only for the studied ice type, geometry, specimen size, ice temperature, and the range of applied load experienced in the experiments. Schapery's model was selected in this study, as it is able to capture the sort of time-dependent behavior known to occur in ice and produces a simple and expedient way to help understand the observed behavior. More thorough analysis with a physically based approach is left to the future.