Theoretical Modelling of the Degradation Processes Induced by Freeze–Thaw Cycles on Bond-Slip Laws of Fibres in High-Performance Fibre-Reinforced Concrete

High-performance fibre-reinforced concrete (HPFRC) is a composite material in which the advantages of fibre-reinforced concrete (FRC) are combined with those of a high-performance concrete (HPC), which mitigates the weaknesses of conventional concrete and improves its overall performance. With the aim to reduce the long-term maintenance costs of structures, such as heavily loaded bridges, HPFRC is highly recommended due to its major durability performance. Specifically, its good antifreezing property makes it suitable for application in cold regions where cyclic freeze–thaw conditions cause the concrete to degrade. In this paper, a numerical simulation of the degradation processes induced by freeze–thaw cycles on bond-slip laws in HPFRC beam specimens has been developed so as to assess their effect on the flexural response of specimens as the fibres’ volume percentage changes. Their cracking strength, postcracking strength, and toughness were predicted, with the present model being able to predict the cracking strength, postcracking strength and toughness of the HPFRC beam element under bending load conditions. Its accuracy was confirmed by comparing the model predictions with experimental results.


Introduction
Over the last few years, high-performance fibre-reinforced concrete (HPFRC) has been widely used to strengthen ageing concrete structures [1][2][3][4][5][6][7] and to control the crack propagation and displacement in concrete slabs and shells, such as industrial floors, while also improving the seismic response of structural elements, such as columns, beams, and walls [8,9]. Moreover, HPFRC is highly recommended in aggressive environments (e.g., marine environments, higher altitudes, northern areas) due to its high durability which is suitable for long-term structures and heavily loaded bridges to reduce any maintenance costs [10,11]. It is also considered a sustainable material for the manufacturing of small thickness elements without steel rebars, resulting in a reduction in the CO 2 footprint [12]. It is well-known how HPFRC is a composite material in which the advantages of fibrereinforced concrete (FRC) are combined with those of a high-performance concrete (HPC), reducing the weaknesses of conventional concrete and improving its durability and mechanical performance. The addition of discontinuous fibres (e.g., steel fibres [13][14][15][16], synthetic fibres [17][18][19][20][21], natural fibres [22][23][24], basalt [25,26]), carbon and glass fibres [27]) in the HPC as well as in the concrete in general is able to significantly reduce its brittle behaviour, thus improving cracking, postcracking strength, and toughness [28][29][30] as well as its durability such as freeze-thaw resistance. Thanks to its dense microstructure, high-performance concrete also has a low permeability, resulting in a good resistance to various external agents such as chloride attacks [31] and carbonation [32] as well as freezing and thawing cycles [33][34][35][36]. Its good antifreezing property makes it suitable for application at both high altitudes and in northern areas where cyclic freeze-thaw conditions are one of the main causes of two types of concrete degradation: surface scaling, which is the loss of cement paste from the exposed surface, and internal crack growth, which makes the concrete crumble and deteriorate. Both phenomena can reduce the quality of concrete throughout its lifetime. Over the last few years, the research on evaluating the freezing resistance of HPFRC has significantly increased, with several relevant achievements having been obtained: in Feo et al. [37], the effects of 75 freeze and thaw cycles on both the dynamic moduli of elasticity, cracking and postcracking strength, as well as the toughness of HPFRC beam specimens reinforced with steel fibres were evaluated; in [38], it was reported how the incorporation of basalt fibres can reduce the influence of freeze-thaw on the damage and failure process of the beam specimen under a bending test; in [39], it was studied how mineral admixtures (e.g., blast furnace slag, fly ash, silica fume, and metakaolin) contained in the HPC matrix possess an excellent frost resistance; in [40], an experimental investigation on the freeze-thaw resistance of HPC containing air-cooled slag (AS) and water-cooled slag (WS) was conducted; in [41], it was explained how adding nanosilica to the concrete makes it durable by enhancing its properties such as impermeability, porosity, and acid resistance. However, based on our knowledge, such studies are experimental and not many predictive models have been proposed that capture the mechanical response of HPFRC, especially under freeze and thaw cycles.
In this study, the previous theoretical model developed by the authors [42], as an extension of a meso-scale formulation of a cracked hinge model implemented in a Matlab code [43], has been improved to predict the effect of freeze and thaw cycles on the flexural behaviour of HPFRC specimens as the fibres' volume percentage changes. This model is intended to take into account explicitly the behaviour of the two typical "phases" of fibre-reinforced cementitious composites (i.e., the cement-based matrix and the spread reinforcement, as well as with their interaction). The kinematics of the proposed model was inspired by the so-called "cracked-hinge" approach [44][45][46][47][48], but both the random spatial distribution and orientation of fibres and the crack-bridging effect of fibres is explicitly simulated. The present model is able to estimate the cracking, postcracking strength, and toughness of a HPFRC beam element under bending load conditions. Its reliability was confirmed by comparing the model predictions with the experimental results obtained in [37].

