DOI: 10.5281/zenodo.8407693
title: Dataset and Software for the Radial Reference Earth Model REM1D

description: What are the Earth's average physical properties (e.g. density, shear/bulk attenuation, velocity, anisotropy)?

Can astronomic-geodetic and full-spectrum seismological data (~1-3200 seconds) be reconciled?

How can strong lateral variations (e.g. crust) be accounted for (i.e. theory, observations) in radial models?

What can we robustly interpret about the Earth's interior based on its radial features?



Is the outer core well mixed and undergoing adiabatic compression?

Can spin transitions in iron-bearing minerals be detected in the lower mantle?

What is the bulk composition of different regions? Is it pyrolitic everywhere?

Is there a pervasive thermo-chemical boundary layer atop the core-mantle boundary?

What is the radial extent and origin of large-scale anisotropy in the mantle?

Where does the most energy dissipation occur and what is its dominant mechanism?




The spherically-averaged profiles of elasticity, density, and attenuation reflect the bulk composition, temperature profile, and dominant processes of the Earth's heterogeneous interior. The two exhaustive manuscripts below describe the new radial (one-dimensional) reference Earth model (REM1D), which serves as the update to the preliminary reference Earth model (PREM; Dziewonski and Anderson, Phys. Earth Planet. Inter., 1981). The first manuscript describes new modeling concepts and reference datasets while features of REM1D and geological interpretations are discussed in the second manuscript. We extend the Full Spectrum Tomography (FST) technique, which is necessary for constructing reference models that account for lateral heterogeneity and geographic bias in diverse data. 

Summary Figure: REM1D_summary.png

Project Website: rem3d.org

Feedback/Questions? Please get in touch with Raj Moulik rajmoulik.com at moulik@caa.columbia.edu 

Reference:



Please cite the following works if you use this data or software.



Moulik P. & G. Ekström (2025) Radial Structure of the Earth: (I) Model Concepts and Data, Phys. Earth Planet. Inter. doi:10.1016/j.pepi.2025.107319

Moulik P. & G. Ekström (2025) Radial Structure of the Earth: (II) Model Features and Interpretations. Phys. Earth Planet. Inter. doi:10.1016/j.pepi.2025.107320



You can also cite the dataset and software from this Zenodo page (Optional).



Moulik, P. (2025) Dataset and Software for the Radial Reference Earth Model REM1D. In Phys. Earth Planet. Inter. (v1.0). Zenodo. doi: 10.5281/zenodo.8407693




Highlights:



Features and Interpretations (Graphic_abstract_I.png)



Radial reference Earth model (REM1D) is introduced for average physical properties

Peak anisotropy (aS=3.90%) & attenuation (Qmu~60-80,Qkappa~386) in the upper mantle

Consistent with pyrolitic mantle (<800 km) and basal thermo-chemical boundary layer

