API Reference¶
Complete API documentation for gdrift. All classes and functions are organized by functionality.
Core Data Loading¶
load_dataset¶
load_dataset ¶
Load a dataset from local cache, downloading if necessary.
Args: dataset_name (str): Dataset name (without .h5 extension) table_names (list, optional): Specific tables to load. Defaults to []. return_metadata (bool, optional): Whether to return file-level metadata.
Returns: dict: dictionary with all the datasets (and optionally metadata tuple)
1D Radial Profiles¶
RadialEarthModel¶
RadialEarthModel ¶
Class representing reference Earth Models such as PREM or AK135 Composite object containing multiple radial profiles representing different Earth properties, such as shear wave velocity (Vs), primary wave velocity (Vp), and density. T
Attributes: depth_profiles (dict): A dictionary of RadialProfile instances.
Initialise the RadialEarthModel with a dictionary of radial profiles instances.
Args: profiles (dict of RadialProfile): Profiles for different properties, keyed by property name.
get_profile ¶
Retrieve a profile by its name.
Args: property_name (str): The name of the property profile to retrieve.
Returns: AbstractProfile: The profile corresponding to the specified property name.
Raises: ValueError: If the specified property name does not exist in the model.
RadialEarthModelFromFile¶
RadialEarthModelFromFile ¶
Bases: RadialEarthModel
A class for loading radial profiles from a dataset.
This class extends SplineProfile to specifically handle loading,
and utilizing available profiles related to profiles in the mantle.
Attributes: model_name (str): The name of the model/dataset from which profiles are loaded. description (str, optional): A brief description of the profile's purpose or characteristics.
Initialize a radial Earth model by loading profiles from an HDF5 dataset.
Loads all property profiles from a registered dataset file. The HDF5 file must contain a "depth" array and one or more property arrays (e.g., "density", "vs", "vp"). Each property is wrapped in a SplineProfile for interpolation.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model_name
|
str
|
Name of the registered dataset (e.g., "1d_prem", "1d_solidus_Andrault_et_al_2011_EPSL"). Must exist in the dataset registry. |
required |
description
|
str
|
Human-readable description of the model. If None, no description is set. Default is None. |
None
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If model_name is not in the dataset registry. |
KeyError
|
If the HDF5 file does not contain a "depth" array. |
Examples:
>>> import gdrift
>>> andrault = gdrift.RadialEarthModelFromFile(
... "1d_solidus_Andrault_et_al_2011_EPSL")
>>> print(andrault.get_profile_names())
['solidus temperature']
>>> solidus = andrault.get_profile("solidus temperature")
>>> T = solidus.at_depth(410e3)
PreliminaryRefEarthModel¶
PreliminaryRefEarthModel ¶
Bases: RadialEarthModelFromFile
Initialises the Preliminary Reference Earth Model (PREM). This model is based on the work by Dziewonski and Anderson (1981) and provides a reference Earth model that can be queried at specific depths for various profiles.
References: Dziewonski, Adam M., and Don L. Anderson. "Preliminary reference Earth model." Physics of the Earth and Planetary Interiors 25.4 (1981): 297-356.
The object is of type RadialEarthModel and is initialized by loading profiles from an existing dataset. Each profile is represented as a SplineProfile object.
Attributes: prem_profiles (list): A list of SplineProfile objects representing different profiles in the PREM dataset.
Initialize the Preliminary Reference Earth Model (PREM).
Loads the PREM dataset (Dziewonski & Anderson, 1981) containing reference profiles for density, seismic velocities, elastic moduli, pressure, and gravity as a function of radius/depth.
The model spans from Earth's center to the surface and includes discontinuities at major boundaries (e.g., core-mantle boundary, 410 km, 660 km discontinuities).
Notes
PREM is the standard 1D reference model for seismology and geodynamics. It represents a spherically symmetric, non-rotating, oceanless Earth.
Examples:
>>> import gdrift
>>> prem = gdrift.PreliminaryRefEarthModel()
>>> print(prem.get_profile_names())
['density', 'vs', 'vp', ...]
>>> rho_profile = prem.get_profile("density")
>>> rho_670km = rho_profile.at_depth(670e3)
References
Dziewonski, A. M., & Anderson, D. L. (1981). Preliminary reference Earth model. Physics of the Earth and Planetary Interiors, 25(4), 297-356. https://doi.org/10.1016/0031-9201(81)90046-7
HirschmannSolidus¶
HirschmannSolidus ¶
Bases: RadialEarthModel
A class representing the solidus model based on the work of Hirschmann (2000).
Attributes: solidus_profile (HirschmannSolidusProfile): The solidus profile based on Hirschmann (2000).
Initialize a Hirschmann solidus radial Earth model.
Creates a RadialEarthModel containing a single profile: the Hirschmann (2000) solidus temperature. This is a convenience wrapper around HirschmannSolidusProfile that conforms to the RadialEarthModel interface.
The solidus can be accessed via get_profile("solidus temperature") or
directly through the solidus_profile attribute.
Examples:
>>> import gdrift
>>> solidus_model = gdrift.HirschmannSolidus()
>>> solidus_profile = solidus_model.get_profile("solidus temperature")
>>> T_410km = solidus_profile.at_depth(410e3)
>>> print(f"Solidus at 410 km: {T_410km:.1f} K")
See Also
HirschmannSolidusProfile : The underlying solidus implementation
3D Seismic Models¶
SeismicModel¶
SeismicModel ¶
SeismicModel(model_name, nearest_neighbours: int = 8, default_max_distance: float = 200000.0, labels=[])
Bases: EarthModel3D
SeismicModel is a class for handling 3D seismic models.
This class inherits from EarthModel3D and is used to load and manage seismic models. It allows for the initialization of a seismic model with a specified name, number of nearest neighbours, and a default maximum distance. The model data is loaded from a dataset file, and quantities and coordinates are set accordingly.
Parameters:
model_name : str Name of the seismic model to load nearest_neighbours : int, optional Number of nearest neighbours to consider for interpolation. Defaults to 8. default_max_distance : float, optional The default maximum distance for the model in meters. Defaults to 200e3. labels : list, optional Specific labels to load from the dataset. If empty, loads all available fields.
Attributes: model_name (str): The name of the seismic model. nearest_neighbours (int): The number of nearest neighbours to consider. default_max_distance (float): The default maximum distance for the model in meters.
The full list of available tomography models is exposed as gdrift.AVAILABLE_SEISMIC_MODELS.
See the Data Catalog for descriptions and citations.
3D Earth Models (not seismic)¶
3D models of a mantle property that is not a seismic wave speed, for example the mantle
viscosity of Lloyd et al. (2024), load directly into an EarthModel3D. Their names in the
registry start with 3d_earthmodel_. Give the name without this prefix:
import gdrift
visc = gdrift.EarthModel3D.from_dataset("Lloyd_et_al_2024_GJI")
xyz = gdrift.geodetic_to_cartesian(-75.0, -100.0, 200e3) # lat, lon (degrees), depth (m)
log10_eta = visc.at("log10_viscosity", xyz) # log10(Pa s)
A field that spans many orders of magnitude is stored and interpolated as its log10, and the
name of the field says so. The list of these models is exposed as
gdrift.AVAILABLE_EARTH_MODELS_3D.
EarthModel3D¶
EarthModel3D ¶
Bases: AbstractEarthModel
Three-dimensional Earth model with KD-tree spatial interpolation.
A container for 3D gridded geophysical data with efficient spatial queries using scipy's cKDTree. Supports multiple interpolation kernels (IDW, Gaussian, Wendland, linear, cubic, nearest neighbor) and handles both Cartesian and geographic coordinate systems.
Typical workflow:
1. Create model: model = EarthModel3D(nearest_neighbours=8)
2. Set coordinates: model.set_coordinates(x, y, z)
3. Add quantities: model.add_quantity("dvs", dvs_array, "label")
4. Query: value = model.at(lat=45, lon=0, depth=500e3, quantity="dvs")
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
nearest_neighbours
|
int
|
Number of nearest neighbors to use for interpolation. Higher values produce smoother interpolation but increase computational cost. Default is 8 (suitable for most 3D tomography models). |
8
|
default_max_distance
|
float
|
Maximum search distance in meters for finding neighbors. Points farther than this from all data points return NaN. Default is 200e3 (200 km), suitable for typical mantle convection resolution. |
200000.0
|
Attributes:
| Name | Type | Description |
|---|---|---|
coordinates |
ndarray or None
|
Nx3 array of (x, y, z) coordinates in normalized units (scaled to
R_earth). Set via |
available_fields |
dict
|
Dictionary mapping quantity names to (data_array, label) tuples.
Populated via |
tree |
cKDTree or None
|
KD-tree for fast spatial queries. Built automatically when
coordinates are set. None until |
nearest_neighbours |
int
|
Number of neighbors for interpolation (from initialization). |
default_max_distance |
float
|
Maximum search distance in meters (from initialization). |
Examples:
>>> import gdrift
>>> import numpy as np
>>> # Create synthetic 3D model
>>> model = gdrift.EarthModel3D(nearest_neighbours=8)
>>> x = np.linspace(-0.5, 0.5, 10) # normalized coords
>>> y = np.linspace(-0.5, 0.5, 10)
>>> z = np.linspace(-0.5, 0.5, 10)
>>> xx, yy, zz = np.meshgrid(x, y, z, indexing='ij')
>>> coords = np.column_stack([xx.ravel(), yy.ravel(), zz.ravel()])
>>> model.set_coordinates(coords[:, 0], coords[:, 1], coords[:, 2])
>>> # Add velocity perturbation
>>> dvs = np.random.randn(len(coords)) * 0.02
>>> model.add_quantity("dvs", dvs, "dln(Vs)")
>>> # Query using geographic coordinates
>>> value = model.at(lat=45.0, lon=10.0, depth=1000e3, quantity="dvs")
See Also
SeismicModel : Subclass for loading seismic tomography models interpolate_to_points : Underlying interpolation engine
Initialize a 3D Earth model with interpolation parameters.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
nearest_neighbours
|
int
|
Number of nearest neighbors for interpolation. Default is 8. |
8
|
default_max_distance
|
float
|
Maximum search distance in meters. Default is 200e3 (200 km). |
200000.0
|
from_dataset
classmethod
¶
Load a 3D model that is not seismic tomography from the dataset registry.
The dataset is "3d_earthmodel_
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model_name
|
str
|
Name without the "3d_earthmodel_" prefix, one of
|
required |
nearest_neighbours
|
int
|
Number of neighbours for the interpolation. Default is 8. |
8
|
default_max_distance
|
float
|
Distance in metres beyond which a query returns NaN. Default is 200e3, which is more than twice the lateral spacing of the 65,341-point layers that the gridded datasets use (about 88 km). Within this distance above the top layer or below the bottom layer, a query returns the value of that layer (no extrapolation). For the viscosity model of Lloyd et al. (2024), this means that points up to 200 km above the surface get the crust value, and points up to 200 km below 2891 km get the value of the deepest layer. |
200000.0
|
labels
|
list of str
|
Load only these quantities. Default is all quantities. |
None
|
Returns:
| Type | Description |
|---|---|
EarthModel3D
|
The model with its coordinates and quantities set. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
Notes
The fields are interpolated as stored. A field that spans many
orders of magnitude is therefore stored as its log10 (for example
"log10_viscosity" in log10(Pa s)), and the caller converts the
interpolated value back with 10 ** value if needed.
The lateral resolution is that of the stored layers (65,341 points, about 88 km apart), which can be coarser than the original model. The description of each dataset in the registry says how it was resampled.
at ¶
at(label: Union[str, List[str]], coordinates: array, kernel=None, extrapolate=False, **kernel_params)
Get the value of a quantity at specified coordinates using various interpolation kernels.
Parameters:
label : str or list
Field label(s) to interpolate
coordinates : np.ndarray
Query coordinates (N x 3 array)
kernel : str, optional
Interpolation kernel. If None, uses the class attribute
default_kernel (falls back to 'idw'). Options:
- 'idw': Inverse distance weighting
- 'gaussian': Gaussian kernel with adaptive or fixed bandwidth
- 'idw_power': IDW with adjustable power parameter
- 'exponential': Exponential decay kernel
- 'wendland': Wendland compactly supported kernel
extrapolate : bool, optional
If False (default), query points whose nearest data point is
farther than default_max_distance return NaN. If True,
the nearest neighbors are used regardless of distance.
**kernel_params : dict
Kernel-specific parameters:
- For 'gaussian': sigma (bandwidth, auto-computed if not provided)
- For 'idw_power': power (default 2.0)
- For 'exponential': decay_length (default 50000m)
- For 'wendland': support_radius (default 100000m)
Returns:
np.ndarray Interpolated values at query coordinates
Thermodynamics¶
ThermodynamicModel¶
ThermodynamicModel ¶
Bases: object
Thermodynamic lookup table for mantle mineral physics properties.
Provides 2D interpolation of pre-computed mineral physics properties (density, seismic velocities, elastic moduli) as a function of depth and temperature. Tables are based on self-consistent thermodynamic databases (Stixrude & Lithgow-Bertelloni 2011, 2021) and stored in HDF5 format with bivariate spline interpolation.
The model supports: - Forward queries: (depth, temperature) → property (e.g., Vs, Vp, rho) - Inverse queries: (depth, velocity) → temperature - Multiple compositions and database versions (see MODELS_AVAIL, COMPOSITIONS_AVAIL) - On-the-fly computation of Vs/Vp from elastic moduli
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model
|
str
|
Thermodynamic database version (e.g. "SLB_21"). See MODELS_AVAIL for all available versions, derived from the datasets.json manifest. |
required |
composition
|
str
|
Mantle composition. Available options depend on the model. See COMPOSITIONS_AVAIL for all available compositions (e.g. "pyroliteCFMAS", "pyroliteNCMAS", "bulk-oceanic-crustCFMAS"). |
required |
temps
|
array_like
|
Temperature grid for subsampling (Kelvin). If None, uses full temperature range from dataset. Default is None. |
None
|
depths
|
array_like
|
Depth grid for subsampling (meters). If None, uses full depth range from dataset. Default is None. |
None
|
extrapolate
|
bool
|
Whether to allow extrapolation outside the table bounds. If False, queries outside the range raise ValueError. Default is False. |
False
|
Attributes:
| Name | Type | Description |
|---|---|---|
model |
str
|
Database version (e.g., "SLB_21"). |
composition |
str
|
Mantle composition (e.g., "pyroliteCFMAS"). |
extrapolate |
bool
|
Extrapolation flag. |
_tables |
dict
|
Dictionary mapping property names to Table objects. Keys are the property names from the HDF5 file (e.g., "rho", "bulk_mod", "shear_mod", "v_s", "v_p"). |
Examples:
>>> import gdrift
>>> # Load SLB 2021 pyrolite model
>>> tm = gdrift.ThermodynamicModel("SLB_21", "pyroliteCFMAS")
>>> # List available properties
>>> print(tm.available_tables())
['rho', 'bulk_mod', 'shear_mod', 'v_s', 'v_p', ...]
>>> # Forward query: temperature to Vs
>>> vs = tm.temperature_to_vs(temperature=1600, depth=500e3)
>>> print(f"Vs = {vs:.1f} m/s")
>>> # Inverse query: Vs to temperature
>>> T = tm.vs_to_temperature(vs=4500, depth=500e3)
>>> print(f"Temperature = {T:.1f} K")
>>> # Generic property lookup
>>> rho = tm.temperature_to_property("rho", temperature=1600, depth=500e3)
Notes
- All depths are in meters from the surface
- All temperatures are in Kelvin
- Velocities are in m/s, densities in kg/m³, moduli in Pa
- Not all model/composition combinations are available
- Use
available_tables()to see which properties are loaded - Vs and Vp are computed from moduli if not present in HDF5 file
See Also
regularise_thermodynamic_table : Smooth phase transitions apply_anelastic_correction : Convert to seismic-frequency velocities compute_swave_speed : Calculate Vs from shear modulus and density compute_pwave_speed : Calculate Vp from bulk/shear moduli and density
References
Stixrude, L., & Lithgow-Bertelloni, C. (2011). Thermodynamics of mantle minerals—II. Phase equilibria. Geophysical Journal International, 184(3), 1180-1213. https://doi.org/10.1111/j.1365-246X.2010.04890.x
Stixrude, L., & Lithgow-Bertelloni, C. (2021). Thermal expansivity, heat capacity and bulk modulus of the mantle. Geophysical Journal International, 228(2), 1119-1149. https://doi.org/10.1093/gji/ggab394
temperature_to_property ¶
Convert temperature and depth to a material property value.
Args: property_name (str): Name of the property (e.g. 'rho', 'alpha', 'Cp', 'vs', 'vp'). temperature: Temperature value(s). depth: Depth value(s).
Returns: Interpolated property value(s) at the given temperature and depth. Returns NaN if extrapolate=False and inputs are out of bounds.
temperature_to_vs ¶
Convert temperature and depth to shear wave velocity (Vs).
Convenience wrapper for temperature_to_property("vs", ...).
Computes Vs from shear modulus and density using the relation:
Vs = sqrt(shear_modulus / density).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
temperature
|
float or array_like
|
Temperature in Kelvin. |
required |
depth
|
float or array_like
|
Depth in meters from the surface. |
required |
Returns:
| Type | Description |
|---|---|
float or ndarray
|
Shear wave velocity in m/s. Returns NaN for out-of-bounds queries if extrapolate=False. |
Examples:
temperature_to_vp ¶
Convert temperature and depth to compressional wave velocity (Vp).
Convenience wrapper for temperature_to_property("vp", ...).
Computes Vp from bulk modulus, shear modulus, and density using:
Vp = sqrt((bulk_modulus + 4/3 * shear_modulus) / density).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
temperature
|
float or array_like
|
Temperature in Kelvin. |
required |
depth
|
float or array_like
|
Depth in meters from the surface. |
required |
Returns:
| Type | Description |
|---|---|
float or ndarray
|
Compressional wave velocity in m/s. Returns NaN for out-of-bounds queries if extrapolate=False. |
Examples:
temperature_to_rho ¶
Convert temperature and depth to density.
Convenience wrapper for temperature_to_property("rho", ...).
Queries the pre-computed density lookup table.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
temperature
|
float or array_like
|
Temperature in Kelvin. |
required |
depth
|
float or array_like
|
Depth in meters from the surface. |
required |
Returns:
| Type | Description |
|---|---|
float or ndarray
|
Density in kg/m³. Returns NaN for out-of-bounds queries if extrapolate=False. |
Examples:
vs_to_temperature ¶
vs_to_temperature(vs: Number, depth: Number, bounds: Optional[Union[Tuple[float, float], Tuple[ndarray, ndarray]]] = (300, 7000)) -> Number
Convert S-wave velocity (vs) to temperature at a given depth.
Parameters:
vs : Number The S-wave velocity. depth : Number The depth at which the temperature is to be calculated. bounds : Optional[Union[Tuple[float, float], Tuple[numpy.ndarray, numpy.ndarray]]], default=(300, 7000) The bounds for the temperature calculation. It can be a tuple of floats or numpy arrays.
Returns:
Number The temperature corresponding to the given S-wave velocity and depth.
vp_to_temperature ¶
vp_to_temperature(vp: Number, depth: Number, bounds: Optional[Union[Tuple[float, float], Tuple[ndarray, ndarray]]] = (300, 7000)) -> Number
Convert P-wave velocity (vp) to temperature at a given depth.
Parameters:
vp : Number The P-wave velocity. depth : Number The depth at which the temperature is to be calculated. bounds : Optional[Union[Tuple[float, float], Tuple[numpy.ndarray, numpy.ndarray]]], default=(300, 7000) The bounds for the temperature calculation. It can be a tuple of floats or numpy arrays.
Returns:
Number The temperature corresponding to the given P-wave velocity and depth.
Phase Transition Regularization¶
regularise_thermodynamic_table ¶
regularise_thermodynamic_table(slb_pyrolite: ThermodynamicModel, temperature_profile: AbstractProfile, regular_range: Dict[str, Tuple] = default_regular_range)
Return a copy of a thermodynamic model with phase-transition jumps in temperature removed.
The density, S-wave and P-wave speed tables are regularised one by one
with derive_then_integrate. In each depth row, every temperature cell
whose slope dV/dT is outside regular_range gets a slope interpolated
from accepted neighbouring cells, and the row is integrated back along
temperature. The integration constant is fixed by the anchor condition
V_reg(depth, T_anchor(depth)) = V_raw(depth, T_anchor(depth)),
with T_anchor = temperature_profile.at_depth(depth). This holds exactly,
also when T_anchor is between two table nodes. In a depth row where no
slope is replaced, the regularised row equals the raw row.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
slb_pyrolite
|
ThermodynamicModel
|
The thermodynamic model to regularise. |
required |
temperature_profile
|
AbstractProfile
|
1D profile of the anchor temperature in kelvin, normally the radial
mean temperature of the convection model. At this temperature the
regularised and the raw tables agree. Keep it inside the temperature
range of the table: outside it, the raw value at the anchor is NaN
unless the model was created with |
required |
regular_range
|
dict
|
Maps "rho", "v_s" and "v_p" to a tuple (lower, upper) of accepted
values of the temperature derivative of that table. Defaults to
|
default_regular_range
|
Returns:
| Type | Description |
|---|---|
RegularisedThermodynamicModel
|
A subclass of ThermodynamicModel that returns the regularised density, S-wave and P-wave speed tables. |
Anelasticity¶
CammaranoAnelasticityModel¶
CammaranoAnelasticityModel ¶
CammaranoAnelasticityModel(B: Callable, g: Callable, a: Callable, solidus: SplineProfile, Q_bulk: Callable = lambda x: 10000, omega: Callable = lambda x: 1.0)
Bases: BaseAnelasticityModel
Anelasticity model that follows the homologous temperature approach of Cammarano et al. (2003).
The shear quality factor is Q_S = B * omega^a * exp(a * g * T_m / T), where
T_m is the solidus (melting) temperature at the given depth. Cammarano et al.
(2003) assume a constant bulk quality factor Q_K of 1000 in the upper mantle
and 10,000 in the lower mantle. The default Q_bulk here is 10,000 at all
depths, which is higher than the paper value in the upper mantle.
Initialize the model with the given parameters.
Args: B (Callable): Scaling factor for the Q model. g (Callable): Activation energy parameter. a (Callable): Frequency dependency parameter. solidus (SplineProfile): Solidus temperature profile for mantle. Q_bulk (Callable): Bulk quality factor (default is 10000). omega (Callable): Seismic frequency (default is 1).
from_q_profile
classmethod
¶
Create a CammaranoAnelasticityModel from a predefined Q-profile.
Uses the six Q-profiles (Q1-Q6) from Cammarano et al. (2003), with upper/lower mantle parameter switching at 660 km depth. The Q factor is computed as:
Q = B * omega^a * exp(a * g * T_solidus / T)
where B controls the overall attenuation amplitude (related to grain size), g is a dimensionless activation energy parameter, a is the frequency exponent, and T_solidus is the solidus temperature at the given depth.
All six profiles and their values are from Table 1 of Cammarano et al. (2003). B is in s^a.
======= ============= ============= Profile B (UM / LM) g (UM / LM) ======= ============= ============= Q1 0.5 / 10 20 / 10 Q2 0.8 / 15 20 / 10 Q3 1.1 / 20 20 / 10 Q4 0.035 / 2.25 30 / 15 Q5 0.056 / 3.6 30 / 15 Q6 0.077 / 4.95 30 / 15 ======= ============= =============
UM = upper mantle (< 660 km), LM = lower mantle (>= 660 km). The profiles test a range of g. Cammarano et al. (2003) take g of about 30 for olivine in the uppermost mantle (Karato, 1993) and g = 10-20 as more appropriate for the lower mantle (Karato & Karki, 2001). Higher g increases the temperature and pressure dependence of Q_S. For each g, B was adjusted so that Q_S along a reference adiabat with a 1300 C potential temperature agrees with radial seismic attenuation models. All profiles use a = 0.2 (frequency exponent) and omega = 1.0 Hz, as in the paper. The paper puts the upper/lower mantle jump in Q_S at the ringwoodite-perovskite transition, which is 660 km here.
Args: q_profile (str): One of "Q1" through "Q6".
Returns: CammaranoAnelasticityModel: The configured model.
Raises: ValueError: If q_profile is not one of Q1-Q6.
GoesAnelasticityModel¶
GoesAnelasticityModel ¶
Bases: BaseAnelasticityModel
Anelasticity model following Goes et al. (2000, JGR, Appendix A).
Uses activation energy and activation volume to compute the shear quality factor Q_mu via:
Q_mu = A * omega^a * exp(a * (H* + P * V*) / (R * T))
where A is a pre-exponential factor, H* is activation energy (J/mol), V* is activation volume (m^3/mol), P is pressure from PREM, R is the gas constant, and T is temperature (K).
Goes et al. (2000) apply it to upper-mantle temperatures at 50-200 km depth. Below 660 km depth, Q is set to a very large value (effectively no attenuation).
Initialize the Goes anelasticity model.
Args: A (float): Pre-exponential scaling factor. H_star (float): Activation energy in J/mol. V_star (float): Activation volume in m^3/mol. a (float): Frequency exponent. omega (float): Seismic frequency in Hz (default 1.0). Q_bulk (float): Constant bulk Q (default 1000). Goes et al. (2000) take Q_K = 1000 as a round value of the upper-mantle Q_K = 943 of model QL6 (Durek & Ekström, 1996). max_depth (float): Depth below which Q is set very high (default 660e3 m).
from_q_profile
classmethod
¶
Create a GoesAnelasticityModel from a predefined Q-profile.
Uses the Goes et al. (2000) parameterization from Table A2 where the Q factor is computed as:
Q_mu = A * omega^a * exp(a * (H* + P * V*) / (R * T))
======= ====== ========== ============= =============== Profile a A H* (kJ/mol) V* (cm^3/mol) ======= ====== ========== ============= =============== Q1 0.15 0.148 500 20 Q2 0.25 2.0e-4 584 21 ======= ====== ========== ============= ===============
Args: q_profile (str): One of "Q1" or "Q2".
Returns: GoesAnelasticityModel: The configured model.
Raises: ValueError: If q_profile is not Q1 or Q2.
apply_anelastic_correction¶
apply_anelastic_correction ¶
Apply anelastic corrections to seismic velocity data using the provided "anelastic_model" within the low attenuation limit. The corrections are based on the equation: \(1 - \frac{V(anelastic)}{V(elastic)} = \frac{1}{2} \cot(\frac{\alpha \pi}{2}) Q^{-1}\) as described by Stixrude & Lithgow-Bertelloni (2005, JGR, doi:10.1029/2004JB002965, Eq. 10).
thermo_model (ThermodynamicModel): The thermodynamic model containing temperature and depth data.
ThermodynamicModel: A new thermodynamic model with anelastically corrected seismic velocities.
The returned model includes the following methods with anelastic corrections:
- compute_swave_speed: Calculates the anelastic effect on shear wave speed.
- compute_pwave_speed: Calculates the anelastic effect on compressional wave speed.
The compute_swave_speed method applies the anelastic correction to the shear wave speed using the provided
anelastic model. The compute_pwave_speed method applies the anelastic correction to the compressional wave speed
using the quality factor derived from the equations provided by Don L. Anderson & R. S. Hart (1978, JGR, eqs 1-3).
The corrections are applied by meshing the depths and temperatures from the thermodynamic model and computing the quality factor matrices for shear and bulk moduli. The corrected seismic velocities are then calculated and returned as new tables with the corrected values.
Utility Functions¶
Coordinate Transformations¶
geodetic_to_cartesian ¶
Convert geographic coordinates to Cartesian coordinates.
Parameters: lat (float or numpy.ndarray): Latitude in degrees. lon (float or numpy.ndarray): Longitude in degrees. depth (float or numpy.ndarray): Depth below Earth's surface in km. earth_radius (float): Radius of the Earth in km. Default is 6371 km.
Returns: tuple: Cartesian coordinates (x, y, z).
cartesian_to_geodetic ¶
Convert Cartesian coordinates to geographic coordinates.
Parameters: x (float or numpy.ndarray): x coordinate in km. y (float or numpy.ndarray): y coordinate in km. z (float or numpy.ndarray): z coordinate in km. earth_radius (float): Radius of the Earth in km. Default is 6371e3 m.
Returns: tuple: Geographic coordinates (lat, lon, depth).
nondimensionalise_coords ¶
Convert dimensional Cartesian coordinates to nondimensional form.
Applies linear radial scaling to map physical coordinates (meters) to a nondimensional reference frame. Commonly used in geodynamic simulations to improve numerical conditioning.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
float or array_like
|
Cartesian coordinates in meters (origin at Earth's center). |
required |
y
|
float or array_like
|
Cartesian coordinates in meters (origin at Earth's center). |
required |
z
|
float or array_like
|
Cartesian coordinates in meters (origin at Earth's center). |
required |
R_nd_earth
|
float
|
Nondimensional radius for Earth's surface. Default is 2.22. |
2.22
|
R_nd_cmb
|
float
|
Nondimensional radius for core-mantle boundary. Default is 1.22. |
1.22
|
Returns:
| Type | Description |
|---|---|
tuple of (float or ndarray)
|
Nondimensionalized (x', y', z') coordinates. |
Notes
Uses linear scaling: r' = a*r + b, where a and b are determined by mapping R_earth → R_nd_earth and R_cmb → R_nd_cmb. Angular coordinates (theta, phi) are preserved.
See Also
dimensionalise_coords : Inverse transformation
dimensionalise_coords ¶
Convert nondimensional Cartesian coordinates back to meters.
Inverse of nondimensionalise_coords. Maps nondimensional coordinates
back to physical units (meters) using linear radial scaling.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
float or array_like
|
Nondimensional Cartesian coordinates. |
required |
y
|
float or array_like
|
Nondimensional Cartesian coordinates. |
required |
z
|
float or array_like
|
Nondimensional Cartesian coordinates. |
required |
R_nd_cmb
|
float
|
Nondimensional radius for core-mantle boundary. Must match the value used in nondimensionalisation. Default is 1.22. |
1.22
|
R_nd_earth
|
float
|
Nondimensional radius for Earth's surface. Must match the value used in nondimensionalisation. Default is 2.22. |
2.22
|
Returns:
| Type | Description |
|---|---|
tuple of (float or ndarray)
|
Dimensional (x, y, z) coordinates in meters. |
Notes
Uses inverse linear scaling: r = a*r' + b, where a and b are determined by mapping R_nd_earth → R_earth and R_nd_cmb → R_cmb.
See Also
nondimensionalise_coords : Forward transformation
Gravity and Pressure¶
compute_gravity ¶
Compute gravitational acceleration at each radius based on the enclosed mass.
Args: radius (numpy.ndarray): Array of radii from the center. mass_enclosed (numpy.ndarray): Array of cumulative mass enclosed up to each radius.
Returns: numpy.ndarray: Array of gravitational acceleration at each radius.
compute_pressure ¶
Calculate the hydrostatic pressure at each radius based on the density and gravitational acceleration.
Args: radius (numpy.ndarray): Array of radii from the center to the surface. density (numpy.ndarray): Array of densities at each radius. gravity (numpy.ndarray): Array of gravitational accelerations at each radius.
Returns: numpy.ndarray: Array of pressures calculated from the surface inward to each radius.
compute_mass ¶
Compute the mass enclosed within each radius using the cumulative trapezoidal rule.
Args: radius (numpy.ndarray): Array of radii from the center of the Earth or other celestial body. density (numpy.ndarray): Array of densities corresponding to each radius.
Returns: numpy.ndarray: Array of cumulative mass enclosed up to each radius.
Constants¶
Usage Examples¶
See the Examples section for Jupyter notebooks demonstrating these APIs in action.