Outline of the Experimental Results
The present study is part of a wider research whose experimental part was already published into details in a previous paper [37]. A brief summary about the obtained results is reported hereafter, for the readers' sake.
Three different HPFRC mixtures CM0, CM1, and CM2, obtained by fixing the HPC matrix and varying the fibre volume fraction, Vf, in the set 0%, 1.25%, and 2.50%, respectively, were examined. Short steel fibres, Dramix OL 13/0.20 [49], with an aspect ratio, l f /d f , equal to 65 were chosen as the reinforcement of the HPC matrix whose mix design was provided by the manufacturer [50]. For each type of HPFRC mixture, eight standard 150 × 150 × 600 mm prismatic specimens (PS) and five standard 150 × 150 × 150 mm cubic samples (CS) were cast. At the end of the curing period, the prismatic specimens for each HPFRC mixture were subjected to 75 freeze-thaw cycles according to UNI 7087-2017 [51].
Subsequently, the prismatic specimens were tested under a four-point bending setup according to UNI 11039-2 [52] in which the vertical load (P) and the corresponding average "Crack Tip Opening Displacement" (CTOD avg ) were monitored during each test. All the cube samples were only tested in compression to evaluate, according to EN 12390-4 [53], the compressive load, F u , as well as the compressive strength, f c , for any changes of the fibres' volume fraction. Table 1 reports the average values of the first crack load, P lf , of the first crack strength, f (lf,avg) , and the equivalent postcracking strengths, f (eq(0-0.6),avg) and f (eq(0.6-3),avg) , before and after the freeze-thaw cycles for each type of HPFRC mixture.  (If,avg) , and of the equivalent postcracking strengths, f (eq(0-0.6),avg) and f (eq(0.6-3),avg) , before (NFT) and after the freeze-thaw cycles (FT) for each type of HPFRC mixture (CM0, CM1, and CM2).

Mix.
U N FT

Theoretical Model
The abovementioned experimental results were here used to improve the cracked-hinge model [42] in which a meso-mechanical approach was adopted with the aim of predicting the bending response of HPFRC beam elements under normal environmental conditions.

