Effect of Microstructure and Relative Humidity on Strength and Creep of Gypsum

The wide range of gypsum facies observed all over the world and the strong heterogeneity that may be present even within a single facies often cause an inhomogeneous mechanical response that, if neglected, may be particularly dangerous in the framework of underground excavations. In addition, gypsum is particularly sensible to the presence of water. The high relative humidity conditions often registered in underground gypsum quarries may imply an additional worsening of mechanical properties. In the present study, the strength and the creep response of a natural gypsum rock facies are investigated, considering the influence of material heterogeneity and relative humidity conditions. The heterogeneity of the material, quantified with MIP and SEM analyses, is observed to strongly affect the mechanical response. To this intrinsic mechanical variability, the influence of an external parameter as the relative humidity is observed to generate an additional reduction of material strength and to increase the creep strain rate in the long-term tests. The effect of all these elements in the underground quarry framework is discussed and a constitutive model of these experimental results is provided.


Introduction
The increasing demand of gypsum as raw material in a wide range of industrial fields (e.g. building engineering, clinical applications, agriculture, food industries and production of paints) implies its extensive exploitation by both open pit and underground quarries all over the world. Following the US Geological Survey on Mineral Commodity (February 2020), approximately 140 million tonnes of gypsum were exploited through mining activities worldwide during 2019.
The stability assessment of these world-wide diffuse mining sites requires a detailed knowledge of mechanical response of gypsum as natural rock. Unlike synthetic gypsum, natural gypsum may present heterogeneous features, even within a single rock facies, that complicate the mechanical characterisation. In the common practice of underground excavation, the lack of a specific focus on the relationship between gypsum geological variability and strength parameters often leads to ignore important information for the mechanical characterisation of the ore deposits, in detriment of both safety and productivity. In addition, the frequent presence of water or elevate humidity in the underground drifts induces a deleterious effect on the mechanical properties that, if neglected, may create serious risk scenarios [e.g. roof collapses or pillar failure (Wang et al. 2008;Xia et al. 2019]. Several research groups have focussed on the influence of weathering and water on the physical and mechanical properties of gypsum rocks. Auvray et al. (2004) proposed a study of the ageing of gypsum in an underground mine, describing the degradation of microstructure and physical properties (i.e. density, porosity, velocity of seismic waves). Yilmaz (2010) described the worsening of Unconfined Compression Strength (UCS) and Young modulus of gypsum samples as a consequence of water saturation. Castellanza et al. (2010) reported on the stability assessment of an abandoned underground gypsum quarry, proposing an evaluation of pillars stability based on the water deteriorating effect. Liang et al. (2012) considered the influence of different temperature conditions and salinity concentrations in the water. The description of microscale structures produced during water dissolution may be found in Yu et al. (2016), in Meng et al. (2016) and in Meng et al. (2018). Recently, Zhu et al. (2019) proposed a multiscale investigation of the phenomenon, suggesting that the deterioration of mechanical properties depends on the change in internal microstructure that, with the continuous increase of immersion time, brings to the destroying of connections between the rock particles through the hydrolysing and weakening of crystal bonds at the tip of microcracks.
The interaction between water and sulphate-rich rocks (e.g. anhydritic clays) was also observed to cause severe expansions due to the precipitation of gypsum from the previous dissolution of anhydrite, with serious damage to the infrastructures Ramon et al. 2017;Ramon et al. 2013;Ramon and Alonso 2018;Alonso et al. 2015).
The underground quarries of gypsum are generally exploited by room-and-pillar methods. This method results in pillars supporting the action of a static load for long times (i.e. the life-time of the active quarry exploitation). If the rock is not affected by a rapid unstable failure under high stresses, a time-dependent deformation or a creep damage under the action of long-term stress can however take place. The creep damage could also be delayed behind the active mining exploitation for several years or even decades.
Despite the importance of the long-term behaviour of gypsum, scarce experimental data in natural gypsum rocks are available in the scientific literature. Hoxha et al. (2005Hoxha et al. ( , 2006 proposed an experimental investigation and a constitutive contribution that relates the long-term creep behaviour of gypsum to the relative humidity. Moirat et al. (2006) presented an experimental investigation on the creep response of saccharoid gypsum in saturated conditions. Their results highlighted the achievement of accelerated creep conditions for applied stresses higher than 35% of the uniaxial shortterm strength. Their experimental data were used as a reference for the formulation of a constitutive model in Nedjar and Le Roy (2013), who proposed a viscoelastic model with a damage-like mechanism active from the very first phases of the creep tests.
With the aim of providing an analysis of the most effective parameters on the mechanical response of gypsum rock in underground quarry environments, this study proposes an experimental investigation of short-term strength and creep behaviour of natural gypsum samples, considering the influence of microstructural heterogeneity, relative humidity conditions and orientation of anisotropy planes. The recently described gypsum facies named branching selenite (Lugli et al. 2010) was chosen as testing material for its large presence and industrial exploitation all over the Mediterranean basin. With the aim of considering the strong heterogeneous and anisotropic structure of this rock facies, the present study provides a microstructural characterisation of the material. The effect of the material microstructure and anisotropy orientation on the mechanical response is experimentally tested, also considering the interaction with different relative humidity conditions. The study then proposes a constitutive contribution that models the experimental results. The relative humidity, included in the model as a variable, is a good parameter to describe the variability of natural conditions in underground drifts. The paper closes with some considerations on the scale effect, contextualising the experimental results in a real quarry framework.

Tested Material and Microstructure
From a geological point of view, the tested material is a Miocene-aged gypsum. It precipitated during the Messinian Salinity Crisis, a period of anomalous salinity conditions that involved the entire Mediterranean Basin between 5.9 and 5.3 Millions of year ago, with the consequent precipitation of enormous volumes of gypsum and salt rocks in Italy, Spain, Greece and in the deep Mediterranean Basin (Dela Pierre et al. 2011Pierre et al. , 2016Lozar et al. 2018;Sabino et al. 2020).
In this study, we considered material from the northwest Italy, and precisely from the Monferrato area (Fig. 1a). In the area, gypsum rocks with different facies (e.g. massive selenite, branching selenite, gyps-rudite and gypsarenite) have been exploited for several decades by both open pit and underground quarries, with development of mining tunnels that can reach lengths of tens of kilometres.
Branching selenite (i.e. the specific object of this study) is a gypsum rock facies with a mean porosity of 5-6%, a mean bulk density of 22 kN/m 3 and a mean gypsum content of 90%. It consists of millimetre gypsum crystals organised in sub-horizontal lenses or nodules (a few centimetres large) surrounded by layers of fine-grained material (mainly gypsum, clay minerals, carbonate minerals, quartz and feldspars). The resulting structure shows a clear anisotropy that is underlined by the preferential orientation of nodules and gypsum crystals (Fig. 1b, c). This material heterogeneity is periodically reproduced at the different observation scales: at the rock-mass scale (Fig. 1d), it describes big branches of coarse gypsum wrapped and separated by layers with finer mean grain size (Caselle et al. 2019b).
The branching selenite samples investigated in this study were obtained from two quarry sites in the Monferrato area (Calliano and Moncalvo quarries) and from a survey campaign of 17 drillings in the same evaporitic succession. Core drillings were performed in vertical direction, i.e. perpendicularly to the sub-horizontal stratification and to the main sedimentary discontinuities. In the quarry sites, the material is nowadays exploited in underground environments that are excavated with road header machines, in a drift configuration characterised by room and rib pillars.
All the samples show features that allow for a classification as branching selenite facies. However, the meso-textural appearance of samples is strongly heterogeneous. Figure 2 shows the range of variability, from samples with evident macroporosities and branching structure clearly identified by colour and differential erosion (in the left of the Figure) to compact and more uniform samples, with branching structure only identified by a change in grain size or slight differences in colour (on the right).
Based on meso-textural features, two material groups were defined (denominated, for sake of simplicity, A-type and B-type, Fig. 3). A-type samples are greyish, with a wellcompacted appearance. Their branching selenite structure is macroscopically only identified by the differences in grain size. B-type samples are more yellow-brownish and the layering of the branching selenite structure is thin and well evident for the brownish colour and the differential erosion.
This two material groups do not include the complete range of variability of the rock: samples with intermediate characteristics cannot be satisfactorily described by any of these categories. "A-type" and "B-type" features describe, indeed, the "extreme" characteristics of the rock (end-members) and were used as a trace in the physico-chemical and mechanical characterisation of the material.
A detailed analysis of the textural and compositional features of these two material groups can be found in Caselle et al. (2019a, b). The results, summarised in Table 1, include the evidences of: -A difference in grain-size distribution, with a coarser and highly sorted size distribution in B-type gypsum. To complete this dataset, the present study proposes a microstructural analysis of the material, to observe and quantify the presence, shape and nature of microcracks and pores.  The first technique (MIP) offers the possibility to analyse, through the injection of mercury in the material with increasing pressure steps, the volume of pores that corresponds to each entrance diameter (i.e. throat size), in a range between 3.5 nm and 500 μm (Giesche 2006;Romero and Simms 2008). The second technique (SEM-SE) provides micro-morphological images of the pores themselves. In addition, SEM-EDS (energy dispersive X-ray spectroscopy) analyses were performed pointwise on the samples to associate chemical information to the morphological images.

Mechanical Tests
As already introduced, the mechanical response of branching selenite gypsum rock was tested in short-term and long-term test configurations, to describe the strength and the creep of the material respectively. All the tests were performed in uniaxial compression.
The uniaxial short-term strength tests were performed with a servo-controlled press applying an axial stress with constant strain rate of 10 −5 s −1 . Both A-type and B-type gypsum samples were tested. The former (A-type) were tested as prismatic samples, with mean sizes of 100 mm × 50 mm × 50 mm, to facilitate the comparison with the creep tests (performed on cubic samples). The latter (B-type) were tested as cylinders due to the low quality of the rock and the difficulties in the cutting of prismatic shape.
To evaluate the influence of the anisotropy orientation on the results, a set of prismatic A-type samples was also prepared and tested with layering orientation angles ranging between 30° and 90°. The angles are defined as showed in Fig. 4. Table 2 summarises the test conditions of all the performed uniaxial short-term strength tests.
For the long-term creep tests, a specific set-up was designed, to guarantee the uniformity of stress conditions over long periods of time (i.e. several months) and to optimise the use of laboratory equipment.
The set-up contemplates a lever system that simultaneously applies the load to three cubic gypsum samples with different sizes (Fig. 5). Since the load distributes on the different areal surfaces of the samples, the resulting stress applied to each of them is different. Metallic plates separate the samples and transfer the load from one sample to the others. Two dial gauges for each plate measure the axial displacements and give additional information about the uniformity of the strain across the section of the sample (i.e. absence of plate rotation) (Fig. 5b). Table 3 summarises the areas of the samples and the correspondent applied stresses. Tests were continued until stable steady state conditions were reached, with the measurement of a constant strain rate over time.
Due to the long times required for the tests, only A-type material (more homogeneous and allowing for more repeatable conditions) was tested. The test repeatability of A-type material allowed for the evaluation of the influence of other parameters: considering the relevance of the anisotropy orientation in the quarry assessment, one set of   In the creep test with the most disadvantageous conditions (i.e. 10 MPa and RH = 90%), a special precaution was taken in view of the possible acceleration and failure of the sample. A metal support was prepared to bear the axial load in case of collapse of the upper sample (Fig. 6). The height of this metal support is 10% lower than the sample. If, during the test, the sample reaches an axial strain of 10% (i.e. creep failure), the higher plate will touch the metal support and the test will continue on the two lower samples.

Control of Relative Humidity
Short-term strength tests and long-term creep tests were performed both in free environment and with relative humidity control.
In the first test modality, samples are free to equilibrate with the relative humidity of the laboratory (i.e. Laboratory of Soil and Rock Mechanics of the Department of Geotechnical Engineering and Geosciences of UPC, Barcelona and Geotechnical Laboratory of the Department of Structural, Geotechnical and Building Engineering of Polytechnic, Turin). In both laboratories the environmental humidity is maintained constant with the help of a system of conditioned air. The measurement of humidity over a period of 6 months in Turin laboratory returned values ranging between 38-40% and 50-52%. In Barcelona laboratory, despite a shorter period of monitoring, a range of humidity between 38-40% and 55-57% was measured. A mean value of relative humidity of 45-50% can be therefore attributed to the free tests, but, for the long-term conditions, no assurances of humidity stability is given.
In the second mode, samples were tested within a close environment (i.e. plastic bags), with controlled relative humidity imposed with specific water-salt solutions. Figure 7 shows the circuit used to create controlled humidity conditions in short-term and creep tests, respectively.
Silicone pipes connect the plastic bags with a close recipient with a water solution oversaturated with a specific salt. At the equilibrium, the salt in solution controls the relative humidity of the environment in contact with the water free surface. Distilled water, with a nominal relative humidity of 100%, was chosen to create conditions of high relative humidity, while an oversaturated lithium chlorite solution, with a nominal relative humidity of 11%, was used to impose the low humidity conditions.
The connection of the recipient with an air pump assures the circulation of the atmosphere with the desired relative humidity and the maintenance of equilibrium conditions between the different parts of the circuit. Depending on the type of pump, the exiting air from the bags directly connects with the air pump itself, closing the circuit (Fig. 7c), or is directed to a new recipient with the same water solution, to avoid the injection of uncontrolled air (Fig. 7a).
Before the starting of the test, the circuit was kept active for a minimum of 3-4 days to reach the equilibrium and allow the humidity to permeate the gypsum samples.
Relative humidity was monitored during the test with probe hygrometers. Considering some air dispersion in the system, regular values of relative humidity of 90 and 25% were obtained in high and low humidity tests respectively. The long permanence of the high humidity conditions in a close environment brought, during the creep tests with 90% of relative humidity, to the precipitation of water drops in the bags. At the end of the test, electrical resistivity of this water precipitation was measured. The value (1.8 mS/cm) is coherent with data about gypsum-saturated water solutions, testifying the occurrence of a process of gypsum dissolution in the plastic bags by the precipitated water.

Microstructural Analysis
The results of the MIP tests are summarised in Fig. 8a and can be described as follows: -< 0.1 μm-no pores were detected neither for sample A and for sample B -0.1-1 μm-both sample A and sample B show an increasing concentration of porosity towards coarser pore sizes, even if in sample B the growth is not uniform, but characterised by a multi-modal distribution. -1-6 μm-high pore concentration for both sample A and sample B. Sample A has its absolute maximum between 1 and 2 μm. Again, sample B pore distribution is multimodal. -6-20 μm-maximum concentration of porosity for sample B, while sample A has almost zero porosity in this size class. -20-100 μm-neither sample A and sample B have significant pore concentrations in this size range -> 100 μm-both sample A and sample B show an increase of measured porosity in this interval, suggesting the presence of big pores, maybe even bigger then higher resolution of the method (500 μm).
The principal differences are, therefore, the porosity peak between 6 and 20 μm in sample B, which does not have an analogous in sample A, and the multi-modal pore-size distribution of sample B.
SEM morphological micrographs in Fig. 8b, c suggest that the first pore-size class (smaller than 1 μm) should correspond to the finest grain-size class of the rock (smaller than 3-4 μm).
This grain-size class mainly consists of non-gypsum minerals. EDS spectra (e.g. Fig. 9b) indicate the predominant presence of Ca and Mg, with the additional presence of K, Si and Al. EDS analysis on these very fine materials cannot be reliably associated with a single particle, but returns a general composition. These materials can be, therefore, considered as a dolomitic plus terrigenous cement of the rock.
An image threshold filter was applied on the micrographs to quantify the size of the pores (Fig. 9a-c). The resulting frequency distribution (Fig. 9d) shows that the mean pore size is compatible with the first pore-size class in the MIP results (0.1-1 μm).
The features of this dolomitic and terrigenous matrix are similar for the A-type and the B-type samples. The grain and pore sizes, in general, are similar, even if slightly coarser in sample B.
The second pore-size range recognised in Fig. 8a (1-6 μm) can be associated with the void spaces at the contact between gypsum crystals. This kind of porosity can be found in both gypsum typologies.
The following pore-size class (6-20 μm) can be explained considering that gypsum crystals in B-type samples are not always well adherent and open spaces and cracks are visible along the grain boundaries. The opening of these contacts was measured on the micrographs (e.g. Fig. 10), returning a range between 10 and 20 μm.
Eventually, big pores (larger than 100 μm) were recognised in SEM images of both samples, even if their presence is more frequent and evident in B-type gypsum. Some of these pores show the presence of groups of acicular crystals of sulphate minerals, probably connected with the circulation of sulphate-oversaturated fluids through the material (Fig. 11).  Table 1, suggest the existence of a primary variability of the material, mainly in terms of texture and gypsum content (i.e. the consequence of the branched structure of the deposit) that is enhanced by a circulation of fluids. The coarser materials, with a higher primary porosity, facilitate a preferential permeation of the circulating fluids. The water, due to the long permanence in the pores, starts a process of gypsum solution, enlarging the pre-existing pores and dissolving the material along the crystal boundaries, creating a rock class with frequent microcracks. Figure 12a, b show the results of uniaxial short-term tests with different humidity conditions in A-type and B-type gypsum, respectively. The comparison between Fig. 12a and Fig. 12b highlights the clear distinction among A-type and B-type samples in terms of both strength and stiffness. Considering the free environment humidity curves (i.e. RH = 50%), a strength decrease from the 18-20 MPa of A-type samples to the 6-8 MPa of the B-type samples can be observed.

Effect of Relative Humidity
The post-peak curves of A-type samples appear to be characterised by a series of instable stress drops. According to Brace and Byerlee (1966) and Brantut et al. (2011), these stress drops can be considered as a release of small fractions of the stress supported by the sample in consequence of the development of a microcrack or a small rock-bridge failure, in analogy with the seismic slip at the scale of natural earthquakes. This behaviour testifies an instable step-wise crack coalescence process, as suggested by the analysis proposed in Caselle et al. (2020). This kind of evidences are not observable, however, in B-type samples that present uniform and regular post-peak curves. Since, as suggested by the microstructural analyses, the material properties of B-type gypsum were affected by circulation of fluids in the pores, that created a rock class with frequent microcracks, the systematic absence of stress drops in this material suggests that the post-peak evolution of the failure occurs through the coalescence of pre-existing microcracks instead that through the creation of new cracks. Figure 13 summarises the values of strength and the correspondent relative humidity. The graph allows for the description of the effect of humidity on the material Fig. 8 a Results of the MIP tests on representative samples of A-type and B-type gypsum. b, c SEM micrographs of the pores in the interval between 0.1 and 1 μm for A-type and B-type samples, respectively. d, e Contact between the grains in sample A and sample B, respectively. f, g Cracks and opening along the grain contacts in B-type sample. h, i Examples of pores bigger than 100 μm in sample A and sample B, respectively response within each gypsum class. In the case of A-type gypsum, the uniaxial strength of samples in free environment (i.e. RH = 50%) and under low relative humidity (i.e. RH = 25%) are similar or even slight higher in the former. The results show, on the other hand, a clear decrease in strength and elastic modulus in samples with high humidity conditions (i.e. RH = 90%). The presence of a humid environment causes a decrease in peak strength and an increase in axial strain in correspondence of the peak and the post-peak phase, describing a longer and less sharp softening.
The high heterogeneity of B-type samples makes difficult to recognise and isolate the influence of an external parameter as the relative humidity. However, despite some outliers, the trend observed for humid A-type samples (i.e. decrease of strength and stiffness with the increase of relative humidity) can be re-proposed for these samples.
The percentage loss of strength generated by the increase in relative humidity is significantly higher in B-type material than in A-type material. This is reasonable, considering the initial higher porosity and presence of microcracks that favours the permeation of water molecules, increasing the rock specific surfaces in contact with the water.

Effect of Anisotropy Orientation
The results of uniaxial short-term tests on A-type gypsum samples with different orientations of the layering are summarised in Fig. 14a. Results are divided in horizontal (0°), vertical (80°-90°) and intermediate (30°-60°) angles. As clearly visible in Fig. 14b, the trend of the uniaxial strength with the anisotropy orientation sees the highest values in correspondence of horizontal and vertical anisotropy planes, with a strength decrease with intermediate angles. The most Fig. 9 a Micrograph representing the finest matrix of the rock. b EDS spectrum, showing the presence of dolomite and clay minerals. c Pore distribution obtained with an image processing of a. d Frequency distribution of pores in c significant variation, however, is the change in behaviour in the stress-strain curves. As already stated, samples with horizontal anisotropy show a long post-peak softening phase with the presence of little stress drops. On the other hand, samples with vertical anisotropy reach the collapse just after the peak. Moreover, the Young Modulus is significantly higher, with lower axial strain in correspondence of the peak point. An intermediate behaviour is manifested by the samples with intermediate anisotropy angles: their curves show a soft deformation until the strength peak and in the initial part of the post-peak. This is however suddenly interrupted with an abrupt loss of strength, in some way corresponding with the first stress drop in the samples with horizontal anisotropy. This behaviour may be considered as the consequence of the parallel orientation between the pre-existing weakness surfaces and the final failure surface. In consequence of this, when the first important microcrack develops, a displacement step of 0.01 mm (i.e. the displacement of the piston for each second considering a mean sample length of 100 mm and the imposed strain-rate of 10 -5 s −1 ) is enough to bring to the complete coalescence of the failure and to the collapse of the sample.

Effect of Relative Humidity
Figure 15a-c show the results of creep tests under high, environmental and low humidity conditions, respectively. Starting point of the curves is after one hour from the application of the load, to nullify the contributions due to the instantaneous deformations. The results correspond to samples with horizontal anisotropy planes (anisotropy angle = 0º).
The difference between Fig. 15a and Fig. 15b, c testifies the strong influence of relative humidity on the results. With high relative humidity (Fig. 15a), the material creeps faster, with a higher strain rate recognisable from the very first data. The curves are regular and well defined, with clear definition of transient and steady state phases of the creep. A clear dependence of the behaviour from the applied stress is observable: samples loaded with 2 and 3 MPa do not show big differences in the measured strain, but samples under 4 MPa and, particularly, 10 MPa clearly show an increase in measured strain and strain rate. However, when the test was interrupted, no evidence of test acceleration was already observed for any of the samples.
Samples deformed without humidity control or under low relative humidity (Fig. 15b, c) show a more irregular trend, with highly noisy strain records. Even if an initial transient phase is maintained, the curves briefly reach a condition of absence of creep, i.e. a strain rate too low to be quantifiable. In free-environmental conditions (Fig. 15b), the curves do not identify a clear relationship between applied stress and final strains and strain rates. Samples under 6 and 10 MPa have similar behaviour to sample under 2 MPa. Without the influence of the humidity, an increase in axial stress seems to not necessarily imply an increase in strain rate. The influence of other parameters (e.g. sample variability, temperature or scale effect that causes an increase in stiffness with the decrease of sample sizes) is probably more effective than the applied stress by itself. It has to be considered, however, that samples under 6 and 10 MPa were tested in Turin laboratory, while the other three samples were tested in Barcelona laboratory. A difference in the environmental humidity (about 40% in Turin and about 60% in Barcelona) can be imagined as a partial explanation of this apparent incompatibility among the results.
Tests performed under low relative humidity (Fig. 15c) registered a further decrease in strain rate. In the long-term, the graph highlights a substantial absence of strain. After 190 days, the relative humidity was manually increased (up to RH = 60%) by changing the water solution in the circuit. As can be seen, this change in the environmental conditions generates a sudden change in the strain path, with an increase in strain rate, at least for the samples with 3 and 4.5 MPa of applied stress.
The gradual decrease and stabilization of strain rate during the evolution of the creep tests is well described by the graph in Fig. 16. In the Figure, the strain rates of samples with RH = 90% and free-environmental humidity conditions are reported against time in a semi-logarithmic plot. The Figure confirms the difference in behaviour as a consequence of the changes in relative humidity and the high noisiness of data recorded in the absence of humidity control.
The decreasing strain-rate trend of the curves in Fig. 16 can be approximated with straight lines in a time-logarithmic scale (Fig. 17a-c). In Fig. 17d, the angular coefficients relating the creep strain and the logarithm of time (defined as tan(λ)) are plotted as a function of the applied stress. The result clearly shows the dependence of creep behaviour on the increase of relative humidity and applied stress. Again, the anomalous trend of the two samples tested in Torino laboratory is evident. Figure 18a shows the strain-time curves of creep tests on samples with anisotropy planes oriented at 45°. The curves of the tests with equivalent conditions (i.e. stresses = 2, 3, 4 MPa and RH = 90%) but horizontal anisotropy are reported on the graph to facilitate the comparison. Results show very low levels of strain and strain rate in samples deformed with applied stresses of 2 and 3 MPa (i.e. largely lower than in the equivalent samples with horizontal anisotropy). On the other hand, sample with the highest applied stress (4 MPa) starts with similar low strain rates, but, after 24 days, suddenly accelerates, reaching the highest strain rates in the dataset.

Effect of Anisotropy Orientation
The anomalous behaviour of this sample may be considered in analogy with the behaviour described in short-term  tests for samples with inclined layering. In both the cases, we observed a sudden change in the material response, with an acceleration of the deformation. Hence, the creep acceleration could actually correspond to an acceleration of viscous deformation due to the parallelism between strain concentrations and material layering. However, it was observed only in one sample, therefore no assurance of repeatability of this result can be given. Figure 18b, showing the strain rate-time data in semilogarithmic scale, confirms the anomalous behaviour of this sample. The initial trend of decreasing strain rate is drastically interrupted by a sudden acceleration at time of about 2 × 10 6 s (~ 23 days), that brings to strain-rate values of 1 × 10 -9 s −1 . A new phase of strain-rate deceleration, however, brings to the new assessment of the general mean values of the dataset (between 1 × 10 -10 and 1 × 10 -11 s −1 ). Nedjar and Le Roy (2013) proposed a constitutive model for the uniaxial short-term strength and long-term creep tests in gypsum rock. The model was created and calibrated on the experimental results reported in Moirat et al. (2006) (i.e. creep tests on saccharoid gypsum samples under uniaxial load and water-saturated conditions). It offers an accurate simulation of creep behaviour of gypsum. However, the formulation is calibrated on saturated conditions, not including the relative humidity (or the saturation) as a variable. To provide a constitutive contribution able to include the variation of humidity, we propose an integration into the Nedjar and Le Roy (2013) model.

Constitutive model
The original model, consisting on the parameters reported in Table 4, can be described, in a one dimension simplified version (used for the here-presented simulations), with the following formulation: where σ and ε are axial stress and axial strain, E is the Young modulus and D and ε v are two internal variables describing the damage and the viscous part of the strain respectively. Since the model relies on the definition of these two internal variables, a brief description of the elements that compose them is provided hereafter.
The first internal variable (damage D) of Eq. (1) is a scalar between 0 and 1 (where 0 is undamaged and 1 completely damaged). Its formulation remains unmodified in the shortterm and long-term behaviour of the material. This means that, once the parameters are defined for the instantaneous behaviour, they can be directly used for the simulation of creep. The definition of damage is controlled by three elements: 1. the function S(ε), describing the source of the damage, that depends on the history of total strains. It is defined by the parameter ζ 1 , following where ⟨ n ⟩ + is the extensional contribution of the n component of the strain tensor. 2. the 'damage exponent' m, that depends on three subparameters (m 1 , m 2 and m 3 ) and is written as a damagedependent variable in the form:  In the short-term uniaxial test curves, the first point where the equality (1 − D) m S( ) = W is verified corresponds to the elastic onset, when the material behaviour departs from the linear elastic curve. For a fixed value of ζ 1 , the parameter W can be therefore determined using the strain value measured in the onset point, where the damage is still 0 and the equality can be written as S( ) = W.
The second internal variable (viscous strain ε v ) is defined as the sum of a number l of contributions following: Each i-th element v i includes an elastic and a viscous contributions, as illustrated in the rheological model (Fig. 19). The elastic and the viscous parts of the i-th elements are defined by a stiffness parameter (E i ) and a viscosity parameter (η i ) respectively. These material parameters appear in the model in terms of relaxation time (τ i = η i /E i ) and stiffness factor (ω i = E/E i , where E is the Young modulus). To simulate the test progress, the increment of the internal variable ε v for a small time interval Δt is iteratively calculated with the following equation: where v i corresponds to each element of ε v (for i = 1,…, l) and n and (n + 1) represent two consecutive time instants, where t n+1 − t n = Δt. Following Eq. (5), the value of ε v in each time step is then calculated as the sum of the l contributions v i . In our simulation, the number l of contributions was fixed to 2. The viscous strain ε v , as well as the relaxation time (τ i ) and the stiffness factor (ω i ), are, therefore, vectors of two elements.
To introduce the humidity variable in the model, a sensitivity analysis was first performed to evaluate the influence of each parameter on the modelled curves. The graphs in Figs. 20, 21 and 22 show the results of a set of simulations performed using as reference the modelled curves obtained with the set of parameters reported in Table 4. These parameters were calibrated using as reference the stress-strain curves of the uniaxial tests with free-environmental relative humidity (with a peak strength of 20 MPa and a correspondent strain of 0.01- Fig. 12a). For the creep test, a constant stress of 10 MPa was considered. The graphs, obtained by varying the 8 parameters, one by one, in a certain range, describe the specific influence of each of them on the final curves.   The Young modulus E (Fig. 20a, b) is the most effective parameter on both short-and long-term curves. An increase in E tends to increase the strength and decrease the creep propensity of the material.
The damage threshold W and the source of damage parameter ζ 1 (Figs. 20c-f) have similar effects, even if inverted. They are, indeed, both involved in the inequality that controls the damage activation ( (1 − D) m S( ) ≤ W ). Both parameters affect the UCS curves modifying the position of the onset point and, consequently, the final strength and strain. They also have a visible effect on the creep deformation. It is, however, possible to identify a 'minimum creep curve' corresponding to the absence of damage (i.e. the strains are, for the entire duration of the creep test, too low to overcome of the damage threshold). Below that curve, the effect of all the damage-related parameters is null.
The parameters related with the definition of the m-exponent (m 1 -m 2 -m 3 - Fig. 21) are only effective in the shortterm regime, changing the shape of the curve. They are, however, ineffective in the long-term regime in the considered stress and strain conditions, since the damage threshold is not overcame during the creep test. Their effect, indeed, being related to the evolution of damage, is subordinate to the parameters of damage threshold: if the threshold is not exceeded, their effect is null.
Since the viscous strain-related parameters (τ and ω) only influence the time-dependent behaviour by definition, the sensitivity analysis on the short-term curves is not reported. Figure 22 shows, on the other hand, the influence of an increase in relaxation time (τ) and stiffness factor (ω) on the creep curves. As can be seen, the relaxation time has the sole effect of changing the shape of the curve (i.e. how fast the transient phase is concluded). The influence of the stiffness factor is more effective, inducing big changes in the measured strains.
Based on these information, the effect of the relative humidity on the experimental results was described by working on the Young modulus and the damage threshold W (i.e. the two most effective parameters on both short-term and long-term curves). Due to the low sensitiveness of shortterm behaviour to the changes in relative humidity below 50%, we assumed the existence of a threshold value of relative humidity (that we fixed at RH = 50%, for the absence of further constraints) that triggers the change in the damage mechanism. Figure 23a, b propose the best-fitting curves for the experimental results of short-term tests. Starting from the parameters reported in Table 4, the Young modulus E and the damage threshold W were modified until the best fitting of the all the experimental data was obtained. The final values of E and W are reported in Fig. 23c, d, as a function of the relative humidity, obtaining two formulations (for A-type and B-type gypsum respectively) that quantify this humidity-related change of damage mechanism. The equations reported in Fig. 23c, d describe, therefore, how the Young modulus and the damage threshold change if the relative humidity exceeds the threshold of 50%.
Unlike the short-term strength tests, the creep tests registered a certain sensitivity also to the changes in relative humidity below 50%, requiring the introduction of an additional mechanism. Based on the sensitivity analysis in Fig. 22, the more effective stiffness factor (ω) was selected to simulate this mechanism.
The curve fitting reported in Fig. 24 follows the relationships proposed in Fig. 23. The good fitting of creep data using the new proposed parameters for the short-term compression curves assures a validation of the modified model. The additional modification of stiffness factor (ω) was adjusted on the basis of to obtain the best possible fitting. In Fig. 24b (i.e. creep curves in free-environmental humidity) values of RH = 40% and RH = 60% were considered for the tests performed in Turin and in Barcelona laboratory respectively. The final result describes well experimental data, even if in Fig. 24b some of the curves (e.g. the creep test under 10 MPa) are completely out of trend. The data fitting in Fig. 24c (i.e. tests with low relative humidity) well describes the initial part of the test, while the fitting decreases in the long term. This is mainly due to the impossibility to reproduce, with the model, the data noise due to the temperature variations in the room. The best fitting values of stiffness factor are reported in Fig. 25 as a function of the relative humidity, allowing for the definition of a logarithmic dependence equation.

Discussions and Conclusions
The heterogeneity of mechanical and creep response of natural gypsum rock as a function of microstructure and relative humidity plays an important role in the stability assessment of mining sites. The presented laboratory outcomes find,  Table 4 indeed, a specific application and a practical meaning in the framework of underground excavations, where the heterogeneity of the material and the humidity conditions may create unexpected risk scenarios.

Microstructural Heterogeneity and Layering Orientation
The presented results in terms of textural, compositional and microstructural characterisation suggest that branching selenite gypsum has a primary variability that is a consequence of the macroscale structure (i.e. big "branches" of coarse gypsum immersed into the finer-grained gypsiferous matrix). The mesoscale features of the rock (i.e. the size distribution of the crystals, the porosity and the mineralogical composition) depend on the position of the sample in this structure (e.g. within the branches of coarser gypsum or in the surrounding material, with higher or lower concentration of coarser lenses). In addition to this primary heterogeneity, the microstructural analyses showed evidences of water dissolution processes. The data suggest, indeed, that the long permanence of water in portions of the rock-mass have been producing changes in the microstructure, opening microcracks, big pores and structures of detachment along the grain boundaries. Figure 26 proposes a schematisation of a hypothetical branching selenite gypsum layer that shows the distribution in the rock-mass of these heterogeneities (both primary and water-induced). The internal organisation of the gypsum layer is controlled by the presence of the elongated gypsum branches that influences both material heterogeneity and layering orientation. The scheme reports four tunnels in different positions that exemplify some typical situations.
Tunnel A is in a peripheral position with respect to the nucleation point of the branches. In this position, the gypsum lenses are less abundant and, hence, the rock is predominantly characterised by the fine matrix, creating a situation similar to Fig. 27a. Tunnel B, on the other hand, is nearer to the branches origin. In this case, the rock is coarser and has higher gypsum content. The rock appearance could be, therefore, more similar to Fig. 27b. In addition to this primary variability, Fig. 26 exemplifies the possibility of a local concentration of water, due, as instance, to karst processes (tunnel D). In these conditions, a microstructural worsening of rock properties has to be expected.
The overall heterogeneity of the material was schematised, in this study, with the definition of two end-members (A-type and B-type gypsum, describing respectively a good quality and a poorer material). The results of short-term strength tests returned strength ranges of 18-20 MPa and 6-8 MPa for A-type and B-type gypsum, respectively.
In the framework of underground excavations, the consideration of the material heterogeneity assumes a fundamental role for the stability assessment. In the drift portions where B-type material is prevalent (e.g. in the presence of karst phenomena, as in tunnel D in Fig. 26), a mean strength value, representative of the entire rock body, would be largely overestimated, creating serious risk scenarios. On the other hand, the size of underground pillars mainly composed by good quality material (A-type gypsum) may be overestimated if a mean value of strength is considered.
The schematisation of the rock body proposed in Fig. 26 exemplifies another potential risk factor: the  Table 4 change of orientation of the anisotropic layering. As described in the Figure, despite the mean horizontal orientation of the gypsum branches in tunnels A and B, the structure defines a moderate to very inclined orientation of the layering in correspondence of the nucleation point (tunnel C). As shown by the experimental results, the orientation of mesoscale anisotropy (layering) introduces an additional element in the control of the mechanical response. Intermediate angles of anisotropy, with weak surfaces orientated almost parallel to the final angle of failure, facilitate the crack coalescence, reducing the bulk strength of the material. These unfavourable layering orientations were also observed to generate unexpected strains and strain rates in the long-term creep deformation, with a not easily predictable behaviour.

Effect of Relative Humidity
The experimental results presented in this paper show the evident worsening of mechanical response of gypsum rock in presence of high humidity conditions. The mechanical tests underlined a clear decrease in short-term strength and stiffness in presence of high relative humidity (i.e. RH = 90%). The water worsening effect was observed to be particularly deleterious on intrinsic low-quality material (i.e. B-type gypsum), creating very disadvantageous mechanical conditions, with peak strength values of only 2 MPa.

Fig. 23
Experimental and modelled curves for A-type (a) and B-type (b) gypsum. c Dependence of the Young modulus on the relative humidity. d Dependence of the damage threshold on the relative humidity High relative humidity strongly influences also the creep response of the material (only tested for A-type samples). With relative humidity of 90%, the creep tests registered strain rates of one order of magnitude higher than the equivalents in free-environmental humidity.
In underground drifts of gypsum quarries, the atmosphere is often very humid even when liquid water is not directly present. Considering the stress configuration of the pillars, the quantification of expected creep deformation is fundamental for the stability assessment. After the end of active exploitation, the conditions of long-term stability may be further compromised by the interruption of the pumping operations, involved in the dewatering of gypsum orebody during the exploitation. The consequent restoring of groundwater level may bring to the new saturation of the rock mass, creating the environmental conditions that maximise the time-dependent deformations in gypsum. Furthermore, the flow of fresh water, potentially under-saturated in sulphates, could trigger important dissolution processes on the free surfaces of the pillars, worsening the mechanical quality of the material. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 26
Schematic representation of the branching selenite structure. The scheme includes 4 hypothetical tunnels in different positions of the structure. Tunnel A is peripheral with respect to the origin of the branches. The rock in that point will be finer and with lower gypsum content. Tunnel B is narrower to the branches origins, with a con-sequent coarser grain size and higher gypsum content in the rock. Tunnel C is in correspondence of the origin of the branches. In this position, the orientation of the layering is inclined. Tunnel D is in correspondence of a karst area, with water-induced modification of the rock microstructure