High-pressure behaviour and elastic constants of 1M and 2M₁ polytypes of phlogopite KMg₃Si₃AlO₁₀(OH)₂

In the present work, the elastic properties of both 1M and 2M₁ phlogopite polytypes, KMg₃Si₃AlO₁₀(OH)₂ (monoclinic crystal system) were investigated from PV equation of state fitting and by analysis of the fourth-rank elastic tensor. The analysis was performed within the Density Functional Theory framework, using all-electron Gaussian-type orbitals basis sets and the B3LYP functional corrected a posteriori to include long-range interactions (B3LYP-D*). In general, the elastic properties of the two polytypes were strongly anisotropic, with the axial moduli ratio M(a) : M(b) : M(c) being close to $\mathrm{4}:\mathrm{4}:\mathrm{1}$. The volume-integrated third-order Birch-Murnaghan equation of state fitting parameters at 0 K were K₀=57.9(2) GPa, $K^{\prime }=\mathrm{8.29}\left(\mathrm{7}\right)$ and V₀=489.82(3) ų for phlogopite-1M, which were very close to those of the 2M₁ polytype, i.e., K₀=58.3(1) GPa, $K^{\prime }=\mathrm{8.71}(\mathrm{8}$) and V₀=978.96(9) ų. The monoclinic elastic tensors obtained for the two polytypes of phlogopite, which have never been experimentally reported for both minerals so far, were in line with the PV behaviour of the mineral, providing further data related to the directional dependence of the elastic properties and seismic wave propagation. The elastic properties from both PV hydrostatic compression and from the elastic moduli tensor were discussed against the available experimental and theoretical data in the scientific literature, extending the knowledge on this important trioctahedral phyllosilicate.

1 Introduction

Phlogopite [commonly abbreviated with Phl, K(Mg,Fe)₃Si₃AlO₁₀(OH)₂, with Mg $/$ Fe ratio greater than 2], is a ubiquitous mica mineral that can be found in different geological environments and pressure-temperature conditions, i.e. igneous, sedimentary, and metamorphic rocks (Icenhower and London, 1995; Trønnes, 2002), and one of the water and potassium reservoirs in the Earth's mantle (Virgo and Popp, 2000; Tutti et al., 2000).

Figure 1 (a) Structural model of the phlogopite-2M₁ model viewed along the [100] direction. In the bottom right corner, a graphical representation of the tetrahedral rotation angle α_r is shown. (b) Schematic view of the stacking of the TOT layers in the 1M (left) and 2M₁ (right) polytypes. The black arrows link two points related by symmetry.

From the crystallographic point of view, like other phyllosilicates, phlogopite is made of a 2:1 layer formed by two tetrahedral sheets (labelled as T) that sandwich an octahedral one (O). These T-O-T (or TOT) layers present strong (covalent/ionic) in-plane bonds and a negative charge because of Al^III/Si^IV substitutions in the T sheets, which is balanced by potassium cations in the interlayer. Thus, the structure is held together along the c crystallographic axis by a mix of Coulombic (electrostatic) and weak (van der Waals) interactions (Ventruti et al., 2009). Phlogopite presents different polytypes due to different arrangement of the TOT layers along the [001] direction, with the 1M and the 2M₁ being the most common in nature both belonging to the monoclinic crystal system (Lacalamita et al., 2012), while a third polytype, 3T (trigonal), is known to be very rare (Gatta et al., 2011). A graphical representation of the most common polytypes is shown in Fig. 1.

Phlogopite is also a widely adopted mineral in the technological field, for instance in the manufacture of ceramics and glasses (King et al., 2000; Ariane et al., 2023) and, more recently, gained attention as a potential dielectric (insulating) material for optical and electronic applications because it presents a band gap of ca. 5–6 eV that can be modulated by microchemical variations, e.g., by Fe^II/Mg^II substitutions (Ulian and Valdrè, 2023a).

Micas, and then phlogopite, could be used as geothermobarometers to explain petrogenetic processes occurring in high-pressure and high-temperature (non-ambient) conditions (Guidotti and Sassi, 1976), being also a possible carriers of water and potassium in the Earth's mantle (Sudo and Tatsumi, 1990). As suggested by Melzer and Wunder (2001), since phlogopite may be formed in the mantle wedge that overlies a subducting slab, it could play a relevant role in the element ratios of large ions (K, Rb, Cs) in the subduction zone. Also, the bulk chemical composition of altered rocks such as mica-amphibole-rutile-ilmenite-diopside (MARID) is very similar to the chemistry of phlogopite in an olivine – quartz – alkali ternary [(Mg,Fe)₂SiO₄–SiO₂(quartz)–K₂O+Na₂O] (alkali) ternary (Sweeney et al., 1993). From the technological perspective, phyllosilicates are considered low-cost insulating materials that could be used in two-dimensional optoelectronic applications, e.g., 2D transistors and heterostructures (Ulian and Valdrè, 2025), thanks to their easy exfoliation, which is related to the perfect cleavage on their (001) planes. Thus, detailed knowledge of the elasticity of phlogopite is of utmost importance for both geological studies, e.g., to explain geophysical observations involving phlogopite-bearing rocks (Geng et al., 2024; Van Reenen et al., 2023; Wang et al., 2023), and to devise new technological uses of this mineral for advanced electronic and optical applications (Cadore et al., 2022). Many experimental studies were focused on the equation of state of phlogopite (Pavese et al., 2003; Comodi et al., 2004; Gatta et al., 2011), but mainly on the 1M polytype. As also observed by Chheda et al. (2014), no measurement of the full monoclinic elastic constant tensor has ever been reported in the scientific literature. The only experimental evidence of this elastic property is related to elastic wave propagation (EWP) experiments that provided just the 5 independent constants, considering the mineral as pseudo-hexagonal (Aleksandrov et al., 1974; Alexandrov and Ryzhova, 1961). The full elastic tensor of phlogopite-1M was only recently provided from ab initio simulations within Density Functional Theory (DFT) at the local density approximation (LDA) and generalized gradient approximation (GGA) levels (Chheda et al., 2014). However, the authors did not employ any correction to include long-range interactions in the physical treatment of phlogopite. Albeit the general consensus is that the TOT layers in micas are held together by the interlayer cations (K⁺), hence the interaction is mostly of Coulombic (electrostatic) nature, previous works on muscovite and phlogopite showed that the contribution from van der Waals interactions is not negligible in determining the crystal structure, elastic and thermodynamic properties of these minerals (Ulian and Valdrè, 2015a; Ulian and Valdrè, 2023a). Furthermore, no studies on the elastic properties of phlogopite-2M₁ were ever carried out and reported. The differences in the elastic behaviour of the two phlogopite polytypes could be indeed subtle, especially in terms of the second-order elastic moduli tensor, where small variations in the directional properties between the 1M and 2M₁ are hypothesised.