Assumptions and Formulation
In this new model, however, a different transition zone length and a modified bondslip law of the fibres were taken into account in order to estimate the effects of freeze-thaw cycles on the flexural behaviour of standard specimens of a length L, width b, depth h, and transversely notched at the midspan section for a depth equal to h 0 ( Figure 1 By denoting x, y, and z, the axes of a Cartesian coordinate system centre of the midspan section, a random distribution of the fibres in prismatic volume was generated (by utilising the standard random nu Matlab) in which x(f,k), y(f,k), and z(f,k) (with between 1 and nf), and α(y By denoting x, y, and z, the axes of a Cartesian coordinate system originating at the centre of the midspan section, a random distribution of the fibres inside the specimen prismatic volume was generated (by utilising the standard random number generator of Matlab) in which x (f,k) , y (f,k) , and z (f,k) (with k between 1 and n f ), and α (y,k) and α (z,k) are the three coordinates of the fibre centroid (G (f,k) ) and the two relevant angles, respectively. In order to consider the bridging effect offered by the fibres, the total number of fibres in the midspan section, n f , was determined as a ratio between the fibres' volume fraction, V f , and the cement matrix volume, V c . Furthermore, the model was based on the following assumptions: (i) The flexibility was distributed in the central part of the specimen for a length equal to "s" while a rigid body behaviour was exhibited by the remaining end parts ( Figure 2). (ii) The midspan cross-section was discretized in n c layers as shown in Figure 3. The average axial strain of the k-th layer, ε k , before crack formation, and the crack-opening displacement, w k , after the crack formation, can be easily expressed for the k-th layer (k = 1, . . . , n c ) as in Equations (1) and (2).
Equations (1) and (2) are typical of the "cracked hinge" model family (after their original formulations by Hillerborg et al. [46] and Olesen [47]). Specifically, Equation (2) rests on the assumption that plane sections remain plane and Equation (1) on the assumption of a characteristic length "s" which can be defined to convert axial displacements (at the numerator of the right-hand side) in axial strains (at the left-hand side).
(iii) Consequently, the average value of the axial stress, σ c,k , at k-th strip can be determined as a function of the axial deformation, ε k , before cracking, or a function of the crackopening displacements, w k , after cracking. The stress-strain and stress-displacement relationships assumed in this paper are reported in Section 3.2.2. (iv) A transition length, l t , was introduced in the notched cross-section (Figure 4), which starts from the top of the notch to the top of the integral part of the section, in order to consider the possible microdamage phenomena produced by the notching process. The mechanical meaning of this quantity is discussed in details in a previous paper [42] and omitted herein for the sake of brevity. Therefore, a reduced value of the width, b k , inside the transition zone was considered which can be evaluated with an exponential law as in Equation (3) where l k and α are the distance of the k-th strip from the top of the notch and the coefficient of the exponential law, respectively: (v) The bridging effect offered by the fibres was taken in to account by introducing the action, F k,j , mobilised at the j-th step of the incremental analysis as in Equation (4): in which z f ,k , z c,j and ϕ j are the position of the k-th fiber in the cross section, the position of the neutral axis and the rotation of the two rigid blocks at the j-th step of the increment analysis ( Figure 2), respectively. In particular, the axial stress of the k-th fiber, σ f ,k depends on the sliding of one of the two parts embedded in the two sides of the crack and, therefore, it was correlated to the bond stresses, τ, mobilized on its lateral surface due to the crack opening displacement w k : Moreover, it should be noted that only a low number of these fibres, n f < n f , will cross the crack which will potentially open in the middle of the beam. With the above assumptions, the position z c,j of the neutral axis for the j-th imposed rotation ϕ j can be determined by solving the following equilibrium equation along the longitudinal axis which can be written as follows: where n c is the number of layers into which the midspan section is discretized ( Figure 3) and n t is the number of layers of reduced width [42].
The solution of Equation (6) brings us to determining z c,j , at the j-th step of the incremental analysis. It is worth highlighting that this solution can only be obtained numerically: the well-known bisection method was employed to determine the actual value of z c,j which should obviously be included in the range between 0 and h-h 0 ( Figure 3).
Then, the corresponding values of the bending moment M j can be obtained as follows: which is obviously related to the applied vertical load P (in 3-point bending): The corresponding CTOD j value of can be obtained from Equation (2) by just replacing the generic value of z k with the position of the crack tip (z k = −h/2 + h 0 from Figue 2).
Therefore, for each value of the imposed rotation ϕ j , a couple (CTOD j , P j ) can be determined, and then, the Force-CTOD graph can be incrementally determined up to failure.
Furthermore, the constitutive laws adopted for the HPC matrix and short steel fibres have the same shape as well as mathematical expressions as those already presented by the authors [42]. However, the unknown parameters in the constitutive laws were calibrated in Section 4 on the experimental results already summarized in Section 2.

Stress-Strain Relationships for Concrete in Compression and in Tension
The stress-strain relationship when the concrete is in compression was described by Equation (7), according to 54, in which can be calculated as the ratio between the strain, , and the strain at maximum compressive stress, 1 , can be evaluated as in Equation (8); whereas is the plasticity number which depends on the elastic modulus in

Stress-Strain Relationships for Concrete in Compression and in Tension
The stress-strain relationship when the concrete is in compression was described b Equation (7), according to 54, in which can be calculated as the ratio between the strain , and the strain at maximum compressive stress, 1 , can be evaluated as in Equatio (8); whereas is the plasticity number which depends on the elastic modulus i

Stress-Strain Relationships for Concrete in Compression and in Tension
The stress-strain relationship when the concrete is in compression was described by Equation (7), according to [54], in which η can be calculated as the ratio between the strain, ε c , and the strain at maximum compressive stress, ε c1 , can be evaluated as in Equation (8); whereas k is the plasticity number which depends on the elastic modulus in compression, E c , and on the secant modulus from the origin to the peak compressive stress, E c1 , as in Equation (9). The last one can be evaluated as a function of ε c1 and the concrete compressive strength ( f cm ). The schematic representation of the above stress-strain relationship in compression is shown in Figure 5.
The constitutive law when the concrete is in tension, evaluated according to the "fictitious crack method" employed in [44], presents a bilinear relation until the tensile strain of the k-th strip, ε k , reaches the conventional value equal to 0.00015 (Figure 6a). A linear elastic behaviour is described by Equation (12a) until the section is uncracked after which the behaviour is expressed by the linear Equation (12b) as follows: in which σ ct , E ct , and ε ct represent the tension stress (in MPa), the elastic modulus under tension load (in MPa), and the tensile strain, respectively, whilst f ctm (in MPa) indicates the tensile strength. Beyond this level, a softening constitutive stress-crack opening law was considered due to the opening of a crack in the k-th strip of the cross-section, Equation (13). Consequentially, the residual tension σ ct must be expressed as a function of the crackopening displacement, w, (Figure 6b) in which w 1 and w c are dependent on the fracture energy as defined in [54]: Materials 2022, 15, x FOR PEER REVIEW compression, , and on the secant modulus from the origin to the pea stress, 1 , as in Equation (9). The last one can be evaluated as a function concrete compressive strength ( ). The schematic representation of the strain relationship in compression is shown in Figure 5. The constitutive law when the concrete is in tension, evaluated accord titious crack method" employed in Error! Reference source not found., pres relation until the tensile strain of the k-th strip, , reaches the convention to 0.00015 (Figure 6a). A linear elastic behaviour is described by Equation section is uncracked after which the behaviour is expressed by the linear as follows: in which , , and represent the tension stress (in MPa), the elastic m  Figure 6. Schematic graph of the stress-strain relation for uniaxial tension: (a) for ≤ 0.00015 stress-strain behaviour is described by a bilinear relation; (b) for > 0.00015, the stress-st behaviour is described by a softening constitutive stress-crack opening law.

Modified Bond-Slip Model for Short Steel Fibres
The mathematical relation of the local − constitutive law adopted in this stud provided in Equation (14) in which was considered equal to the crack-opening placement , at the level of the k-th fibre while the six unknown parameters ( , , , , , ) have to be calibrated using the experimental data of Section 2: • ≤ ( ) Figure 6. Schematic graph of the stress-strain relation for uniaxial tension: (a) for ε k ≤ 0.00015, the stress-strain behaviour is described by a bilinear relation; (b) for ε k > 0.00015, the stress-strain behaviour is described by a softening constitutive stress-crack opening law.

Modified Bond-Slip Model for Short Steel Fibres
The mathematical relation of the local τ − s constitutive law adopted in this study is provided in Equation (14) in which s was considered equal to the crack-opening displace- ment w f ,k at the level of the k-th fibre while the six unknown parameters (i.e., s el , s R , s u , τ el , τ R , τ u ) have to be calibrated using the experimental data of Section 2: Consequently, the nonlinear bond-slip law is divided in three different branches, as expressed in Equation (12): • a linear-elastic behaviour up to the stress level corresponding to matrix tensile strength, identified by the two parameters s el , τ el ; • a hardening behaviour, characterized by the formation of many microcracks in the HPFRC mix, identified by the two parameters s R and τ R ; • a constant behaviour defined by the two parameters s u and τ u .
The curve presents an adequate "shape" for describing the global pull-out response of short steel fibres embedded within the HPC matrices. For the sake of simplicity, the current assumption includes the effect of the fibre orientation in space (from 0 • to 45 • ) with respect to the matrix surfaces, as demonstrated in a previous study [55]. However, more accurate assumptions could be formulated with the aim to take into account both the axial deformation of fibres (which is neglected in this study, as it is focused on "short" fibres) and the aforementioned effect of fibre orientation with respect to the transverse section of the specimen at midspan.

Inverse Identification of the Relevant Material Laws
An inverse identification procedure was carried out with the aim to determine the values of the model parameters that lead to minimizing the difference between the measured and predicted Force-CTOD curves. Specifically, the cylindrical compression strength, f cm , the transition length l t , the exponential parameter, α, and the six parameters of the bondslip law (i.e., s el , s R , s u , τ el , τ R , τ u ) were considered as variable with some quantitative restrictions (e.g., s el < s R < s u , τ el < τ R and τ R = τ u ) intended at respecting the mechanical consistency of the model. In order to evaluate how the values of these parameters depend on the effect of freeze-thaw cycles, several numerical simulations were carried out. In particular, three groups of 100 simulations each, assuming n c = 50 and s = 300 mm, were run as described below.

•
In the first one, the cylindrical compression strength, f cm , the transition length l t , and the exponential parameter, α, were calibrated on the flexural response of the conditioned CM0 specimens (labelled CM0-FT). Experimentally, a 21% reduction in the cylindrical compression strength, f cm , was observed on the conditioned specimens compared to not conditioned ones. This reduction was taken in account to calibrate the value of the transition length, l t , whose value, in the present model, was assumed equal to 85 mm (with an increase of 21% compared to that used in the previous model [42] in which the flexural behaviour of unconditioned CM0 specimens (labeled CM0-NFT) was predicted with a transition length, l t , equal to 70 mm). In both models, the coefficient of the exponential law, α, was considered constant and equal to 0.40 (Table 3). Figure 7 shows both the average experimental P − CTOD ,avg curve (lightblue line) and the average numerical P − CTOD ,avg curve (pink line) obtained with the present model employed for the CM0-FT specimens.

•
In the second one, the six parameters of the bond-slip law (i.e., s el , s R , s u , τ el , τ R , τ u ) were calibrated on the flexural response of conditioned CM1 specimens (labelled CM1-FT). A 13% reduction in the parameter τ el was adopted in the calibration of the conditioned specimens compared to the unconditioned ones.

•
In the last one, only the parameter, τ el , was calibrated again on the flexural response of conditioned CM2 specimens (labelled CM2-FT) while all the other parameters were considered constant. A 19% reduction in the parameter τ el was adopted for the conditioned specimens compared to the unconditioned ones. Moreover, as in [42], a 20% reduction in the fibres' volume fraction, V f , was considered in order to take into account the nonuniform fibre distribution. Table 3. Calibration of the input data for unconditioned CM0 specimens (labelled CM0-NFT) and for conditioned CM0 specimens (labelled CM0-FT) used in the previous model of Ref. [42] and in the last one, respectively. CM1-FT). A 13% reduction in the parameter was adopted in the calibration of the conditioned specimens compared to the unconditioned ones.

•
In the last one, only the parameter, , was calibrated again on the flexural response of conditioned CM2 specimens (labelled CM2-FT) while all the other parameters were considered constant. A 19% reduction in the parameter was adopted for the conditioned specimens compared to the unconditioned ones. Moreover, as in Error! Reference source not found., a 20% reduction in the fibres' volume fraction, , was considered in order to take into account the nonuniform fibre distribution. Table 3. Calibration of the input data for unconditioned CM0 specimens (labelled CM0-NFT) and for conditioned CM0 specimens (labelled CM0-FT) used in the previous model of Ref. [42] and in the last one, respectively.

Specimen Designation Model [MPa]
[mm] [-] The results of the last two calibrations were listed in Figure 8 and Table 4, while Figures 9a and 10a show the comparison between the average experimental P − CTOD ,avg curve (violet line) and the average numerical P − CTOD ,avg curve (green line) obtained with the previous model for the CM1-NFT and CM2-NFT specimens, respectively. Whereas Figures 9b and 10b show the comparison between the average experimental P − CTOD ,avg curve (light-blue line) and the average numerical P − CTOD ,avg curve (pink line) used with the present model for the CM1-FT and CM2-FT specimens, respectively. The results of the last two calibrations were listed in Figure 8 and Table 4

Results
The model developed by the authors Error! Reference source not foun used to predict the postcracking response of HPFRC beam elements under cycles.
In order to assess the model accuracy, the theoretical and experimental CM1 and CM2 specimens were compared: first, the average values of the tw postcracking strengths after the freeze-thaw cycles,  Table 5, for the CM1 and CM2 mixtures, in which the values were co those obtained before the freeze-thaw cycles of Ref. [42],

Results
The model developed by the authors [42] was, here, used to predict the postcracking response of HPFRC beam elements under freeze-thaw cycles.
In order to assess the model accuracy, the theoretical and experimental results of the CM1 and CM2 specimens were compared: first, the average values of the two equivalent postcracking strengths after the freeze-thaw cycles, f FT eq(0−0.6),avg and f FT eq(0.6−3),avg , are listed in Table 5, for the CM1 and CM2 mixtures, in which the values were compared with those obtained before the freeze-thaw cycles of Ref. [42], f NFT eq(0−0.6),avg and f NFT eq(0.6−3),avg ; second, the average values of the two working capacity indices after the freeze-thaw cycles,

U FT
1,avg and U FT 2,avg , are summarised in Table 6, for the CM1 and CM2 mixtures, in which the values were compared with those obtained before the freeze-thaw cycles of Ref. [42], U NFT 1,avg and U NFT 2,avg . Table 5. Comparison between the experimental and theoretical average values of the two equivalent postcracking strengths after the freeze-thaw cycles, f FT eq(0-0.6),avg and f FT eq(0.6-3),avg , and those before the freeze-thaw cycles of Ref. [42], f NFT eq(0-0.6),avg and f NFT eq(0.6-3),avg , for the CM1 and CM2 mixtures.   The correlation graphs of the two equivalent postcracking strengths were plotted in Figures 11 and 12, respectively, while Figures 13 and 14 show the correlation graphs of the average values of the two working capacity indices U 1 and U 2 , respectively. For each point, the red bars represent the standard deviation of the experimental and predicted values. It was noted how the convergence with the present model was good so that a lower standard deviation was observed. the average values of the two working capacity indices 1 and 2 , respectively. For each point, the red bars represent the standard deviation of the experimental and predicted values. It was noted how the convergence with the present model was good so that a lower standard deviation was observed.  , and those before the freeze-thaw cycles of Ref. [42],   Table 6. Comparison between the experimental and theoretical average values of the two working capacity indices after the freeze-thaw cycles, 1, and 2, , and those before the freeze-thaw cycles of Ref. [42], 1, and 2, , for the CM1 and CM2 mixtures.

Conclusions
The present paper was aimed at scrutinizing the effects of freeze-thaw cycles on the fundamental mechanical behaviour of the "components" controlling the structural behaviour of HPFRC members in bending. Specifically, based on a series of experimental results obtained by the first two authors in a previous study Error! Reference source not found.,  Figure 14. Correlation between the experimental and theoretical average values of the working capacity indices after the freeze-thaw cycles, U FT 2,avg , and those before the freeze-thaw cycles of Ref. [42], U NFT 2,avg , for the CM1 and CM2 mixtures.

Conclusions
The present paper was aimed at scrutinizing the effects of freeze-thaw cycles on the fundamental mechanical behaviour of the "components" controlling the structural behaviour of HPFRC members in bending. Specifically, based on a series of experimental results obtained by the first two authors in a previous study [37], a simplified "cracked-hinge" model was considered with the aim to identify both the concrete constitutive relationship and the bond-slip law of fibres in "unconditioned" and "conditioned specimens".
Based on the results obtained from the aforementioned inverse identification procedure, the following main considerations can be drawn out:

•
The freeze-thaw cycles effect the cylindrical compression strength, f cm , the transition length l t , and the bond-slip law of fibres, which confirms their significance as relevant parameters controlling the resulting response of HPFRC specimens; • Table 3 shows that the compressive strength f cm undergoes a substantial reduction (in the order of 20%) as a result of the degradation processes induced by the FT cycles; • as for the transition zone, which is a peculiar aspect of the considered model, a moderate increase in its the depth (from 70 mm to 85 mm) can be identified after the FT cycles, whereas its shape (controlled by the exponent α) does not change; • Table 4 points out that the bond-slip law of fibres is also affected by the FT cycles, as, specifically, the elastic limit stress, τ el , (and, consequently, the initial elastic stiffness of the same law) reduces by about 15%, with no changes in the other parameters; • however, under the designers' standpoint (and besides the specific values obtained in the present study), it should be noted that this change affects both serviceability and ultimate limit states in the structural response.
Finally, the generally good agreement between the experimental data and the values obtained from the identified model confirms the mechanical consistency of the latter and its potential accuracy. However, further experimental results are needed to calibrate general relationships between the main parameters controlling the bond-slip law of fibres and the actual number of freeze-thaw cycles: this will be part of the future developments of the present research.