Manuscript draft in the format of Geophysical Journal International · not submitted

Local-earthquake tomography of the Mount Rainier crust as a perturbation of a geology-based three-dimensional model

rainier3d working draft, 27 October 2026. Stage S37, branch paper-review-fixes, commit given in the repository log for docs/tomography/.

Numbers come from the run of 27 October 2026 on a 250 m eikonal grid with the critical zone left as built, from the smoothed and recalibrated starting model with the lateral cutoff growing from zero at the base of the critical zone (127.1 min on a 48-core Linux server), unless another run is named.

1 Introduction

Moran et al. (1999) used local earthquakes to image the P-wave structure of the crust around Mount Rainier, including a low-velocity body 5–18 km beneath the summit. The rainier3d model takes a different route. It builds the upper crust from mapped geology and a rock-physics law, adds a critical-zone model for the top 150 m, fuses the result with regional tomography at long wavelengths, and calibrates a small number of global parameters against PNSN travel times (stage S13). That calibration adjusts depth-dependent factors and therefore cannot correct lateral errors in the geological rules, such as where map units are extended to depth.

Here we ask how much lateral structure the existing PNSN arrival times can add to that model, and with what resolution. We treat the geology-based model as the reference and invert only for smooth perturbations, so that the result can be read as a correction to a prior rather than as an independent image. We choose the regularisation and the stopping point with earthquakes held out by region, which tests whether an update obtained in one part of the network improves predictions elsewhere.

2 Data

The picks are PNSN analyst P and S arrival times from the ComCat phase-data products of all events of magnitude 0.5 or larger located in the 70 × 75 km model domain between 2015 January 1 and the data freeze of 2026 September 23 (3867 events, 117 609 picks), merged with the 91 validation events of stage S13. Only stations inside the domain are used, because the travel-time grid covers the domain only; this leaves 50 stations. Events need at least six picks, four of them P. The resulting data set has 3770 events, 57 659 P and 42 057 S picks, and a median of 25 picks per event (Table 1, Fig. 1).

The seismicity has two dominant clusters (Fig. 1a, c). About 1720 events lie less than 1 km below sea level beneath the summit, spread over the whole period with 620 of them in 2025, and the western Rainier seismic zone extends from about 2 to 15 km depth some 15 km west of the summit. Events are assigned to a fitting and a held-out set by 10 km squares of their catalogue epicentres. Squares are allocated in decreasing order of pick count to whichever set has fewer picks, which gives 1921 fitting and 1849 held-out events with 49 859 and 49 857 picks. The summit square falls in the fitting set, so 81 per cent of the shallow summit events are fitted and most of the held-out events belong to the western seismic zone and the outer domain. The held-out score therefore measures transfer of the update to regions away from the summit cluster.

3 Method

3.1 Starting model and forward problem

The starting model is the rainier3d model of this branch: rule-based geology with the critical zone in its top 150 m, a tapered spherical low-velocity body centred 11.5 km below sea level, geological contacts and unit boundaries low-passed below the critical zone with a lateral cutoff that grows from zero at its base to 750 m 150 m deeper and then with depth, fusion with the USGS Cascadia community velocity model and CRESCENT Gen0 at long wavelengths, the two blended over 8–16 km below the ground, and the S13 calibration redone on that model (run of 2026 October 27). It is resampled on a 250 m travel-time grid (281 × 301 × 139 nodes, 11.8 million, from 4.5 km above to 30 km below sea level). Within the depth range of the top model level, whose cells are 50 m thick, each grid node takes the average slowness of the cells within half a grid spacing of it, so that thin slow layers such as the critical zone keep their vertical travel time; nearest-node sampling would keep or miss them depending on where a node falls. At the first node at or below the ground this lowers the median Vs from 2608 to 2220 m s⁻¹. Grid cells more than one cell above the ground carry 343 m s⁻¹, so rays cannot cross valleys through air at rock speed. Travel-time fields are computed for every station and phase with the fast-marching solver of White et al. (2020), using reciprocity; one field takes about 20 s.

