Applications of a micro-structured brittle damage model to laboratory tests on rocks

A multiscale microstructured brittle damage model is used to describe the behavior of confined rock materials. Plane strain and triaxial tests conducted at the laboratory scale are simulated in terms of boundary value problems. Simulations reveal good predictive qualities of the model to describe the macroscopic features of specimens at failure. The microstructures, oriented in different directions, allow the localization of the macroscopic strain along straight lines, emerging at the macroscale in the form of shear bands. The microstructured material model, characterized by recursive equidistant parallel cohesive-frictional faults, is fully defined by six elastic and inelastic material constants. The model was originally developed in a finite kinematics framework to simulate the dynamic behavior of confined brittle materials (Pandolfi et al. in J Mech Phys Solids 54:1972–2003, 2006). In linearized form, it has been extended and used for the simulation of in-field excavations (De Bellis et al. in: Eng Geol 215:10–24, 2016). The performance of the model in predicting the behavior of small scale rock tests in the laboratory, the object of the present study, has never been investigated. Numerical simulations show that the model is able to capture several important features observed in rocks, in particular the reduction of the overall stiffness for increasing deterioration of the material, fragile to ductile transition, strain localization, shear band formation, and more general size-effect.


Introduction
Rocks are natural, inorganic solids, composed of two or more minerals, that over geologic time undergo severe transformation processes involved with the rock cycle (ranging from igneous, sedimentary to metamorphic rocks). Depending on the mechanical and environmental boundary conditions evolving in time, rocks can be intact or characterized by the presence of discontinuities (fractures, joints, cracks, bedding or foliation planes) caused by the combined action of pressure, temperature, strain localization, and fluid flow. Rock discontinuities can be induced by geological processes or human activities, i.e. drilling, hydraulic fracturing, tunnelling and/or core handling. In either case, the resulting fractured rock masses often exhibit sets of nearparallel, near-equally spaced fractures that are responsible for extremely complex behaviours, whose study is made more challenging by the coupling between mechanics and fluid flows. The onset and development of discontinuities (that evolve according to, e. g., an accumulation process) determine substantial changes both in terms of stiffness and strength degradation (Amitrano 2006) and in terms of porosity and permeability changes. An ensuing feature of natural rocks is the marked mechanical anisotropy associated with the evolution of the microstructure, responsible for the direct effects at the macroscopic scale (Hoek 1983;Shao et al. 2005).
Both brittle and ductile behaviours are observed in natural rock embankments (Paterson and Wong 2005).
In laboratory experiments, the transition from brittle to ductile in the same rock may be obtained by increasing the level of confinement (Horii and Nemat-Nasser 1986) or by inducing with pre-cracks the onset of localized mechanisms that develop within narrow zones (Nikolic et al. 2015).
All these behaviors suggest the definition of advanced numerical models, able to retain memory of the evolving microstructure and to predict the overall non-linear response. In this respect, conventional mechanical approaches, where the problem is modelled at the structural scale adopting phenomenological models for the material, may fail to be predictive under different loading conditions, see, e. g., classical damage and plasticity models (De Borst et al. 2012). Phenomenological approaches are unable to describe the exact mechanisms of degradation, since they are based on strong simplifying assumptions. Thus, they capture the correct mechanical response only with reference to specific loading conditions (Bazant 1986;Cho et al. 1991;di Prisco et al. 2003).
On the other hand, fully micro-mechanical models (Jing and Hudson 2002;Zhu and Shao 2017), possibly relying on discrete (Lisjak and Grasselli 2014;Ingraffea and Heuze 1980), continuous (Pande et al. 1979;Riahi et al. 2010) or hybrid (Munjiza 2004) approaches, provide an accurate description of mechanical and geometrical features of the microstructure and are able to predict its evolution, but, due to often prohibitively high computational costs, their use is, in general, limited to very small domains, of the order of the size of the microstructure.
An alternative approach consists in making use of multiscale models, which represent a fair compromise between the contradictory requirements of describing the real complexity of the microstructure and of keeping acceptable computational costs. Among different strategies successfully applied to rock mechan-ics (Abou-Chakra Guery et al. 2009;Guo and Zhao 2014;Zeng et al. 2014;Guo and Zhao 2016;Liu et al. 2016;Wang and Sun 2018), here we adopt a multiscale model of brittle damage, based on the explicit description of damaging microstructures consisting of recursively nested parallel cohesive-frictional faults embedded in an overall elastic matrix. The brittle damage model has been developed originally within a finite kinematics framework (Pandolfi et al. 2006), and it has been applied to grasp the transition from brittle to ductile behaviors of confined brittle ceramics under dynamic loading. Numerical analyses demonstrated that the model is particularly accurate in simulating confined deviatoric stress states, which is suggested to address applications concerning soil and rock embankments. An obvious extension to confined poro-mechanical problems was proposed later to model hydraulic fracking phenomena in low permeability shale gas and oil reservoirs (De Bellis et al. 2017;Caramiello et al. 2018). A key advantage of the model is the possibility of evaluating, analytically, the changes in permeability, directly related to the onset and evolution of interconnected nested parallel faults. The linearized version of the coupled model, presented in De Bellis et al. (2016), was indeed used at the material point level to reproduce hydro-mechanical triaxial tests, considering different rocks available in the literature, and, in a finite element implementation, to simulate large scale field problems, e.g., deep excavations in dry soils.
In this paper, we investigate the performance of the linearized version of the model in reproducing laboratory experiments involving rocks in a plane or triaxial strain state. We apply the model to plane strain and triaxial tests on hexahedrons made of brittle geological materials (Yumlu and Ozbay 1995a). We demonstrate that the model is particularly efficient in capturing the global response in terms of load versus axial displacement and in reproducing the failure patterns of confined dry brittle materials under severe loading conditions, albeit it is fed with a reduced number of material parameters (six, two of which are elastic) and a rather coarse discretization have been used.
The paper is organized as follows. In Sect. 2 we briefly recall the mechanical features of the dry model. In Sect. 3 we present numerical examples of plane strain and triaxial loading of rock specimens and show the performance of the brittle damage model. In Sect. 4 we highlight the advantages in using the multi-scale brittle damage model and indicate possible applications in geo-mechanics.