For all these reasons, the present work is focused on the elastic properties of phlogopite, considering both 1M and 2M₁ polytypes. To this aim, ab initio DFT simulations were performed (i) to understand the compressional behaviour of the mineral through an equation of state fitting, (ii) to obtain the full monoclinic elastic moduli and (iii) to assess the role of polytypism on the elastic behaviour of phlogopite. All the data were compared and discussed with both experimental and theoretical results found in the scientific literature to extend the knowledge on this important phyllosilicate.

2 Computational Methods

As previously introduced, the simulations reported in the present work were carried out within the Density Functional Theory framework. In detail, the hybrid B3LYP functional (Becke, 1993; Lee et al., 1988), which includes 20 % of exact Hartree-Fock exchange and some non-local contribution to the exchange-correlation terms, was employed throughout the present study. B3LYP is considered a suitable choice when dealing with the simulation of the structural and elastic properties of minerals due to its high accuracy, and it was already adopted with success for the description of several other minerals (Prencipe et al., 2009; Ottonello et al., 2010; Belmonte, 2017; Ulian and Valdrè, 2023b, 2019) and phlogopite (Ulian and Valdrè, 2023a). However, since most generalized-gradient approximation functionals, including their hybrid counterparts, do not adequately include van der Waals (long-range) interactions, the DFT-D2 correction (Grimme, 2006) was adopted, employing the parameters of Civalleri and co-workers (2008) for the B3LYP functional (B3LYP-D* method). According to the proposed methodology, the original B3LYP-D2 parameters were modified as s₆=1, the van der Waals radius of the hydrogen atom, R_vdw(H), was set to 1.30, and the R_vdW of the heavier atoms were scaled by a factor of 1.05, while the C₆ values were the same proposed in the DFT-D2 scheme (see Table S1 in the Supplement). This approach was used in previous studies on the high-pressure behaviour of talc (Ulian et al., 2014, Ulian and Valdrè, 2015c), pyrophillite (Ulian and Valdrè, 2015b), muscovite (Ulian and Valdrè, 2015a), resulting in data in very good agreement with experimental data. The all-electron Gaussian-type orbitals (GTO) used to construct the Kohn-Sham orbitals within the linear combination of atomic orbitals (LCAO) approach were the same used in a recent work on the 1M polytype of phlogopite (Ulian and Valdrè, 2023a). Si, Al, Mg, K, O and H atoms were described with 88-31G* (Nada et al., 1996), 85-11G* (Catti et al., 1994), 8-511d1G (Valenzano et al., 2007), 86-511G (Dovesi et al., 1991), 8-411d11G (Valenzano et al., 2006), and 3-1p1G (Gatti et al., 1994), respectively.

The accuracy of the Coulomb and exchange series was controlled by five thresholds set to 10⁻⁸ (ITOL1 to ITOL4) and 10⁻¹⁶ (ITOL5) for the structural relaxation procedures, whereas the convergence on the total energy was set to 10⁻⁸ Hartree. The reciprocal space was sampled with a $\mathrm{5}\times \mathrm{5}\times \mathrm{5}$ Monkhorst–Pack mesh (Monkhorst and Pack, 1976), corresponding to 39 independent k points.

The lattice and internal geometry were optimised within the same run using a numerical gradient for the unit cell parameters and an analytical gradient for the atomic positions. The Hessian matrix was upgraded with the BFGS algorithm (Broyden, 1970b, a; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970). The tolerances for the maximum allowed gradient and the maximum atomic displacement were set to 10⁻⁵ Ha bohr⁻¹ and $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{5}}$ bohr, respectively.

All the simulations were performed with the CRYSTAL17 code (Dovesi et al., 2018), whereas the QUANTAS software (Ulian and Valdrè, 2022, 2024a) was employed to post-process the elastic data. Graphical representations have been carried out with the molecular graphics program VESTA (Momma and Izumi, 2008).

The initial model of the crystal structure of 2M₁- phlogopite was created from the single-crystal X-ray diffraction data of Lacalamita et al. (2012), which was subsequently optimised within the B3LYP-D* level of theory. The mineral belongs to the monoclinic system (space group C2/c), presenting two independent tetrahedral sites (labelled as T) and two octahedral sites (M), with partial occupancies of the cationic sites with several possible atoms. For example, the TO₄ tetrahedral sites in the unit cell are generally occupied by about 70 % by Si and 30 % by Al (Laurora et al., 2007; Lacalamita et al., 2012), which can be translated in $\mathrm{1}/\mathrm{4}$ sites being randomly occupied by aluminium. However, as also explained in a previous study on the 1M polytype (Ulian and Valdrè, 2023a), theoretical modelling requires each site to be deterministically occupied by a single kind of atom; in the example cited above, a T site could present either Al or Si but not both at the same time. Thus, to preserve the monoclinic lattice of phlogopite when the Al^III/Si^IV substitutions are included in the structure, the 2M₁ polytype with formula KMg₃(AlSi₃)O₁₀(OH)₂ was described with the P2/c space group (23 inequivalent atoms, 88 atoms in the unit cell, Z=4 unit formulas), which is a sub-group of the C2/c symmetry. Within the P2/c space group, phlogopite-2M₁ presented four inequivalent T sites (Si^IV in T1–T3, Al^III in T4) and three non-equivalent O sites (M1–M3). Hence, each T site was related to four SiO₄ or AlO₄ tetrahedra, whereas all octahedral sites are related to two MgO₆ octahedra.

3 Results and discussion

3.1 2M₁-Phlogopite crystal structure

Table 1 reports the lattice parameters, tetrahedral (T_thick) and octahedral (M_thick) sheet thicknesses, interlayer distance (I_thick), and selected polyhedral properties of the 2M₁ phlogopite polytype, obtained from DFT simulations with (B3LYP-D*) and without (B3LYP) the correction to include the effects of long-range interactions. Previous XRD refinements are also shown in Table 1 for a direct comparison.

