5.1.1 Rock under confining pressure

Each rock unit receives a P-wave speed that increases with effective pressure as cracks close: $$V_P(P) = V_\infty - (V_\infty - V_0)\, e^{-P/P^{*}}, \qquad P = \rho_{b}\, g\, d - \rho_{w}\, g\, \max(d - z_w,\ 0),\qquad{(1)}$$ with bulk density $\rho_b$ = 2500 kg m⁻³, water density $\rho_w$ = 1000 kg m⁻³, $d$ the depth below the local ground and $z_w$ the depth of the water table of Sect. 5.1.2 (hydrostatic pore pressure below it, none above).

• Vp/Vs and density. Ice, open water and the three unconsolidated deposits carry their own Vp/Vs and density in configs/petrophysics.csv, partly because the Brocher relations are not valid below Vp = 1.5 km s⁻¹. Where a unit has no ratio or density of its own, both follow Brocher (2005): Vs from his regression fit (his Eq. 6), $V_S = 0.7858 - 1.2344V_P + 0.7949V_P^2 - 0.1238V_P^3 + 0.0064V_P^4$, and density from the Nafe–Drake fit (his Eq. 1), $\rho = 1.6612V_P - 0.4721V_P^2 + 0.0671V_P^3 - 0.0043V_P^4 + 0.000106V_P^5$, with V in km s⁻¹ and ρ in g cm⁻³ (valid for 1.5 < Vp < 8.5 km s⁻¹; S4, rainier3d.petro.relations, rainier3d.properties.assign).
• Attenuation. $Q_S = 0.05\,V_S$ (Vs in m s⁻¹) and $Q_P = 2\,Q_S$ (q block of configs/perturbations.yaml, m1_placeholder); S5 recomputes both from the fused Vs. Because $Q_P = 2\,Q_S$ implies a negative bulk quality factor where Vp/Vs < √(8/3) ≈ 1.63, the 19,476 of 2.24 million non-ice fused cells below that ratio carry an inconsistent Q. The rule has the form proposed by Olsen (2003) for the Los Angeles basin.
• Parameters. Table 9 lists the unit parameters, chosen within published ranges. Every row of configs/petrophysics.csv has the source key m1_placeholder and a one-line basis (for example, granodiorite and quartz diorite for the plutons). The calibration of Sect. 5.3.4 scales them with three global multipliers that act on the eleven rock units, the magma body included, and leave ice and the unconsolidated deposits unchanged (block geology of configs/velocity_calibration.yaml). The contrasts between units therefore come from the table, and their overall level comes from the travel times. [TODO: the published range each table value was taken from, per unit]
• Magma body. Inside it, the middle-crust parameters apply (V₀ = 6.0 km s⁻¹, V∞ = 6.5 km s⁻¹, P* = 50 MPa; row magma_mush of configs/petrophysics.csv). Its anomaly is not a property of the rock law: Vp, Vs and density of the fused model (Sect. 5.3.3) are reduced by up to 10%, 15% and 3% (magma_mush block of configs/perturbations.yaml, m1_placeholder), which raises Vp/Vs by about 6%, with a weight of 1 in the core of the body and a cosine taper to 0 at its edge (Sect. 5.3.1). It is applied after the fusion: applied before it, as in earlier versions of the model, the anomaly was largely removed by the regional low-pass and kept only at the edges of the body.
• Alteration. The alteration intensity $a$ (0 to 1, from S22 in Sect. 5.3.2; zero in air and ice) scales Vp by $(1 - 0.30a)$, Vs by $(1 - 0.35a)$ and density by $(1 - 0.12a)$, within the ranges reported for altered volcanic rock (Heap and Violay, 2021). Fully altered rock is thus 30% slower in Vp, 35% slower in Vs and 12% less dense, with Vp/Vs about 8% higher. S4 applies the factors after the crack-closure law (alteration block of configs/perturbations.yaml, m1_placeholder; the file notes the values are still to be sourced).

5.1.2 Granular and partly saturated media

Soil, Quaternary cover and weathered rock are described by a granular frame whose stiffness depends on effective stress and water content. It is weighted against the rock law by the weathering index $W$ of the critical-zone columns (Sect. 5.2).

Water and effective stress. Pressure head $h$ is hydrostatic below the water table of Fan et al. (2017), which is placed at the glacier bed under ice. (The estimate of Ma et al. fits the wells better, Sect. 7, but its licence forbids derivatives.) Above the water table the suction head $\psi = -h$ sets the effective saturation by the retention curve of Genuchten (1980), $$S_e = \left[1 + (\alpha\psi)^n\right]^{-(1 - 1/n)},\qquad{(2)}$$ with $\alpha$ and $n$ from POLARIS in soil and placeholder classes below. The effective stress is $$\sigma' = \sigma_v - u, \qquad u = \begin{cases} \rho_w g h, & h > 0,\\ -\rho_w g \psi S_e\, W, & h \le 0, \end{cases}\qquad{(3)}$$ where $\sigma_v$ is the overburden of the column, including the weight of glacier ice, and the second line is the suction stress of Lu et al. (2010) with $\chi = S_e$ (Lu and Likos, 2006), which stiffens the unsaturated granular frame. The crack-closure pressure $P$ of Eq. 1 in every L1 to L3 cell uses the same water table.

Granular frame. The dry frame of soil, cover and weathered rock follows the soft-sand model of Dvorkin and Nur (1996): at the critical porosity $\phi_c$ = 0.40 a random pack of spheres has the Hertz–Mindlin moduli (Mindlin, 1949; Walton, 1987) $$K_{\mathrm{HM}} = \left[\frac{C^2(1-\phi_c)^2\mu_s^2}{18\pi^2(1-\nu_s)^2}\,\sigma'\right]^{1/3}, \qquad \mu_{\mathrm{HM}} = \frac{2 + 3f_s - \nu_s(1 + 3f_s)}{5(2 - \nu_s)} \left[\frac{3C^2(1-\phi_c)^2\mu_s^2}{2\pi^2(1-\nu_s)^2}\,\sigma'\right]^{1/3},\qquad{(4)}$$ with coordination number $C$ = 9 and slip fraction $f_s$ = 0.5, which reduces the shear stiffness that the pure Hertz–Mindlin model overpredicts (Makse et al., 1999). Below $\phi_c$ the frame joins the mineral along the modified Hashin–Shtrikman lower bound (Dvorkin and Nur, 1996; Mavko et al., 2020). The grain moduli $K_s$, $\mu_s$ and Poisson’s ratio $\nu_s$ are the Hill average (Hill, 1952) of clay (21 and 7 GPa) and quartz–feldspar (40 and 25 GPa) for the SOLUS100 clay fraction. The pore fluid has the patchy-mixing modulus of Brie et al. (1995), $K_f = (K_w - K_a)S_w^{b} + K_a$, with water saturation $S_w$ and $b$ = 3, and the saturated bulk modulus follows Gassmann’s equation (Mavko et al., 2020), $$K_g = K_d + \frac{(1 - K_d/K_s)^2}{\phi_g/K_f + (1-\phi_g)/K_s - K_d/K_s^2}, \qquad \mu_g = \mu_d,\qquad{(5)}$$ valid at seismic frequencies, which lie in the low-frequency limit of Biot’s theory (Biot, 1956). All rock-physics parameters are m1_placeholder values within the ranges of Mavko et al. (2020).