Mechanical features of the brittle damage model
The linearized brittle damage model adopted here is presented first in De Bellis et al. (2016) in an extended version addressing porous materials. The dry model here considered is characterized by planar microstructures, consisting of nested families of equi-spaced frictional-cohesive faults bounding an elastic matrix, see Fig. 1a. In the case of a single family of faults, with orientation N and spacing L, a displacement jump Δ applied to each fault of the family leads to the discontinuous or singular deformation component where u f is the average displacement due to the fault activity. If, additionally, the matrix undergoes a uniform deformation ε m = sym∇u m , we assume that the total deformation decomposes additively into matrix and discontinuous components At the early stages of damage following the inception before the complete decohesion, faults are assumed to show a cohesive behavior. The residual inelastic bridging between the flanks of the crack is expressed through a cohesive law T = T (Δ) that provides the resisting tractions. The underlying thermodynamic framework postulates the existence of a cohesive energy density per unit surface, Φ(Δ, q), defined in terms of Δ and dependent on an appropriate set of internal variables q describing the current state of the faults. A simple class of mixed-mode scalar cohesive laws Φ = Φ(Δ, q) is obtained by introducing an effective opening displacement Δ defined as Ortiz and Pandolfi (1999) where the material constant β assigns different weights to the normal and tangential components of Δ and |Δ| is the magnitude of Δ, see Fig. 1b. In a cohesive-frictional material, β corresponds to the friction coefficient μ = tan φ, where φ is the friction angle (De Bellis et al. 2016). In scalar terms, the effective cohesive traction follows as and the cohesive tractions derive as Pandolfi et al. (2006) The irreversible nature of damage is enforced by assuming unloading to the origin, which allows to introduce only one internal variable storing the maximum attained effective opening displacement q. The kinetic equation necessary for the evolution of q iṡ q = Δ , if Δ = q andΔ ≥ 0 first loading 0, unloading/reloading .
A typical cohesive law is the linear monotonic loading cohesive envelope combined with unloading to the origin, see Fig. 2 (Ortiz and Pandolfi 1999). The cohesive envelope is a two parameter law, characterized by the cohesive strength T c and the critical energy release rate, G c , represented by the area enclosed by the envelope. Thus the cohesive model is characterized by the three parameters T c , G c and μ.
Upon the attainment of a critical opening displacement Δ c = 2G c /T c , faults lose cohesion completely and cohesive tractions vanish, thus friction associated with fault sliding may remain the sole dissipation mechanism for the material. According to the classical failure criteria that generally apply to geomaterials, friction operates at the fault level concurrently with cohesion during closure, when the normal opening displacement  Δ N = Δ · N is null and only the sliding displacement Δ S = Δ − Δ N N may be nonzero. Following the variational approach described in Pandolfi et al. (2002), friction on the faults is accounted for by introducing a dual kinetic potential per unit area Ψ * , and so the frictional forces can be derived as

