Calculating the cohesion and internal friction angle of volcanic rocks and rock masses

Rock failure criteria are key input parameters for models designed to better understand the stability of volcanic rock masses. Cohesion and friction angle are the two deﬁning material variables for the Mohr-Coulomb failure criterion. Although these can be determined from laboratory deformation experiments, they are rarely reported. Tabulated data for volcanic rocks, calculated using published triaxial results, show that cohesion and friction angle decrease with increasing porosity. If porosity is known, these empirical ﬁts can provide laboratory-scale cohesion and friction angle estimations. We present a method to upscale these parameters using the generalised Hoek-Brown failure criterion, discuss the considerations and assumptions associated with the upscaling, and provide recommendations for potential end-users. A spreadsheet is provided so that modellers can (1) estimate cohesion and friction angle and (2) upscale these values for use in large-scale volcano modelling. Better constrained input parameters will increase the accuracy of large-scale volcano stability models.

Volcanoes are inherently unstable structures. First, they are built by a combination of endogenous growth and the exogenous accumulation of the products of effusive and explosive eruptions, and are often characterised by materials with disparate physical and mechanical properties [Heap and Violay 2021]. Second, their complex construction, coupled with endogenous forcing resulting from, for example, dyke emplacement and magma accumulation (e.g. cryptodome formation), results in oversteepened and unstable slopes. Third, processes acting within the volcano, such as hydrothermal alteration, can reduce the strength of volcanic rocks [del Potro and Hürlimann 2009;Julia et al. 2014;Wyering et al. 2014] and create weakened layers that can promote instability, sliding, and collapse [López and Williams 1993;van Wyk de Vries et al. 2000;Reid et al. 2001;Ball et al. 2018]. As a result, partial flank collapse is not uncommon at volcanoes worldwide [Siebert et al. 2010]. The partial collapse of a volcanic flank can result in debris avalanches and flows, laterally-directed explosions, and high-energy pyroclastic density currents that pose a substantial hazard to those communities living close to the volcano [Siebert et al. 2010;Cole et al. 2015]. Due to the hazard posed, flank deformation is often closely monitored by volcano observatories [e.g. Moretti et al. 2020].
Largescale open-source models have been developed and Corresponding author: marlene.villeneuve@unileoben.ac.at used to assess flank instability scenarios [e.g. Scoops3D, used by Reid et al. 2015;Ball et al. 2018;Finn et al. 2018;Schaefer et al. 2019;Kereszturi et al. 2021] and LaMEM, [used by Kaus et al. 2016;Heap et al. 2021a;b], and there are a variety of commercial software packages available (e.g. FLAC by Itasca Consulting Group , RS2 by Rocscience † , and PLAXIS 2D by Virtuosity ‡ ) that have also been used to assess volcano stability [e.g. Apuani et al. 2005a;Moon et al. 2009;Schaefer et al. 2013;Wallace et al. 2021]. Rock physical and mechanical properties are important input parameters in slope stability models. However, while data compilations for properties such as Young's modulus [Heap et al. 2020b] and uniaxial compressive strength [Heap and Violay 2021] are available for volcanic rocks, other input parameters are less readily available. For example, the required input parameters for defining failure criteria, such as the cohesion and internal friction angle for the Mohr-Coulomb failure criterion, of volcanic rocks are rarely reported, as noted by Ball et al. [2018]. To assist volcano modellers, we first detail a simple method to determine the cohesion and internal friction angle using data from laboratory deformation experiments, which we then use to create a table of values for volcanic rocks using published triaxial deformation data. Empirical relationships between porosity and cohesion and internal friction angle are provided, which can be used to esti-mate these values in the absence of laboratory data. We then outline in detail a method to upscale these values to account for macrofractures in the rock mass, using the generalised Hoek-Brown failure criterion. We show how uniaxial compressive strength and Hoek-Brown fitting parameter, (m i ), values required to upscale cohesion and internal friction angle to the scale of the rock mass, are derived from laboratory triaxial data. We also provide empirical relationships between porosity and uniaxial compressive strength and Hoek-Brown fitting parameter m i , which can be used in the absence of laboratory data. Finally, we discuss the various considerations and limitations of the upscaling approach presented. A Microsoft Excel® spreadsheet is provided as Supplementary Material to help end-users calculate input parameters for failure criteria for their modelling.

