Tesfahiwet Yemane

Seismic Methods for Imaging Geothermal Systems

Tesfahiwet Yemane

Mike Kendall, Jon Blundy, Thomas Hudson, Petros Bogiatzis, Sacha Lapins

James Hammond, Luca De Siena, Simona Gabriele, Ghebrebrhan Ogubazghi, Tobermory C. Mackay-Champion and others

Department of Earth Sciences, Oxford, UK

GGM | 24-02-2026

Outline

Motivation

This research is motivated by the increasing need for clean energy and critical metals (e.g., Lithium) required for the energy transition.

Lithium Demand Projection (2020-2035) (Statista, 2024).

Motivation

The dual benefit of geothermal brines (geothermal brines are hot, saline solutions from geothermal areas enriched with metals, such as Lithium).

Geothermal brines can be a source of both geothermal energy and critical metals (Yemane et al., 2025, GRL).

Motivation

The annual lithium flux in the Taupo Volcanic Zone (TVZ) could supply approximately 3% of global annual lithium demand while simultaneously producing geothermal energy.

Eight geothermal wells at TVZ, New Zealand (Lucjan et al., 2023, MDPI Resources).

Objectives

This research aims to image geothermal systems (B), including brine reservoirs (A), located above the melt zone (E), and to infer temperature/composition using seismic methods.

Schematic diagram of Geothermal Systems (Hudson et al., 2023, GRL).

Seismicity

The area of interest is the Aluto-Langano geothermal field in Ethiopia, in the MER.

Aluto-Langano Geothermal Field (Yemane et al., 2025, JGR: Solid Earth).

Seismicity

The drilled geothermal wells at Aluto reach temperatures of around 360 °C.

Geothermal wells at Aluto-Langano (Yemane et al., 2025, GRL).

Seismicity

Seismic activity at Aluto and around Aluto volcano from January 2012 to January 2014.

Seismicity at Aluto volcano (Yemane et al., 2025, GRL).

Traditional (Smith et al., 2020, GJI) vs Deep Learning (Lapins et al., 2021 (UGPD), JGR: Solid Earth).

Seismicity

Transfer learning reuses knowledge from one trained model to improve learning on a related task. We applied transfer learning similar to UGPD using a pre-trained model (PhaseNet).

P- and S-wave arrivals picked by the trained model, which performs well for high-magnitude events (Yemane et al). Events detected by the traditional method are used for further processing.

Seismicity

High b-values of 1.97 ± 0.10 at Aluto indicate the presence of fluids beneath the volcano.