∂Δ .
A convenient rate-independent kinetic potential for Coulomb friction can be taken of the form Pandolfi et al. (2006) where |Δ| is the norm of the displacement jump rate and c the contact tractions, leading to frictional forces tangential to the fault plane. The matrix is assumed to be linear elastic and isotropic, characterized by the young modulus E and the Poisson coefficient ν.

Variational update
By time discretization, it is possible to introduce a functional that represents the incremental work-ofdeformation, useful to characterize variationally the behavior of irreversible and frictional materials at every time step (Pandolfi et al. 2006;De Bellis et al. 2016). An incremental process requires the evaluation of the state variables at discrete times t 0 , . . ., t n+1 = t n + Δt. At the time t n the state of the material is assumed to be fully known from the previous steps. At the next time t n+1 , the total deformation ε n+1 is assigned and the new state of the material, given in particular by the variables Δ n+1 and q n+1 , is sought. Note that, for the sake of readability, in the following equations the subscript n + 1 for the current variables is omitted.
Following the approach in Pandolfi et al. (2006), we define an extended incremental work-of-deformation E n (ε m , Δ, q), constrained with the impenetrability of the closing faults and the irreversibility of the cohesive damage. The functional is defined over the time interval Δt as The index n in E n emphasizes the dependence on the initial conditions. E n splits additively into three terms: (i) the strain energy density per unit of volume of the matrix, W m (ε m ), depends only on the strain in the matrix and it is assumed to be convex and positivedefinite; (ii) the cohesive energy density per unit of surface of the faults, Φ(Δ, q); (iii) the frictional dissipation expended in the time interval Δt, expressed through the time discretized potential Ψ * (Δ; ε m , Δ) (Pandolfi et al. 2002). By minimizing the constrained functional, we obtain an incremental strain energy density W n (ε n+1 ), cf. (Ortiz and Pandolfi 1999) The equilibrium stress then follows from the optimized potential as The current tangent stiffness of the material follows directly as Operatively, the constrained optimization is achieved by introducing Lagrange multipliers as explained in De Bellis et al. (2017).