Table 1 Equilibrium geometry of 2M₁-phlogopite as obtained from DFT/B3LYP and B3LYP-D*, compared to previous experimental refinements.

^a X-ray diffraction data of Laurora et al. (2007). ^b single-crystal XRD refinement of Lacalamita et al. (2012) performed on sample BU1_14. T_thick, M_thick, and I_thick are tetrahedral sheet thickness (calculated from the z coordinates of the basal and apical oxygen atoms), the octahedral sheet thickness (calculated from the z coordinates of the apical and hydroxyl O atoms) and interlayer thickness (calculated from the z coordinates of the basal oxygen atoms), respectively. TQE and OQE are the tetrahedral and octahedral quadratic elongations, respectively (Robinson et al., 1971). T and M represent the generic tetrahedral and octahedral sites, respectively.

The simulations at the DFT/B3LYP-D* level of theory are in good agreement with the experimental crystallographic data reported in literature, showing a small underestimation of the unit cell volume (ΔV $=-\mathrm{1.0}\mathit{\%}$), which is the result of a slight expansion of the tetrahedral and octahedral sheet thicknesses (ΔT_thick=1.5 % and ΔO_thick =1.8 %, respectively) and a contraction of the interlayer distance (ΔI${}_{\mathrm{thick}}=-\mathrm{2.5}\mathit{\%}$). These results for the 2M₁-phlogopite are in agreement with previous simulations of the 1M polytype (Ulian and Valdrè, 2023a) and other phyllosilicates and layered minerals (Ulian et al., 2013; Ulian and Valdrè, 2015a, b, 2019), where the inclusion of long-range interactions in the physical treatment is of utmost importance to properly simulate the crystal-chemistry and properties of this kind of structure. When the DFT-D2 correction is not implemented in the simulation framework, a dramatic increase in the unit cell volume can be noted (+6 % with respect to the DFT/B3LYP-D* results, +5 % from the XRD data), which is the result of the increase of the lattice parameters a (+0.9 %), b (+0.9 %), and c (+3.6 %) and the shrinking of the β angle (−2.1 %). Considering the absence of any thermal effects in the Density Functional Theory simulations, the overestimation of the lattice constants and unit cell volume is a further sign of the relevant role played by van der Waals interactions in micas.

3.2 Compressional behaviour of phlogopite 1M and 2M₁ polytypes

The compressional behaviour of phlogopite was modelled considering 10 different unit-cell volumes, both smaller (compressed, 7 volumes) and larger (expanded, 3 volumes), within 88 % and 108 % of the equilibrium geometry volume (V_eq). Then, the internal coordinates and lattice parameters were optimized keeping the selected volume fixed during the procedure, constraining the space group of 1M-phlogopite (P2) and 2M₁-phlogopite (P2/c), a procedure known as “symmetry preserving, variable cell-shape structure relaxation” that was proposed by Pfrommer et al. (1997). The equilibrium geometry of 1M-phlogopite was retrieved from the work of Ulian and Valdrè (2023a) that was obtained with the same computational settings here adopted for the 2M₁ polytype.

Figure 2 Evolution of normalised (a) unit cell volume $V/V_{\mathrm{0}}$, (b) $a/a_{\mathrm{0}}$, (c), $b/b_{\mathrm{0}}$ and (d) $c/c_{\mathrm{0}}$ as a function of pressure for the phlogopite polytypes 1M and 2M₁. The experimental data of Pavese et al. (2003) and Comodi et al. (2004) are reported for a comparison.

The dependence of the total energy E of the mineral at the different volumes V, i.e., the E(V) curve, was described in terms of a volume-integrated 3rd-order Birch-Murnaghan equation of state, BM3 (Birch, 1947), as reported by Hebbache and Zemzemi (2004):

$$\begin{array}{}\text{(1)}& {\displaystyle }E\left(V\right)=E_{\mathrm{0}}+{\displaystyle \frac{\mathrm{9}}{\mathrm{16}}}{K}_{\mathrm{0}}{V}_{\mathrm{0}}\left\{K{}^{\prime }{\left({\mathit{\eta }}^{\mathrm{2}}-\mathrm{1}\right)}^{\mathrm{3}}+\left[{\left({\mathit{\eta }}^{\mathrm{2}}-\mathrm{1}\right)}^{\mathrm{2}}\left(\mathrm{6}-\mathrm{4}{\mathit{\eta }}^{\mathrm{2}}\right)\right]\right\}\text{(2)}& {\displaystyle }\mathit{\eta }={\left({\displaystyle \frac{V_{\mathrm{0}}}{V}}\right)}^{\mathrm{1}/\mathrm{3}}\end{array}$$

where η is a dimensionless parameter, K₀ is the bulk modulus at 0 K, K^′ its pressure first derivative and V₀ the volume at zero pressure. The pressure values at each unit cell volume (reported in Tables S2 and S3 for the 1M and 2M1 polytypes, respectively) were calculated using the fitting parameters K₀, K^′ and V₀ obtained from the previous equation and employed in the well-known P-V formulation of the BM3 (Birch, 1947):

