Elastic-viscoplastic characterization of S2 columnar freshwater ice

This work addresses the time-dependent response of 3m x 6m floating edge-cracked rectangular plates of columnar freshwater S2 ice by conducting load control (LC) mode I fracture tests at -2 °C in the Ice Tank of Aalto University. The loading profile 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 5 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 elastic-viscoplastic; no viscoelasticity or major recovery were 10 detected. In addition, there was no clear effect of the creep loading on the fracture properties: the apparent fracture toughness, failure load, and crack opening displacements.

1.2 kN, 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 load up was applied in approximately 10 seconds and released in approximately 10 seconds. The slopes of the wave on load up and load release were 0.04kN/s, 0.08 kN/s, and 0.12 kN/s for the 0.4kN, 0.8 kN, and 1.2 kN load 95 levels, respectively. Once at the desired hold level, the load was kept constant for a predetermined time interval. The load intervals were multiple of the hold interval for the 0.4 kN load level, ∆t 1 = 126 seconds. For the 0.8 kN and 1.2 kN load levels, the time interval was doubled and quadrupled: 2∆t 1 = 252 seconds and 4∆t 1 = 504 seconds, respectively. The four zero load recovery periods, separating the creep load periods, were also function of ∆t 1 . Three recovery periods were held at zero load level for 5∆t 1 = 100 630 seconds, but the last recovery period was maintained for a longer interval of 10∆t 1 = 1260 seconds.
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 1 .

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 105 cyclic-recovery loading consisted of 3 sequences, each being composed of four fluctuating loads, at the levels of 0.4 kN, 0.8kN, and 1.2 kN. Each cyclic sequence continued for a constant time interval ∆t 2 = 480 seconds. The slopes of the wave on the load up and load release were 1/150 kN/s, 1/75 kN/s, and 1/50 kN/s for the 0.4kN, 0.8kN, and 1.2 kN load levels, respectively. The 0.4kN, 0.8kN, and 1.2 kN cyclic load periods were followed by zero load recovery periods of 1.25∆t 2 = 600 seconds, 1.25∆t 2 = 600 seconds, and 2.5∆t 2 = 1200 seconds, respectively. 110 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 2 = 120 seconds.

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 115 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 120 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 4 https://doi.org/10.5194/tc-2020-237 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License. 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 dis-125 placement output resulting from the whole load input. If the number of steps tends to infinity, the total displacement is given as: where C 0 is the initial, time-independent compliance component and ∆C(t) is the transient, time-dependent component of compliance.

130
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: where C 0 and ∆C are the previously defined components of Boltzmann principle, ψ and ψ are the so-called reduced times 135 defined by: and g 0 , g 1 , g 2 , and a P are nonlinear functions of the load. Each of these functions represents a different nonlinear influence on the compliance: g 0 models the elastic response, g 1 the transient response. g 2 the loading rate, and a P is a time scale shift factor.
These load-dependent properties have a thermodynamic origin. Changes in g 0 , g 1 , and, g 2 reflect third and higher order stress-140 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 0 = g 1 = g 2 = a P = 1, and Schapery's constitutive equation (2) reduces to Boltzmann's equation (1).
Equation (2) contains one time-dependent compliance property, from linear viscoelasticity theory, ∆C and four nonlinear load-dependent functions g 0 , g 1 , g 2 , and a P , which reflect the deviation from the linear viscoelastic response, that need to be 145 evaluated. Schapery's model uses experimental data to evaluate the material property functions in (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 150 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 155 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 160 theory represents the displacement at the crack mouth (δ CMOD ) as the sum of elastic, viscoelastic, and viscoplastic components: In the above equations, ψ and ψ are defined in (3). g 0 , g 1 , g 2 , g 3 , and a P are nonlinear load functions to be determined. The

170
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 experimental problem at hand is not precisely uniaxial, but it is approximated as so, and 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 0 = 1. Schapery stated that g 1 = a P = 1 if the instantaneous jump and drop in the displacement are equal (Schapery, 1969).

175
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: The viscoelastic compliance is assumed to follow a power law in time with a fractional exponent n. This gives: Incorporating each of these conditions, the total displacement is expressed as where δ CMOD , P and t are in m, N, and seconds, respectively. It follows from (10) 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 formulated as a least-squares problem and optimized through the N-M technique, by minimizing the objective 185 function F given by the difference between the model and data, as shown in (11). The components of the total displacement were computed and optimized using MATLAB. 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.
where M i and D i refer to the CMOD values given by the model (10) and the experimental data, respectively. . 2 is the Euclidean norm of a vector. N is the number of data points (≈ 2e6 points).
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, ...). In addition, the analysis does not account for the formation of fracture process zone in the vicinity of the crack tip. Schapery's formulation 195 models the experimental response until crack growth initiation and does not account for crack propagation.
4 Results and discussions

Experimental 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 (3m x 6m) (Gharamti et al., in-press). The  Table 1 shows the measured and computed parameters for the LC experiments. t f represents the time to failure, computed from the fracture ramp. From the failure load and dimensions, an apparent fracture toughness (K Q ) is computed using the weight function procedure outlined in Section 4 of (Gharamti et al., in-press). CMOD is computed at crack growth initiation.ĊMOD 205 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.ṄCOD1 indicates the displacement rate in the vicinity of the tip and is obtained by dividing NCOD1 by the failure time. Fig. 4 gives the results of the apparent fracture toughness K Q , crack mouth opening displacement CMOD, and near crack-tip opening displacement NCOD1 as a function of the loading rate for the DC tests (Gharamti et al., in-press) Table 1 in (Gharamti et al., in-press) presents the elastic modulus (E CMOD ) calculated at the crack mouth for the DC tests;

