Hydrothermal Fluids and Where to Find Them: Using Seismic Attenuation and Anisotropy to Map Fluids Beneath Uturuncu Volcano, Bolivia

Mapping fluid accumulation in the crust is pertinent for numerous applications including volcanic hazard assessment, geothermal energy generation, and mineral exploration. Here, we use seismic attenuation tomography to map the distribution of fluids in the crust below Uturuncu volcano, Bolivia. Seismic P wave and S wave attenuation, as well as their ratio (QP/QS), constrain where the crust is partially and fully fluid‐saturated. Seismic anisotropy observations further constrain the mechanism by which the fluids accumulate, predominantly along aligned faults and fractures in this case. Furthermore, subsurface pressure‐temperature profiles and conductivity data allow us to identify the most likely fluid composition. We identify shallow regions of both dry and H2O/brine‐saturated crust, as well as a deeper supercritical H2O/brine column directly beneath Uturuncu. Our observations provide a greater understanding of Uturuncu's transcrustal hydrothermal system, and act as an example of how such methods could be applied to map crustal fluid pathways and hydrothermal/geothermal systems elsewhere.

Here, we use microseismicity at Uturuncu  to perform seismic attenuation tomography and anisotropy analysis. Seismic attenuation is sensitive to the presence of fluids, with tomography enabling any fluids to be mapped (e.g., Hauksson & Shearer, 2006). Attenuation studies of volcanoes have previously been undertaken (Bohm et al., 2013;Caudron et al., 2019;De Siena et al., 2014;Gudmundsson et al., 2004;Lanza et al., 2020;Lees, 2007;O'Brien & Bean, 2009;Sanders et al., 1995;Zhao, 2001). Seismic anisotropy in volcanic or tectonic regions is predominantly sensitive to fault structures (e.g., Baird et al., 2015), enabling one to identify whether fluids migrate and accumulate along volcanic and/or fault structures (Bacon et al., 2021;Baird et al., 2015;Gerst & Savage, 2004;Johnson et al., 2011;Maher & Kendall, 2018;Nowacki et al., 2018). The novelty of this work lies in using seismic attenuation and anisotropy to identify the dominant attenuation mechanism and fluid saturation level. We then combine these results with auxiliary pressure, temperature, and electrical conductivity profiles (Comeau et al., 2016;Pritchard et al., 2018), to infer fluid composition and ascertain how fluids migrate and accumulate in the crust.

Earthquake Catalog
The Uturuncu earthquake catalog used in this study is from Hudson et al. (2022). Seismicity is detected using the PLUTONS network (Kukarina et al., 2017), which was deployed from the 13 April 2010 to the 27 October 2012. The PLUTONS network comprised of 33 Guralp CMG-3T 120 s seismometers recording at 100 Hz. Earthquakes are detected using QuakeMigrate (Hudson et al., 2019;Smith et al., 2020) and relocated using NonLinLoc (Lomax & Virieux, 2000). The velocity model is the same as in Hudson et al. (2022). This produces a catalog of 1,356 earthquakes, which we use in this study (see Figure S1 in Supporting Information S1). Catalog statistics and filtering parameters are shown in Figure S6 in Supporting Information S1.

Seismic Attenuation
Seismic attenuation is the loss of energy of a seismic wave into the medium during propagation. For this study, we assume that the amplitude of a seismic wave, ( , ) , at a position with a frequency , is approximately described by Aki and Richards (2002): where 0( ) is the source amplitude in the direction from vertical and from North, ( ) is a geometrical spreading factor, ( ) is the raypath length, and ( ) is the attenuation factor for a particular frequency. Therefore, the greater the value of , the greater the attenuation.
In seismology, the convention is to quantify attenuation using the seismic quality factor, Q, the inverse of , which is given by Aki and Richards (2002): where is the velocity of phase (P or S). Here, we approximate , and therefore Q, to be frequency-independent over the bandwidth of the earthquakes studied (0.5-20 Hz ), with the majority of earthquake energy lying between 2 and 10 Hz (see Figure 2, Hudson et al., 2022). The attenuation is measured by fitting a Brune-model (Brune, 1970) to the spectra of the data. We use a spectral-ratio method to isolate source (corner-frequency) and path (Q) effects to fit the Brune-model, removing the trade-off between corner-frequency and Q (see in Supporting Information S1; Lapins, 2021). The windows used for both P and S waves are 0.1 s before and 0.6 s after the relevant phase arrivals, approximately isolating direct-wave arrivals only. A sensitivity analysis of window-length is included in Supporting Information S1.
Assuming Equation 1 is valid for the study site implies that intrinsic attenuation and simple scattering off impedance contrasts dominate, while low-frequency resonant scattering and surface-wave amplification are negligible. Although this assumption might not be valid in highly topographically varying volcanic settings (Cormier, 2021), we attempt to minimize these low-frequency effects (Gabrielli et al., 2020) by applying a high-pass filter at 0.5 Hz. Such surface-wave phases will also typically arrive after our observation window (0.6 s). However, if this is not the case, then the above effects would result in an underestimation of S wave attenuation. As such, any high-attenuation regions that we do observe would only be enhanced if we accounted for the above mechanisms, providing us with confidence in our current interpretations of high-attenuation regions.

Attenuation Tomography
Attenuation tomography uses path-averaged attenuation observations, Q path , between earthquake sources and receivers, to map the attenuation structure. The 3D attenuation tomography performed in this study uses the earthquakes and receivers shown in Figure S1 in Supporting Information S1, with the raypaths shown in Figures  S2 and S3 in Supporting Information S1. Our method is based upon that described in Wei and Wiens (2018), with some alterations due to the local nature of our study area. A full description of the attenuation tomography methodology can be found in Text S1 in Supporting Information S1 (with t* window-length sensitivity analysis in Text S3 in Supporting Information S1).

Q P /Q S Ratios
P wave attenuation is sensitive to scattering from faults and other velocity contrasts, as well as the intrinsic attenuation of the host rock and other compressible media, such as supercritical fluids and gases. S wave attenuation is also sensitive to the scattering from fault structures and the presence of fluids, but in contrast to P waves is insensitive to the compressibility of fluids (M. Chapman, 2003;S. Chapman et al., 2021;Klimentos, 1995). Attenuation tomography therefore not only allows for mapping the presence of any fluids, but Q P /Q S ratios are also diagnostic of fluid saturation, with Q P /Q S < 1 indicating that rock is only partially saturated and Q P /Q S > 1 suggesting fully saturated rock (Amalokwu et al., 2014;Hauksson & Shearer, 2006;Winkler & Nur, 1979). Here, we assume the concept of partially saturated rock in its broadest sense, in that it simply has to have some fraction of compressible fluid or gas present (see Figure 4).

Seismic Anisotropy
Seismic anisotropy in this study refers to S wave velocity anisotropy, measured from shear-wave splitting. Such anisotropy can be broadly attributed to two factors: crystallographic preferred orientation, where individual crystals of the medium preferentially align; and shape-preferred orientation anisotropy, typically caused by the preferential alignment of fractures (Kendall, 2000). As an S wave propagates through an anisotropic region, the energy will be partitioned into two orthogonal components, one oriented in the plane of the fast-direction, , and the other in the plane of the slow-direction. This phenomenon is called shear-wave splitting (Crampin, 1981;Silver & Chan, 1991), measured here using SWSPy (Hudson, 2022). is controlled by the orientation of the fabric and/ or fractures. The arrival-time difference between the fast and slow S waves, , is controlled by the strength and/ or spatial extent of the anisotropy.

Results
The overall attenuation tomography results for Q P , Q S, and Q P /Q S are shown in Figure 1. Tomographic resolution test results, based on the point-spread-function method (Rawlinson & Spakman, 2016), are shown in Figures S2 and S3 in Supporting Information S1. Masks are applied to the results for regions that cannot be interpreted due to insufficient resolution (see Discussion for details). Regional Q P values are ∼200-300, while regional Q S values range from 200 to 400. These values fall within the range observed at similar crustal depths elsewhere (Del Pezzo et al., 1995;Hauksson & Shearer, 2006;Lanza et al., 2020). Regional Q P and Q S both generally increase with depth, as expected for denser crust. A region of lower Q P (high attenuation) extends in the SE-NW direction, intersecting Uturuncu. This high attenuation region extends to depths of 8 km bsl. Additionally, we tentatively suggest that there may be a narrow region of high Q P (low attenuation) extending SW-NE ∼20 km North of Uturuncu. Q S follows the SE-NW trend observed in Q P , except with more isolated pockets of lower Q S (high attenuation) rather than a continuous band. These isolated pockets occur on lateral length scales of ∼5-10 km and extend to depths of up to 3 km bsl. One exception is a highly attenuating region ∼15 km North of Uturuncu, at a depth of ∼12 km bsl. However, it should be noted that this feature is at the limit of the resolvable region (see Figures S2 and S3 in Supporting Information S1).
The majority of the region has Q P /Q S < 1 (see Figures 1g-1i). A significant feature is a column of Q P /Q S << 1 directly beneath Uturuncu's summit, extending from 1 to 2 km bsl to at least 8 km bsl. This feature is shown in greater detail in Figures 2a-2c (labeled A), where Q P /Q S is compared to resistivity tomography results from Comeau et al. (2016). These resistivity data are primarily sensitive to the conductivity of fluids in the pore space and/or fractures, and are therefore valuable for interrogating fluid composition. Feature A corresponds to a region of low resistivity (<10 Ω m ), indicative of high-salinity brines (Afanasyev et al., 2018). Isolated features with Q P /Q S > 1 also exist. One such feature, beneath Cerro Loromayo (Fracchia, 2009;Soler et al., 2007), is shown in detail in Figures 2d-2f (labeled C). This feature also corresponds to a region of low resistivity, as do other Q P / Q S > 1 regions. Although the strikes of the Q P /Q S results are less clear than for Q P or Q S , there are features that appear to have strikes SE-NW (e.g., through Uturuncu) and SW-NE (e.g., 20 km North of Uturuncu).
Seismic anisotropy observations at Uturuncu are shown in Figure 3. The rose diagrams in Figure 3a compare the orientations of the fast S wave polarisation (yellow) to the source polarizations derived during the S wave splitting analysis (blue), as well as fault strikes from Hudson et al. (2022) (gray) and previously published seismic anisotropy results (red) (Maher & Kendall, 2018). The fault strikes and source polarizations, which theoretically should be oriented fault-parallel, are primarily oriented SW-NE and NNW-SSE. However, the fast-directions are predominantly oriented NNW-SSE, with a minority of fast-directions oriented SW-NE lying to the north of Uturuncu (see Figure 3b). Data showing the orientation and strength of anisotropy with depth are plotted in Figure  S5 in Supporting Information S1, which suggest that the strength of anisotropy decreases as depth increases.

Faults and Fluids
Seismic attenuation can be caused by both absorption and scattering mechanisms. Here, we do not attempt to quantify the contribution of each mechanism. We instead identify key trends in the orientation of attenuation features and compare them to anisotropy and fault strike information, so as to infer what the overall attenuation shows.
S wave source polarizations obtained from the shear-wave splitting analysis (see Figure 3) show two dominant orientations, NW-SE and NE-SW. S wave source polarizations are interpreted to correspond to fault orientations, Masked areas correspond to regions where features of at least 4 km size cannot be resolved (see Figures S1 and S2 in Supporting Information S1 for resolution tests). Additionally, Q P /Q S results masked where Q P and/or Q S > 500 (see Discussion for justification on masks). The time-window for all Q measurements is 0.7 s.
since the polarisation will align with fault slip for a double-couple source. Fault strikes obtained from principal component analysis in Hudson et al. (2022) agree with S wave polarizations observed in this study, showing the same fault orientations. However, fast S wave velocity anisotropy directions primarily align with the NW-SE striking faults. At depths shallower than 8 km bsl, high S wave attenuation is confined to small pockets and high P and S wave attenuation are generally oriented in a narrow band with a strike NW-SE (see Figure 1). The attenuation structure deeper than 8 km bsl transitions to become approximately homogeneous, with anisotropy also decreasing with depth. This is likely an artifact of data resolution limits (see Figure S2 in Supporting Information S1).
We suggest that the observation of fault strikes in two directions yet attenuation, and to some extent S wave anisotropy, aligning with only the NW-SE direction is caused by the presence of fluids within the NW-SE faults but not the NE-SW faults. This is likely controlled by the crustal stress fabric, with NW-SE faults preferentially  Figure S1 in Supporting Information S1 for location on a map). Profiels CC' and DD' intersect the summit of Cerro Loromayo. Hashed regions are where resistivity >10 Ω m, for which any brines would be low salinity or not present (Afanasyev et al., 2018). Labels A-C are as referred to in Figure 4 and the Discussion. Results are masked as in Figures 1g-1i. opening and entraining fluids, while approximately perpendicular NE-SW faults are closed with little/no fluids present. The crustal stress state is presumably governed locally (<10 km from Uturuncu) by deformation at the volcano Pritchard et al., 2018) and regionally by tectonic stresses. Both fractures can scatter seismic energy due to reflection off fracture interface (Zhu et al., 2007). However, fluids in the pore space and along fractures can also cause attenuation due to mechanisms such as squirt flow (M. Chapman, 2003;M. Chapman et al., 2002). We infer that fluids play an important role in facilitating the observed attenuation structure. Furthermore, seismicity has b-values >1, consistent with fluids reducing fault normal stresses . Additionally, the location of a column of P wave attenuation is in agreement with a hypothesized route for fluid ascent (Del Potro et al., 2013;Gottsmann et al., 2017;Pritchard et al., 2018). Based on the above evidence, we therefore suggest that the observed attenuation structures are predominantly caused by faults with fluids ascending and/or accumulating along them.