Fault inception and orientation
Given the deformation ε n+1 at time t n+1 , in the case of an undamaged material, the variational formulation will establish whether the insertion of faults is energetically favorable or not. In the numerical algorithm, we compare two end states of the material, one with newly inserted faults, and the other with no damage. The correct end state will be associated with the lowest incremental energy density W n (ε n+1 ). Clearly, the insertion of faults corresponds to the attainment of an irreversible condition, and the subsequent analysis will always include the faults. The optimal energy is characterized by the optimal orientation N of the faults, obtained at once with the other state variables at time t n+1 , by considering the additional unit length constraint for the fault orientation |N| = 1. We recall here the main results, while the complete discussion can be found in De Bellis et al. (2016). The time-discretized extended potential (1) is modified as By optimizing the constrained functional, the resulting set of nonlinear equations reflect two possible sit-uations. The first situation is characterized by normal opening of the faults, i. e., Δ = Δ N = Δ · N > 0 and Δ S = 0, see Fig. 3a showing the Mohr's circles for the critical stress state. The situation corresponds to the attainment of the Rankine criterion for failure, and it cannot occur if the stress in the matrix is compressive in all directions. The optimality equations reduce to the symmetric eigenvalue problem The eigenvalue Λ is the maximum tensile principal stress experienced by the matrix, the most energetically favorable since it leads to the largest effective traction T and to the least expenditure of cohesive energy. The faults orient themselves so that N is aligned with the direction of the opening.
The second situation is characterized by shear sliding, and it occurs when the matrix is compressed in all directions. In this case, the incipient faults are necessarily closed and deform by sliding, thus Δ = Δ S and Δ N = 0. The resulting optimality equations imply N to be normal to the plane where the matrix shear stress satisfies the Mohr-Coulomb failure criterion, see the orientation α I in Fig. 3b showing the Mohr's circles for the critical stress state. The optimality equations can be restated in the classical scalar form where M = Δ/|Δ| is the direction of the displacement jump and μT c defines the material cohesion (shear resistance at null normal stress). In the sliding mode, the incipient faults orient themselves along planes of maximum shear stress for the matrix. The sliding mode can be activated when the maximum shear planes are under compression. Finally, in cases in which both the opening and sliding modes can operate, they are evaluated in turn and the operative mode is chosen to be the energyminimizing one. The model described in the previous sections depends on an assigned length L, which is the second geometrical feature of the fault microstructure. The parameter L can be computed variationally as part of the incremental update in Eq. 1 if a suitable energetic contribution is added. The additional energy can be thought as the one necessary to accommodate the fault family in order to satisfy the boundary conditions, see Fig. 4. The satisfaction of the compatibility between the faults and their container demands the presence of boundary layers, located where the faults meet a confining boundary, which have a finite thickness and penetrate into the faulted region to a certain depth, see Fig. 4. The corresponding misfit energy, accounting for dislocations piling up at the borders, is of nonlocal type and stored in the boundary layers (Pandolfi et al. 2006). The length scale L of the microstructure results from the balancing of the elastic-cohesive-frictional non-convex energy functional and of the new scale-dependent nonlocal term. The competition allows for the choice of the most energetically favorable size of the microstructures. The resulting relaxed energy is tractable numerically, since it eliminates any degeneracy of the nonconvex model. In particular, for the linearly decreasing cohesive model chosen here, the fault separation L results are defined analytically and are independent of the opening displacement Δ: where L 1 is the size of the container (e. g., the size of the domain in field applications), G the shear modulus and the length L 0 represents an additional feature with the potential to distinguish between ductile and brittle behaviors (Pandolfi et al. 2006). Specifically, a small value of the ratio L 0 /Δ c corresponds to the formation of numerous faults, that may reproduce a diffused displacement jump leading to a ductile behavior, while a large L 0 /Δ c corresponds to the formation of a single fault where the displacement jump is localized and describes a brittle behavior, see Fig. 5. This modification introduces a length scale in the model. The coefficient L 0 is the sixth parameter of the brittle damage material.
More complex microstructures (recursive faulting) can effectively be generated by applying the previous construction recursively. In the first level of recursion, we simply replace in Eq. 1 the elastic strainenergy density W m (ε m ) of the matrix by the effective strain-energy density W n (ε) of the first faulting pattern. The substitution can be iterated, leading to a recursive definition of W n (ε n+1 ). The recursion will end when the matrix between the faults remains elastic, thus W n (ε) = W (ε). Note that the small deformation tensor will result in the sum of the various contributions: The level of recursion defines the rank of the microstructure.