One medium from soil to rock. At each depth the rock end member is the crack-closure law of the unit beneath, evaluated at the same effective stress $\sigma'$. The bulk and shear moduli of the medium are the Hill average of the granular and rock moduli, weighted by $W$, $$M = \tfrac12\left[W M_g + (1-W)M_r\right] + \tfrac12\left[\frac{W}{M_g} + \frac{1-W}{M_r}\right]^{-1}, \qquad \rho = W\rho_g + (1-W)\rho_r,\qquad{(6)}$$ for $M = K, \mu$, so that velocity varies continuously with effective stress and saturation from soil to fresh rock, and where $W$ = 0 the medium is exactly the rock law of Sect. 5.1.1. Permeability blends in the same way, from the conductivity of POLARIS and the cover and weathered classes to the crustal permeability–depth curve of Ingebritsen and Manning (1999) in rock, capped at 10⁻¹² m².

5.1.3 Shared inputs of the two laws

The two laws are built to agree where they meet, and the parts of the model share their inputs:

• One pore pressure. Both laws use the same water table, placed at the glacier bed under ice, so the crack-closure pressure $P$ of Eq. 1 in every L1–L3 cell and the effective stress $\sigma'$ of Eq. 3 in the columns share one pore pressure. The overburden is not yet shared: the columns integrate their own densities, ice at 917 kg m⁻³ included, while Eq. 1 uses a constant $\rho_b$ = 2500 kg m⁻³ from the ground or ice surface, so ice is loaded as rock in the crustal cells.
• One rock law. Where $W$ = 0 the medium is the rock law of Sect. 5.1.1, unchanged. The calibration multipliers of Sect. 5.3.4 act on that law, so they change the crustal cells and the rock end member of every column together.
• Continuity. At the base of the columns (150 m), Vs differs by a median 0.5% (90th percentile 1.5%) from the L1 cell centred 25 m below (cz.base_continuity in outputs/paper/numbers.json). The lateral low-pass of the geology starts from zero at the column base (Sect. 5.3.3), so the fusion does not open a step there.
• What is not yet consistent. The overburden of the crustal cells (above). The quality factors follow a rule of their own (Sect. 5.1.1). Cell averages over a saturated weathered layer reach Vp/Vs above 3 in the top L1 cell (Sect. 5.2). The water table is that of Fan et al. (2017), not a solution of a groundwater model on these columns.

5.3.1 Three-dimensional geology

Units are extended to depth by explicit column rules, in the manner of the San Francisco Bay region model (Aagaard and Hirakawa, 2021): the geology comes first and the velocity rules second. For a cell at elevation $z$ and depth $d$ below the ground, the rules of Table 11 are applied from the top down.

• Edifice. The edifice footprint covers the Rainier volcanics and ice mapped within 12 km of the summit (Columbia Crest, configs/domain.yaml) on the 100 m surface grid, closed morphologically over 500 m (a 5 × 5 cell element) with holes filled, keeping the connected region that contains the summit (S3, rainier3d.geomodel.rules.edifice_footprint; geometry.edifice in configs/units.yaml, m1_placeholder values). Its base is the pre-volcanic surface: the DEM elevations of the plutonic and supracrustal outcrops in a two-cell (200 m) ring around the footprint, interpolated linearly across it (nearest value where the linear interpolant is undefined) and kept at or below the glacier bed. The base lies at 804–2234 m NAVD88 over a 252 km² footprint (/surface node of data/processed/geomodel.zarr).
• Basement. Beneath the deposits and the edifice, the basement unit is taken from the nearest outcrop: every map cell takes the plutonic or supracrustal unit mapped there, or else that of the nearest such outcrop on the 100 m surface grid (nearest-neighbour fill, rainier3d.geomodel.rules.bedrock_units, S3). Map symbols of the 1:100,000 state geologic map (Washington Geological Survey, 2026a) reach the model units only through the ordered regular-expression rules of configs/units.yaml (crosswalk_rules).
  • The supracrustal units (the Ohanapecosh, Fifes Peak and Stevens Ridge formations, the Eocene volcanic rocks and the Puget Group, and also the pre-Tertiary Russell Ranch Formation) extend to 4 km below sea level (supracrustal_base_z = −4000 m NAVD88 in configs/units.yaml, an m1_placeholder); the middle crust fills the column below.
  • The Miocene plutons (Tatoosh, White River, Carbon River and Nisqually; map symbols starting Mi, MOi, PLOi, PLMi) extend to 10 km below sea level (pluton_base_z = −10000 m NAVD88 in configs/units.yaml, an m1_placeholder). They therefore continue through L2 into the top four 1 km layers of L3, where they cover about 700 km² in plan view (2,706 L3 cells).
  • The Mashel Formation (Miocene sedimentary rocks, map symbols starting Mc) is capped at 300 m thickness below the ice and deposits (max_thickness_m: 300 in configs/units.yaml, an m1_placeholder); below that, the column takes the nearest mapped basement unit other than the Mashel (rainier3d.geomodel.rules.build_level).
• Unconsolidated deposits. Their thickness beneath any ice is set per map symbol where the Description of Map Units of the 1:100,000 map states one (geometry.symbol_thickness_m in configs/units.yaml), taking the midpoint of the stated range, from zero when only a maximum is given (an m1_placeholder rule). Eight symbols state a thickness (Table 12); they cover 7.0% of the box and 26.3% of its unconsolidated cells. Elsewhere the thickness is the default of the model unit (geometry.unconsolidated_thickness_m, all m1_placeholder): glacial drift 15 m, alluvium and colluvium 10 m, lahar deposits 15 m, and 5 m for mapped open water, which is given the properties of saturated alluvium. S1 writes the result as unconsolidated_thickness on the surface grid. Most deposits are thinner than half a 50 m L1 cell, so they appear only where a cell centre falls inside them: 6618 of the 7.39 million L1 cells (6189 with the unit defaults alone).