3.2 Location

Every event is relocated in every trial model by a grid search below the ground on every second node (500 m), followed by a Huber-weighted Gauss–Newton refinement on travel times interpolated on the full grid (Huber 1964; the grid-search-then-refine scheme follows Lomax et al. 2000). Pick uncertainties are 0.14 s for P and 0.23 s for S, the root-mean-square residuals of the PNSN locations of the S13 events; picks beyond five standard deviations receive zero weight. Hypocentres may lie at most 50 m above the ground.

3.3 Parameterisation and Fréchet derivatives

The model is

V(x) = V₀(x) exp[w(d) m(x)],

where d is the depth below the ground or the glacier bed and w is 0 in the top 150 m, the critical-zone columns of the starting model, and rises linearly to 1 at 500 m. The forward model, the kernels and the written model use the same w, so the critical zone keeps its rock physics and joins the updated rock below without a step. The perturbation m is interpolated trilinearly between nodes spaced 3 km horizontally (24 × 26) and placed at nine depths, −3, 0, 2, 4, 6, 9, 13, 18 and 25 km below sea level, separately for P and S (5616 nodes per phase). The perturbation applies to rock only. The derivative of a travel time with respect to node k is −∫ s w w_k dl along the ray, where s is the slowness and w_k the trilinear weight; rays are traced down the gradient of the station's travel-time field in steps of a quarter grid spacing. For the best-sampled S node, the ray derivatives of 18 253 picks agree with finite differences through the eikonal solver with a slope of 0.9995 and a correlation of 0.9993 (Fig. 2a).

3.4 Separation of hypocentral parameters

For each event, the rows of the system are scaled by the square root of the location weight divided by the pick uncertainty, and projected onto the orthogonal complement of the partial derivatives with respect to the hypocentre and the origin time (Pavlis & Booker 1980). The projected system relates the residuals that relocation cannot remove to the velocity perturbation.

3.5 Regularisation and held-out choice of the smoothing

Each step minimises the projected misfit of the fitting events plus λ_s ‖Dm‖² + λ_d ‖m‖², where D stacks second differences along the two horizontal axes and depth (with the coefficients of uneven spacing in depth), both weights are relative to the mean squared column norm of the projected matrix, and λ_d = 0.03. The system is solved with LSQR (Paige & Saunders 1982), in the manner of Thurber (1983). At the first iteration the linearised held-out misfit is computed for λ_s between 0.1 and 100 (Fig. 2b). It decreases monotonically towards small weights and flattens below 0.3 (0.4885 at 0.3, 0.4860 at 0.1), where the damping rather than the smoothing controls the solution. We take the largest weight whose score is within 1 per cent of the best, λ_s = 0.3, as the smoothest model that the held-out events do not penalise.

3.6 Stopping rule and safeguards

After each step all events are relocated in the new model and the held-out misfit, the mean of (residual/σ)² over P and S picks, is recomputed. A step is kept if it lowers this misfit by more than 0.5 per cent. Otherwise half the step is tried, and if that also fails the last kept model is final. The first three steps lowered the held-out misfit from 0.502 to 0.282, 0.272 and 0.2701. The fourth step raised it to 0.2707, by 0.2 per cent, and its half lowered it to 0.2700, by less than 0.1 per cent, so the inversion stopped with the model of the third step (Fig. 2c). Every trial model is saved before its travel times are computed and checked for non-finite or non-positive velocities, and its travel-time fields must be computed within three times the time taken for the starting model, otherwise the trial counts as rejected. A first run with λ_s = 0.1 stalled in the eikonal solves of its second iteration for 1.5 hours before these safeguards were added; we have not identified the cause.

4 Ray coverage