$$\begin{array}{}\text{(3)}& P\left(V\right)=P_{\mathrm{0}}+{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{K}_{\mathrm{0}}\left({\mathit{\eta }}^{\mathrm{7}}-{\mathit{\eta }}^{\mathrm{5}}\right)\left\{\mathrm{1}-{\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}\left(\mathrm{4}-K{}^{\prime }\right)\left({\mathit{\eta }}^{\mathrm{2}}-\mathrm{1}\right)\right\}\end{array}$$

with P₀ the reference pressure (0 GPa). In Table 2, the BM3 equation of state parameters for both 1M and 2M₁ phlogopite are reported. A graphical representation of the normalised volume and axis lengths as a function of calculated pressure is presented in Fig. 2., alongside the experimental results from X-ray diffraction (Comodi et al., 2004; Pavese et al., 2003) and the theoretical calculations of Chheda et al. (2014) carried out at the LDA and GGA level of theory.

A direct comparison between the different experimental and computational approaches is however possible on the phlogopite-1M, which is better depicted in Fig. S1 (Supplement) considering the absolute values of the cell volume and lattice vectors. In general, compared to the high-pressure XRD refinements, the crystal structure from the DFT/B3LYP-D* approach well agrees with the results of Pavese et al. (2003) at low pressure (P<5 GPa), however there is a slight overestimation of the unit cell volume at high pressure (see Figs. 2a and S1a). Analysing the lattice vectors, it can be noted that the a and b lengths are very close to the experimental ones at each pressure state, whereas the c-axis shows a stiffer behaviour and an overestimation that increases with pressure. By considering the results of Comodi et al. (2004), the present B3LYP-D* data are slightly underestimated in terms of unit cell volume, which is due to the lower values of the a- and b-axes, whereas the evolution of the clattice parameters is correctly described by the proposed computational approach (see Fig. S2b–d). The stiffer behaviour of both polytypes at higher pressures elucidated by the $V/V_{\mathrm{0}}$ trend is related to the static conditions of the simulations, i.e., the zero-point and thermal effects were neglected. It is suggested that the use of the quasi-harmonic approximation, which is however beyond the scope of the present work, could reduce the mismatch between the theoretical and experimental data.

It is interesting to note that if the long-range interactions are not included in the theoretical framework, the B3LYP functional leads to a general overestimation of the unit cell volume and lattice parameters. The a and b vectors are closer to the results of Comodi and collaborators (2004) than those of Pavese et al. (2003), however, the c lattice parameter is severely overestimated when compared to both XRD results. All these considerations are reflected in the PV fitting using the BM3 equation of state, where the B3LYP-D* approach resulted in a small overestimation of the bulk modulus K₀ when compared to the experimental data (see Table 2), and the theoretical K' and V0 values between the ranges obtained from the XRD refinements. Conversely, the B3LYP data shows a large underestimation of the bulk modulus of about 18 GPa, and an overestimation of the cell volume.

Table 2 Equation of state parameters and associated standard deviations (σ) obtained for the 1M and 2M₁ polytypes of phlogopite.

Notes: GTO – Gaussian-type orbitals; PAW – projector-augmented wave; SCXRD – single-crystal X-ray diffraction; SXRPD – synchrotron X-ray powder diffraction.

Compared to previous theoretical results, the crystal structure of the phlogopite 1M polytype and the related PV fitting obtained from the uncorrected B3LYP functional are almost identical to the GGA values of Chheda et al. (2014). Instead, the B3LYP-D* approach, which includes the modified DFT-D2 correction, provided crystallographic results that are between those calculated using the LDA (lower bound in Fig. S1) and the GGA functionals (upper bound). Interestingly, the LDA approach led to a correct description of the c lattice parameter with pressure, whereas the generalised-gradient approximation functional resulted in the same overestimation noted with the B3LYP approach. Thus, GGA and hybrid B3LYP functionals provide an adequate description, i.e., within the experimental information, of the a and b lattice parameters, where strong ionic/covalent bonds are responsible of the TOT structure; though, both approaches fail in properly describing the c-axis behaviour. This is a further confirmation of the hypothesis that the forces acting along the [001] direction of micas are not only electrostatic/ionic, but there is also a contribution from weak van der Waals interactions (Ulian and Valdrè, 2015c, 2023a), as demonstrated using the DFT-D2 correction. The authors are aware that the B3LYP-D* approach, whose DFT-D2 parameters were empirically recalibrated on molecular crystals (Civalleri et al., 2008), is just one of the possible methods to properly treat long-range interactions, and other more ab initio methods could be employed, e.g., the vdW functionals developed of Dion et al. (2004). However, the proposed theoretical framework was already suitable to provide new insights into this fundamental topic regarding the bonding nature of mica phyllosilicates.

The axial compressibility, M₀, was obtained using a linearized form of the third-order Birch-Murnaghan equation of state, as described by Angel et al. (2014), which was also employed to re-calculate the experimental compressibility reported by Pavese et al. (2003) and Comodi et al. (2004) including the uncertainties on pressure and lattice parameters. The results are shown in Table 3, alongside the theoretical results of Chheda and collaborators (2014). The results are consistent for all the simulations, as the bulk modulus $K_{\mathrm{0}}=\left[M_{\mathrm{0}}^{-\mathrm{1}}\left(a\right)+M_{\mathrm{0}}^{-\mathrm{1}}\left(c\right)+M_{\mathrm{0}}^{-\mathrm{1}}\left(c\right)\right]$ calculated from the axial compressibility is within 2.64 % (1M) and 1.67 % (2M₁) from the K₀ value obtained from the equation of state fit. In general, it can be noted that the K₀ values are almost the same between the two polytypes. However, the 1M and 2M₁ phlogopite models showed a slightly different behaviour in terms of axial compression. In fact, while the M₀(a) and M₀(b) values are the same between the two phlogopite models, the M₀(c) modulus is significantly higher in the 2M₁ polytype, with a stiffness increase of about 10 % with respect to the 1M one. Similar results were obtained by Chheda et al. (2014) from DFT simulations and by Pavese et al. (2003) from XRD experiments. Also, the recalculated values of Comodi et al. (2004) with the linearized formulation lead to a bulk modulus within 0.57 % the value obtained from the PV EoS fitting.

Table 3 Axial compressibility parameters obtained from linearized equation of state fit.

^a Present work. ^b Pavese et al. (2003). ^c Comodi et al. (2004). ^d Chheda et al. (2014). GTO – Gaussian-type orbitals; PAW – projector-augmented wave; SCXRD – single-crystal X-ray diffraction; SXRPD – synchrotron X-ray powder diffraction.

Crystal structure data of both 1M and 2M₁ polytypes are reported in Tables S2 and S3, respectively, providing further insights into the pressure effects on the phlogopite internal geometry.

Figure 3 (a) Variation as a function of pressure of the octahedral flattening angle ψ (red symbols) and tetrahedral rotation angle α_r (blue symbols) in the 1M and 2M₁ phlogopite models. (b) Evolution of the tetrahedral TO₄ and octahedral MO₆ polyhedral volume with pressure.

The TO₄ tetrahedra showed a stiff behaviour, with zero-pressure bulk moduli equal to 282(3) and 285(2) GPa for the 1M and 2M₁ polytypes, respectively, whereas the octahedral MO₆ units were more compressible, presenting a bulk modulus of 143(1) GPa for both phlogopite polytypes (Fig. 3a). This means that the TOT layer of the two minerals presented similar elastic properties. As also explained in previous works on phyllosilicates (Comodi et al., 2004; Chheda et al., 2014; Ulian and Valdrè, 2015c), the variation of the tetrahedral rotation angle α_r (Donnay et al., 1964b, a) is the main mechanism that accommodated the compression of the mineral. The α_r value is greater than 0 for both phlogopite polytypes at equilibrium, which means that the TO₄ mesh is not hexagonal (ideal) but forms a di-trigonal ring (see Fig. 1), and α_r increases with pressure as reported in Fig. 3b. This trend is in line with the simulations performed by Chheda and collaborators (2014), and with the general compressional behaviour of trioctahedral micas reported by several authors (Comodi and Zanazzi, 1995; Gatta et al., 2015, 2011; Pavese et al., 1999, 2003; Pawley et al., 2002). In addition, it is interesting to note that the octahedral flattening angle ψ became smaller by increasing pressure in the 1M polytype, whereas a direct proportionality was observed for the 2M₁ one.

3.3 Elastic constants and derived quantities

To provide a complete analysis of the elastic behaviour of phlogopite polytypes, the second-order elastic tensor of the two models was also calculated. The approach relies on stress-strain relationships based on total energy E(V, ε) calculations through a Taylor expansion in terms of the strain components truncated at the second order:

$$\begin{array}{}\text{(4)}& E\left(V,\mathit{\epsilon }\right)=E\left(V_{\mathrm{0}}\right)+V\sum _{\mathit{\alpha }}{\mathit{\sigma }}_{\mathit{\alpha }}{\mathit{\epsilon }}_{\mathit{\alpha }}+{\displaystyle \frac{V}{\mathrm{2}}}\sum _{\mathit{\alpha }\mathit{\beta }}{C}_{\mathit{\alpha }\mathit{\beta }}{\mathit{\epsilon }}_{\mathit{\alpha }}{\mathit{\epsilon }}_{\mathit{\beta }}+\mathrm{\dots }\end{array}$$

where ε is the strain, σ is the stress, C is the 6×6 matrix representation of the elastic tensor (Voigt's notation), V is the volume of the strained unit cell and V₀ is the volume at equilibrium. In this expression, α, β=1, 2, 3, … 6. The adiabatic second-order elastic moduli are related to the second derivatives of the total energy E on the strain according to:

$$\begin{array}{}\text{(5)}& C_{\mathit{\alpha }\mathit{\beta }}={\displaystyle \frac{\mathrm{1}}{V}}{\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}E}{\partial {\mathit{\epsilon }}_{\mathit{\alpha }}\partial {\mathit{\epsilon }}_{\mathit{\beta }}}}\right|}_{\mathrm{0}}\end{array}$$