Numerical results
The brittle damage model has been implemented into an in-house finite element code and has been used to reproduce hydro-mechanical tests conducted in the laboratory at the material point level, i.e., without modelling the boundary value problem of the specimen and only considering the constitutive response in terms of stress-strain curves (De Bellis et al. 2016). While the quality of the mechanical response was good, there was no possibility to capture the failure mechanism of the specimen.
We are interested in verifying the ability of the brittle damage model to provide a predictive response in terms of failure mechanism with reference to laboratory tests conducted on samples of soil and rock materials. Specifically, we consider the experiments documented The same displacement jump Δ may be spread over multiple faults, leading to a diffused deformation, or applied to a single fault, leading to a sharp separation between the two cohesive surfaces in the technical note (Yumlu and Ozbay 1995b), consisting of plane strain and triaxial tests. The experimental campaign referred to has been focused on capturing both global responses, in terms of stress-strain curves and local responses, i.e. failure mechanisms, of different rock specimens including coal, sandstone, norite and quartzite. Special attention has been devoted to investigating the influence of confinement, which may result in significantly different responses ranging from the brittle to the ductile regime. The triaxial experiments were carried out in agreement with the standards prescribed by the International Society for Rock Mechanics ISRM.
In the experiments, plane strain and triaxial conditions were imposed using a stiff frame test rig, sandwiching the hexahedral specimen between two rigid plates connected by bolts, parallel to the plane y-z. Plates provided the confinement necessary to avoid out-of-plane displacements in the x direction, in the case of plane strain tests, and the desired pressure for the fully triaxial tests. The confinement pressure on the directions y and z is provided by means of two additional plates; the rig was able to provide a confinement pressure up to 10 MPa. As specified in Yumlu and Ozbay (1995b), the experimental setup was only able to approach plane strain conditions, which are indeed very difficult to assign in experiments. In fact, the displacement in the direction normal to the plane strain can be determined only by considering the relative stiffness of the frame and the rock specimen, so that, depending on the rock type, only "quasi-plane strain" conditions can be guaranteed. The experiments evinced that the results obtained for norite and quartzite were strongly affected by the setup, while a better constraint was achieved for the sandstone and coal rock specimens, stiffer than norite and quartzite. On these grounds, in the following we will simulate only the sandstone experiments.
Material parameters for sandstone rock used in the simulations are listed in Table 1, including the elastic and frictional parameters available from Yumlu and Ozbay (1995b) and the three additional parameters for the brittle damage model, the latter calibrated with preliminary calculations at the material point level. Table 1 Material properties adopted in the numerical tests for the simulation of the sandstone experiments documented in Yumlu and Ozbay (1995b) 0.26 68 128 900 120 The first three parameters are taken from the experimental data, the last three parameters have been calibrated with preliminary calculations at the material point level Fig. 6 Schematic of the boundary conditions applied in the numerical simulation in the plane x-z, for the a plane strain case and b the triaxial case. In both cases, on the planes normal to the y direction, a confinement pressure σ c was applied The laboratory experiments revealed significant differences in the response of the rocks to plane strain and triaxial conditions. The main discrepancies are: (i) higher compressive strength in the z direction observed in plane strain conditions; (ii) thoroughly distinct responses between plane strain and triaxial conditions in the residual stage; (iii) earlier triggering of the yielding stage, and ensuing deviation from the linear behaviour, observed in triaxial tests (Yumlu and Ozbay 1995b). The experimentalists justify these differences with the fact that, in plane strain conditions, the principal stress in x direction increases because of the constrained displacement. Other observed differences are related to the formation of shear bands. In triaxial tests, localized deformations were noticed to begin simultaneously with the onset of yielding, well before the peak stress. By the way of contrast, in plane strain tests, shear bands appeared abruptly in correspondence to the attainment of the peak stress.
We are interested in verifying the performance of the brittle damage model in the two cases. Hexahedral specimens with dimensions 30 × 30 × 10 mm are used for all the tests, and in all test three different confinement stresses, σ c = 3, 5, 8 MPa, are applied. Figure 6 visualizes the boundary conditions adopted in the xz plane for the simulations of plane strain and triaxial tests, respectively. Results are reported in terms of deviatoric stress versus axial strain, defined as where R z is the total nodal reaction in direction z, A 0 is the area of the specimen normal to the z direction and σ c is the confinement pressure, applied to the specimen surfaces normal to the y direction, and where u a is the average relative displacement of the two specimen surfaces in the z direction and H the length of the specimen in direction z (3 mm). We also report volumetric deformations, which are computed as where u x2 , u x1 , and u y2 , u y1 , are the average displacements on the two faces of the specimen normal to the x and y directions, and H measures the dimension of the specimen in y direction. Numerical results are obtained using the mesh shown in Fig. 7, selected after preliminary calculations on coarser and finer meshes that confirmed the meshindependence of the multiscale model. The computational finite element mesh comprises 40,817 nodes and 27,648 10-noded tetrahedra. MPa. In keeping with the experimental results, numerical simulations show that an increase in the confinement level induces a stiffer response, accompanied by a higher peak value and post-peak behaviors. In all cases the numerical model is able to reproduce the initial linear response and the value of the peak stress in a very satisfactory manner. The post-peak softening branch is in general too steep, while the hardening branch is well captured. Note that the model is not able to describe the initial nonlinear response that is observed in the experimental curve for a low confinement. This discrepancy can be explained by the presence of irregularities on the interface between the plates and the specimen that are not described in the numerical model (Klerck et al. 2004).
In every plane strain specimen, experimental failure patterns testify the formation of a well marked shear fracture on the x-z plane, running across the full thickness in diagonal direction. As the confinement increases, the fracture becomes straighter and the ends reach the vertices of the x-z surface, so that the specimen splits into two equal halves. Figure 9 compares the numerical and experimental damage patterns for the Experimental results show the formation of a sharp crack in the same direction. It is interesting to observe that the orientation of the faults at the microscopic level is very different from the macroscopic orientation of the shear band. In fact, the shear band emerges as the results of small sliding along the faults, see Fig. 10.
Numerical simulations were able to point out the overall important phenomenon of dilatancy that was not measured in the experimental work. Figure 11 shows the evolution of the volumetric strain with the axial strain as observed in the numerical outcomes. Regardless of the confinement level, in a fully elastic compressive regime the model experiences volumetric compression up to the peak stress. At the onset of fault inception and in the subsequent phases, the model exhibits dilatancy, an increase of the volumetric deformation progressing until the cohesive stage is fully developed. Afterwards, when no further cohesive dissipation is observed, the model keeps expanding at a constant rate (which can be null). In keeping with common findings in rock mechanics, dilatancy is markedly influenced by the level of the confinement: high confinement pressure reduces the compressive deformation but prevents the development of positive volumetric deformations, so that dilatancy is observed only in the case of lower confinements.