• Magma body. A slow body below the summit is a sphere centred beneath Columbia Crest at 11.5 km below sea level (center_z = −11500 m NAVD88) with a radius of 6.5 km, so it spans 5–18 km below sea level (geometry.magma_body in configs/units.yaml, source key moran_1999). As unit magma_mush it replaces every non-air unit it contains: 305 L2 cells and 1168 L3 cells. The weight of its anomaly (magma_weight) is 1 within 4.5 km of the centre (389 L3 cells) and falls to 0 at the radius as a squared cosine over 2 km (radius_m, taper_m, m1_placeholder). It follows the low-Vp body imaged 5–18 km beneath the summit by Moran et al. (1999). [TODO: confirm whether the 5–18 km of Moran et al. (1999) is below sea level or below the summit; configs/units.yaml line 59 marks the depth reference as still to verify]

5.3.2 Hydrothermal alteration from the helicopter electromagnetic survey

Clay-bearing altered rock is electrically conductive, while fresh lava and ice are resistive. The 1996 helicopter electromagnetic (EM) and magnetic survey of Rainier (Rystrom et al., 2000) measured apparent resistivity at four frequencies (33 kHz, 4737 Hz, 4341 Hz and 837 Hz) on a 50 m grid. Finn et al. (2001) used it to map collapse-prone altered zones. Script S22 turns these grids into an alteration field in three steps.

1. Registration. The grids are converted from NAD27 to the model datum. Their orientation is checked against the data, not assumed: the magnetic anomaly correlates with high-passed topography at 0.60 as read, and at 0.08 if flipped north–south.
2. Intensity. The survey saturates at a different resistivity at each frequency, so each frequency $f$ is referred to its own fresh level $L_f$. $L_f$ is the median log resistivity over the surveyed edifice lavas: 4.35, 3.89, 3.74 and 3.00 (log₁₀ Ω m). The intensity is $$a_f = \mathrm{clip}\!\left(\frac{L_f - 0.3 - \log_{10}\rho_a}{0.7},\,0,\,1\right),\qquad{(8)}$$ which starts at half the fresh resistivity and reaches 1 at a tenth of it.
3. Depth. Each frequency senses to a depth of investigation of half its skin depth, $503\sqrt{\rho_a/f}$ m, capped at 150 m. The three-dimensional field is the maximum over frequencies of $a_f$ tapered below that depth. Depth is measured below the glacier bed, and the field is restricted to the edifice lavas and young volcanic rocks.

Over the 126 km² of surveyed edifice lava, 10.6 km² reaches an intensity of 0.5 or more (Fig. 9). That corresponds to 1.9 km³ of equivalent fully altered rock in the model grid. Table 13 gives the distribution around the summit.

The concentration east of and around the summit agrees with the exposed east–west belt of altered rock described by Finn et al. (2001) and John et al. (2008). The EM field sees only the top ~150 m. It cannot see the altered rock buried beneath fresh cover on the upper west flank, which Finn et al. (2001) inferred from magnetic modelling. The survey’s terrain-correlated apparent magnetisation is published as a separate map.

5.3.3 Fusion with the regional models

Regional models. The regional model has two parts:

• Down to 9.9 km below the ground: the USGS Cascadia velocity model v1.7 (Wirth et al., 2025), for Vp and Vs, from its levels L01 (0–100 m every 10 m, then every 100 m to 1.2 km) and L2 (1.5–9.9 km every 300 m) on a 200 m horizontal grid, with L2 interpolated onto the L01 grid; its level L3 (10.8–59.4 km) is not used (rainier3d.fusion.regional.load_cvm). Its depth axis is below the ground surface. For the shallow level this is confirmed by a median top-sample Vs of 194 m s⁻¹ over the model domain (data/raw/regional/cvm17_domain.nc; 206 m s⁻¹ over the larger area of configs/sources.yaml), and for the deeper level it is inferred from continuity: at the summit, Vs is 2841 m s⁻¹ at 1.2 km in L01 and 2870 m s⁻¹ at 1.5 km in L2 (configs/sources.yaml, key cvm17).
• Below 9.9 km: CRESCENT Gen0 Vs (He et al., 2026), with Vp from Brocher (2005) through his regression $V_P = 0.9409 + 2.0947V_S - 0.8206V_S^2 + 0.2683V_S^3 - 0.0251V_S^4$ (km s⁻¹). The CRESCENT depth axis (−4 to 100 km) is taken as kilometres below sea level; the model is resampled onto a 2 km UTM grid padded by 5 km around the box and interpolated linearly to the cells, and its Vs uncertainty is carried as vs_unc_regional.

From one regional model to the other. At 9.9 km below the ground, the base of the Cascadia data, CRESCENT is faster than the Cascadia model: the jump Δ = ln(V_CRESCENT / V_CVM) of Vs has a median of 0.062 over the columns of the domain (5th–95th percentiles 0.017–0.085; regional.cvm_crescent_jump in outputs/paper/numbers.json). Δ of each column is spread over 8–16 km below the ground with a raised-cosine weight W(d), 0 at 8 km and 1 at 16 km: ln V = ln V_CVM + WΔ above 9.9 km and ln V = ln V_CRESCENT − (1 − W)Δ below (fusion.regional_blend_m in configs/domain.yaml, m1_placeholder; rainier3d.fusion.regional). The regional model is continuous at 9.9 km, each model is used only where it has data, and the window follows the topography. A switch at 9.9 km, as in earlier versions of the model, left this step in the model: the static correction of the calibration is continuous in depth, and the update of the local-earthquake tomography (Sect. 5.3.7) is too smooth to remove it.

Below 1 km, both are multiplied by a static, depth-dependent correction fitted to the PNSN arrivals (Sect. 5.3.4). The corrected regional Vs is capped at Vp/1.6, so that its Vp/Vs is at least 1.6 (Christensen, 1996); the geology side is not bounded this way.

Strategy. Neither model is right at every wavelength. The regional models resolve the crust at wavelengths of several kilometres and are constrained by teleseismic, ambient-noise and local-earthquake data; they carry no information on the mapped units. The geology model carries the mapped units and their contrasts, but its absolute level rests on laboratory-style rules. The fusion keeps each where it is informative: the regional model supplies the long wavelengths and the geology model the short ones. The combination is done in the logarithm of velocity, so contrasts are relative: $$\ln V = \mathrm{LP}_{\lambda_c}\!\left(\ln V_{\mathrm{reg}}\right) + \left[\ln V_{\mathrm{geo}} - \mathrm{LP}_{\lambda_c}\!\left(\ln V_{\mathrm{geo}}\right)\right],\qquad{(9)}$$ where LP is a horizontal Gaussian low-pass filter with half amplitude (a quarter of the power) at the cutoff wavelength $\lambda_c$, that is $\sigma = \sqrt{2\ln 2}\,\lambda_c / 2\pi$.

Implementation details.