Calculating cohesion and internal friction angle from laboratory data
The Mohr-Coulomb criterion, τ = c +σ tan φ f , where τ is the shear strength at effective confining pressure σ , c is the cohesion, and φ f is the internal friction angle, is a commonly used failure criterion for soil and rock. This criterion is derived from laboratory triaxial deformation experimental data (σ 1 > σ 2 , and σ 2 = σ 3 , where σ 1 , σ 2 , and σ 3 are the maximum, intermediate, and minimum principal stresses, respectively) [Labuz and Zang 2012]. Cylindrical samples (typically 10-62 mm in diameter), prepared from a single block of material, are deformed until macroscopic failure (the formation of a throughgoing shear fracture) at a constant stress or strain rate under different confining pressures, P c . These experiments can either be performed on dry samples or samples saturated with water (with or without a pore fluid pressure). If a pore fluid pressure is used, it is commonplace to assume that the effective pressure, P ef f , is the confining pressure, P c , minus the pore fluid pressure, P p . To calculate cohesion and internal friction angle, at least three experiments are required (at different confining or effective pressures) and all the experiments should be performed in the brittle regime [Labuz and Zang 2012]. The formation of a shear fracture is the hallmark of brittle deformation; high pressures, which can promote ductility, should be avoided in experiments on porous volcanic rocks; [see Heap et al. 2015a]. One of the experiments can be performed under uniaxial conditions (where σ 1 > σ 2 and σ 2 = σ 3 = 0). What is required from each of the experiments is the value of σ 1 at macroscopic sample failure and σ 3 (the confining pressure for a dry experiment); when a pore fluid pressure is used, σ 1 and σ 3 are assumed to be equal to the effective pressures, σ 1 and σ 3 ).
To derive the Mohr-Coulomb failure criterion, the triaxial results can be either plotted as Mohr circles in shear stress versus effective normal stress space, or as stress points in shear stress versus mean effective normal stress space [Craig 1997]. The latter method is easiest to implement analytically, for example using a spreadsheet: the series of triaxial results are first converted to maximum shear stress, τ m , where τ m = (σ 1 − σ 3 )/2, and mean effective normal stress, σ m , where σ m = (σ 1 + σ 3 )/2 . A best-fit linear regression to these data, τ m = a · σ m + f , provides two fitting parameters, a and f . The standard error of these coefficients is calculated and provided as a "±" after the value. The cohesion and internal friction angle are then computed using the following transformation from maximum shear stress versus mean effective normal stress space to maximum shear stress versus effective normal stress space [Craig 1997]: An example, using three datapoints for Volvic trachyandesite (France; data from Heap and Violay [2021]), is shown in Figure 1. For these data, a = 0.636 ± 0.047 and f = 16.8 ± 3.48 ( Figure 1). The errors on a and f were propagated through Equations 1 and 2, and are also provided as a "±" after the value. The internal friction angle, φ f , and the cohesion, c , are then 39.5 ± 3.56 and 21.8 ± 5.97 MPa, respectively.
Using this technique, we provide here laboratoryscale values of cohesion and internal friction angle derived for a range of volcanic rocks and textures (andesite, basalt, dacite, and tuff), using data from the . For Volvic trachyandesite, the linear regression fitting parameters required to calculate internal friction angle and cohesion (using Equations and , respectively) are a' = .6 6 ± . 6 8 and f' = 6.8 ± . .
published literature (Table 1). This is based on a compilation of well-documented triaxial experimental data that are conducted at room temperature, similar strain rates, and also include porosity measurements. The derived data given in Table 1 are plotted in Figure 2 as a function of average connected porosity (henceforth referred to simply as "porosity") for the dataset. Figure 2 shows that cohesion and internal friction angle both decrease as a function of increasing porosity. In the absence of triaxial data with which to calculate the cohesion and internal friction angle parameters, they could be estimated using the porosity, ϕ, with the following empirical regressions to the data from Table 1: As with all empirical relationships, the regressions presented in Equations 3 and 4 should only be used within the porosity range of the dataset on which they are based (0.04 to 0.3). Measures of the correlation between porosity and φ f and c for these datasets are given as Spearman's rank correlation, , because it is applicable to monotonic datasets regardless of the regressions ultimately selected, whether linear or nonlinear, unlike the goodness of fit tests Pearson's r or r 2 , which are only applicable to linear regressions. Spearman's provides a measure of the correlation between two variables where a value of 1 and −1 are perfect positive and negative correlations, respectively, and 0 represents no correlation (as described in ). Spearman's for the relationship between porosity and cohesion is −0.72 (strong correlation) and between porosity and internal friction angle is −0.24 (weak correlation). Since the Spearman's does not provide a measure of goodness for tested regressions (linear, power, logarithmic, exponential), the standard error of the estimate (SEE; Lane [2003]) was used as the measure of the goodness of fit for different regression types. The fitting function with the lowest SEE was then selected as the preferred function. The SEE is 6.0°for Equation 3 and 8.2 MPa for Equation 4. The SEE is slightly higher at porosity > 0.15 than at porosity < 0.15 (11°versus 5.8°, 14.7 MPa and 8.1 MPa, for Equations 3 and 4, respectively), highlighting the higher variability in the data at higher porosity.