Triaxial conditions
We consider only the case with 5 MPa confinement. Fig. 12 compares the experimental and numerical deviatoric stress versus the axial strain curve. A global reduction of the tangent stiffness attributable to yield-  ing is observed, both in experiment and simulation, well before the attainment of the peak stress. The peak is correctly captured as well as the post-peak softening behavior. In contrast with the plane strain condition results, triaxial tests show a final plateau with constant deviatoric stress for increasing axial strain, although, with respect to the experimental results, the numerical plateau is characterized by a lower value of the residual stress. In terms of mechanical response, the brittle damage model captures the complex behavior of the sandstone. Figure 13 visualizes the contour map of the sliding displacement component at the time of the formation of a neat shear band, corresponding to the point A in Fig. 12. The failure patterns resemble the ones observed in the plane strain condition, in keeping with the experimental results (Yumlu and Ozbay 1995b). It is interesting to notice that, in the absence of the plane strain boundary conditions, a slight out of plane non-symmetric deformation is observed in the numerical results, see the arrow in Fig. 13, which was not observed under plane strain conditions.

Conclusions
In this paper, we presented new applications of a dry brittle damage model in linearized version, by simulating laboratory tests on rocks documented in Yumlu and Ozbay (1995b). The study has been carried out on simple and well-documented cases in order to reveal the MPa. Sliding displacement map after the formation of a marked diagonal shear band, i.e., point A in Fig. 12. The arrow points to a out-of-plane deformation due to the triaxial conditions mechanical features of the brittle damage model and to clarify the relevance of the underlying microstructure in the predictability of the model. We simulated laboratory plane strain and triaxial tests on rock samples that were conceived with the goal of understanding the effects of loading conditions and of the confining pressure on peak and residual strengths (Peters et al. 1988).
Interestingly, the brittle damage model has been able to capture all the distinctive features of the material behavior in terms of stiffness, peak load, and post peak behaviour under different values of confinement. Furthermore the model was able to predict typical crack patterns, characterized by well defined pseudo-planar shear bands, observed in the experiments. As pointed out in Labuz et al. (1996), in the case of plane strain conditions, the failure patterns are marked and well defined. For triaxial tests, the onset of fracture is anticipated by more diffuse phenomena that can be ascribed to yielding, which tends to localize in a shear band leading to a localized rupture, see also (Besuelle 2001) where Vosges sandstone was tested.
The model, in all the numerical tests, has been able to describe the different responses of the sandstone, reproducing both the overall mechanical response and the failure patterns, in keeping with the typical brittleto-ductile transition manifested by geological materials under confinement. Remarkably, all the results have been obtained by using the same set of material parameters, with no need to tune them according to the particular loading condition examined. This property is a natural outcome of the microstructured nature of the model that can be characterized by several length scales and does not suffer of mesh dependency when used in discretized domains.
We conclude that the brittle damage model in the linearized version is a very promising material model for geomechanical problems, especially considering the very small numbers of characterizing material parameters.
Possible applications in the fields of geomechanics include fracking processes (Caramiello et al. 2018), excavations, slope instability, geothermal heating and cooling, and others.