Identifying and Mapping Fluids and Their Composition
Assuming that the observed attenuation is primarily caused by fluids within faults, then relevant questions are what the fluid composition is and does it vary spatially. For simplicity, we consider only three representative fluids: CO 2 ; H 2 O (or analogous brines); and molten rock. Diagnostic data for addressing these questions are Q P /Q S ratios, which indicate whether the crust contains a proportion of fluids or gases that are compressible (Q P /Q S < 1) or the crust is fully saturated with incompressible fluids (Q P /Q S > 1) (see Section 2.2.2 and Amalokwu et al., 2014;S. Chapman et al., 2021;Hauksson & Shearer, 2006;Tisato et al., 2015). Our interpretation of Q P /Q S ratios is based on the assumption that both Q P and Q S results are of adequate quality to obtain their ratio, within their resolution limits. Additionally, pressure and temperature profiles, as well as resistivity data, can be used in combination with Q P /Q S to constrain fluid composition spatially. Figure 4 shows a schematic diagram summarizing an inferred map of different fluids at Uturuncu based on the following discussion. Figure 4 also includes crustal pressure and temperature profiles (Pritchard et al., 2018 and references therein) and schematic compressibility and conductivity profiles that are used to further inform this discussion.
Below sea-level, features are less well resolved by the attenuation tomography than at shallower depths. However, one feature we are confident of is the Q P /Q S << 1 column directly beneath Uturuncu, labeled A in Figure 4 (also see Figures 2a-2c), which extends from ∼2 to 8 km bsl. This region of crust is highly faulted (Hudson et al., 2022), with Q P /Q S implying that these faults contain at least some compressible fluid component. Given the lithostatic pressures, it is unlikely that any gaseous phase is present at these depths. Intriguingly, H 2 O becomes supercritical at these depths (see Figure 4a), driving the compressibility up and therefore Q P /Q S down (see Figure 4b). This feature has resistivities 10 Ω m (Comeau et al., 2016), consistent with high salinity brines (Afanasyev  Figure 4c) and not molten rock (Laumonier et al., 2017). This is consistent with (Gottsmann et al., 2017), which suggests the shallowest possible depth of partial melt might be 6 km bsl. These conductivities also rule out this feature being solely comprised of supercritical CO 2 (see Figure 4c). Furthermore, the sudden transition from a weak Q P /Q S < 1 signal to a strong Q P /Q S << 1 signal at ∼2 km bsl, where H 2 O/brines become supercritical, implies that this region contains at least some supercritical, compressible brine phase. However, this region of crust could also contain CO 2 , or other volatiles, in combination with brines.
At depths shallower than sea-level, Q P /Q S ratios vary <1 (red regions, Figure 2, labeled B in Figure 4) and >1 (blue regions, Figure 2, labeled C in Figure 4), indicating that some regions of the shallow crust are partially saturated and others are fully saturated. This observation is expected <1 km from the surface, since this region of Bolivia is arid, yet surface snow and H 2 O are present in places, which likely percolate into the shallow crust. However, saturated regions of the crust extend to depths at sea-level (>4 km below surface), too deep for groundwater aquifers. We interpret the regions in the vicinity of Uturuncu's summit with Q P /Q S < 1 and corresponding low resistivities (see label B, Figures 4b and 4d) as a 2-phase brine-CO 2 mixture, which could be referred to as "sparkling CO 2 ". Similar observations of bubbles within fluids have been suggested to cause significant attenuation elsewhere (Tisato et al., 2015). These brines and CO 2 likely ascend from the APMB, via feature A, along faults in the crust.
Regions of Q P /Q S > 1, implying crust saturated with incompressible fluids, typically lie away from Uturuncu's summit. One particularly obvious example is that focused on in Figures 2d-2f, labeled C in Figure 4d. Again, the feature corresponds with low resistivities, indicative of high-salinity brines. This feature is therefore interpreted to be brine-saturated crust. Furthermore, this feature is particularly interesting as it lies beneath Cerro Loromayo (Fracchia, 2009;Soler et al., 2007; Cerro Loromayo summit at intersection of the profiles CC' and DD' in Figure 2), a volcano that could have exhibited activity as late as 0.9 Ma (Fracchia, 2009). Based on our Q P /Q S results, this feature appears isolated, with no vertical column-like feature below it. Furthermore, it corresponds to a low-density gravity anomaly (MacQueen et al., 2021), indicative of fluid accumulation. We therefore suggest that it is a brine-lens (Afanasyev et al., 2018) that was formed when Cerro Loromayo had an active hydrothermal system, but has subsequently cooled and is no longer fed by CO 2 or other volatiles from the APMB. Features such as this are exciting because they could potentially be a well-endowed accumulation of metal-rich brines (Blundy et al., 2021).
One further feature that we tentatively identify is a band of low Q P attenuation and a correspondingly Q P /Q S > 1, 20 km North of Uturuncu with a NE-SW strike. This feature is parallel to a subset of S wave source polarizations and fault strikes, and is therefore likely associated with faulting. We suggest that attenuation is lower here due to the prevailing stress-field acting to close similarly oriented faults, resulting in a smaller volume of fluids at higher pressures and therefore less likely to contain compressible gases, resulting in Q P /Q S > 1. These observations, if correctly interpreted, have implications for imaging crustal stress-state in similar fluid-rich systems.