Upscaling cohesion and internal friction angle
Laboratory-measured rock physical properties require upscaling to be used in large-scale modelling, to take into consideration the fact that macrofractures at the field scale, which are not present in laboratory samples, decrease the strength of the rock mass. There is   Table ). Empirical fits to the data are given as dashed lines.
no method for directly upscaling cohesion and internal friction angle, however in rock engineering practice the generalised Hoek-Brown failure criterion is first derived, then transformed to the Mohr-Coulomb failure criterion to obtain rock mass cohesion and internal friction angle [Hoek and Brown 1997].
The first step is to derive the generalised Hoek-Brown failure criterion (Equation 5) for fractured rock masses [see Eberhardt 2012;Hoek et al. 2013;Heap et al. 2020b]. The variables in Equation 5 are a combination of laboratory strength measurements (σ ci is uniaxial compressive strength, in MPa) and dimensionless empirical fitting parameters, m b (Equation 6), s (Equation 7), and a (Equation 8).
where σ 1 and σ 3 are the effective maximum and minimum principal stresses, respectively, m i is the dimensionless Hoek-Brown constant, and GSI is the Geological Strength Index. Finally, the dimensionless parameter D describes the extent of blasting damage during large mining/tunnel excavations. In the absence of anthropogenic blasting disturbance, we suggest here that D is zero (as discussed in Heap et al. [2020b]) and we include it here for completeness. The m i and GSI are explained in more detail below. The generalised form of the Hoek-Brown failure criterion defines the compressive strength, σ 1 , of the rock mass for any confining pressure, σ 3 . In the case of volcano modelling, the lithostatic pressure acting on the rock mass at a depth h, can be calculated using σ 3 = ρgh, where ρ is the density of the rock mass (typically taken as 2500-2900 kg m −3 ) and g is the acceleration due to gravity (= 9.81 m 2 s −1 ). σ ci is the uniaxial compressive strength, the peak stress achieved by a rock sample prior to macroscopic failure under unconfined conditions, and can be determined from laboratory testing (uniaxial compressive strength data for volcanic rocks, including a downloadable Ex-cel® spreadsheet containing the data, are available in Heap and Violay [2021]). The dimensionless Hoek-Brown constant, m i , as well as σ ci , can be calculated by fitting the intact rock version (Equation 9) of the generalised Hoek-Brown failure criterion from Equation 5 to laboratory triaxial deformation data (following the mathematics detailed in Appendix 2 of Hoek and Brown [1997]).
As stated in Section 2 for deriving laboratory scale cohesion and internal friction angle, a minimum of three points are needed for data fitting, and five points are needed for reliable data fitting [Hoek and Brown 1997], and we have found that the regression residuals are lower when one of these points is from a uniaxial compressive strength test. By reorganising Equation 9 into the linear form in Equation 10, where x = σ 3 and y = (σ 1 − σ 3 ) 2 , it is possible to solve for m i and σ ci for n specimens using Equations 11 and 12 [Hoek and Brown 1997].
We have calculated m i for the same datasets as in Table 1 (excluding datasets that contain only two datapoints). These values are given in Table 2 and plotted as a function of porosity in Figure 3. The Spearman's correlation between porosity and m i for this dataset is −0.42 (moderate correlation). Therefore, in the absence of triaxial data with which to calculate m i , this parameter can be estimated if the porosity, ϕ, is known using the following empirical fit to the data from Villeneuve et al. [2021]:  Table ). The standard error of the estimate is 8. .
If data are not available for a specific volcano, or laboratory testing is not possible, the uniaxial compressive strength can also be estimated using empirical fits to the compiled datasets given in Figure 4, if the porosity, ϕ, is known. The empirical fit in Equation 14 can be used to provide estimates for the uniaxial compressive strength of volcanic rocks based on the data in Heap and Violay [2021] (Figure 4) using the fitting parameters given in Table 3.
The relationships in Figure 4 exhibit high variability, especially for samples with low porosity, as highlighted in Table 3 where the SEE is given for the full porosity range, and separately for two porosity ranges, < 0.1 and > 0.1. Nevertheless, the Spearman's correlations are strong to very strong ( Table 3), supporting that porosity is a suitable parameter with which to estimate uniaxial compressive strength.
The GSI parameter is a dimensionless parameter that describes the rock mass structure (i.e. fracture density) and the quality of the fracture surfaces (i.e. whether they are weathered and the presence/nature of fracture filling material) [Marinos et al. 2005]. The GSI is a number from 0 to 100, where 100 represents a pristine (fracture-free) rock mass and 0 represents an essentially cohesionless rock mass. Increasing GSI, representing an increase rock mass quality, also increases both rock mass cohesion and internal friction angle ( Figure 5). This is particularly pronounced for rock masses with a GSI greater than 60, where the impact of the fractures is far less than for poorer quality rock masses. Assessing GSI directly in the field is preferable, but in the likely event that this is not feasible, Heap et al. [2020b] concluded, using a dataset of GSI estimates for volcanic rock masses, that a GSI of 55 represents an average volcanic rock mass, which are typically "blocky". Therefore, if GSI is not known, we recommend that end-users also use a value of 55 in the analysis described herein.
The next step is to transform the Hoek-Brown failure criterion to the Mohr-Coulomb failure criterion according to Hoek et al. [2002]. Using this transformation, the rock mass (denoted as subscript "rm") friction angle, φ rm , (Equation 15) and rock mass cohesion, c rm , (Equation 16) can be calculated.
where m b , s, and a are from Equations 6, 7 and 8, respectively, and: The input parameters for Equations 17, 18 and 19 are calculated from laboratory values, unit weight (γ in kN m −3 ), and the empirical parameters from Equations 6, 7 and 8. However, the σ 3max parameter must be linked to the geometry of the problem: h, the height of the slope or volcano in question [Hoek et al. 2002]. In the event of computing rock mass cohesion and friction angle for subsurface excavations, a different version of σ 3max must be used, as described and given in     [Hoek et al. 2002]. This parameter is necessary because the derivation of rock mass cohesion and friction angle is based on linear fits to the curved Hoek-Brown failure envelope, which have different slopes and intercepts depending on over which confining stress range the fit is made. The selection of h will affect the computed rock mass cohesion and internal friction angle, whereby cohesion will be higher and friction angle will be lower for higher slopes, and vice versa. Figure 6 demonstrates this for rock mass cohesion and friction angles for example volcanic rock masses with 0.05, 0.18, and 0.3 (denoted as low, moderate, and high) porosity, a constant GSI of 55, and density of 2650 kg m −3 at slope heights from 0-3000 m.
While it is possible to compute rock mass cohesion and friction angle for any slope height, when modelling a large and/or complex topography, such as an entire volcano, the selection of an appropriate single slope height can be difficult. In such cases it may be preferable, if possible, to use the generalised Hoek-Brown failure criterion directly in the model, rather than using it to derive rock mass cohesion and friction angle for use in the Mohr-Coulomb failure criterion. To illustrate this, we computed three different two-dimensional finite element models (using the commercial software RS2 ) of a hypothetical slope in volcanic rock using the Hoek-Brown and Mohr-Coulomb parameters for the low, moderate, and high porosity volcanic rock masses from Figure 6, using a GSI = 55 and a slope height = 250 m. The use of a two-dimensional representation of the slope neglects edge effects and possible contributions from the out-of-plane geometry, however this approach is sufficient for this demonstration. We used the shear strength reduction method [Dawson et al. 1999;Hammah 2005;Diederichs et al. 2007]. The shear strength reduction method is an automated numerical modelling solution method for slope stability analysis in which the material properties are decreased or increased, the model is solved, and the convergence of the model provides the indication for whether or not the slope is stable. If the model does not converge to a solution, the slope is considered unstable. The amounts by which the original material parameters must be increased or decreased to ensure model convergence are compared to the original material parameters to compute the strength reduction factor (SRF). This is analogous to a factor of safety against slope failure, whereby  Table . a value less than unity denotes instability, and a value greater than unity denotes stability.
In all three porosity scenarios, the strength reduction factor was above 2, suggesting the modelled edifice is stable for the given input parameters and boundary conditions (Table 4). While the strength reduction factor and maximum shear strain are similar between the failure criterion models, suggesting that both failure criteria provide consistent results, the differences that exist are not consistent (e.g. Hoek-Brown models are not consistently more stable than Mohr-Coulomb models; Table 4).
Another measure of the impact of failure criterion selection on model results is the strength factor, which is calculated by comparing rock mass strength at every element in the finite element mesh with the principal and shear stress invariants considering both in-plane stresses and the out-of-plane stress. A stress factor > 1 represents material strength greater than the stresses,  Table . whereas a stress factor < 1 represents material strength less than the stresses. The Hoek-Brown failure criterion allows strength to be directly computed for each element at the exact stresses in the element. By contrast, the Mohr-Coulomb failure criterion computes the strength for a curve fitted at a particular slope height (in this case 250 m), using the stresses in the element. This means that it can overpredict the strength factor at slope heights (or depths within the slope) greater and less than the slope height used to derive rock mass cohesion and friction angle. This is demonstrated in Figure 7, where the models that used the Mohr-Coulomb failure criterion have smaller areas with SRF ≤ 3. Areas at the upper (shallower than 250 m) and lower (deeper than 250 m) portions of the edifice have higher SRF than in the Hoek-Brown models because the strength is overpredicted here. There are also very small concentrations with low strength factor along the topographic boundary which are likely numerical artefacts associ- ated with mesh shape in these locations. This demonstrates that stress-induced failure (strength factor <1) could be under-or overpredicted in areas at different stress conditions when using the Mohr-Coulomb failure criterion for rock masses.
Where the model code (or software) only accommodates a single value of cohesion and friction angle, we suggest that a sensitivity analysis approach be first undertaken using upper, middle, and lower ranges for slope height to capture the various effects that parameter selection have on the modelled volcanic processes. Care must be taken in the interpretation of these results, however, since they will not provide upper (best case scenario), middle, and lower (worst-case scenario) end-member results. Because of the divergent variation of rock mass cohesion and friction angle with slope height, the models will provide variable results with potential impacts on deformation magnitude and style, as well as stability. We would certainly encourage any model developers, if possible, to implement the generalised Hoek-Brown failure criterion directly into modelling codes to eliminate these effects.
Indeed, sensitivity analysis of all of the parameters in a model should be considered best practice, in particular considering the high variability of rocks and rock masses and the inherent uncertainty in the resulting input parameters, as demonstrated by the propagated errors given in Table 1 and the standard error of the estimate (SEE) for the porosity-based relationships developed for the datasets we present here. The propagated errors on The SEE can be used to calculate upper and lower bound values with which to conduct sensitivity analyses, and if the errors are assumed to be normally distributed, then probabilistic analyses could also be conducted.