1. Fused quantities. Vs and Vp/Vs are fused rather than Vp and Vs separately, so the fused Vp/Vs always lies between the two input ratios.
2. Cutoff wavelength. $\lambda_c$ grows with depth below the ground: 6 km at the surface, 10 km at 2 km, and 20 km from 10 km down. The values approximate the lateral resolution of the regional models.
3. Air cells. Cells above the ground are excluded from the filters by normalised convolution.
4. Overshoot. A Gaussian high-pass overshoots at sharp contrasts and puts fast rims around slow bodies. Six alternating projections remove them: each cell is clamped between its geology and regional values, then the regional low-pass is restored.
5. Top kilometre. The top 300 m keep the geology model unchanged, and a linear taper reaches the fused model at 1 km. The regional models do not resolve this depth range; their shallowest layer is a generic near-surface model.
6. Rule contacts. The column rules of Sect. 5.3.1 make flat-lying contacts (the edifice base, the base of the cap and of the Mashel Formation) and vertical unit boundaries that would reach the model as steps of up to a factor of two between neighbouring cells; the horizontal low-pass of Eq. 9 does not act on flat contacts, and the geology is kept as it is in the top 300 m. Below the 150 m critical-zone columns, the geology is therefore low-passed before the fusion, in the vertical with a cutoff that grows to 1 km 250 m below the column base and horizontally with a cutoff that grows from zero at the column base to 750 m 150 m below it and then with the depth below it, to 2 km (lambda_z_m, lambda_h_ramp, lambda_h_start_m, lambda_h_m in configs/domain.yaml, m1_placeholder; rainier3d.fusion.blend.smooth_geology). The critical zone is neither smoothed nor smoothed into the cells below, and mapped contacts stay sharp in it and are softened more with depth below it, where the extrusion rule is least constrained. The horizontal cutoff starts from zero because the critical zone above is not smoothed: a cutoff of 750 m in the first cell below it would put a vertical step at the column base next to every unit boundary (in L1 of the geology, the 95th percentile of the step of ln Vs at the base rises from 0.045 to 0.081); the ramp keeps it at 0.047. The rate of the ramp trades this step against the lateral steps just below the base (99th percentile 150–1000 m below the ground 0.148 with the ramp, 0.124 without it, 0.256 with no lateral smoothing, in the same geology).
7. Magma body. Its anomaly (Sect. 5.1.1) is applied to the fused model, so neither the low-pass nor the clamping of item 4 removes it.
8. Density. Density follows Vp through the ratio of Nafe–Drake densities, which preserves the density contrasts between units.

Interfaces that are not geology. Three construction steps can leave steps in the model that no rock contact explains: the boundary between the levels L2 and L3 at −6 km, the change from the Cascadia model to CRESCENT, and the vertical walls of map units extended to depth. Table 14 measures them as the absolute step of ln Vs between neighbouring cells of a uniform 250 m grid (rainier3d.report.roughness), before and after the changes described above. “Before” is the model built from the code and calibration of 14 October 2026 (a switch from the Cascadia model to CRESCENT at 9.9 km, a horizontal cutoff growing slowly from zero below the critical zone, and each level clamped to its own cell centres when resampled); “after” is the published model, recalibrated. Both are measured with the same definitions. The step at the L2/L3 boundary is not a disagreement of the levels: in the published model, the median change of ln Vs from the bottom cell centre of L2 to the top cell centre of L3, 625 m lower, is +0.0077, and the vertical gradient of L2 predicts +0.0073 over that gap (model.l2_l3_continuity). Clamping each level to its own cell centres put the whole change over the gap into one 250 m cell of the uniform grid; interpolating across the gap spreads it over the cells it spans. The uniform grids of S9 and S26 and of the client now interpolate linearly across the gap between levels (rainier3d.export.grids.uniform). The step from the Cascadia model to CRESCENT is measured on the steps inside the blend window, 8–16 km below the ground, without the pair of cells across a level boundary, which has its own row; a switch at 9.9 km and the steepest part of a blend fall in the same window. Its reference is the steps 16.5–20.5 km below the ground, where only CRESCENT is used. Percentiles over all steps are compared with the reference; the largest step per column is given too but has no reference of the same length, because the largest of 32 steps exceeds the largest of the 16 that the reference window holds. The blend lowers the 99th percentile from 0.017 to 0.0054, 2.5 times the reference: the largest step per column, 0.0032 at the median, is the steepest slope of a raised cosine that spreads a 6% jump over 8 km (0.003 per 250 m), so a smaller step needs a wider window, which would shift CRESCENT further. The all-step reference at 7–9 km below the ground now overlaps the top of the blend window and the steeper depth profile of the recalibrated regional correction (Table 18), so it is no longer a quiet reference. The lateral filter lowers the largest steps (95th and 99th percentiles) and raises the median, because a wall is spread over several cells.

Fig. 10 shows the three models along a section through the summit. The fusion is checked by an invariant (Table 15): below 1 km depth, the low-passed fused model must equal the low-passed regional model within an RMS of 0.03 in ln V.

5.3.4 Calibration on P and S−P travel times

Data. The data are the analyst picks of the 88 PNSN earthquakes with M ≥ 2 in the box between 2015 and 2026 that have at least six picks, four of them P. They comprise 1823 P and 1280 S picks at 50 stations, with origins and picks from the USGS ComCat catalogue (U.S. Geological Survey, 2026d). The event list is committed (configs/validation_events.csv), so every run uses the same set. Picks are weighted by the RMS of PNSN’s own location residuals: 0.14 s for P and 0.23 s for S.

Travel times. They are computed on a 500 m grid with pykonal’s point-source solver (White et al., 2020), one solve per station and phase by reciprocity.

• Solver accuracy. Against analytic travel times on a grid of the same size (70 × 75 × 24 km) and spacing (0.5 km), for an off-grid source and 360 receivers, its RMS error is 8 ms in a uniform 6 km s⁻¹ model and 11 ms with the gradient 4 + 0.1z km s⁻¹; fteikpy (second-order fast sweeping) gives 22 and 46 ms. At 0.25 km spacing the pykonal errors fall to 4.7 and 5.7 ms at about nine times the cost (scripts/bench_eikonal.py, docs/eikonal_benchmark.md, pykonal 0.4.1, run on 24 September 2026).
• Topography. The travel-time grid is the S6 grid (validation block of configs/domain.yaml): 141 × 151 × 70 nodes at 500 m, from 4.5 km above to 30 km below sea level. Cells more than one cell (500 m) above the ground carry the speed of sound in air, 343 m s⁻¹, so rays follow the rock and cannot cut across valleys (rainier3d.validate.locate.with_topography). The one-cell skin keeps every station in rock; station elevations match the DEM within ±70 m (median +5 m). Against rock-filled air, the air handling changes the station-mean residuals by a median of 0.2 ms.
• Cross-check. For station OBSR, NonLinLoc’s Grid2Time P times on the exported NonLinLoc grids (S9, rainier3d_nll_500m.tar.gz) agree with the pykonal times to 12 ms on average (RMS 31 ms, 85 events). This number comes from a separate session’s record and has not been rerun.