The C_αβ values were calculated by imposing a certain amount of strain ε along the crystallographic direction corresponding to the component of the dynamical matrix. This procedure is automated in the CRYSTAL code (Perger et al., 2009), and in the present simulations the default values that control the calculation of the elastic moduli were employed. To calculate the elastic tensor, the crystal was oriented with the b and c crystallographic axes parallel to the y and z Cartesian axes, respectively. Within this crystal orientation and considering the monoclinic lattice, the 6×6 matrix C of the elastic moduli in Voigt's notation is given by:

$$\begin{array}{}\text{(6)}& \mathbf{C}=\left(\begin{array}{llllll}{C}_{\mathrm{11}}& C_{\mathrm{12}}& C_{\mathrm{13}}& \mathrm{0}& C_{\mathrm{15}}& \mathrm{0}\\ & C_{\mathrm{22}}& C_{\mathrm{23}}& \mathrm{0}& C_{\mathrm{25}}& \mathrm{0}\\ & & C_{\mathrm{33}}& \mathrm{0}& C_{\mathrm{35}}& \mathrm{0}\\ & & & C_{\mathrm{44}}& \mathrm{0}& C_{\mathrm{46}}\\ & & & & C_{\mathrm{55}}& \mathrm{0}\\ & & & & & C_{\mathrm{66}}\end{array}\right),\end{array}$$

whose values are reported in Table 4, considering both B3LYP and B3LYP-D* approaches. As expected, the inclusion of long-range interactions leads to a stiffer behaviour of phlogopite, especially for the diagonal terms C_αα.

Table 4 Elastic moduli C_αβ (in GPa) of phlogopite polytypes obtained in the present DFT simulations, compared to previous theoretical and experimental data.

Notes: GTO – Gaussian-type orbitals; PAW – projector-augmented wave; EWP – Elastic wave propagation. ^a Present work. ^b Chheda et al. (2014). ^c Aleksandrov et al. (1974).

Very few experimental data were reported in the scientific literature because suitable single crystals of phlogopite are rare, and many of these specimens are required for a complete characterisation of the elastic tensor of low-symmetry phases. Nevertheless, the present simulations are in good agreement with the experimental results of Alexandrov and Ryzhova (1961) and Aleksandrov et al. (1974), where the authors measured the elastic moduli from two crystals, one from the Slyudyanka and one from the Aldan region, near Lake Baikal (Russia). It is worth noting that only 5 independent elastic tensor components related to the hexagonal symmetry were reported instead of the 13 moduli required for the actual monoclinic symmetry of phlogopite. Also, two pairs of elastic moduli were equal, i.e., C₁₁= C₂₂ and C₄₄=C₅₅, and the off-diagonal components C₁₅, C₂₅, C₃₅ and C₄₆ were zero due to symmetry constraints. The present theoretical results at the B3LYP-D* level of theory indeed showed that C₁₁≈C₂₂ and C₄₄≈C₅₅, suggesting the near-isotropic elastic behaviour of the basal plane. Vaughan and Guggenheim (1986) measured the the full monoclinic elastic tensor for the dioctahedral equivalent of phlogopite, i.e., muscovite KAl₂Si₃AlO₁₀(OH)₂, finding a similar elastic behaviour with a small, but not insignificant, deviation from the hexagonal symmetry. Within this computational framework, the calculated C_αα moduli of phlogopite are overestimated by about 31 GPa (α=1, 2, 3) and 13 GPa (α=4, 5, 6) on average, but the trends of the diagonal components $C_{\mathrm{11}}\approx C_{\mathrm{22}}>C_{\mathrm{33}}$ and $C_{\mathrm{44}}\approx C_{\mathrm{55}}<C_{\mathrm{66}}$ are the same as those experimentally obtained from elastic wave propagation experiments (Aleksandrov et al., 1974; Alexandrov and Ryzhova, 1961). The uncorrected B3LYP functional resulted smaller absolute deviations from the experimental data (ca. 4 and 2 GPa for the uniaxial and biaxial tensor components, respectively) and the same $C_{\mathrm{11}}\approx C_{\mathrm{22}}>C_{\mathrm{33}}$ relationship; however, the C₄₄ modulus was smaller than the C₅₅ one, and also much smaller than the one calculated at the B3LYP-D* level. Regarding the off-diagonal components, both approaches resulted in C₁₂>C₁₃, although the uncorrected functional showed C₁₃>C₂₃, while the inclusion of the van der Waals interactions provided C₁₃≈C₂₃, which is in better agreement with the experimental evidence.

