A hyperelastic model for soils with stress-induced and inherent anisotropy

In this paper, modelling of the superposition of stress-induced and inherent anisotropy of soil small strain stiffness is presented in the framework of hyperelasticity. A simple hyperelastic model, capable of reproducing variable stress-induced anisotropy of stiffness, is extended by replacement of the stress invariant with mixed stress–microstructure invariant to introduce constant inherent cross-anisotropic component. A convenient feature of the new model is low number of material constants directly related to the parameters commonly used in the literature. The proposed description can be incorporated as a small strain elastic core in the development of some more sophisticated hyperelastic-plastic models of overconsolidated soils. It can also be used as an independent model in analyses involving small strain problems, such as dynamic simulations of the elastic wave propagation. Various options and features of the proposed anisotropic hyperelastic model are investigated. The directional model response is compared with experimental data available in the literature.


Introduction
An elastic stress-strain relation is the core of all the elastoplastic constitutive models. It describes the initial stiffness and influences substantially the modelled pre-failure behaviour of soils. In the case of natural soils, there are two important features of the initial elastic stiffness, namely barotropy and anisotropy. Barotropy is typically taken into account as the dependency of stiffness on the actual stress level that is represented by the mean stress p or by one of the principal stress components r i . Models used in the geotechnical practice, e.g. [49], mostly incorporate the exponential empirical relation in the form proposed by Ohde [45] or Janbu [31]. In this relation, stress dependency of stiffness concerns elastic moduli in the isotropic Hooke's law and the Poisson's ratio m is assumed constant. This results in an isotropic hypoelastic relation which, of course, is dependent neither on the stress obliquity (e.g. g ¼ q=p or K ¼ r 3 =r 1 ) nor on the direction of loading.
Influence of the stress obliquity is often described as stress-induced anisotropy. In the case of stress-induced anisotropy, a directional distribution of stiffness changes depending on the current stress position in the principal stress space. Stress-induced anisotropy vanishes under isotropic stress conditions. However, stiffness of natural soils remains anisotropic even under isotropic stress conditions. This is due to the anisotropic microstructure of soil fabric that was developed as a result of deposition and diagenesis [23,30,39]. This component of anisotropy is described as inherent anisotropy. Its properties can be investigated only under isotropic stress conditions. In natural soils, the actual directional dependency of elastic stiffness is a superposition of both stress-induced and inherent anisotropy. Such situation is illustrated schematically in Fig. 1. It is an exemplary problem of the slope stability, in which the principal axes of stress n ri and microstructure x mi are generally not collinear. They are also different from the applied geometrical coordinate system x i which usually corresponds to the direction of gravity.
A superposed anisotropy of the actual elastic stiffness is usually complex. Principal stress directions may rotate relative to the soil microstructure as a result of geotechnical processes (e.g. excavation, tunnel drilling). Hence, the directional dependency of stiffness changes with the angle between the principal axes of stress and of microstructure. In the standard laboratory tests, such as triaxial or oedometric compression, these systems are collinear and rotation of the principal stress directions is not possible. Contrarily to stress-induced anisotropy, inherent anisotropy of the elastic stiffness remains rather constant even after significant stress changes [32]. This seems to be evident in natural heavily overconsolidated clays. In reconstituted or moderately overconsolidated clays, inherent anisotropy may evolve due to irreversible straining as it is shown in experimental investigations by Mitaritonna et al. [40].
In the recent years, the importance of stiffness nonlinearity and anisotropy in the small and intermediate strain ranges has been shown in many analyses of practical cases in geotechnics [1,3,11,16,35,42,47]. Anisotropy of the elastic stiffness is also known to influence the yielding characteristics of soils in the undrained conditions [44,50,56].
A common type of inherent anisotropy is the crossanisotropy, also called the transverse isotropy. In this case, the elastic material properties are symmetric about axis that is normal to the plane of isotropy. In soils, the plane of isotropy is identified with the prevailed particle orientation or bedding planes. The symmetry axis, here indicated by the unit vector v, is collinear with the deposition direction. This direction is usually vertical; however, it may be inclined due to geological processes, especially in the older overconsolidated clayey deposits, as depicted in Fig. 1.
In the classical linear elastic description, a cross-anisotropic stiffness requires the definition of five independent material constants [46]. Different sets of five independent cross-anisotropic parameters may be selected. Using notation with indices h and v to denote the horizontal and the vertical directions, respectively, these parameters are usually chosen as: E v , Young's modulus in vertical direction; E h , Young's modulus in horizontal direction; m vh , Poisson's ratio for horizontal strain due to vertical strain; m hh , Poisson's ratio for horizontal strain due to complementary horizontal strain; G vh , shear modulus in vertical plane. This naming convention assumes the horizontal plane of isotropy. However, if vector v deviates from vertical direction, plane h À h is also rotated relative to the horizon. Hence, one needs to note that indices h and v are associated with the cross-anisotropic microstructure.
Based on the experimental observations, Graham and Houlsby [21] proposed a simplification of the five-parameter cross-anisotropic stiffness. They reduced the number of material constants to three independent parameters leaving E v and m hh , denoted by E Ã and m Ã ; respectively, and introducing the anisotropy coefficient a that imposes the following condition: where G hh is the shear modulus in the plane of isotropy: Values of the anisotropy coefficient a in the overconsolidated soil deposits are generally higher than 1.0, which indicates higher stiffness in the plane of isotropy, i.e. G hh [ G vh and E h [ E v . In normally consolidated or lightly overconsolidated soils, stiffness may be higher in the direction of symmetry axis [37,55]; however, the opposite situation is also reported [9,36]. Intensity of inherent cross-anisotropy, expressed as a deviation of a from 1.0, is also observed to be higher in clayey deposits [8,26] than in sands. This is due to the fact that the total Fig. 1 Geometrical x i , stress n ri and microstructural x mi principal axes in a representative soil element and in an exemplary boundary value problem of the slope stability analysis anisotropy in sands is mostly dominated by stress-induced component [27,34]. Mašín and Rott [38] refined the definition of a into three coefficients corresponding to the Young's and shear moduli, and the Poisson's ratios, respectively: Additionally, to compare the anisotropy coefficients, they proposed two anisotropy exponents x GE and x Gm : It follows that in the cross-anisotropic stiffness proposed in [21], the imposed values of anisotropy exponents are: x GE ¼ 0:5 and x Gm ¼ 1:0. However, after careful inspection of the recent experimental data for natural clays, this is shown in [38] not to be necessarily valid and approximate average value of x GE ¼ 0:8 is suggested. Due to difficulties in the measurement of the anisotropic Poisson's ratios, the values of a m and x Gm are characterised by large scatter; hence, further experimental investigations are needed. It is suggested in [38] that, for practical purposes, x Gm ¼ 1:0 should be adopted. The classical linear cross-anisotropic stiffness is not stress-dependent. One should be careful when introducing an empirical barotropy of cross-anisotropic elastic parameters into this model. This is because of the constraints imposed on these parameters by the elastic energy conservation requirement [46].
The objective of this work is to present an elastic model capable of reproducing the stiffness barotropy, stress-induced anisotropy and inherent cross-anisotropy. For this purpose, we use the framework of hyperelasticity [25,28,43]. The proposed model is based on an improvement of the elastic potential function by Vermeer [51]. The new refined potential function is not dependent on the stress invariant allowing stress-induced anisotropy only but instead on the new mixed stress-microstructure invariant. The mixed invariant introduces inherent anisotropy component in a way proposed by Boehler and coworkers [5][6][7]. This method allows to control the values of anisotropy coefficients and the exponents defined in [38]. The initial works on the proposed refinement of Vermeer's elastic potential are reported in monograph [10] as a part of a hyperelastic-plastic model of overconsolidated fine-grained soils.
The general framework of incorporating inherent anisotropy based on the microstructure within the hyperelasticity has been recently presented by Houlsby et al. [29] and validated with the experimental evidence by Amorosi et al. [2]. Another notable works on anisotropic hyperelastic models for soils have been reported in the literature by Gajo and Bigoni [17,18] and Xiao et al. [54]. An interesting method of incorporating inherent cross-anisotropy into hyperelastic model by scaling of stiffness has been shown by Niemunis et al. [44].

Formulation of the model
In the elasto-plastic modelling, the initial stiffness is considered elastic within a small region in the stress or strain space. A constitutive model of elastic behaviour should provide a closed stress or strain loop within the elastic locus when material undergoes a closed strain or stress loop, respectively. Additionally, neither the dissipation nor the generation of the elastic energy is allowed in such process. It means that the following conditions must be fulfilled: where C t and D t are the fourth-order tangent elastic compliance and stiffness tensors, respectively. If the thermodynamic requirements from Eq. 8 are satisfied, the model is considered hyperelastic. The most convenient way to formulate a truly conservative hyperelastic model is to derive it from the elastic potential function W [25,28,43]. It is a scalar function of stress-negative Gibbs free energy WðrÞ or elastic strain-Helmholtz free energy Wðe e Þ. An elastic potential function is not a subject to any direct measurements; hence, it should be elaborated, for example, by inspection and analysis of its derivatives constituting the compliance or stiffness operators. In the case of first stress derivative of WðrÞ, we obtain a secant compliance tensor C s relating the elastic strain and stress: To obtain a tangent compliance tensor C t that relates the elastic strain and stress rates, the second stress derivative of WðrÞ should be calculated: If an elastic potential is expressed as a function of elastic strain Wðe e Þ, we obtain secant and tangent stiffness tensors, D s and D t , by analogous differentiations with respect to the elastic strain.
A comparison of different stiffness formulations and investigation of the influence of their parameters may be achieved by confronting the so-called response envelopes [22]. Response envelope is a polar diagram of the tangent stiffness or compliance tensor. It is usually shown in the triaxial stress or strain plane ( ffiffi ffi 2 p r 3 À r 1 ; ffiffi ffi 2 p e e 3 À e e 1 ) as a closed curve representing the stiffness or compliance response to a circular stress or strain probe, respectively. Schematic description of obtaining the stiffness response envelope due to strain probe is shown in Fig. 2. Response envelopes can be obtained also experimentally. However, it is a complicated technique that still remains under development. Experimental response envelopes for sands have been reported by Danne and Hettler [13] and Knittel et al. [33].

Basic hyperelastic model with stress-induced anisotropy
Among different proposals of elastic potential functions available in the literature, we decided to make use of the stress-dependent formulation by Vermeer [51]. This function will be later modified and extended to account for inherent anisotropy. In our experience, it is an optimal solution that combines the practical functionality and simplicity associated with a small number of material parameters. The elastic potential has the following form: where G ref 0 is the small strain shear modulus measured at the isotropic reference stress p ¼ p ref , b is a material parameter controlling the order of stress dependency of stiffness, and Q is the stress invariant defined as: After differentiation (Eq. 10), the tangent compliance tensor is obtained: where G 0 is the stress-dependent shear modulus: After analytical inversion of the compliance tensor C t , we obtain the tangent stiffness tensor: Tangent stiffness D t is a homogeneous function of stress of order m ¼ 1 À b. However, parameter b also influences the directional distribution of stiffness in the principal stress space. If the considered stiffness response is compared to that of the linear isotropic Hooke's law under isotropic stress conditions, parameter b may be related to the Poisson's ratio as follows: It is important to note some limitations of the model applicability related to the coupling between m and m which is controlled by the model parameter b. It is obvious that condition b 6 ¼ 0 must be satisfied; otherwise, it implies m ¼ 0:5 and the infinite stiffness, see Eq. 15. However, condition b 6 ¼ 0 eliminates also the possibility of linear stress dependency of stiffness (m ¼ 1:0) as observed in normally consolidated clays. Hence, application of the model and its refinement presented in this paper is developed to simulate the elastic behaviour of overconsolidated fine-grained soils Fig. 2 Response envelope in the triaxial plane. Version with a circular strain probe and the ellipsoidal response in the triaxial stress plane (r 2 ¼ r 3 ). Unlabelled dots correspond to purely deviatoric (De e v ¼ 0) and purely volumetric (De e q ¼ 0) strain increments or granular deposits, for which the values of m are reported within the range from 0.3 to 0.7, e.g. [4,8,14]. The coupling between m and m controlled by the value of b is shown in Fig. 3. Additionally, response envelopes that illustrate the characteristic properties of the employed hyperelastic stiffness are presented. Size of the envelopes increases with the stress level which indicates the barotropy. Stress-induced anisotropy is evident by rotation of the envelopes resulting from the stress obliquity, i.e. K 6 ¼ 1:0. Elongation of the response envelopes, which shows a higher stiffness in compression along the stress paths K ¼ const , is directly controlled by the values of b. This effect may be compared to influence of the Poisson's ratio on stiffness in the isotropic Hooke's law. In Fig. 3, response envelopes of the hyperelastic stiffness are equivalent to those of the isotropic Hooke's law only for p ¼ p ref and K ¼ 1:0 with the shear modulus and Poisson's ratio calculated with Eqs. 14 and 16, respectively.

Definition of inherent cross-anisotropy
To define a cross-anisotropic microstructure, we use second-order tensor M proposed in [7]. Its components are calculated from the following dyadic product: where the unit vector v identifies the symmetry axis normal to the plane of isotropy. Using spherical coordinates /; h in the geometrical framework shown in Fig. 4, the Cartesian coordinates of v are: After operation (Eq. 17), the general form of tensor M representing the cross-anisotropic microstructure is as follows: In this paper, the horizontally oriented plane of isotropy with the vertical symmetry axis is usually assumed, i.e. v ¼ ½1; 0; 0 T and h ¼ 0, which leads to:

Formulation based on the mixed stressmicrostructure invariant
The modelling aim is to incorporate information on the cross-anisotropic microstructure into the framework of hyperelasticity using tensor M. We start from the basic Fig. 3 Parameter b coupling the order of stress dependency of stiffness m and the Poisson's ratio m as compared to the Hooke's elasticity at isotropic reference stress conditions. Left hand side: relation between b; m and m. Right hand side: response envelopes of Vermeer's basic hyperelastic model [51] for different stress ratios K, stress levels p and b values. Response envelopes are scaled for visualisation purposes function used to derive the stiffness or compliance operators, i.e. the elastic potential function W. According to [5,7], a general scalar function of two symmetric tensorial arguments, here f ðr; MÞ, can be expressed as a function of the following three groups of scalar invariants: They represent the stress, microstructure and jointed stress-microstructure invariants, respectively. In the general case of cross-anisotropy (Eq. 19), all three invariants of microstructure in Eq. 22 are equal to unity and M ¼ M 2 . Hence, influence of the microstructure is significant only in the case of jointed invariants tr ðr Á MÞ and tr ðr 2 Á MÞ. To denote functions of both stress and jointed invariants, we apply the overline symbol ðÞ. An elastic potential function may be expressed in the following form: with parameters n i as the additional scalar multipliers. The basic elastic potential without inherent cross-anisotropy in Eq. 11 is a function of only one stress invariant Q. An analogous joint stress-microstructure invariant Q M may be defined as: As in Eq. 24, we can group the invariants Q and Q M to obtain the following mixed invariant: where c 1 and c 2 are the new material constants. Finally, the incorporation of inherent cross-anisotropy into the elastic potential is attained by simply replacing the stress invariant Q in Eq. 11 with its mixed equivalent Q: Constant c 1 controls the intensity of pure stress dependency of stiffness. Constant c 2 introduces inherent cross-anisotropic stiffness component with a fixed orientation of the symmetry axis, collinear with vector v. In the case of c 1 ¼ 1:0 and c 2 ¼ 0:0, inherent cross-anisotropic component is deactivated and the basic potential function from Eq. 11 is recovered. As a rule, when c 1 is kept constant and c 2 [ 0:0, inherent stiffness in the plane of isotropy is higher than in the direction of symmetry axis, i.e. a G ; a E [ 1:0. Otherwise, when c 1 ¼ const and c 2 \0:0, inherent cross-anisotropy coefficients are a G ; a E \1:0. Practically, the influence of parameter c 1 is a proportional scaling of the stiffness. The influence of c 1 on the directional distribution of stiffness is minor and its intensity depends on the current value of c 2 . In order to simplify the presentation of the new model, we prefer to keep c 1 ¼ 1:0 in the simulations and to control the directional properties of stiffness only by parameter c 2 . We address the problem of parameter identification later in this paper.
In derivation of the secant hyperelastic strain-stress relation and the tangent compliance tensor, we need the partial derivative of the mixed invariant Q: The secant strain-stress relation is obtained as: with The influence of c 2 value controlling inherent cross-anisotropic stiffness component may be observed by a qualitative analysis of plots showing the elastic strain distribution from Eq. 29. It is presented in Fig. 5. The influence of c 2 value is clear when comparing stress Another important feature of the proposed model, which can be illustrated on the basis of secant strain-stress relation (Eq. 29), is the shape of the constant volumetric strain e e V and the shear strain e e q contours. They are usually presented in the p À q triaxial plane, e.g. [2,28]. Such volumetric and shear strain contours for the proposed model are presented in Fig. 6.
It should be noted that e e V ¼ const contours correspond to the undrained triaxial compression stress paths for elastic behaviour. Hence, the imposed inherent cross-anisotropy using parameter c 2 influences direction of the undrained stress paths. Consequently, if the proposed hyperelastic description is a part of an elasto-plastic model, it will impact the undrained shear strength.
Finally, the tangent compliance of the hyperelastic model with stress-induced anisotropy and inherent crossanisotropy is obtained from the following differentiation: where and The tangent stiffness matrix can be derived analytically or numerically by inversion of the tangent compliance tensor in the Voigt notation. For details of such operations, the reader is referred to [41]. It should be noted that the derivation of hyperelastic stiffness via the inversion of compliance matrix, applied in numerical calculations presented in this paper, is an alternative way to the classic and more elegant method presented in, for example, [29] or [44]. Elastic potential expressed as a function of elastic strain may be obtained from stress functions (11) or (27) by Legendre transform. Finally, the hyperelastic stiffness may be obtained directly from differentiation of the transformed potential with respect to the elastic strain.
3 Features of the model and material respectively. The initial isotropic stress condition is r 0 ¼ Àp 0 d and the cross-anisotropic microstructure is represented by tensor M as defined in Eq. 20. Components of the incremental elastic strain response to the stress increments Eqs. 34 and 35 are calculated using tangent compliance tensor (Eq. 31): De e ¼ C t ðr 0 ; MÞ : Dr: The following formulas for the inherent cross-anisotropic stiffness moduli are obtained: Using the elastic strain responses to the stress increments (Eq. 35), we define the cross-anisotropic Poisson's ratios: Now, the anisotropy coefficients (Eqs. [3][4][5] can be derived for the proposed model: In practice, it is feasible to determine the following pairs of the cross-anisotropic parameters in triaxial tests at a reference isotropic stress state represented by p 0 ¼ p ref : • G ref vh , a G , using vertical and horizontal pairs of bender elements; • E ref v , a G , using bender elements as above in drained compression with high accuracy axial strain measurements; • E ref uv , a G , using bender elements as above in undrained compression with high accuracy axial strain measurements.
In order to simplify estimation of the model parameters, it is proposed to fix the constant value of c 1 ¼ 1:0 and to adjust c 2 with a G using Eq. 44: In undrained triaxial compression tests, the small strain undrained Young's modulus in direction of the symmetry axis E ref uv can be measured [20]. In order to derive E ref uv , we have to inspect the model response to the following elastic strain increment:

5: ð49Þ
Components of the incremental elastic stress response under isotropic stress condition r 0 ¼ Àp 0 d are calculated using the tangent stiffness tensor: and the undrained Young's modulus can be calculated from: After inversion of the above equation, assuming p 0 ¼ p ref , can be obtained as follows:

Relations between inherent anisotropy coefficients
Combining the relations for anisotropy coefficients (Eqs. [43][44][45], it is possible to define a E and a m as functions of parameters a G and b: The obtained relations, unlike in Eqs. 6 and 7, are not exponential and the corresponding anisotropy exponents x GE and x Gm depend on the values of a G and b. In Fig. 7, a graphical comparison of relations (Eqs. 53, 54) with exponential functions is presented for b ¼ 0:5, which implies order of stress dependency m ¼ 0:5 often assumed in overconsolidated clays and sands. Another graphical comparison concerns relations between anisotropy coefficients for different values of b as presented in Fig. 8. Since a G is independent of b, we have examined changes of a E and a m with b for two opposite cases of inherent crossanisotropy, i.e. higher stiffness along the symmetry axis (a G ¼ 0:6) and higher stiffness in the plane of isotropy (a G ¼ 2:0). Keeping the same range of b, we have also shown changes of exponent x GE for seven selected values of a G between 0.6 and 3.0 as presented in Fig. 9.

Mixed anisotropy at axisymmetric and true triaxial stress conditions
Changes of stiffness resulting from the superposition of variable stress-induced anisotropy and constant inherent cross-anisotropy may be illustrated qualitatively by comparison of response envelopes elaborated for different initial stress ratios K and selected values of anisotropy coefficient a G . Such response envelopes at axisymmetric stress conditions are presented in Fig. 10.  Fig. 10 that for a G \1:0, the radial stress response to purely Fig. 7 Relations between anisotropy coefficients a E ; a m and a G reproduced by the model for b ¼ 0:5. In the case of relation a E ða G Þ, it is close to findings presented by Mašín and Rott [38]. The relation a m ða G Þ deviates from a m ¼ a G as assumed in the cross-anisotropic model by Graham and Houlsby [21]. However, experimental evidence in this case is still very scattered [38] Fig. 8 Relations between anisotropy coefficients a E ; a m and a G reproduced by the model for two opposite cases of inherent cross-anisotropy (a G ¼ 0:6 and a G ¼ 2:0) and different values of parameter b volumetric strain probes rotates anticlockwise and the radial stress response to purely deviatoric strain probes rotates clockwise relative to the response obtained for a G ¼ 1:0. As a result of such behaviour, the response envelopes seem to rotate anticlockwise, if a G \1:0. If a G [ 1:0, the rotations occur in the opposite direction.
Definitions of the anisotropy coefficients may be extended to account for axisymmetric stress conditions. This allows to investigate the mixed anisotropy quantitatively and to show how the stress obliquity K influences the resulting stiffness. For this purpose, one needs to inspect the model response to stress increments (Eqs. 34,35) starting from the following initial axisymmetric stress state:  . 36), we obtain the same relation for anisotropy coefficient a G as in Eq. 44. It is an important feature of the model that a G is constant and independent of the stress obliquity at axisymmetric stress conditions. Analogously, it can be shown that a G is also constant at a true triaxial stress state, i.e. for three different normal stress components and zero shear stress components. Hence, in the extended definitions of cross-anisotropic moduli. we use a G assuming c 1 ¼ 1:0 and (Eq. 46). The following formulas for the extended cross-anisotropic stiffness moduli are obtained: ð56Þ where Regarding the cross-anisotropic Poisson's ratios at axisymmetric stress conditions, the following formulas are obtained: where a m ða G ; b; KÞ ¼ The above extended definitions of anisotropy coefficients can be used to show changes of the mixed anisotropy with the stress ratio K and coefficient a G at axisymmetric stress states. This is presented for b ¼ 0:5 in Fig. 11. At the isotropic stress (K ¼ 1:0), anisotropy coefficients represent the inherent values from Eq. 53 and 54. Deviation from the isotropic stress (K 6 ¼ 1:0) results in a proportional increase or decrease in the value of anisotropy coefficients.
In real soils, the mixed anisotropy of initial stiffness may be more complex. All principal stress components may have different values with the principal axes rotated with respect to the geometrical axes. Moreover, the plane of isotropy may be non-horizontal. To visualise stiffness obtained with the model for such conditions, an orientation distribution function . À1C related to the fourth-order tangent compliance tensor can be used [52]: . À1C ðC t ; nÞ ¼ ðn i n j C t ijkl n k n l Þ À1 : It reflects a scalar value which is equal to the Young's modulus in the direction indicated by the unit vector n.
After parametrisation of n in terms of spherical coordinates, as in Eq. 18, spherical plots of . À1C can be produced to illustrate the mixed anisotropy of the proposed hyperelastic stiffness. This is shown in Fig. 12 where cases (a)-(d) present the growing complexity of stress and microstructural conditions in sequence. In Fig. 12a, both stress and microstructure are isotropic and Young's modulus distribution forms perfect sphere like in the isotropic Hooke's stiffness. In Fig. 12b, example true triaxial stress state is introduced; however, principal stress axes are collinear with the geometrical axes and the microstructure is still isotropic. The distribution forms ellipsoid reflecting stress-induced anisotropy from the basic hyperelastic Vermeer's stiffness in Eq. 15. The axis of symmetry v is shown vertical, but its orientation has no influence because a G ¼ 1:0. Inherent cross-anisotropy is introduced in Fig. 12c by imposing a G ¼ 2:0. This results in an increase of Young's modulus in all horizontal directions and the distribution diverges from a regular ellipsoid. Note that geometrical, microstructural and stress axes are still collinear and the plane of isotropy is horizontal. Rotation of the plane of isotropy representing the inherent microstructure is introduced in Fig. 12d. In all presented spherical plots, The obtained results are presented for sands and clays separately. In all referred element tests, the plane of isotropy in Fig. 11 Changes of anisotropy coefficients a E and a m with the stress ratio K at axisymmetric stress conditions the cross-anisotropic description is horizontal excluding the case of Opalinus Clay. For details on the considered tests, the reader is referred to the original papers.
As the main calibration tool to obtain the values of model parameters, FindMinimum function was used that is available in the algebra program MATHEMATICA [53]. The minimum of the sum of squared differences between the experimental and the model results was searched.

Toyoura Sand
Data on Toyoura Sand [27] involve values of E v and E h . Small strain static cyclic loading in drained triaxial tests was conducted on prismatic samples. Values of the shear moduli are not provided.
Triaxial compression with r h ¼ const was performed on the considered reconstituted sample of Toyoura Sand, starting from the isotropic stress state at p ¼ 100 kPa up to the stress ratio K ¼ 0:29. The sample was unloaded to the initial isotropic state along the same stress path and triaxial extension with r v ¼ const up to K ¼ 3:5 followed. Cyclic loading was performed in the vertical and horizontal directions at isotropic and various anisotropic stress states during both triaxial compression and extension.
Experimental results and simulations with the calibrated model parameters are shown in Fig. 13. Values of E v and E h in the left hand side of Fig. 13 were divided by the void ratio function f ðeÞ ¼ ð2:17 À eÞ 2 =ð1 þ eÞ [24]. In the cases of K 0:4 and K ! 2:5, a decrease in the experimental values of the Young's moduli can be seen, which corresponds to the excess of the elastic range. Hence, these data were not considered during calibration of the model parameters. Values of the model parameters are shown in the right hand side plot.

Ham River Sand
Static cyclic loading and bender element triaxial tests were conducted on reconstituted Ham River Sand samples [34]. The results provide values of E v , E h , G vh and G hh for various testing conditions.
Considered are the results of drained tests. Samples were normally consolidated under K ¼ 0:45. The experimental data obtained for different values of p and the

London Clay
Data on the small strain stiffness parameters were obtained with the cyclic loading and bender element triaxial tests performed on the intact samples of London Clay [19].
Samples taken from different depths were consolidated to the in situ stress states. Both the drained Young's moduli and the shear moduli investigated at three different levels of p, representing the in situ stress, are regarded in this study. Each of the three different values of p corresponds to different value of the stress ratio K.
The experimental values of E v , E h , G vh and G hh and the simulations are shown in Fig. 15. Stress ratios K corresponding to the three different stress levels and values of the model parameters are all given in the right hand side plot. Points corresponding to each pair of p and K are denoted in the simulated plots with the bigger indicators.

Gault, Kimmeridge, Oxford and London Clays
Experimental results referred in [8] provide data on E v , E h , G vh and G hh for four different clays. Elastic stiffness parameters were investigated on both intact and reconstituted samples of Gault, Kimmeridge, Oxford and London Clays with the drained static probes and bender element triaxial tests.
All samples were tested under nearly in situ stress conditions. Values of the small strain stiffness parameters obtained for each clay and the simulated curves are plotted in Fig. 16. Stress ratios K corresponding to different stress

Opalinus Clay
Anisotropic behaviour of Opalinus Clay was investigated on intact samples tested perpendicular (S-specimens) and parallel (P-specimens) to bedding plane in the triaxial apparatus [14,48]. Additionally, data on samples with bedding orientation inclined at 45 with respect to the loading direction (Z-specimens) are provided in [48]. The resolution of drained triaxial compression curves (e 1 À q) presented in the referenced reports is not sufficient to focus on the details of small strain behaviour. However, the qualitative simulation of initial stiffness for P, Z and S specimens could be conducted and only the initial linear segments of compression curves are considered.
In the cross-anisotropic material description of P-specimens and Z-specimens, the symmetry axis and the plane of isotropy is not collinear with vertical and horizontal directions, respectively. Hence, one needs to distinguish between the principal axes of stress and of microstructure.
Samples reported in [14] were consolidated isotropically to three different levels of p, i.e. 2, 5, 10 MPa and 2, 5, 12 MPa, in the case of P-and S-specimens, respectively. This was followed by the drained triaxial compression. Experimental and simulated relations between deviatoric stress q and axial strain e 1 for both P-and S-specimens are shown in Fig. 17. Solely the initial small strain elastic range is considered. Values of the model parameters are shown in the upper right hand side plot. In Fig. 17, only one Pspecimen is considered, on which a multistage test was conducted. Due to lack of the experimental data for p ¼ 10 MPa and p ¼ 12 MPa in the case of S-and P-specimen, respectively, only the simulated relation is plotted.
In Fig. 18, experimental and simulated values of the corresponding Young's moduli and a E are presented. Since the samples were tested parallel and perpendicular to bedding, vertical Young's moduli E 1 can be interpreted as E v and E h in the case of S-specimens and P-specimens, respectively. The modelled data were calculated with the values of parameters shown in the upper right hand side of Fig. 17. As already mentioned, the simulations were conducted only for the small strain elastic region. However, the experimental Young's moduli in Fig. 18 concern Fig. 17 Experimental [14] and simulated relation between q and e 1 for Opalinus Clay. Data are presented for P-specimen and S-specimens under different values of r 2 ¼ r 3 . Only the initial linear segments of compression curves are shown due to limited accuracy of the tests within the small strain region Fig. 18 Experimental [14] and simulated elastic stiffness parameters for Opalinus Clay. Left hand side: E 1 plotted for different values of p equal to confining stress. E 1 obtained for P-specimens and S-specimens corresponds to E h and E v , respectively. Right hand side: a E plotted for different values of p different ranges of strain than shown in Fig. 17. These were derived from the linear part of q À e 1 unloading-reloading loops for S-specimens and from the initial part of q À e 1 plot up to 25% of the peak shear strength for P-specimens.
Results analogous to these in Fig. 17 but regarding also Z-specimens are shown in Fig. 19 for selected data from [48]. Same set of model parameters as in Fig. 17 was used in all simulations considering Opalinus Clay.

Conclusions
A simple anisotropic hyperelastic model is presented in this paper. It is a refined version of the basic hyperelastic model proposed by Vermeer [51]. The refinement concerns incorporation of inherent cross-anisotropy with the mixed stress-microstructure invariant. Properties of the obtained anisotropic hyperelastic stiffness are investigated and compared with laboratory results from literature. The current experimental evidence is still limited in many aspects of small strain anisotropic behaviour of soils, and further investigations are needed. However, the proposed model is capable of a robust description of the pure inherent crossanisotropy and the mixed anisotropy observed in various tests conducted on both sands and clays. Due to the small number of parameters that can be easily related to laboratory results, the model proves to be a simple tool allowing reproduction of the small strain mixed anisotropy, being the superposition of stress-induced and inherent anisotropy as observed in natural soils.
The main application of the model can be its incorporation into some more advanced hyperelastic-plastic models to properly simulate the initial small strain stiffness and the pre-failure undrained behaviour. The further works related to the proposed model are considered, and they will concern its FE implementation within existing elastoplastic models by substitution of the popular isotropic hypoelastic formulation, e.g. [12].
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://creativecommons. s; s ij : Deviatoric part of the stress tensor, s ij ¼ r ij þ p d ij ; v; v i : Unit vector defining the symmetry axis in a cross-anisotropic material; W: Elastic potential function, Wðe e Þ-Helmholtz free energy function WðrÞ-negative Gibbs free energy function; x GE ; x Gm : Anisotropy exponents, a G ¼ a xGE E , a G ¼ a xGm m ; a: Anisotropy coefficient in Graham and Houlsby model, Functions of both stress and jointed stress-microstructure invariants; Á: Single contraction, a Á b ¼ a i b i or A Á b ¼ A ij b j ; : : Double contraction, D : r ¼ D ijkl r kl or r : r ¼ r ij r ij ; : Dyadic or outer product, A ¼ a b ¼ a i b j ¼ A ij ; k k: Euclidean norm, kxk ¼ ffiffiffiffiffiffiffi x i x i p or krk ¼ ffiffiffiffiffiffiffiffiffiffi r ij r ij p ; tr ðÞ: Trace of a tensor, tr r ¼ r ii ; ðÞ n : Power of a tensor, r n ¼ r Á r Á . . . Á r zfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflffl{ n or tr r 3 ¼ tr ðr ik r kl r lj Þ ¼ r ik r kl r li Fig. 19 Experimental [48] and simulated relation between q and e 1 for Opalinus Clay. Data are presented for P-specimens, Z-specimens and S-specimens. Only the initial linear segments of compression curves are shown due to limited accuracy of the tests within the small strain region