The Gutenberg-Richter distribution is represented as: \( \color{#ff9900}{\log_{10} N = a - bM} \tag{1} \) It characterizes the magnitude-frequency distribution of earthquakes in a specific region (Gutenberg & Richter, 1944), where N is the number of earthquakes equal to or exceeding magnitude M, while a and b denote constants that describe the rate of seismic activity and the relationship between the rate of smaller and larger earthquakes, respectively.

b-values at Aluto-Langano geothermal field (Yemane et al., 2025, JGR: Solid Earth).

Seismicity

Focal mechanisms show normal faulting in the direction of rift extension and full-moment tensor inversions suggest shear-failure with fluids potentially activating existing faults.

Focal mechanisms and lune plots illustrating the source PDF (Yemane et al., 2025, JGR: Solid Earth).

Seismicity

The melt is situated at a shallow depth in the northwest of Aluto below SDFZ. The magmatic and hydrothermal systems are connected through faults and fracture systems.

The schematic 3D cross-section of the Aluto-Langano geothermal field (Yemane et al., 2025, JGR: Solid Earth).

Attenuation tomography

We use direct wave attenuation (using coda normalisation method), coda wave attenuation (using coda wave decay) and peak delay time (envelope broadening).

An event recorded by one of the stations (A01E) along with its signal envelope (Yemane et al., 2025, GRL).

Attenuation tomography

Coda wave attenuation measures absorption and effectively detects fluid reservoirs, high temperature areas, and melt.

We use the energy envelope decay of coda, described by Aki and Chouet (1975) as a function of lapse time from the event origin time (t) and frequency (f), which is given by: \(E(t,f)=S(f){t}^{-\alpha }\,\mathrm{exp}\left(\frac{-2\pi ft}{{Q}_{c}(f)}\right),\) where E(t,f) is the power spectral energy density, S(f) includes frequency-dependent source and/or site term, t is the lapse-time from the origin of the event, α is a constant related to geometrical spreading. The equation can be linearized by taking the logarithm and rearranged so that \[\frac{\ln\left[E(t,f)\cdot {t}^{\alpha }\right]}{2\pi f}=\frac{\ln\left[S(f)\right]}{2\pi f}-\frac{1}{{Q}_{c}(f)}t.\]

Frequency dependence of coda quality factor (Qc) (Yemane et al., 2025, GRL).

Attenuation tomography

We use the MLTWA (Fehler, 1992) method to calculate seismic albedo (\(B_{0}\)) and inverse extinction length (\(L_{e}^{-1}\)), which are input parameters to the sensitivity kernel.

\[ E_k(r) = \int_{t_k}^{t_k + \Delta t} E(r,t) \, dt, \quad k = 1, \dots, 3, \label{eq:sieq2} \]

where \(E_k(r)\) is the time energy integral, \(E(r,t)\) is the seismic energy density at distance \(r\) and time \(t\), \(t_k\) is the start time of the k-th window, and \(\Delta t\) is the duration of the window.

\[ E^{\text{obs}}_k(r) = \log_{10} \left( 4 \pi r^2 \frac{E_k(r)}{\int_{t_{\text{coda}}}^{t_{\text{coda}} + \Delta t} E_k(r, t) \, dt} \right), \label{eq:sieq3} \]

\[ \begin{split} E_{i,j}^{3D} \left[ r_{ij}, t \right] &\approx \frac{W_0 \exp \left( -L_e^{-1} v t \right)}{4 \pi r_{ij}^2 v} \delta \left[ t - \frac{r_{ij}}{v} \right] + W_0 H \left( t - \frac{r_{ij}}{v} \right) \\ &\quad \frac{\left( 1 - \frac{r_{ij}^2}{v^2 t^2} \right)^{\frac{1}{8}}}{\left( \frac{4 \pi v t}{3 B_0 L_e^{-1}} \right)^{\frac{3}{2}}} \cdot \exp \left( -L_e^{-1} v t \right) F \left[ v t B_0 L_e^{-1} \left( 1 - \frac{r_{ij}^2}{v^2 t^2} \right)^{\frac{3}{4}} \right] , \label{eq:sieq4} \end{split} \]

where \(E_{i,j}^{3D} \left[ r_{ij}, t \right]\) is the approximation of the energy transport equation (estimate the energy envelope of a seismogram), \(W_0\) is the source energy, \(v\) is the seismic velocity (3.5 km/s), \(\delta\) is the Dirac's delta function, \(H\) is the Heaviside step function, and \(F\) is given by \(F[x] = e^x \sqrt{1 + \frac{2.026}{x}}\).

A plot of the logarithm of the normalized energy integral versus hypocentral distance (in km) (Yemane et al., 2025, GRL).

Attenuation tomography

The best fitting values of \(B_{0}\) and \(L_{e}^{-1}\) are obtained by minimizing the total misfit between the observed and theoretical energies.

\[ M(L_e^{-1}, B_0) = \sum_{i=1}^{N} \sum_{k=1}^{3} \left[ E^{\text{obs}}_k(r_i) - E^{\text{theo}}_k(r_i, L_e^{-1}, B_0) \right]^2 , \label{eq:sieq6} \]

\[ K^{3D}_{i,j} \left[ \rho, t, B_0, L_e^{-1}, v \right] = \int_0^T E^{3D} \left[ r_{s\rho}, \tau, B_0, L_e^{-1}, v \right] \times E^{3D} \left[ r_{\rho r}, T - \tau, B_0, L_e^{-1}, v \right] \, d\tau, \label{eq:sieq7} \]

where \(\rho\) is a point in space with coordinates \(\{i,j,z\}\), \(r_{\rho r}\) and \(r_{s\rho}\) are point-to-receiver and source distance respectively, and \(\tau\) is the integration variable. The kernel functions are used to construct the rows of the inversion matrix of the grid which is given by:

d = Gm \[ \begin{bmatrix} Q_{c_1} \\ Q_{c_2} \\ \vdots \\ Q_{c_N} \\ \end{bmatrix} = \begin{bmatrix} k_{11} & k_{12} & \cdots & k_{1M} \\ k_{21} & k_{22} & \cdots & k_{2M} \\ \vdots & \vdots & \ddots & \vdots \\ k_{N1} & k_{N2} & \cdots & k_{NM} \\ \end{bmatrix} \begin{bmatrix} Q_{c}^{1} \\ Q_{c}^{2} \\ \vdots \\ Q_{c}^{M} \\ \end{bmatrix} \]

Here, \(K_{ij}\) represents the source-receiver pair, with \(i = 1, \ldots, N\), and \(j=1, \ldots, M\). \(Q_{c_N}\) is the \(Q_c\) measurement for each source-receiver pair and \(Q_{c}^{M}\) is \(Q_c\) associated with each block of the grid.

Best fits for seismic albedo \(B_{0}\) and inverse extinction length \(L_{e}^{-1}\) (Yemane et al., 2025, GRL).

Attenuation tomography

Peak delay time measures scattering and is sensitive to lithlogical variations and structural features (faults and fracture systems).

For each frequency band \(f\), the dependence of peak-delay on hypocentral distance \((R, \text{in km})\) is given by: \( \log_{10}(T_{pd}^t(f)) = A_r(f) + B_r(f) \log_{10}(R) \) where \(A_r(f)\) and \(B_r(f)\) are the coefficients of a linear fit. The variation in peak delay is calculated as the difference between the measured peak delay time \(T_{Pd}^{obs}\) and the theoretical peak delay \(T_{Pd}^t\) at the corresponding hypocentral distance: \[ \Delta \log_{10}(T_{pd}(f)) = \log_{10}(T_{pd}^{obs}(f)) - \log_{10}(T_{pd}^t(f)) \] The variation in peak delay time is mapped using a regionalisation approach, which uses the peak delay time measurements.

Logarithmic plot of peak-delay (in seconds) as a function of travel time (in seconds) (Yemane et al., 2025, GRL).

Attenuation tomography

The weighted average value of peak delay is allocated to each block (blocks intersected by at least five rays are considered) (De Siena et al., 2016, EPSL) .

Ray path density map for peak delay at 0 km bsl at 3 Hz, showing the number of rays intersecting each block (Yemane et al., 2025, GRL).

Attenuation tomography

High Qc^-1 anomalies align with high-temperature areas and areas of elevated CO₂. High T_pd anomalies correspond to structural features which act as fluid pathways.

Absorption (Yemane et al., 2025, GRL)

Scattering (Yemane et al., 2025, GRL)

Resistivity (Hochstein et al., 2017, GT)

Attenuation tomography

High enthalpy wells correlate with high Qc^-1 and mass flow rate and circulation loss correlate with high T_pd. LA-3 and LA-6 are the hottest wells, while LA-4 and LA-7 are the most permeable.

Location of geothermal wells (Yemane et al., 2025, GRL)

Productive vs Non-Productive areas at Aluto-Langano

Absorption vs Scattering (Yemane et al., 2025, GRL)

Travel-time Tomography

We use final frequency kernels (Bogiatzis et al., 2022) instead of ray paths with a direct sparse solver for 445 seismic events (2176 P and 1165 S arrivals).

Seismic events used for travel-time tomography (Yemane et al., 2026, GJI).

Travel-time Tomography

The initial velocity model for travel-time tomography and the ray path coverage. We have a very good ray coverage of both P and S waves.

Velocity model (Yemane et al., 2026, GJI).

P and S ray path coverage (Yemane et al., 2026, GJI).

Travel-time Tomography

Elevated Vp/Vs ratios at around 0 km bsl near productive geothermal wells and across faults. indicating high fluid content and/or high temperature.

Vp model of the last iteration

Vs model of the last iteration

Vp/Vs model of the last iteration (Yemane et al., 2026, GJI)

Travel-time Tomography

The Derivative Weight Sum (DWS) shows good ray coverage. across different depth slices particularly for P-waves, with sensitivity decreasing at greater depths.

DWS for Vp

DWS for Vs

Vp model of the last iteration (Yemane et al., 2026, GJI)

Travel-time Tomography

The full resolution matrix (RDE) from the final iteration indicates generally high resolution. in the centre of the caldera, particularly for P-waves.

Resolution for Vp

Resolution for Vs

Vp model of the last iteration (Yemane et al., 2026, GJI)

Travel-time Tomography

The off-diagonal elements of the resolution matrix show low smearing in both horizontal and vertical sections at different depths. We have good resolution because of good ray coverage.

Horizontal smearing (2 km)

Horizontal smearing (at different location)

Vertical smearing (2 km) (Yemane et al., 2026, GJI)

Travel-time Tomography

The methods reveal a partial melt body at around 7 km bsl M , feeding shallower reservoirs; supercritical fluid (S), brine (H), fluid-rich fracture (F), and steam (C), migrating along faults.

Vp, Vs, Vp/Vs, absorption, scattering and conductivity models beneath Aluto-Langano (Yemane et al., 2026, GJI).

Travel-time Tomography

Seismic tomography and other methods revelal a partial melt M, supplying fluids/volatiles to S, H, F, and C, with migration along faults and fractures.

Schematic model of the geothermal system at Aluto-Langano (Yemane et al., 2026, GJI).

Ambient Noise Tomography

We perform cross-correlation and stacking using data from all stations at Aluto-Langano to derive phase velocities and shear-wave velocities.

Seismic stations used for ambient noise tomography (in Prep).

Ambient Noise Tomography

Vertical cross-correlations of all station pairs computed on five-minute-long windows and stacked over two years sorted according to the interstation distance.

Stacked vertical cross-correlations (in Prep).

Ambient Noise Tomography

The average phase velocity is estimated by minimizing the misfit between the observed wrapped phase from all station pairs and the theoretical phase of the Bessel function.

Observed amplitude decay compared to zero-order Bessel function (in Prep).

Ambient Noise Tomography

S-wave velocity models are obtained by extracting and inverting one-dimensional dispersion curves at all target locations (from phase-velocity tomography maps).

Vs depth profiles for each chain compared to the reference model (in Prep).

Posterior Vs models with the reference model (in Prep).

Ambient Noise Tomography

We compute Rayleigh-wave sensitivity kernels to assess the depth sensitivity of the inverted surface-wave models.

Rayleigh-wave phase velocity sensitivity kernels (in Prep).

Ambient Noise Tomography

Shear wave velocity models show low shear velocities near productive geothermal wells at approximately 2 km below sea level.

Shear wave velocity models at 0 and 2 km bsl and across longitude and latitude along the center of the caldera (in Prep).

Future Work

We use all the geophysical observations to estimate temperature of the reservoir using Gaussian Processes.

Temperature estimation using Gaussian Processes.

Future work

The temperature along the latitude and longitude cross-sections shows high temperatures around the productive geothermal wells.

Temperature cross-sections of the Aluto–Langano Geothermal Field.

Future work

Uncertainity of temperature measurements at different depth slices of the Aluto–Langano Geothermal Field.

Uncertainty of temperature measurements.

Future work

The geothermal power potential is linearly related to both temperature and permeability estimated from scattering.

Prospectivity index of the Aluto–Langano Geothermal Field.

Future work

Bulk conductivity of 5% porous granite saturated with H₂O–NaCl fluids under supercritical geothermal conditions (10 wt% NaCl) superimposed on H₂O–NaCl phase diagrams.

Bulk conductivity estimation to compare with MT measurements. We also plan to do Rock Physics modeling to estimate the fluid and gas saturations.

Conclusion

Aluto is seismically active, focal mechanisms are normal and some with left lateral strike-slip. High b-values suggest the presence of fluids.
The magmatic and hydrothermal systems are connected through pre-existing faults and fractures.
Absorption mapping delineates areas of high fluid content and high-temperature regions.
Peak delay time mapping delineates structural (faults and fracture systems) and lithological features.
Seismic tomography and other methods reveal a partial melt (M) supplying fluids/volatiles to S, H, F and C, with migration along faults and fractures.
Seismic methods can be used as tools for geothermal exploration especially if integrated with other geophysical and geological data.
In the future we plan to do a multiparameter geophysics to estimate temperature and if possible to differentiate melt from brine.