Calculation spreadsheet
The methods for obtaining the input parameters for the Mohr-Coulomb and Hoek-Brown failure criteria of intact rock and rock masses detailed in Sections 2 and 3 have been implemented in a Microsoft Excel® workbook provided in the Supplementary Material. The workbook has two worksheets, one for intact rock and one for rock masses. The intact sheet contains the mathematics described in Section 2 and takes confining pressure (σ 3 ) and maximum axial stress (σ 1 ) as input. It provides intact cohesion (c ), internal friction angle (φ f ), computed uniaxial compressive strength (σ ci ) and the Hoek-Brown constant mi. It provides the propagated standard error for cohesion and internal friction angle. It also provides a plot of the Mohr circles in shear stress versus effective normal stress space, as well as the linear fit to the stress data in maximum shear stress versus mean effective normal stress space. The rock mass sheet contains the mathematics described in Section 3 and takes the geological strength index (GSI), damage parameter (D), depth of interest (this is the depth parameter h in Equation 18), unit weight (γ), as well as intact elastic modulus and intact Poisson's ratio (if desired). It provides rock mass cohesion (c rm ), rock mass internal friction angle (φ rm ), for both slopes (the mathematics for which are given in Equation 18) and subsurface excavations (the mathematics for which can be found in Hoek and Brown [1997]). It also provides the rock mass uniaxial compressive strength (σ c ; solved by setting σ 3 to 0 in Equation 5), and the Hoek-Brown constants m b , s and a. The rock mass elastic modulus and rock mass Poisson's ratio (if desired) are computed according to the methods described in Heap et al. [2020b]. The calculation and output cells are locked, to ensure integrity of the spreadsheet, however the equations are visible and can be interrogated or used to implement into other software or programming languages.