We trace the ray of every pick in the final model and record, in each 250 m cell, the ray length and the number of distinct rays. P rays pass through 1.06 million of the 11.8 million cells and S rays through 0.89 million; 116 614 cells are crossed by ten or more P rays, and the densest cell, 11 km south and 3.5 km west of the summit at 0.75 km elevation, by 3590. The total ray length is 1.11 million km for P and 0.82 million km for S. Aggregated in 1 km cubes (Fig. 3), coverage is densest beneath the summit in the top 2 km and along the western seismic zone between 2 and 12 km depth. It thins rapidly below 12 km and is nearly absent below 15 km. Directly beneath the summit, few rays pass below about 5 km depth, because the summit earthquakes are shallow and the deeper earthquakes lie to the west (Fig. 3b). The sampled rays (Figs 4 and 5) show the two families: short, steep paths from the summit cluster to the stations on the edifice, and longer paths from the western seismic zone that cross the volume between 2 and 12 km depth at moderate incidence.

5 Resolution

We test resolution with checkerboards of ±5 per cent and cell sizes of 12 and 6 km in all three directions, with P and S of opposite sign so that leakage between phases is visible (Figs 6–8). Synthetic data are computed with the projected kernels of the final model, Gaussian noise at the pick uncertainty is added, and the patterns are inverted with the same weights as the data. This is a linear test around the final model and does not include changes of ray paths or hypocentres. Checkerboard tests describe the recovery of one pattern and can overstate resolution of other structures (Lévêque et al. 1993), so we read them as an upper bound.

The 12 km pattern is recovered with correlations of 0.66–0.74 between sea level and 6 km depth, and with 45–54 per cent of its amplitude (Table 2). Recovery falls to 0.49–0.56 at 9 km and to 0.29–0.44 at 13 km. The node layer 3 km above sea level recovers nothing (|r| ≤ 0.09), and the layer at 25 km recovers the pattern with the wrong sign; both are controlled by the smoothing operator. The 6 km pattern is recovered with correlations of 0.37–0.43 and 18–20 per cent of its amplitude in the top 2 km, where the summit earthquakes and the dense stations on the edifice provide crossing rays, and with 0.21–0.23 between 4 and 9 km. The recovered 6 km fields are recognisable near the summit and noisy elsewhere (Fig. 7).

6 Results

6.1 Misfit

The three kept steps reduce the residuals of both sets (Table 3, Fig. 9). For held-out events the RMS falls by 20 per cent for P and 33 per cent for S, and the median absolute residual by 34 and 42 per cent. The reduction for fitting events is larger, 38 and 46 per cent; the difference between the two sets is the part of the update that does not transfer across regions. The RMS of held-out residuals averaged per station (stations with at least 20 picks) falls from 0.052 to 0.023 s for P and from 0.152 to 0.066 s for S, and the largest S station mean from 0.39 to 0.18 s (Fig. 10). Since the critical zone is not changed, these delays are explained by the rock between 150 m and a few kilometres below the stations.

6.2 Velocity perturbations

Nodes with a derivative weight sum above 5 km make up 52 per cent of the P nodes and 46 per cent of the S nodes (Fig. 11, bottom row). Within them the perturbation has a mean of +0.001 and an RMS of 0.047 for ln Vp, and a mean of −0.005 and an RMS of 0.057 for ln Vs; the 5th to 95th percentiles are −0.082 to +0.066 for P and −0.093 to +0.061 for S. Between sea level and 13 km the mean ranges from −0.012 at sea level to +0.008 at 6–9 km for P, and from −0.029 at sea level to +0.015 at 13 km for S.