Limitations
A fundamental limitation of any tomographic method is the spatial resolution. This spatial resolution is governed by the density of raypath coverage, which is shown in Figure S1 in Supporting Information S1. The P wave coverage is significantly better than the S wave coverage, owing to the greater number of picked P wave arrivals. The Q S and Q P /Q S tomography results are therefore limited, both in overall spatial extent and the minimum size of feature that can be resolved. The Q P /Q S results in Figure 2 are masked for regions where we are not confident that we can resolve features greater than 2 km in size. Significantly, we cannot observe the APMB with any confidence, primarily because only a small number of earthquakes from below the APMB are of sufficient quality to use in the attenuation tomography inversions. Although one could exclude events deeper than the APMB from the tomography inversion, it is likely that this would have minimal impact on the shallow resolution while providing no resolution deeper than 8 km bsl. The resolution of the tomography could be improved by detecting more S waves and/or increasing the number of stations in the network.
Another potential limitation is that we invert for Q P −1 and Q S −1 separately, producing our Q P /Q S tomographic map by division of each unit of the tomographic model individually, as in other studies (e.g., Hauksson & Shearer, 2006). However, this division method has two potential issues. First, both Q P −1 and Q S −1 can take near-zero values for low-attenuation regions, which could lead to large variations in Q P /Q S . Second, even if the regularization and smoothing parameters happen to be equal for both the Q P and Q S inversions, the resolution and associated uncertainty in the tomography results may differ. These issues can be mitigated by performing a direct inversion for Q P /Q S . We do not perform a direct inversion for Q P /Q S , as the complex velocity structure means that one cannot assume the same raypath for a given P wave and S wave from the same source to the same receiver, a required condition if inverting directly for Q P /Q S using our method. Others have avoided this limitation by inverting for Q P /Q S by fixing Q P using the Q P tomography solution and then only varying Q S in the inversion (Wei & Wiens, 2020). However, for our data set, absolute values of Q P /Q S are highly dependent on the starting model. Instead, to maximize confidence in our direct division-derived Q P /Q S results, we mask out regions that have small Q P −1 and/or Q S −1 (Q > 500) to minimize any near-zero division issues. We also minimize the effect of differing Q P −1 and Q S −1 resolution and uncertainties by masking out regions where either 4 km size Q P or Q S features cannot be resolved (see Figures S1 and S2 in Supporting Information S1). Furthermore, the poorer Q S resolution results in greater smearing than in Q P (see resolution tests in Supporting Information S1), leads to potential underestimation in the deviation of Q P /Q S . We therefore generally only interpret features with Q P /Q S significantly less or significantly greater than one, for which we are confident in our binary saturated versus compressible interpretation. Overall, we therefore have confidence in the remaining Q P /Q S results presented in Figures 2 and 4.
A final limitation of note is the ambiguity of Q P /Q S ratios for identifying fluid type and the exact ratio of gas to liquid. Here, we require temperature and pressure profiles to identify the most likely fluid/s associated with a given Q P /Q S ratio. However, pressure and temperature profiles may well be inaccurate, especially as they only describe variations in one-dimension, depth. V P , V S and V P /V S results could provide additional constraint, although velocity is less sensitive to fluid saturation than attenuation (Winkler & Nur, 1982). Other geophysical measurements such as density (Del Potro et al., 2013;Ward et al., 2014) would provide additional constraint, which in combination with other parameters might constrain porosity. However, one particularly useful observation would be resistivity tomography (Comeau et al., 2016), which would aid constraint of the presence of conductive brines versus resistive gases. Ideally, one would perform a joint inversion, including all these parameters to simultaneously constrain additional parameters, such as porosity, conductivity, and density, in order to identify fluid composition and prevalence better.