The different performance of the B3LYP and B3LYP-D* approaches in calculating the diagonal and off-diagonal components of the stiffness tensor depends on the underlying physics of the elastic response. The diagonal elastic constants are dominated by the response of the strongly bonded T–O–T layers. When dispersion is neglected (B3LYP), long-range interactions are underestimated, leading to an overexpanded interlayer region. In Gaussian basis-set calculations, this partly compensates the overestimation of forces caused by Pulay stress (as discussed below), producing an accidental error cancellation and a seemingly better agreement for the diagonal terms. In contrast, the off-diagonal terms describe shear-normal coupling and interlayer-intralayer stress transfer, which are very sensitive to the correct description of weak interactions. Without dispersion (B3LYP), the layers are too decoupled, and the coupling terms C_αβ (α≠β) are underestimated or show incorrect relative magnitudes. Including dispersion (B3LYP-D*) restores realistic interlayer distances and mechanical coupling, improving the prediction of the off-diagonal elastic constants. This further confirms that the proper treatment of long-range interactions is relevant for the fundamental analysis of the elastic properties of layered silicates, especially when using Gaussian-type orbitals basis sets.

The present results at the B3LYP(-D*) level of theory are also in line with the simulations performed by Chheda and collaborators (2014) using LDA and GGA functionals and plane-waves basis sets, albeit some notable differences are present (see Table 4). Considering the cited work, it is known that the local density approximation generally underestimates lattice vectors and overestimates the forces, resulting in a stiffer behaviour of simulated phases. Conversely, GGA functionals are generally underbinding, hence the calculated elastic properties can be slightly underestimated. Then, hybrid functionals like B3LYP, which include a small fraction of the exact Hartree-Fock exchange to the total energy of the system, are expected to provide results that are between the LDA/GGA extremes if, and only if, other computational parameters (basis set type and quality, k-point sampling, converge criteria, etc.) are the same. Another source of discrepancy between the present elastic results and those of Chheda et al. (2014) is indeed the basis set type. Localised Gaussian-type orbitals generally provide slightly larger C_αβ values than those obtained with plane-waves basis sets in density functional theory simulation because of the Pulay forces that arise from the Hellmann-Feynman theorem and the use of atom-centred GTOs. As previously explained (Ulian et al., 2021; Ulian and Valdrè, 2018), the Pulay forces are associated with the basis set superposition error (BSSE), which is due to the incompleteness of GTOs and that artificially increases the calculated stress/force within a solid phase. The authors are aware that an approach like the geometrical counterpoise (gCP) method (Brandenburg et al., 2013; Kruse and Grimme, 2012) devised for molecule and molecular crystals could at least partially remove this BSSE-related error. However, the current implementation is limited to some general-purpose GTO basis sets that do not include the one employed in the present work, and a thorough analysis of its suitability for inorganic crystals and minerals is still ongoing (Ulian and Valdrè, 2024b). Besides the Pulay stress, the general overestimation of the theoretical results, both the present one and those of Chheda et al. (2014), could be also due to the athermal conditions, i.e., the DFT simulations were carried out at 0 K without zero-point energy and any thermal contribution, while the elastic wave propagation measurements are typically carried out at room temperature (about 300 K). Regarding the present work, the use of gCP and including temperature effect, e.g., with the quasi-harmonic approximation, can provide results that are in more agreement with the experimental findings.

For practical applications, it is commonly preferred the use of the polycrystalline averages, namely the bulk (K) and shear (μ) moduli, calculated with the Voigt–Reuss–Hill approach (Hill, 1952; Nye, 1957) according to the formulas:

$$\begin{array}{}\text{(7)}& {\displaystyle }{K}_{\mathrm{V}}=\left(\mathrm{1}/\mathrm{9}\right)\left[C_{\mathrm{11}}+C_{\mathrm{22}}+C_{\mathrm{33}}+\mathrm{2}\left(C_{\mathrm{12}}+C_{\mathrm{13}}+C_{\mathrm{23}}\right)\right]\text{(8)}& {\displaystyle }{K}_{\mathrm{R}}={\left[S_{\mathrm{11}}+S_{\mathrm{22}}+S_{\mathrm{33}}+\mathrm{2}\left(S_{\mathrm{12}}+S_{\mathrm{13}}+S_{\mathrm{23}}\right)\right]}^{-\mathrm{1}}\text{(9)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{V}}={\displaystyle \frac{\left[C_{\mathrm{11}}+C_{\mathrm{22}}+C_{\mathrm{33}}+\mathrm{3}\left(C_{\mathrm{44}}+C_{\mathrm{55}}+C_{\mathrm{66}}\right)-\left(C_{\mathrm{12}}+C_{\mathrm{13}}+C_{\mathrm{23}}\right)\right]}{\mathrm{15}}}\text{(10)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{R}}={\displaystyle \frac{\mathrm{15}}{\mathrm{4}\left[S_{\mathrm{11}}+S_{\mathrm{22}}+S_{\mathrm{33}}-\left(S_{\mathrm{12}}+S_{\mathrm{13}}+S_{\mathrm{23}}\right)\right]+\mathrm{3}\left(S_{\mathrm{44}}+S_{\mathrm{55}}+S_{\mathrm{66}}\right)}}\end{array}$$

where $S_{ij}=C_{ij}^{-\mathrm{1}}$ are the component of the compliance tensor, i.e., the inverse of elastic moduli tensor C, and the subscripts V and R indicate the Voigt (upper) and Reuss (lower) bound of the isotropic elastic properties. The Young's modulus E and the Poisson's ratio υ were calculated from the following relations:

$$\begin{array}{}\text{(11)}& E={\displaystyle \frac{\mathrm{9}K\mathit{\mu }}{\mathrm{3}K+\mathit{\mu }}}\end{array}$$

and

$$\begin{array}{}\text{(12)}& \mathit{\upsilon }={\displaystyle \frac{\mathrm{3}K-\mathrm{2}\mathit{\mu }}{\mathrm{2}\left(\mathrm{3}K+\mathit{\mu }\right)}}\end{array}$$

The polycrystalline mechanical properties were reported in Table 5. The K_R values obtained for both phlogopite 1M and 2M₁ polytypes were 58.8 and 59.4 GPa, respectively, which are consistent with those calculated through the BM3 fitting procedure (57.9 and 58.3 GPa, respectively).

The elastic anisotropy of phlogopite was described using the universal anisotropic index, A^U, as proposed by Ranganathan and Ostoja-Starzewski (2008):

$$\begin{array}{}\text{(13)}& A^{\mathrm{U}}=\mathrm{5}{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{V}}}{{\mathit{\mu }}_{\mathrm{R}}}}+{\displaystyle \frac{K_{\mathrm{V}}}{K_{\mathrm{R}}}}-\mathrm{6}\ge \mathrm{0}.\end{array}$$

For all but isotropic systems, for which A^U=0, the ratio between the Voigt and Reuss moduli is not equal to unity, thus A^U>0. The percentage anisotropy in compression (A_K) and shear (A_μ) was also calculated according to the following formulas (Ranganathan and Ostoja-Starzewski, 2008):

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& A_K=\mathrm{100}{\displaystyle \frac{K_{\mathrm{V}}-K_{\mathrm{R}}}{K_{\mathrm{V}}+K_{\mathrm{R}}}}\\ & A_{\mathit{\mu }}=\mathrm{100}{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{V}}-{\mathit{\mu }}_{\mathrm{R}}}{{\mathit{\mu }}_{\mathrm{V}}+{\mathit{\mu }}_{\mathrm{R}}}}\end{array}\end{array}$$

The anisotropies obtained via these expressions were reported in Table 5, showing that the A^U values of the two phlogopite polytypes are very similar (about 3.8 at the B3LYP-D* case), and between the values calculated from the theoretical results of Chheda and collaborators (2014) using the LDA (2.9) and GGA (4.4) functionals. There is a slight discrepancy between the different DFT approaches and the experimental results (A^U=11.8), which is related to the high anisotropy of the shear moduli. It is suggested that the high ${\mathit{\mu }}_{\mathrm{V}}/{\mathit{\mu }}_{\mathrm{R}}$ ratio is probably due to the very small values of the C₄₄=C₅₅ elastic tensor components that were experimentally measured by Aleksandrov et al. (1974) with elastic wave propagation methods. Similar figures were found for the bulk and shear anisotropies, with A_K≈14 % and A_μ≈26 % calculated at the B3LYP-D* level of theory, intermediate values between those obtained using LDA and GGA functionals by Chheda and co-workers (2014). Conversely, the uncorrected B3LYP functional led to an anomalously high anisotropy index (35.2) due to the very large difference between the Voigt and Reuss shear moduli (A_μ≈77 %).

Table 5 Isotropic elastic properties (bulk K, shear μ, and Young's E moduli, and Poisson's ratio υ) of phlogopite, alongside the universal anisotropy index A^U and the percentage of anisotropy of the bulk (A_K) and shear (A_μ) moduli.

^a Present work. ^b Chheda et al. (2014). ^c Aleksandrov et al. (1974). Subscripts V and R indicate the Voigt (upper) and Reuss (lower) bounds, respectively. ^d Values calculated in the present work from the second-order elastic moduli reported by Aleksandrov et al. (1974).

To further assess the validity of the linear model used to calculate the axial compressibilities, the axial moduli M(a), M(b) and M(c) were calculated from the elastic compliance tensor components according to the expressions reported by Mookherjee et al. (2016):

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& M\left(a\right)={\mathit{\beta }}_a^{-\mathrm{1}}={\left(S_{\mathrm{11}}+S_{\mathrm{12}}+S_{\mathrm{13}}\right)}^{-\mathrm{1}}\\ & M\left(b\right)={\mathit{\beta }}_c^{-\mathrm{1}}={\left(S_{\mathrm{12}}+S_{\mathrm{22}}+S_{\mathrm{23}}\right)}^{-\mathrm{1}}\\ & M\left(c\right)={\mathit{\beta }}_c^{-\mathrm{1}}={\left(S_{\mathrm{13}}+S_{\mathrm{23}}+S_{\mathrm{33}}\right)}^{-\mathrm{1}}\end{array}\end{array}$$

For instance, the calculated moduli were M(a)=344.3, M(b)=365.2 and M(c)=88.0 GPa for phlogopite-1M, and M(a)=360.5, M(b)=353.9 and M(c)=89.1 GPa for the 2M₁ polytype, values consistent with those reported in Table 3 obtained from the equation of state fit and that provide another measure of the anisotropy of the mineral, with a ratio $M\left(a\right):M\left(b\right):M\left(c\right)\approx \mathrm{4}:\mathrm{4}:\mathrm{1}$ for both polytypes.

Directional single-crystal elastic properties (Young's modulus E, linear compressibility β= K⁻¹, shear modulus μ and Poisson's ratio υ) were calculated in the Cartesian space using the QUANTAS code (Gaillac et al., 2016; Ulian and Valdrè, 2022). The interested reader can find in the cited works all the formulations employed to calculate the spatial dependence of the cited properties. Differently from the other elastic properties reported above, where no significant variation between the 1M and 2M₁ phlogopite polytypes was evinced, the directional properties showed some effect due to the different TOT stacking (see Fig. S1). For example, the shape of the Young's modulus on the (xy) plane (Fig. S1a) is slightly more compressed along NE-SW and NO-SE directions in the polytype-1M, while it appeared more isotropic in phlogopite-2M₁. The shear moduli on the (xz) plane of the latter polytype (see Fig. S1c) is slightly canted, with the oblate maxima direction almost parallel to the east-west direction.