Between sea level and 6 km depth, velocities increase by 5–8 per cent in an arc north-east, east and south of the summit and decrease by similar amounts to the west and south-west, around and above the western Rainier seismic zone (Fig. 11). Within 4.5 km of the summit, the node layer 3 km above sea level, inside the edifice, is slower by 6.3 per cent in Vp and 4.4 per cent in Vs on average, and the layer at 2 km depth is faster by 3.9 per cent in Vp and 1.9 per cent in Vs. The 12 km checkerboard is not recovered at the layer above sea level (Table 2), so the edifice anomaly cannot be interpreted with the present parameterisation. The vertical sections (Fig. 12) show a fast upper crust immediately east of the western seismic zone and a slow column beneath the summit from about 5 to 20 km depth, which extends through the low-velocity body of the starting model. Ray coverage beneath the summit is nearly absent below 5 km (Fig. 3b), so this column is the downward continuation of the edifice anomaly by the smoothing operator and carries no information from the data. Outside the sampled nodes the perturbation reaches 0.25 for P and 0.38 for S, mostly in the corners of the domain; the second-difference regulariser does not penalise linear gradients and extends them into unsampled volume.

The critical zone is unchanged by construction. In the L1 cells that average the critical-zone columns the perturbation is zero, except in the cells that straddle the 150 m base, where the taper begins and the change is at most 2.3 per cent in Vp and 3.0 per cent in Vs. The invariants of the model hold in model_let.zarr: the vertical S travel time of the top two rock cells equals that of the critical-zone columns to 10⁻³ in all 81 388 ice-free columns, Vp/Vs stays at most 4.70 in cells that contain critical zone and between 1.57 and 2.16 in rock below, and the fused Vp/Vs beneath the tomography factor lies between its geology and regional inputs. An otherwise identical run without the taper changed Vs in the critical-zone cells by a median of −2 per cent and by more than 5 per cent in 30 per cent of them, and broke the travel-time agreement with the columns in 79 per cent of the columns; it fitted the held-out data slightly better (Table 4), by absorbing near-station delays into a layer that the data do not resolve.

6.3 Hypocentres

Relocation in the final model moves epicentres by a median of 462 m and hypocentres shallower by a median of 654 m (Fig. 13). The depth change is larger beneath the summit (median −766 m for 1665 events within 5 km of it) than elsewhere (−476 m), and 136 events move by more than 3 km in depth. The median formal depth uncertainty falls from 422 to 266 m. Since the mean velocity change is close to zero, the upward shift beneath the summit is unlikely to be a uniform trade-off between velocity and depth; it may be related to the slower edifice above sea level, which lies at a node layer that the checkerboard test does not resolve. The present station geometry does not constrain the absolute depth of the summit cluster independently.

6.4 Comparison of runs

Four earlier runs used the same data and splits (Table 4). The first used the 500 m grid of the S13 calibration, nearest-node sampling of the starting model and the smoothing weight of its own scan, λ_s = 10, and no critical-zone taper. It kept one step, reduced held-out residuals less, resolved only 12 km structure, and produced a positive mean perturbation of +2.2 per cent in Vp and +1.6 per cent in Vs, coupled to a larger upward shift of hypocentres. The second used the 250 m grid with slowness averaging but no taper. Between these two the grid spacing, the sampling of the starting model and the smoothing weight all change, so the table does not attribute the differences to any one of them; the stronger resolution on the finer grid follows in large part from the weaker smoothing, which the held-out events chose, and the finer grid made the kernels more accurate (slope 0.998 against 1.02). The taper alone was tested on the starting model of the second run (run of 13 October 2026): leaving the critical zone as built raised the held-out RMS by 0.0014 s for P and 0.0053 s for S, kept the resolution, and made the result consistent with the petrophysical model of the near surface. The third run (14 October 2026) added the taper on a starting model recalibrated after the geology of the fusion was smoothed. The fourth run (21 October 2026) started from a smoother model: the two regional models are blended over 8–16 km below the ground instead of switching at 9.9 km, the lateral smoothing of the geology starts at a 750 m cutoff just below the critical zone, and the calibration was redone on that model. With the same data, splits and parameters it kept three steps instead of two, and its residuals, resolution and perturbation pattern were close to those of the third run. The run of this report (27 October 2026) starts from the model in which the lateral cutoff grows from zero at the base of the critical zone, recalibrated again; it kept three steps, with residuals and resolution close to those of the fourth run.