The Bigger Picture
Mapping what fluids accumulate where in the subsurface is important for numerous applications, including: volcanic hazard assessment; efficiently exploiting geothermal systems; and in searching for brines rich in metals critical for a green energy transition (Blundy et al., 2021). For volcanic hazard assessment, it is important to discriminate between hydrothermal and partial melt storage regions in order to accurately discriminate the volume of melt at a given volcano. Q P /Q S combined with temperature and pressure profiles makes this possible. Theoretically, our approach may also allow one to measure the proportion of melt, if densities are adequately constrained. For geothermal system characterization, Q P /Q S can elucidate and map fluid-saturated crust versus potentially dry crust, providing improved spatial constraint when targeting geothermal prospects. This constraint will reduce geothermal exploration risk. Similarly, Q P /Q S can also benefit the endeavor to find new, sustainable mineral resources in the form of metal-rich brines. Q P /Q S , in combination with temperature, pressure, density, and conductivity profiles (Comeau et al., 2016;Del Potro et al., 2013;Ward et al., 2014), not only allows for the spatial constraint of metal-rich brine accumulation but also potentially the volume fraction of these brines within the crust. Brines with higher-than-average metal contents will have higher conductivities, with such regions accumulating when at supercritical depths (Blundy et al., 2021). This information is critical for maximizing the success of extracting metal-rich brines from geothermal systems going forward.