225
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 post-peak response as can be distinguished from Fig. 5b which displays more decay behavior than Fig. 5a. The gradual decay of the load portrays the time dependency in the behavior of freshwater ice.  Fig. 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 Fig. 3a and Section 2.2) and ≥ 1.25∆t 2 (Cyclic test, Fig. 3b  Figs. 6a and 6b 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 245 loading cycle.
Interestingly, the ice behavior in the current study differs from previous experimental creep and cyclic work on freshwater ice. Large delayed elastic or recoverable component has been previously observed. Several researchers performed creep experiments at lower temperatures (Mellor and Cole, 1981Cole, , 1982Cole, , 1983Cole, 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).

250
When unloaded, the ice exhibited creep recovery. According to his analysis: during loading, the internal stresses opposing the dislocation motion increases; upon unloading, the movement of dislocations produced the reversible deformation and is caused by the relaxation of internal stresses. Sinha (Sinha, 1978;Sinha et al., 1979) concluded that high-temperature creep of polycrystalline ice 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.

255
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 columnar freshwater ice tested at -2°C did not show a delayed elastic effect, and the microstructural changes were mainly irreversible upon unloading?
An influencing factor is the mechanisms taking place in the process zone. There is possibility that the loading conditions 260 produced dislocations which would ordinarily generate some viscoelastic deformation upon unloading. However, local damage in the process zone relieved the internal stresses that are needed to drive the dislocation recovery. 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., in-press) that scale had an effect at the tested loading rates. It is probable that the specimen size influenced the time-dependent 265 deformation of freshwater ice. When the specimen dimensions are several meters, apparently viscoelasticity is not an important deformation component.
The ice response indicates that the combined effects of the geometry, the applied loading profiles, the warm temperature (-2°C), and the testing conditions triggered an instantaneous transformation from the primary stage to the steady-state stage, resulting in permanent irreversible deformation that accumulated after each creep/cyclic-recovery sequence. This concludes 270 that the response of columnar freshwater S2 ice in these tests is elastic-viscoplastic.

Schapery-optimization modelling analysis
The nonlinear theory, outlined in Section 3, was used to analyze the experiments. The results of the initial optimization trials confirmed the previous analysis; the viscoelastic component δ ve CMOD had no effect on the final fit between the data and the model. The variation of the constants corresponding to δ ve CMOD , κ, a, b, d, and, n didn't affect the final converged values of 275 C e , C vp , and c.
The final optimization runs were carried out by considering the elastic and viscoplastic components (first and last terms of Eq. (10)) only. This resulted in 2 parameters, C e and C vp , and one exponent c, that need to be optimized. The optimization converged results are given in Table 2 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 recovery displacements were fitted accurately by the model in the three recovery periods. The experimental response for the creep-recovery test RP16 appeared to conform to the model results, but the model under estimated the recovery displacement in the first two cycles with a maximum misfit of ≈ 2e −5 µm. Schapery's model has been tested for creep-recovery sequences 290 of saline ice with an increasing load profile (Schapery, 1997;Adamson and Dempsey, 1998;LeClair et al., 1999LeClair et al., , 1996. This is the first application of the model with a load profile of increasing and decreasing load levels (Fig. 3a). The model succeeded to follow the increasing and decreasing load levels and the corresponding recovery phases. The model generated the peak displacement values for all load levels with a misfit of 1e −5 µm for the last two creep cycles. The observed misfit (1-2 e −5 µm) is small and should be ignored. It is related to the accuracy of the measurement line (LVDT + amplifier + data processing 295 unit) which is affected by many environmental and technical factors. Thus, the implemented model provided a good fit with the data over the creep-recovery and cyclic-recovery sequences.
Considering the fracture ramp, Schapery's nonlinear equation succeeded to model 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 The applicability of the proposed model and the fitted parameters are 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.

Conclusions
In the present work, five 3m x 6m freshwater S2 ice specimens were tested under creep/cyclic-recovery sequences followed by a monotonic ramp at -2°C. The tests were load controlled and led to complete fracture of the specimen. The purpose of this 315 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, and fracture parameters were computed and analyzed. The results showed that the creep and cyclic sequences 320 had no clear effect on the apparent fracture toughness, the failure load, and the crack opening displacements. The ice response at the testing conditions was elastic-viscoplastic, and the steady-state regime dominated the loading phases. 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 (Duval, 1978), no viscoelasticity was detected in this study. The primary (transient) creep stage was absent. This indicates that the collective effects of the testing conditions and plate configuration 325 triggered an instantaneous transformation to the steady-state regime resulting in permanent (unrecoverable) displacement.
In the modeling part, Schapery's nonlinear constitutive model was applied for the displacement response at the crack mouth.
The elastic-viscoplastic formulation succeeded to predict 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.

330
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.
Code and data availability. The code used for material modeling is written in Matlab. Scripts used for analysis and more detailed information   is the part of that total that is due to the viscoelastic mechanism and the rest is due to viscous processes.