Figure 4 Lamber equal area (stereographic) projection of the upper hemisphere of the seismic anisotropy calculated for (a) 1M and (b) 2M₁ polytypes of phlogopite. In each panel, the left picture is related to the primary velocity, the central one is the S-wave anisotropy (related to the fast and slow shear waves) and the rightmost stereoplot is the polarization of the v_S1 (fast secondary) shear velocity. All the quantities were calculated at the DFT/B3LYP-D* level of theory.

Finally, the phase (seismic) velocities were calculated by solving Christoffel's equation (Musgrave, 1970) using the routines implemented in the QUANTAS code (Jaeken and Cottenier, 2016; Ulian and Valdrè, 2022). The results are shown in Fig. 4, where the compressional wave velocity (v_P, P-wave, primary), the anisotropy of the shear velocities (S-wave, secondary), and the polarization of the fast shear wave (v_S1) are reported for both phlogopite polytypes. From the stereographic projections, particularly the S-wave anisotropy, it is immediately visible an east-west mirror plane normal to the b-axis that is related to the monoclinic crystal system of the mineral. This figure is in very good agreement with the theoretical analysis of Chheda et al. (2014), with some degree of deviation of the absolute values of the seismic velocities and anisotropy due to the stiffer elastic moduli obtained with the combination of the Gaussian-type orbitals basis sets and hybrid B3LYP-D* functional.

Table 6 reports the maximum and minimum v_P values, the percentage of v_P anisotropy (Av_P) and the range of S-wave anisotropy (Av_S), as calculated from the present GTO/B3LYP(-D*) simulations for both phlogopite polytypes, together with previous experimental (Aleksandrov et al., 1974) and theoretical (Chheda et al., 2014) data related to the phlogopite 1M. The results are in general good agreement, with a slight overestimation of the maximum P-wave velocity and a smaller minimum primary wave velocity obtained from the uncorrected B3LYP functional. The shear wave anisotropy at the B3LYP-D* level of theory is within the values obtained by Chheda and collaborators (2014), with the maximum Av_S value being underestimated. This latter datum is instead more than double when calculated from the B3LYP results, and very overestimated compared to the experimental value reported by Aleksandrov et al. (1974).

As mentioned in the introduction, the knowledge of the elastic properties could help explain and interpret some geophysical observations reported in previous analyses. For instance, a large delay time (δt) of about 1–2 s was measured in certain subduction zones, like the Ryukyu and Calabrian arcs, in addition to the trench-parallel shear wave splitting (Long and Silver, 2008). Chheda and collaborators (2014) suggested that peridotitic rocks with olivine as the major component would require a very large thickness ranging from 100 to 150 km to justify this large δt value. However, phyllosilicates like talc (Mainprice et al., 2008; Ulian et al., 2014), muscovite (Ulian and Valdrè, 2015a), antigorite (Mookherjee and Capitani, 2011) and phlogopite (Pavese et al., 2003; Chheda et al., 2014) present significantly higher anisotropic behaviours than olivine. Thus, their presence in the subduction layer would require smaller thickness (∼10–15 km) to explain both the trench parallel shear wave polarization anisotropy and the large delay time (Chheda et al., 2014). Indeed, according to previous investigations (Trønnes, 2002), phlogopite has a wide PT stability, thus it can be found in different settings, e.g., in the upper mantle and subduction zone, the latter being also characterised by hydrating conditions that could further stabilise this trioctahedral mica and other hydrous phases.

Table 6 Maximum and minimum compressional wave velocity (v_P, km s⁻¹), P-wave anisotropy (Av_P, %), and maximum and minimum S-wave anisotropy (Av_S, %) of phlogopite.

^a Present work. ^b Chheda et al. (2014). ^c Aleksandrov et al. (1974). ^d Values calculated in the present work from the second-order elastic moduli reported by Aleksandrov et al. (1974).

4 Conclusions

In the present work, using ab initio Density Functional Theory simulations, the elastic properties of two phlogopite polytypes, namely the 1M and the 2M₁, were investigated considering hydrostatic compression up to about 15 GPa and the monoclinic elastic tensor. The following conclusions can be drawn:

1. despite the different stacking of the TOT layers, the elastic properties of phlogopite were not severely affected by polytypism. Considering the numerical errors of the fitting procedure, the bulk modulus K₀ obtained from the third-order Birch-Murnaghan equation of state was in the range 55–60 GPa for both polytypes.

2. the compressional behaviour of phlogopite is highly anisotropic, with the a and b crystallographic axes being about four times less compressible than the c one, as noticed from the linear moduli. This agrees with the experimental observations made on dioctahedral and trioctahedral micas, and it is due to several factors, i.e., the high compressibility of the interlayer region, the low compressibility of the tetrahedral and octahedral polyhedral, and the accommodation of the stress within the TOT layers via variation of the tetrahedral rotation angle α_r.

3. the elastic tensors of phlogopite-1M and phlogopite-2M₁ have comparable C_ij components (in Voigt's notation), with a difference in sign and absolute value in the C₂₅ one. This is probably the reason of the small differences observed in the directionally dependent elastic properties.

4. the DFT simulations were able to correctly describe the behaviour of the known 1M polytype of phlogopite, extending the knowledge to the 2M₁ one. In particular, the use of the DFT-D2 correction provided a robust description of the physical forces acting along [001], namely the stacking direction of the TOT +i layers.

It is also known that the elastic properties of micas can be deeply affected by the presence of cationic/anionic substitutions in the T and O sheets, e.g., Fe^II substituting Mg^II. These are currently under investigation and will be the subject of a future work.

From the mineralogical, petrological, and geophysical perspectives, the knowledge of the physical properties of phlogopite is of utmost importance for deepening the understanding of the thermodynamic and kinetic stability of altered rock assemblages, such as the mica–amphibole–rutile–ilmenite–diopside system. Indeed, the phlogopite chemical composition is similar to that of the ternary system olivine–quartz–alkali, which is also close to the bulk composition of the cited altered rocks. Also, the elasticity of the trioctahedral mica is relevant both for bulk and 2D applications of phlogopite in materials science, owing to the dielectric properties of the mineral that can be tuned by the exfoliation in the basal plane and cationic substitutions, particularly in the octahedral sheet. First principles simulations can help in shedding light on these topics, providing useful resources for the interpretation of geophysical phenomena and the development of more, sustainable materials for advanced applications.