Lower mantle (mu/kappa, Poisson's ratio) consistent with spin transitions in Fe-bearing minerals

Outer core is neutrally stable, well mixed and undergoing adiabatic compression






Concepts and Data (Graphic_abstract_II.png)



New framework for radial reference models as the spherical average of heterogeneity

Reference data of astronomic-geodetic constants, body/surface waves & normal modes

Theoretical limitations, geographic biases & non-linear crustal effects are quantified

Parameterization is adjusted to satisfy datasets & expectations from mineral physics

Rapid convergence is achieved using analytical kernels & recent a priori constraints




Summary Abstracts:



Features and Interpretations: A new reference model is presented for the spherically-averaged profiles of elasticity, density and attenuation, which reflect the bulk composition, temperature profile and dominant processes of the Earth's heterogeneous interior. This study discusses the features of REM1D and geological interpretations while the underlying modeling concept and reference datasets are described in a companion manuscript. All physical parameters in REM1D vary smoothly between the Mohorovicic and 410-km discontinuity, thereby excluding the 220-km discontinuity in earlier models. REM1D predicts arrival times of major body-wave phases in agreement (±0.8 s, normalized misfit ψpb ≤ 0.25 s) with widely used but theoretically incomplete isotropic models optimized for earthquake location. Substantial radial anisotropy is present only in the shallowest mantle (∼250 km) with peak values of shear-wave (aS = 3.90% , ξ = 1.08) and compressional-wave anisotropy (aP = 3.78% , φ = 0.93) between ∼125–150 km, consistent with textures that can form by the alignment of intrinsically anisotropic minerals in this deforming region. The upper mantle (24.4-410 km) is the most dissipative region with a finite bulk attenuation (Qκ ~ 386) and strong shear attenuation (Qµ ~ 60–80) that peaks at a depth of ~150–175 km in the mechanically weak asthenosphere. An olivine-rich pyrolitic composition is broadly consistent with REM1D structure in the upper mantle and extended transition zone (<800 km) with step changes across the 410-km and 650-km discontinuities. Features of the lower mantle can be reconciled with: (i) effects of thermally driven convection throughout the central lower mantle (771–2741 km) leading to an apparent subadiabaticity in the stratification parameter ηB, (ii) effects of spin transitions in iron-bearing minerals that manifest as distinct linear segments in modulus and Poisson’s ratios (µ/κ, σP) on either side of a complex transition region (~1300–1700 km, 52–73 GPa), (iii) a thermal boundary layer with steeper superadiabatic gradients than near the surface, which ultimately exceed the critical gradients for both vP and vS (but not for density ρ) at a depth of 2741 km, and (iv) chemical stratification in the bottom ~500–750 km of the mantle that acts to suppress the thermal effects. Signatures of this thermo-chemical boundary layer are: (i) a gradual increase of density and steep positive gradients with depth (dρ/dz) in the lower mantle, (ii) large values of the stratification parameter (ηB > 1.03) followed by an abrupt reduction to values below one near the core-mantle boundary (CMB), (iii) variations in bulk modulus with pressure κl = dκ/dp that are inconsistent with Equations of state (EoS) expectations of a uniform composition, (iv) very steep negative vP and vS gradients that form a low-velocity zone in the Dll region. The vP and ρ variations in the outer core have steep and the derivative properties are consistent with a neutrally stable region comprising a well-mixed iron alloy undergoing adiabatic compression (ηB ~1, N 2~0, negative κll). REM1D is readily extendable due to its modular construction and represents the average physical properties, features essential for geological interpretations and the construction of a three-dimensional reference Earth model.

Concepts and Data: A framework is introduced for developing a radial reference model that incorporates diverse observations and techniques for improving the constraints on bulk Earth structure. This study describes new modeling concepts and reference datasets while features of the reference Earth model REM1D and geological interpretations are discussed in a companion manuscript. Recent measurements from various techniques have improved in precision and are broadly consistent, and are summarized as best estimates with uncertainties. We construct a reference dataset comprising normal-mode eigenfrequencies and quality factors, surface-wave dispersion curves, impedance constraints and travel-time curves from body waves, and astronomic-geodetic observations. Classical radial reference models do not account for the theoretical effects and observational biases resulting from heterogeneity in the crust and mantle. We address three issues that account for lateral variations in the modeling of average elastic, anelastic and density structure. First, current ray coverage of traveling waves is biased towards structure in the northern hemisphere, leading to faster velocities especially in the lower mantle. Second, horizontal wavelength of the heterogeneity that a traveling wave encounters is assumed to be much greater than that of the corresponding normal mode in most ray-theoretical and finite-frequency formulations of wave propagation. Effects of the full volumetric sensitivity on local eigenfrequencies and phase velocities that are ignored with this approximation exceed the data uncertainty for both fundamental spheroidal (Rayleigh waves, T ≥ 220 s) and toroidal modes (Love waves, T ≥ 120 s); waves at these longer periods cannot be modeled solely in terms of radial variations along the ray path. Third, non-linear effects from the strongly heterogeneous crustal structure are substantial for shorter-period waves (T ≤ 100 s) and need to be accounted for while deriving radial models. After accounting for these issues on heterogeneity, rapid convergence for average structure is facilitated by utilizing a priori constraints from recent literature, analytical sensitivity kernels that account for physical dispersion, and a flexible parameterization comprising polynomial functions and cubic B-splines. By adopting a higher-order polynomial for density than the elastic structure, artifacts that imply strong inhomogeneity and non-adiabaticity are avoided in potentially well-mixed regions like the outer core. Derivative properties like the gradient of bulk modulus with pressure (κl = dκ/dp) and the Bullen’s stratification parameter (ηB) are adjusted in the core to match expectations from mineral physics without deteriorating the fits to reference datasets. A cubic polynomial parameterization in the lower mantle is adequate to capture possible changes in the gradients of the modulus ratio (µ/κ) associated with spin transitions in iron-bearing minerals. Radial reference models need to account for lateral heterogeneity and prior geological information in their construction to accurately represent the bulk average properties of a heterogeneous Earth.


Data Products:



REM1D Model Coefficients



Elastic/density structure (REM1D.elas), anelastic structure (REM1D.anelas), and corresponding card deck file listing properties in 750 concentric layers/shells (REM1D.cards). These files are appropriate for various normal-mode, surface-wave and astronomic-geodetic calculations. A summary highlight figure is also provided for the REM1D model (REM1D_summary.png).

Same as above but modified for calculating body-wave arrivals (REM1D_TT.elas, REM1D_TT.anelas, REM1D_TT.cards). Values in the upper crust are extended to the surface and ocean is removed. These files are appropriate only for the body-wave calculations.

All physical parameters (VPH,VSH,VPV,VSV,VP,VS,Eta,A,C,N,L,F,Eta,Rho,Qmu,Qkappa) and derivative properties (Poisson's ratio, Pressure, dK/dP, Bullen, Gravity) from REM1D evaluated in 750 concentric layers/shells inside the Earth (REM1D_layer_properties.txt). These values are evaluated at the reference period of 1 second from the model coefficients listed in the files above (REM1D.elas, REM1D.anelas) using the equations in Paper II. Elastic properties at other periods after accounting for physical dispersion can be calculated following equation 8 in Paper II and using the Python script below.



Python Script



A pure Python 3 script (evaluate_REM1D.py) to evaluate REM1D physical parameters at any depth and period accounting for physical dispersion. This script needs coefficient files (REM1D.elas, REM1D.anelas) in the same directory. Some examples from the command line:



Help: evaluate_REM1D.py -husage: evaluate_REM1D.py [-h] [-m MODEL] [-T PERIOD] [-p PARAMETER [PARAMETER ...]] [-b BOUNDARY] -d DEPTH_IN_KM                         [DEPTH_IN_KM ...]Evaluates REM1D physical parameters at the reference period of the default 1 second e.g. to calculate all parameters at150 km, enter: evaluate_REM1D.py -d 150options:  -h, --help            show this help message and exit  -m MODEL, --model MODEL                        Radial reference Earth model defined at the reference period (default: 1 second)  -T PERIOD, --period PERIOD                        Reference period (default: 1 second)  -p PARAMETER [PARAMETER ...], --parameter PARAMETER [PARAMETER ...]                        Physical parameter. Can give multiple values e.g. rho vpv vsv qkappa qmu vph vsh eta  -b BOUNDARY, --boundary BOUNDARY                        If the depth queried is an internal boundary, + evaluates the value above while - evaluates the value                        below the discontinuity. Ignored and assumed to be + if depths queried have repeated values.  -d DEPTH_IN_KM [DEPTH_IN_KM ...], --depth_in_km DEPTH_IN_KM [DEPTH_IN_KM ...]                        Depth queried in km. Can give multiple values e.g. 24.4 75. 150. 225. 300. 410. 530. 650. 650.                        820. 1320. 1820. 2320. 2550. 2791. 2891. If depth is repeated, then values above and bottom are                        queried assuming the user has queried a discontinuity.

Calculate all parameters at 150 km: evaluate_REM1D.py -d 150

Calculate density at 100 km: evaluate_REM1D.py -d 150 -p rho

Calculate all parameters at 150 km and 200 seconds accounting for physical dispersion following equation 8 in Paper II: evaluate_REM1D.py -d 150 -T 200

Calculate all parameters at 100, 150, 650, 1500, and 2891 km. Note that the discontinuities are repeated to evaluate above (+) and below (-) the internal boundary; if not repeated, the code will only provide values above (+) by default: evaluate_REM1D.py -d 100 150 650 650 1500 2891 2891

Calculate all parameters explicitly below (-) the 650-km discontinuity: evaluate_REM1D.py -d 650 -b -





Reference Normal Mode/Surface Wave Data



Reference datasets of phase-velocity perturbations w.r.t. PREM predictions between 4-40 mHz



Love waves (avgdcbyc_Love.ref)

Rayleigh Waves (avgdcbyc_Rayleigh.ref)



Reference normal mode datasets 



Radial Modes - Eigenfrequencies (RADIAL_F.ref), Quality Factors (RADIAL_Q.ref)

Spheroidal Fundamental Modes - Eigenfrequencies (SFUND_F.ref), Quality Factors (SFUND_Q.ref)

Spheroidal Overtones - Eigenfrequencies (SOVER_F.ref), Quality Factors (SOVER_Q.ref)

Toroidal Fundamental Modes - Eigenfrequencies (TFUND_F.ref), Quality Factors (TFUND_Q.ref)

Toroidal Overtones - Eigenfrequencies (TOVER_F.ref), Quality Factors (TOVER_Q.ref)





REM1D Data Fits and Predictions



Fits to reference datasets



Radial Modes (FITS_RADIAL.txt)

Spheroidal Fundamental Modes, including Rayleigh waves (FITS_SFUND.txt)

Spheroidal Overtones (FITS_SOVER.txt)

Toroidal Fundamental Modes, including Love waves (FITS_TFUND.txt)

Toroidal Overtones (FITS_TOVER.txt)



Normal mode catalogs for periods down to 24 seconds in binary (REM1D_024.rts) and HDF formats (REM1D_024.h5). These calculations use the new reference astronomic-geodetic constants (e.g. revised Gravitational constant G) from Paper I. These can be read using standard HDF5 modules (e.g. h5py) or the open-source AVNI package. Conversion from the binary file was done using the following commands:



from avni.data import NM

NM.write_modes_hdf('REM1D_024.rts','REM1D_024.h5')



Non-linear crustal corrections derived from REM1D overlain with CRUST2.0 heterogeneity for Love (nonlinear_CRUST2_Love.csv) and Rayleigh waves (nonlinear_CRUST2_Rayl.csv). These files contain domega and dc/c that need to be added to the eigenfrequency and dispersion predictions from REM1D, respectively, following equations 13 and 14 in Paper I.