What is estimated. Fifteen parameters $\boldsymbol\theta$ scale the two models that are fused, and the hypocentres and origin times of the 88 events are re-estimated in every model tried. The travel-time data therefore constrain the velocity parameters only through the part of the misfit that relocation cannot absorb. This avoids the circularity of scoring a model with catalogue hypocentres that were located in a 1D model with station corrections.

Parameters.

• Rock physics (3). Log-multipliers on three quantities shared by all rock units: the zero-pressure velocity V₀ and the crack-closure pressure P* of Eq. 1, and Vs. The Vs multiplier acts at fixed Vp, so it is a change of Vp/Vs. Ice and unconsolidated deposits are not scaled.
• Regional correction (12). A static depth profile $m(d)$ of the log-factor that multiplies the regional velocity, $V_{\mathrm{reg}} \to V_{\mathrm{reg}}\,e^{m(d)}$, one profile for Vp and one for Vs. Here $d$ is depth below the ground. Each profile is piecewise linear between knots at 0, 1, 2, 4, 7, 11, 16 and 25 km and constant below 25 km. The values at 0 and 1 km are fixed at zero, because the fused model takes the top kilometre from the geology model (Sect. 5.3.3); the six deeper values are free.

Forward problem. For a trial $\boldsymbol\theta$, the model is rebuilt from its inputs by the same code that builds the published model. At $\boldsymbol\theta = 0$ it returns the uncalibrated model cell for cell, which checks that the calibration fits the model that is delivered and not an approximation of it. The steps are:

1. Vp of every rock-unit cell from Eq. 1 with the scaled V₀ and P*, and Vs from that Vp through the unit’s Vp/Vs ratio, or the regression of Brocher (2005) where the unit has none, times the Vs multiplier;
2. the regional Vp and Vs multiplied by $e^{m(d)}$;
3. the fusion of Eq. 9, including the geology-only top 300 m and the taper to 1 km;
4. P and S travel times from every station on the 500 m grid, with air above the ground, as described under Travel times;
5. every event relocated in that model (see Relocation below).

Objective. With $t_{ij}$ the pick of phase at station $j$ for event $i$, $T_j$ the travel-time field, and $(\mathbf{x}_i, \tau_i)$ the hypocentre and origin time, the normalised residual is $$r_{ij} = \frac{t_{ij} - \tau_i - T_j(\mathbf{x}_i;\boldsymbol\theta)}{\sigma_{\mathrm{ph}}},\qquad{(10)}$$ where $\sigma_{\mathrm{ph}}$ is the pick uncertainty of the pick’s phase: $\sigma_P = 0.14$ s for P picks and $\sigma_S = 0.23$ s for S picks. The calibration minimises $$\Phi(\boldsymbol\theta) = \sum_{ij} w_{ij}\, r_{ij}^2 + \boldsymbol\theta_g^{\mathsf T} \mathbf C_g^{-1} \boldsymbol\theta_g + \lambda_s \lVert \mathbf D \mathbf m \rVert^2 + \lambda_d \lVert \mathbf m \rVert^2,\qquad{(11)}$$ where the hypocentres minimise the first term for each $\boldsymbol\theta$. The weights $w_{ij}$ are the Huber weights of the relocation (threshold 1.5σ), so outlying picks count less. $\boldsymbol\theta_g$ holds the three rock-physics log-multipliers, with prior standard deviations of 0.3, 0.7 and 0.1 on the diagonal of $\mathbf C_g$ (factors of 1.35, 2 and 1.1). $\mathbf m$ holds the regional log-factors at all eight knots, with the two fixed at zero, and $\mathbf D$ takes second differences along depth. The smoothing weight $\lambda_s$ and the damping weight $\lambda_d = 0.01$ are relative: both are multiplied by the mean diagonal of the data term for the regional parameters.

Update. Each iteration linearises the residuals about the current model and hypocentres, $\mathbf r(\boldsymbol\theta + \delta\boldsymbol\theta) \approx \mathbf r - \mathbf G\,\delta\boldsymbol\theta$, with $\mathbf G = \partial \mathbf T / \partial \boldsymbol\theta$ scaled like the residuals.

• Separating hypocentres from velocity. For each event, the columns $\mathbf H_i$ of partials with respect to $(\mathbf x_i, \tau_i)$ are formed at its current location. The event’s rows of $\mathbf G$ and $\mathbf r$ are projected onto the orthogonal complement of $\mathbf H_i$, through a QR factorisation of $\mathbf H_i$ (Pavlis and Booker, 1980). What remains is the part of the misfit that no shift of that hypocentre can remove. Events with five or fewer usable picks keep at most one degree of freedom after the projection and are left out of the update.
• Step. With the projected $\tilde{\mathbf G}$, $\tilde{\mathbf r}$ and the regularisation matrix $\mathbf L$ of Eq. 11, $$\boldsymbol\theta \leftarrow \boldsymbol\theta + \left(\tilde{\mathbf G}^{\mathsf T}\tilde{\mathbf G} + \mathbf L\right)^{-1}\left(\tilde{\mathbf G}^{\mathsf T}\tilde{\mathbf r} - \mathbf L\,\boldsymbol\theta\right).\qquad{(12)}$$
• Relocation. Every event is then relocated in the updated model before the next linearisation, so velocity and hypocentres are updated in alternation, as in the minimum one-dimensional model of Kissling et al. (1994). Relocation is a grid search over nodes at or below the ground, within 15 km of the catalogue epicentre, on the sum of |r| with the median origin time, followed by Huber Gauss–Newton (Huber, 1964) on interpolated travel times with a penalty on hypocentres more than 50 m above the ground. It follows the grid-search-then-refine design of NonLinLoc (Lomax et al., 2000).

Derivatives.

• Rock physics. A rock-physics multiplier changes every rock-unit cell, and its effect passes through the fusion filters, so its column of $\mathbf G$ is a one-sided finite difference ($h = 0.05$ in log) through the whole forward problem, with the hypocentres held fixed. That costs one full rebuild and eikonal solution per parameter, and it is recomputed at every iteration on the fitting half.
• Regional correction. For knot $k$ the derivative is the ray integral $$\frac{\partial T}{\partial m_k} = -\int_{\mathrm{ray}} \beta(d)\, \phi_k(d)\, s\, \mathrm d\ell,\qquad{(13)}$$ where $s$ is slowness, $\phi_k$ the linear interpolation weight of knot $k$ at depth $d$, and $\beta(d)$ the share of the regional model in the fused model: 0 above 300 m, rising linearly to 1 at 1 km. Rays are traced down the gradient of each station’s travel-time field (Thurber, 1983). The integral assumes that a depth-only factor on the regional model passes through the fusion unchanged. The low-pass filter of Eq. 9 passes a factor that varies slowly across the box, but the clamping of fused values between the two input models is not linear, and is the likely cause of the difference below. Against a finite difference through the full forward problem, for the Vs knot at 4 km over all 1280 S picks, the ray derivatives correlate at 0.97 with a slope of 0.81: they underestimate the sensitivity by about 20%. The residuals are always evaluated with the full forward problem, so this error does not bias the data fit. At convergence it acts, to first order, like regularising the regional correction about 20% more weakly than $\lambda_s$ and $\lambda_d$ state.