Conclusions
We present seismic attenuation and anisotropy results at Uturuncu volcano, Bolivia. 3D attenuation tomography shows higher-than-background P wave attenuating structures that align parallel to the orientation of fractures (NNW-SSE). The presence and orientation of these fractures are evidenced by seismic anisotropy. Higher-than-background S wave attenuation is isolated to smaller localities, sporadically located along the strike of the high P wave attenuation regions. Q P /Q S ratios, indicative of crustal fluid saturation, show that most of the crust is partially saturated, with only a few pockets of fully saturated crust. We interpret saturated crust above sea-level to be brine saturated crust in its normal state, likely a brine lens that developed during active volcanism that has now ceased. A column of particularly low Q P /Q S directly beneath Uturuncu is interpreted to comprise of supercritical fluids, most likely H 2 O or metal-rich brines, potentially containing volatiles. This likely feeds CO 2 and other volatiles into a shallow hydrothermal system.
We show that high seismic attenuation features can be attributed to fluid accumulation along faults, if constrained by other observations such as seismic anisotropy and fault orientations. Furthermore, if one has constraint over pressure and temperature profiles through the crust, as well as data constraining crustal resistivity, then it is possible to identify the most likely fluid compositions responsible for each seismic attenuation signature. We therefore conclude that fluid accumulation and composition can be mapped seismically. Such observations at other locations could be used for a range of applications, from volcanic hazard assessment to the exploration of metal-rich brine deposits that could potentially be exploited to facilitate the green-energy transition.

Data Availability Statement
All the seismic data used in this study are available on the publicly accessible IRIS data repository (YS and XP network codes [see Pritchard, 2009;West & Christensen, 2010]) (https://ds.iris.edu/ds/nodes/dmc/). The earthquake catalog used is detailed in Hudson et al. (2022) and the gridded tomography data are available on Zenodo (Hudson, 2023b). The software used to measure path-average attenuation for each earthquake is SeisSrc-Moment (Hudson, 2020), available open source. The shear-wave splitting analysis was undertaken using SWSPy (Hudson, 2022), another open source package. The code used to perform the attenuation tomography is available via Zenodo (Hudson, 2023a).