7 Discussion

The held-out reduction of 20 per cent for P and 33 per cent for S shows that the arrival times contain lateral information that the geology-based model, its fusion with regional tomography and its depth-dependent calibration do not capture. The update transfers across regions, since every held-out event lies in a square that contains no fitting event. It does not transfer completely: the fitting events gain more, by 38 and 46 per cent.

The station-mean residuals fall by more than half without station terms and without any change to the critical zone. The top 150 m beneath each station are carried by the petrophysical critical-zone model and the topography, and the delays that remain after the starting model are explained by the rock between 150 m and a few kilometres. A run that also perturbed the critical zone fitted the held-out data only slightly better while breaking its consistency with the rock physics, which indicates that the critical-zone model is not the main source of the station delays. Where it is wrong, the remedy is to calibrate its parameters (weathering depth, water table) against observations that resolve the top 150 m, not to scale its velocities in a tomography that does not.

Resolution is set by the station geometry, the depth distribution of the earthquakes and the regularisation. The 12 km checkerboard is recovered with half of its amplitude between sea level and 6 km, and the 6 km checkerboard only in the top 2 km near the summit. Below 12 km, and beneath the summit below 5 km, rays are too few for the perturbations to be interpreted, which includes the depth range of the low-velocity body of the starting model. The weak smoothing chosen on the held-out events also makes the update rougher than on the 500 m grid, and the roughness at scales below about 6 km is not resolved.

Before the update replaces the model product, we consider five steps necessary. The travel-time grid should be extended beyond the model domain so that stations outside it contribute; of the 6876 picks of the S13 events, 3108 remain after restriction to stations and events inside the domain. Damping should increase where the derivative weight sum is small, or the update should be applied only there, so that unsampled volume returns to the starting model. A test with structures of realistic shape should be added, together with a nonlinear version of the resolution test that includes relocation. The two grids should be compared at the same smoothing weight, to separate the effect of the grid from that of the regularisation. More S arrivals, for example from automatic picking of the 2025 nodal deployment, would reduce the larger S uncertainties.

8 Conclusions

A regularised local-earthquake tomography on a 250 m eikonal grid, with relocation in every model and region-based held-out validation, reduces held-out P and S residuals at Mount Rainier by 20 and 33 per cent relative to a geology-based, calibrated three-dimensional model, and halves the mean residual at stations without station terms, while leaving the petrophysical critical-zone model of the top 150 m unchanged. The data resolve structure of about 12 km between sea level and 9 km depth, and about 6 km in the top 2 km near the summit. Within the sampled volume the update is a pattern of 5–8 per cent contrasts with a mean close to zero; below 12 km and beneath the summit below 5 km it is not constrained by rays.

Data and code availability

Picks are from the USGS ComCat phase-data products for PNSN (network code uw), queried up to the data freeze of 2026 September 23 and cached in data/raw/pnsn/catalog_2015-01-01_2026-09-23_m0.5. The inversion is scripts/37_tomography.py with rainier3d.validate.tomo; the figures are drawn by scripts/37_tomography_figures.py from outputs/tomography/, which also holds the ray density (ray_density.nc), the sampled ray paths (rays_sample.json) and every trial model. All parameters are in the tomography block of configs/domain.yaml; the node spacing, depths, damping, block size, stopping tolerance and grid spacing are placeholders chosen by the authors. The result is data/processed/model_let.zarr; the 500 m run is kept in outputs/tomography_3770ev_500m/ and data/processed/model_let_500m_20261013.zarr, and the 250 m run without the critical-zone taper in outputs/tomography_250m_notaper/ and data/processed/model_let_250m_notaper.zarr. How the whole model is built is recorded in docs/model_construction.md.

References