Validation and choice of smoothing. Events alternate between a fitting half and a held-out half in origin-time order, and held-out events are relocated in every iterate, so their misfit is an independent score. At the first iteration, $\lambda_s$ is chosen from 0.01, 0.1, 1 and 10 by the normalised RMS of the held-out separated residuals predicted by the linearised step: 0.713, 0.714, 0.723 and 0.760. The adopted value is 0.01, the smallest tried and so at the edge of the scan; it differs from 0.1 by 0.0009 in the score, so the score does not resolve the weight below 0.1.

Iterations and uncertainties. Three iterations on the fitting half are followed by two on all events, which reuse the last rock-physics derivatives (Table 16); the run takes 40 min with 48 workers on a Linux server. Posterior standard deviations are the square roots of the diagonal of $(\tilde{\mathbf G}^{\mathsf T}\tilde{\mathbf G} + \mathbf L)^{-1}$ at the last iteration, with the covariance multiplied by the mean reduced χ² of the two phases, 0.562, the mean of 0.438 for P and 0.687 for S (scripts/13_joint_calibration.py). They are linearised. They do not include the error of the ray derivatives, or the trade-off with hypocentres beyond what the projection removes.

The fit is resolved as follows (Table 17, Table 18, Fig. 11, Fig. 12).

• V₀. It is the parameter the data require: ×1.45 ± 0.03 (×1.31 without the critical zone, in a calibration that also predates the smoothing of the geology in the fusion). The table rock is too slow near the surface. With the multiplier, the zero-pressure Vp of the rock law is 4.1 km s⁻¹ for the Rainier volcanics and 4.9 km s⁻¹ for the Ohanapecosh Formation. Because the multiplier is global, the Miocene plutons reach the cap of 0.98 V∞ = 6.08 km s⁻¹ at zero pressure (4.5 × 1.45 = 6.5), so their rock law starts at the cap; the velocity drop toward the surface comes from the critical zone alone.
• P*. It is not resolved: its log-multiplier ranged from −0.29 to +0.31 across the iterations after the first and was still rising at the last one (final ×1.36, posterior s.d. 0.15 in ln; the linearised s.d. does not capture this drift), while the held-out RMS changed by less than 2 ms after the second iteration. It trades off against V₀: over the last three iterations both rose (ln V₀ from 0.240 to 0.373, ln P* from 0.008 to 0.305), and the calibration of 21 October 2026, with a lateral smoothing that started at 750 m directly below the critical zone, ended at V₀ ×1.32 and P* ×0.56.
• Rock Vs. It falls by 4.2 ± 1.7% (ln multiplier −0.043 with posterior s.d. 0.017, against a prior s.d. of 0.1). Because it acts at fixed Vp, it raises the Vp/Vs of every rock unit by about 4%.
• Regional correction. It changes Vp by at most 7% (+6.9% at 2 km, −4.8% at 7 km). It makes Vs 5–9% faster at 2–7 km and 3–7% slower at 16–25 km; with the rock multipliers, the median Vp/Vs at 2–4 km below the ground falls from 1.83 in the regional model to 1.76 in the fused model.
• Agreement between the two models. After calibration, the geology model and the corrected regional model agree to within 5.1% between 0.3 and 4 km depth (largest difference +0.050 in ln Vs at 0.3–1 km). Uncalibrated, the geology is up to 20% slower (−0.227 in ln Vp at 0.3–1 km), and the fusion invariant exceeds its tolerance in L1 (0.031 for Vp, 0.032 for Vs). These are mean ln(V_geology / V_regional) by depth below ground over L1–L3 without air and ice, from model.zarr and from the uncalibrated model that the validation builds with the same code and the critical zone (model_nocal.zarr, pixi run nocal; outputs/paper/numbers.json). In the top 300 m, which the fused model takes from the geology alone, the calibrated geology, critical zone included, remains 19% faster in Vp and 25% faster in Vs than the regional model (mean ln ratios 0.17 and 0.22).

Why this parameterisation. Two simpler parameterisations were tested against the same data (Table 19).

• Vs-only depth factor, catalogue hypocentres fixed. It fits S as well as the adopted model, but pushes Vp/Vs to 1.60–1.67 at 2–4 km. That is low for crustal rock and lower still than expected beneath a volcano with a hydrothermal system: the S−P misfit of fixed hypocentres is partly a location error. This is the earlier S12 calibration (scripts/12_calibrate_vs.py, retired and kept in the git history; configs/vs_calibration.yaml), which raised Vs by up to 10% at 2–7 km; it left 30,086 fused cells below Vp/Vs = √(8/3) in the model of that time, against 19,476 for the adopted model.
• Depth factors on the regional Vp and Vs at all depths, with relocation. It fits as well as the adopted model, but places the correction in the regional model where the fusion does not use it. It requires Vp 7–19% faster in the top 2 km, where the fused model follows the geology, and it breaks the fusion invariant (0.050). This is the first S13 run (configs/velocity_calibration_v1.yaml): Vp factors of 1.19, 1.14 and 1.07 at 0, 1 and 2 km, and an L1 invariant of 0.050 for Vp and 0.044 for Vs (outputs/relocation/s5_joint.log).
• Adopted parameterisation. Fitting the geology’s rock physics and correcting the regional model only below 1 km fits equally well, satisfies the invariant (0.011 for Vp and 0.012 for Vs, L1 in outputs/fusion_report.csv) and keeps the unit contrasts in physical parameters.

Alteration and the calibration. The first calibration without the critical zone was run with a conduit-centred alteration field; replacing it with the EM-based field of Sect. 5.3.2 changed the relocated RMS from 0.095 to 0.093 s for P and left S at 0.190 s. The calibration with the critical zone was run directly on the EM-based field.

5.3.5 Validation

Held-out events. The 44 events of the held-out half are relocated in every iterate but enter no fit until the last two iterations (Table 16). After the three iterations on the fitting half, their RMS has fallen from 0.125 to 0.094 s for P and from 0.355 to 0.188 s for S, while the fitting half ends at 0.096 and 0.191 s; the held-out misfit is not larger than the fitted one. These are the independent scores of the calibration. The smoothing weight $\lambda_s$ is chosen on the same half (see Validation and choice of smoothing), so they are not fully blind to it. The two final iterations use all 88 events, and the published model (critical zone and EM-based alteration, data/processed/model.zarr) relocated with S14 gives 0.093 s for P (949 picks) and 0.188 s for S (655 picks) on the same 44 events, against 0.131 and 0.263 s in the PNSN 1D model; these events are no longer independent of that model. The split is in time, not in space: held-out and fitting events share source regions and ray paths (outputs/paper/numbers.json, scored from outputs/relocation/<model>/picks.csv with the split of S13: event identifiers in origin-time order, every second one held out).