Summary
We provide a workflow guidance for determining cohesion and internal friction angle of intact and fractured volcanic rock using a combination of published experimental data and failure criteria, presented visually in Figure 8. The workflow is useful for deriving a wide range of input parameters for modelling volcanic processes at a variety of scales. By linking the intact rock parameters to porosity, which is easily measured and widely reported, it is possible to estimate intact cohesion, internal friction angle, uniaxial compressive strength and Hoek-Brown constant m i in Figure : Two-dimensional finite element modelling results of slope stability for low, moderate, and high porosity volcanic rocks with a Geological Strength Index (GSI) of and density of 6 kg m − , calculated using the m i values from Table , using generalised Hoek-Brown failure criterion and Mohr-Coulomb failure criterion derived for a slope height of m. The blue stress block in the lower left shows that vertical and horizontal gravitational stresses are equal. Only strength factors ≤ are displayed to highlight differences between the modelling methods.
the absence of laboratory mechanical data. The rock mass quality can be determined via field observations or by using an overall average value, such as GSI = 55. The intact parameters, uniaxial compressive strength (σ ci ) and the Hoek-Brown constant m i , are then combined with the rock mass quality to determine rock mass strength, defined either directly using the generalised Hoek-Brown failure criterion parameters m b , s, and a, or via a slope-height related transformation to obtain the Mohr-Coulomb failure criterion parameters: rock mass cohesion and rock mass internal friction angle. We provide a Microsoft Excel® workbook with two worksheets, one to calculate intact cohesion and internal friction angle, and a second to calculate rock mass cohesion and friction angle, as well as the generalised Hoek-Brown parameters m b , s, and a. We demonstrate the use of standard errors and propagation of these errors when deriving cohesion and internal friction angle from triaxial data as a way to highlight the uncertainty associated with deriving failure criterion parameters from real rock samples. These errors could be further propagated to the empirical relations between porosity and the failure criterion parameters. We also provide assessments of the quality of the correlation and the standard error of the estimate for the empirical relations, again to highlight the uncertainty arising from experiments on highly variable rocks. These errors could be used for sensitivity or probabilistic analysis in stability modelling.
We end with some recommendations for conducting a sensitivity analysis to establish the effect of slope height selection on model results when using the fitted Mohr-Coulomb failure parameters. Finally, where possible, we encourage the implementation of the generalised Hoek-Brown failure criterion directly into numerical software, as it is in commercial geotechnical modelling packages such as FLAC (by Itasca Consulting Group ), RS2 (by Rocscience † ), and PLAXIS 2D (by Virtuosity ‡ ).

Data availability
All data and equations used are either previously published, presented in the manuscript, or are available in the Microsoft Excel® spreadsheet that accompanies this contribution as Supplementary Material.

Copyright notice
© The Author(s) 2021. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.