All models scored the same way. Table 19 relocates all 88 events in each model with the same locator and topography. Relocation alone does not rescue the uncalibrated three-dimensional model: its S residuals stay larger than those of the 1D model. The S−P misfit therefore lies in the velocity model, not in the catalogue hypocentres. With the ground bound, no event is placed at the top of the grid, and one of the 88 events ends more than 50 m above the ground in the published model, at the limit of the ground penalty (outputs/relocation/fused_joint/stats.json, above_ground_gt_50m = 1; three in the conduit-alteration run).

Independent tomography. The iMUSH project imaged Vp and Vs around Mount St. Helens from local earthquakes and explosions recorded by a 70-station array (Incorporated Research Institutions for Seismology, 2020; Ulberg et al., 2020).

• Coverage. Its grid spans the whole model box. Its Vp is resolved over 82–90% of the box. Its Vs is resolved over 83–93% of the southern strip (south of 46.65° N) but only 11–19% of the northern half. Script S21 samples our models at its 51,300 resolved nodes below ground. The model is the EarthScope EMC file iMUSH-localEQ-Ulberg-2020.r0.0-n4c.nc, with 1.2 km nodes masked where checkerboard tests do not recover 20 km features; its depth axis is read as kilometres below sea level, as the EMC page states, although the netCDF long_name says below the surface. Coverage is counted over 0–14 km below sea level, and north of 46.75° N between 188 and 1354 Vs nodes fall in each depth bin.
• Beneath Rainier, north of 46.75° N (Table 20). The iMUSH model finds the regional model’s Vs 5–8% too slow at 1–8 km and up to 3% too fast at 11–20 km, the shape of Table 18. With the correction, the difference in Vs is 2% or less below 2 km. Vp/Vs moves from 1.83–1.85 towards the iMUSH value of 1.75–1.77; at 4–8 km our Vp/Vs (1.69) is 0.06–0.07 below it.
• South of 46.65° N, inside the iMUSH array. The iMUSH Vs is within 1% of the Cascadia model at 1–6 km, and the correction makes our Vs about 7% too fast at 2–6 km. This suggests that the correction is local to Rainier, which a laterally varying correction would test.
• Independence. Both studies use PNSN arrivals from 2015–2016. The iMUSH-only stations and explosions make the comparison largely independent in the south, less so in the north.

Catalogue hypocentres. At the ComCat hypocentres, which were located in a 1D model, the three-dimensional model predicts later arrivals than the 1D model at every station: by 0.27 s for P (mean over 50 stations of the station-mean difference) and 0.45 s for S (45 stations with at least five picks); per pick, the means are 0.23 and 0.41 s (outputs/pnsn_residuals.csv). With each event’s mean residual removed, the RMS is 0.130 s (P) and 0.227 s (S) for the 3D model, against 0.159 and 0.287 s for the 1D model. The per-station 3D − 1D delay correlates with the mean 1D residual at 0.72 for P and 0.68 for S (Fig. 15). The 3D structure therefore explains part of what station corrections absorb.

5.3.6 The fused model

Figs. 16, 17 show the fused model along two sections through the summit, and Fig. 18 compares vertical profiles.

• Level L1 (the ground to sea level). L1 is on average 6% faster in Vp and 8% faster in Vs than the regional model. The calibrated rock is stiffer than the shallowest layer of the Cascadia model, a generic near-surface model with a median top-sample Vs of 194 m s⁻¹ over the model domain. The travel times constrain this difference only beneath the stations. The averages are the mean ln(V_fused / V_regional) over all non-air L1 cells, from 0 to 4.4 km NAVD88: +0.057 for Vp and +0.079 for Vs (outputs/fusion_report.csv). The difference is largest in the top 300 m, where the calibrated geology, critical zone included, is 19% faster in Vp and 25% faster in Vs than the regional model.
• Plutons. The Miocene plutons stand out as fast columns to 10 km below sea level, the base set by pluton_base_z. With the calibrated V₀ their rock law starts at the cap of 6.08 km s⁻¹, and their geology Vp is 6.09 km s⁻¹ just below the critical zone (200–300 m) and 6.13 km s⁻¹ at 2 km depth (L1 medians in data/processed/properties_geology.zarr).
• Puget Group. The Puget Group block west of the summit is slow down to 4 km below sea level. It is unit eocene_sedimentary (V₀ = 2.6 km s⁻¹, V∞ = 5.2 km s⁻¹ in Table 9), whose base is the supracrustal base at −4 km NAVD88, and on the summit row it extends from 571.6 to 581.9 km easting (UTM 10N, cell centres).
• Seismicity. Summit earthquakes form a column from the edifice to about 3 km below sea level. WRSZ earthquakes concentrate 4–12 km below sea level, 12–18 km west of the summit. The sections show 4369 PNSN events of M 0.5–3.4 from 2015 to 2026 (web/atlas/data/events.geojson, depths taken as below sea level). Of the 1559 events within 2 km of A–A′ and 3 km of the summit, 90% lie between 2.3 km above and 1.9 km below sea level and 99% above 2.6 km below. Of the 329 West Rainier Seismic Zone (WRSZ) events within 2 km of A–A′ and 8–25 km west of the summit, 90% lie at 4.5–11.8 km below sea level and 10.8–18.2 km west.
• Magma body. The slow body (5–18 km below sea level, Sect. 5.3.1) is applied to the fused model after the fusion, so it is the one structure of the model that the regional models do not constrain at its own wavelength. In L3, where it lies, its core (weight 1) is 15.1% slower in Vs than the middle crust at the same elevations, and the whole body, taper included, 9.4% slower (mean ln ratios −0.163 and −0.099; outputs/paper/numbers.json); L2 holds only cells of its taper. Whether a body of the size imaged by Moran et al. (1999) and Pang et al. (2025) belongs in the model is a question for data that resolve it: the tomography of Sect. 5.3.7 has almost no rays below 12 km.

5.3.7 Local-earthquake tomography

The calibration adjusts global factors and a depth profile, so it cannot correct lateral errors of the geological rules between the wavelengths that the geology encodes and those that the regional models resolve. A local-earthquake tomography (S37, rainier3d.validate.tomo; Supplement S1) perturbs ln Vp and ln Vs of the calibrated model on nodes 3 km apart, using 99,716 analyst picks of 3770 PNSN earthquakes from 2015 to the data freeze, with travel times on a 250 m eikonal grid, every event relocated in every trial model, and earthquakes in 10 km squares withheld from the fit to choose the smoothing and stop the iterations. The perturbation is zero in the top 150 m below the ground or glacier bed, so the critical zone keeps its rock physics, and all invariant tests of the model hold in the result. Three steps reduce the RMS residual of the held-out earthquakes from 0.0904 to 0.0722 s for P and from 0.1802 to 0.1206 s for S, and the RMS of the held-out residuals averaged per station from 0.052 to 0.023 s for P and from 0.152 to 0.066 s for S, without station terms. Checkerboard tests recover 12 km structure between sea level and 6 km depth with correlations of 0.66–0.74 and about half of its amplitude; below 12 km, and beneath the summit below 5 km, ray coverage is nearly absent and the perturbation is set by the regularisation. Within the sampled volume the update has a mean close to zero and an RMS of about 5%. The result is distributed as model_let.zarr beside the model product; it does not replace model.zarr until the resolution tests and the station coverage of Supplement S1 have been extended.

6.5.1 How the strain fields are computed

Both fields are first-order estimates. Neither is computed with Green’s functions of the heterogeneous model or with a rheology other than linear elasticity. This subsection gives the calculation as implemented (rainier3d.geodesy.load, rainier3d.geodesy.volume, script S18) and then lists what it leaves out.

Tectonic strain rate. The horizontal strain-rate tensor $\dot\varepsilon_{\alpha\beta}(x, y)$, with $\alpha, \beta \in \{x, y\}$, is the velocity-gradient fit of Sect. 6 at the surface. At depth it is not calculated but assigned: $$\dot\varepsilon_{\alpha\beta}(x, y, z) = \dot\varepsilon_{\alpha\beta}(x, y, 0), \qquad \dot\varepsilon_{zz} = -\frac{\nu}{1-\nu}\left(\dot\varepsilon_{xx}+\dot\varepsilon_{yy}\right), \qquad \dot\varepsilon_{xz}=\dot\varepsilon_{yz}=0.\qquad{(14)}$$ The vertical component follows from a vanishing vertical stress rate (plane stress). This is the thin-plate view of the upper crust: it holds when the sources of deformation are far or deep compared with the depth of interest and the crust is laterally uniform. No source (megathrust locking, fault slip, block rotation) is inverted from the GNSS velocities.

The shear on the WRSZ is resolved as $$\dot\gamma = -\,\mathbf{s}\cdot\dot{\boldsymbol\varepsilon}\cdot\mathbf{n},\qquad{(15)}$$ with $\mathbf{s}$ the horizontal unit vector along the strike and $\mathbf{n}$ the horizontal normal. The sign is positive for right-lateral motion. The strike is the long axis of the covariance of the WRSZ epicentres.

Edifice-load stress. The load is the rock above the pre-volcanic surface of the geology model. Its columns of height $h$ are summed on 500 m cells into vertical point forces $P_k = \rho g h_k A_k$ ($\rho$ = 2500 kg m⁻³). The forces act on a flat reference plane at the mean elevation of the edifice base (1789 m). The stress at a point is the sum over $k$ of the Boussinesq solution for a point force on the surface of a homogeneous, isotropic, linear-elastic half-space (Johnson, 1985). At horizontal distance $r$ and depth $z$ below the plane, with $R = \sqrt{r^2+z^2}$ and tension positive: $$\sigma_{zz} = -\frac{3Pz^3}{2\pi R^5}, \quad \sigma_{rz} = -\frac{3Prz^2}{2\pi R^5}, \quad \sigma_{rr} = \frac{P}{2\pi}\left[\frac{1-2\nu}{r^2}\left(1-\frac{z}{R}\right) - \frac{3r^2z}{R^5}\right], \quad \sigma_{\theta\theta} = -\frac{P(1-2\nu)}{2\pi}\left[\frac{1}{r^2}\left(1-\frac{z}{R}\right)-\frac{z}{R^3}\right].\qquad{(16)}$$ These stresses depend on Poisson’s ratio ($\nu = 0.25$) but not on the elastic moduli. Points within 250 m below the plane are left out, where the point-force field is singular.

Edifice-load strain. The stress is converted to strain with the isotropic compliance of each cell: $$\varepsilon_{ij} = \frac{\sigma_{ij}}{2\mu} - \frac{\lambda\,\sigma_{kk}}{2\mu(3\lambda+2\mu)}\,\delta_{ij}, \qquad \mu = \rho V_S^2, \quad \lambda = \rho V_P^2 - 2\mu,\qquad{(17)}$$ with $V_P$, $V_S$ and $\rho$ taken from the fused model.

• Cone interior. Above the half-space (higher than 1539 m), the stress is the laterally confined overburden, $\sigma_{zz} = -\rho g d$ and $\sigma_{xx} = \sigma_{yy} = \nu\,\sigma_{zz}/(1-\nu)$, with $d$ the depth below the local ground.
• SHmax. It is the eigenvector of the horizontal stress tensor $(\sigma_{xx}, \sigma_{yy}, \sigma_{xy})$ with the most compressive eigenvalue. It is undefined in the cone interior, where the horizontal stress is isotropic. Azimuths are in degrees east of north, from 0 to 180°.

What the calculation does not include.

• Green’s functions of the 3D model. The stress is that of a uniform half-space, so the stiffness contrasts of the model neither concentrate nor divert it. Applying the local compliance (Eq. 17) makes soft rock strain more under the same stress. The stress field, however, is not in equilibrium in the heterogeneous medium, and its axes do not rotate around stiff plutons or the soft edifice.
• Topography as part of the elastic body. The load sits on a flat plane. The cone above it is not part of the elastic domain, and its interior is approximated in one dimension.
• Rheology beyond linear elasticity. There is no viscoelastic relaxation of the warm lower crust or of the magma body, no plasticity, no pore-pressure coupling and no thermal stress.
• A tectonic source model. The GNSS rates are not explained by locking, fault slip or block motion, so their variation with depth (Eq. 14) is assumed rather than computed.
• Absolute stress. The load gives a static stress and GNSS a strain rate. Without the magnitude of the tectonic background stress the two cannot be summed, so their orientations are reported separately.

A consistent model would compute the same quantities in three steps:

1. a static finite-element solution on the model mesh, with gravity acting on the real topography and heterogeneous $\lambda$ and $\mu$ from the fused velocities;
2. Green’s functions computed in the same 3D elastic model, used to invert the GNSS velocities for their sources before computing strain rates at depth;
3. a viscoelastic lower crust for time-dependent loading.

The load pattern of Fig. 23, tangential beneath the summit and radial at depth, is the qualitative signature of any load on a cone. Heterogeneity would change its magnitudes and rotate its axes near stiffness contrasts. The depth extension of the GNSS axes is the least certain part of the product.