Tomography · Part B: the 2016 Geysers picks

Part A built tomography in a 2D world whose answer we planted. It found four things that decide what a tomographic image may be trusted for:

1. bent rays and the adjoint give the gradient for every cell at once, and Gauss–Newton turns it into a model;
2. the data alone do not fix the model: smoothing and damping fill what the rays leave open, and their strength changes the picture;
3. the null test (a planted truth with no structure) measures what the method invents, and a planted anomaly measures what it recovers;
4. velocity and earthquake locations trade off, and solving for both with parameter separation keeps the velocity from absorbing location errors.

This notebook takes the same method to a real field: 453 earthquakes under The Geysers in 2016, recorded by 41 stations, with 23,191 P and S arrival times picked by a neural network (Session 5) and located with ADLoc. The method runs in 3D now, through ADTomo, the group's tomography code: an eikonal solver with its adjoint, written as a PyTorch function.

Every earthquake sends two kinds of waves: P (compressional) waves, which arrive first, and S (shear) waves, which travel slower. So the model has two unknowns at every place, the P speed Vp and the S speed Vs. Their ratio Vp/Vs (about 1.6–1.8 in crustal rock) is worth watching on its own: what fills the pores changes Vp and Vs differently, so Vp/Vs is sensitive to fluid. At The Geysers, rock whose pores hold steam rather than water has been seen as low Vp/Vs (Julian, Ross, Foulger & Evans 1996, Geophysical Research Letters, doi:10.1029/95GL03321). Part A used one wave type to keep the method clear; here both are inverted.

There is no truth here. So the notebook works in three steps:

• the real data raise the question (B.1–B.3): what does a 1D model leave unexplained, and what 3D model explains it?
• a synthetic twin answers how far to trust the answer (B.4–B.5): the same stations, the same earthquakes, the same picks and the same program, but with planted truths;
• the real result is read in the light of the twin (B.6–B.7).

The questions this notebook answers

1. What do the picks tell us that a 1D velocity model cannot explain? (B.2)
2. How does the 3D fit run, and how long does it take? (B.3)
3. In this geometry, what does the method invent, and what can it recover? (B.4)
4. How strong should the regularisation be, and how was that decided? (B.5)
5. What 3D structure do the real data show, and which parts are stable? (B.6, B.7)

B.1 · The data and the grid

Three tables: the stations, the earthquakes (with the locations ADLoc gave them in a 1D model), and the picks. A pick is one arrival time of one phase (P or S) of one earthquake at one station.

41 stations, 453 earthquakes, 23,191 picks (12,141 P, 11,050 S)

	event_id	station_id	phase_type	phase_score	phase_time
0	2016.357.4439	NC.GAC..EH	P	0.992	2016-12-22 22:57:30.570
1	2016.017.4331	NC.GAC..EH	P	0.991	2016-01-17 14:18:41.890
2	2016.094.5301	BG.MCL..DP	P	0.991	2016-04-03 15:24:15.770

tomo.geometry turns the three tables into what the solver needs:

• a grid of nodes 1 km apart, with x east, y north and z down (km). It covers every station and earthquake with a node to spare, and its top is at or above the highest station;
• for every pick, the station, the phase and the earthquake it belongs to;
• the observed time of every pick, in seconds after the earthquake's catalogue origin time.

One eikonal solve per station and phase gives the travel time from that station to every node (Part A, A.1: by reciprocity, a station can act as the source). So the cost of one forward model is set by the number of station–phase pairs, not by the number of picks.

grid 27 x 22 x 11 nodes, h = 1 km; depth -2 to 8 km below sea level
82 travel-time tables (station-phase pairs) serve 23,191 picks
unknowns of a 3D model: 13,068 (Vp and Vs at every node)

A 3D model is hard to see at once, so every model in this notebook is drawn the same way: maps at fixed depths, and two vertical sections through the middle of the earthquakes. A–A′ runs along the field from north-west to south-east, B–B′ across it from south-west to north-east.

We draw the lateral part of each field: its change, in per cent, from the mean of its own depth layer. The 1D part of the model (how speed grows with depth) is shown separately as a profile. The lateral part is what a 1D model cannot have, and it is the question here.

Colours follow one rule throughout: for Vp, red is slow and blue is fast; for Vp/Vs, red is high and blue is low.

Where the stations and earthquakes are. The earthquakes fill a band about 20 km long along A–A′, most of them between sea level and 4 km depth in the catalogue; the stations sit on the hills above them, up to about 1 km above sea level. The grey lines on the right are the depths of the grid's nodes: 1 km apart, so the earthquakes span only four or five of them.

B.2 · The 1D model, and what it cannot explain

A 3D fit starts from a 1D model: one Vp and one Vs per depth, the same everywhere. Part A started from a known background; here the background must be found too, together with the earthquake locations, since the two trade off (Part A, A.8). The program does it in two fits:

1. 1D velocity and locations. Start from a uniform model (Vp 5.4 km/s, Vs 3.0 km/s) with every earthquake at the centre of the grid, locate each earthquake in it, then alternate Gauss–Newton velocity steps and relocation (parameter separation). The profile is smoothed in depth, and Vp and Vs may not decrease with depth. The fit is run from three uniform starts and the one with the lowest loss is kept, in case one start stops in a local minimum.
2. Cut, then add station terms. Picks whose residual is beyond 5 robust standard deviations are removed from the velocity fit (they still locate their earthquake). The robust standard deviation is 1.4826 times the median absolute deviation from the median: for Gaussian errors it equals the standard deviation, but a few wild picks cannot inflate it. Then the fit is repeated with one station term per station and phase: a constant delay, for the rock just under the station that no 1D model can hold (Session 2, A.10).

Every fit uses the Huber loss (Session 2, A.6), with its threshold at one pick standard deviation, so a bad pick weighs like its distance, not its square.

These two fits take too long to run in class, so we load their result. B.3 runs the 3D fit live.

 depth (km)  Vp (km/s)  Vs (km/s)  Vp/Vs
       -2.0       4.03       2.43  1.658
       -1.0       4.03       2.43  1.658
        0.0       4.19       2.53  1.653
        1.0       4.24       2.58  1.644
        2.0       5.05       3.10  1.631
        3.0       5.20       3.21  1.620
        4.0       5.64       3.50  1.610
        5.0       5.66       3.52  1.609
        6.0       5.67       3.52  1.608
        7.0       5.67       3.52  1.608
        8.0       5.67       3.52  1.608

Vp grows from 4.0 km/s above sea level to 5.6 km/s at 4 km and then levels off; Vp/Vs falls from 1.66 to 1.61. Below about 4 km few rays travel, and the profile there is held by the smoothing and the rule that speed may not decrease, not by data.

The fit also moved the earthquakes. Compare their depths with the catalogue's, which ADLoc computed in a different 1D model.

median depth below sea level: catalogue 1.53 km, after the 1D fit 3.08 km; 10-90 %: 0.4-3.2 km and 2.0-4.5 km

The earthquakes moved about 1.5 km deeper. The fit is better in the new model, but B.7 shows that the data fix absolute depth poorly: a delay common to all S picks trades against depth. Keep in mind that every section from here on draws the earthquakes deeper than the catalogue does.

Now the question that motivates tomography. Put the 1D model and its earthquake locations into the forward model, without the station terms, and look at what is left over: the residual of each pick, observed minus predicted time.

One correction comes first. The eikonal solver works on a grid, and on a 1 km grid its times are too long by about a tenth of a second, more than the pick errors. We compute every time twice, on the model grid and on a grid four times finer, and subtract the difference from the data (invert.grid_correction). Part A could skip this because its data came from the same solver.

grid correction: mean -121 ms, rms 134 ms (1 s, including the first compilation)
residuals in the 1D model: robust standard deviation P 65 ms, S 158 ms

The residuals spread by about 65 ms for P and more than twice that for S. Pick errors alone would give a spread too, but pick errors do not care where the station is. So average the residuals of each station: if a station's mean residual is far from zero, every ray to it is late or early, and that is structure, not noise.

P: median station residual, stations north of 38.80° +6 ms, south -39 ms
S: median station residual, stations north of 38.80° +23 ms, south -99 ms

The pattern is not random. Stations in the south-east of the field see arrivals early, by up to a tenth of a second for P and more for S: the rays there cross faster rock than the 1D model holds. Several stations at the edge of the field, west and east, see them late. Neighbouring stations agree with each other, which pick errors would not do.

A station term would absorb each station's delay, but it cannot say where along the rays the fast or slow rock is. A 3D model can, because rays from many earthquakes cross each place at different angles.

B.3 · The 3D fit

The 3D fit starts from the 1D model and its earthquake locations. Every node of the grid now has its own Vp and Vs, about 13,000 unknowns (B.1), against 23,191 picks. Many nodes are crossed by no ray, so the regularisation of Part A, A.6, decides them. It has three parts, each with a scale (a larger scale is a weaker penalty):

• smoothing of log Vp between neighbouring nodes: ten times stronger sideways than with depth, because the earth changes faster with depth than sideways;
• the same smoothing of Vp/Vs, three times stronger, because Vp/Vs comes from the small difference between S and P times and is noisier than Vp;
• damping of every node toward the 1D model, which holds the nodes no ray reaches.

One difference from Part A: the smoothing penalty grows like the square of a small difference between neighbours but only like the difference itself for a large one (an L1 penalty above a corner l1_eps). A sharp step then costs less than in Part A, so the model can keep a sharp boundary instead of blurring it. The price: on a 1 km grid such a step falls on a row of nodes, and straight edges along the grid rows appear.

The smoothing and damping act on the change from the 1D model, not on the model itself. So the 1D layering costs nothing, and where no ray samples, the model simply stays 1D.

There are no station terms in the 3D fit: the 3D model must explain each station's delay itself. Whether that is right is a real question, and B.4 tests it. The values are tomo.CFG; B.5 shows how they were chosen.

   smooth = 0.99
 smooth_h = 0.099
 smooth_k = 0.297
smooth_kh = 0.0297
     damp = 0.546
   damp_k = 0.546
   l1_eps = 0.01
    delta = 1.0
    niter = 60
   rounds = 3

The full pipeline runs up to 60 Gauss–Newton iterations per round and recomputes the grid correction three times between rounds. We run one round of 15 iterations (a few minutes), then compare with the full result. Each iteration does what Part A, A.8, did: relocate every earthquake in the current model, linearise, take the velocity step with each earthquake's location directions projected out, and keep the step only if the loss drops. The program prints one line per iteration: the loss (misfit plus regularisation, per pick), the misfit alone, the Levenberg–Marquardt damping lam and the largest change of any node.

[fit] 3D eval 0 loss 3.1659e-01 misfit 3.1659e-01

[fit] 3D eval 1 loss 2.6924e-01 misfit 2.6395e-01 lam 3.3e-02 |dv| 8.14e-02 t 9s

[fit] 3D eval 2 loss 2.3697e-01 misfit 2.2531e-01 lam 1.1e-02 |dv| 1.66e-01 t 15s

[fit] 3D eval 3 loss 2.2164e-01 misfit 2.0592e-01 lam 3.7e-03 |dv| 1.71e-01 t 20s

[fit] 3D eval 4 loss 2.1297e-01 misfit 1.9487e-01 lam 1.2e-03 |dv| 9.73e-02 t 26s

[fit] 3D eval 5 loss 2.0843e-01 misfit 1.8947e-01 lam 4.1e-04 |dv| 6.45e-02 t 33s

[fit] 3D eval 6 loss 2.0482e-01 misfit 1.8484e-01 lam 1.4e-04 |dv| 8.29e-02 t 39s

[fit] 3D eval 7 loss 2.0267e-01 misfit 1.8243e-01 lam 4.6e-05 |dv| 4.75e-02 t 45s

[fit] 3D eval 8 loss 2.0070e-01 misfit 1.7991e-01 lam 1.5e-05 |dv| 6.44e-02 t 51s

[fit] 3D eval 9 loss 1.9883e-01 misfit 1.7746e-01 lam 5.1e-06 |dv| 8.52e-02 t 57s

[fit] 3D eval 10 loss 1.9774e-01 misfit 1.7616e-01 lam 1.7e-06 |dv| 3.58e-02 t 63s

[fit] 3D eval 11 loss 1.9684e-01 misfit 1.7503e-01 lam 1.0e-06 |dv| 2.92e-02 t 70s

[fit] 3D eval 12 loss 1.9617e-01 misfit 1.7423e-01 lam 1.0e-06 |dv| 2.68e-02 t 76s

[fit] 3D eval 13 loss 1.9558e-01 misfit 1.7344e-01 lam 1.0e-06 |dv| 4.66e-02 t 82s

[fit] 3D eval 14 loss 1.9500e-01 misfit 1.7265e-01 lam 1.0e-06 |dv| 5.07e-02 t 88s

[fit] 3D eval 15 loss 1.9464e-01 misfit 1.7228e-01 lam 1.0e-06 |dv| 1.60e-02 t 94s

15 iterations in 94 s; loss 0.3166 -> 0.1946

Now the lateral part of the live model, next to the full pipeline's, at 2 km depth and along A–A′. The comparison is only over the nodes the rays sample. The colour scale stops at ±10 %; a few shallow nodes go beyond it.

live vs full, lateral Vp on the sampled nodes: correlation 0.95; rms 5.1 % vs 5.9 %

Fifteen iterations already draw the main pattern: slow rock in the north-west of the field at 1–3 km depth, and fast rock in the south-east. The remaining iterations and the repeated grid correction strengthen it but do not change its shape.

Did the 3D model explain what the 1D model left? Compare the residuals of three models: the 1D model alone (B.2), the 1D model with its station terms, and the full 3D model.

                             robust sd P (ms)  robust sd S (ms)  station medians, sd P (ms)  station medians, sd S (ms)
1D model                                 65.0             158.0                        82.0                       181.0
1D model + 82 station terms              38.0              69.0                         5.0                        30.0
3D model, no station terms               38.0              71.0                        17.0                        56.0

Read the last two rows together. The 3D model fits the picks about as well as the 1D model with one constant delay per station and phase, not better. Both explain the station pattern of B.2.

So the data alone cannot tell two explanations apart:

1. each station sits on its own rock, which delays every ray to it by the same amount: a station term, or equally a thin layer right under the station;
2. the rock between the stations and the earthquakes, down to a few kilometres, varies from place to place: a 3D model.

A station delay is shallow 3D structure too, so these are not opposites. The question is where the slow and fast rock is: right under each station, or spread through the 1–3 km where the rays from many earthquakes cross. Fit cannot answer it. The twin in B.4 can at least say what each choice does when the truth is known.

B.4 · The synthetic twin

B.3 found a 3D model that fits the data much better than the 1D model. Part A taught us not to believe it yet: a regularised inversion draws structure that is not there, smears what is there, and loses amplitude, and how much depends on the geometry of stations and earthquakes. So we build a twin of this data set, with a truth we choose, and run the identical pipeline on it.

A twin is only useful if it is as hard as the real data. This one keeps:

• the same stations, earthquakes and picks: every real pick has a synthetic twin, with the earthquake where the real 1D fit put it;
• the real 1D model as the background of every truth;
• times computed on a grid four times finer than the model grid, so the data carry the solver's grid error and structure between the nodes, as real times do (Part A could not);
• pick errors like the real ones: Gaussian noise, plus a few picks badly wrong and some moderately wrong (about as many as the real residuals show);
• the same program and settings, from the uniform start: 1D fits, cut, 3D fit.

Each truth is run with two random draws of the noise; the scores below are their mean.

           grid: 4x finer than the model
    noise_sigma: 0.4 x (P 0.10 s, S 0.15 s)
      bad_picks: 3 % off by 0.3-1.5 s
 moderately_bad: 8 % off by 0.1-0.5 s, 55 % of them late
     background: the real 1D model
    earthquakes: where the real 1D fit put them
          seeds: 2

Four truths, each asking one question:

truth	what is planted	question
null	the 1D background only	what does the method invent?
anomaly	a slow, low-Vp/Vs body at 1–4 km (a reservoir); a NW–SE gradient; shallow patches under half the stations	are realistic structures recovered, and how strongly?
columns 6 km	±5 % Vp, ±0.04 Vp/Vs in 6 × 6 km columns, all depths	can a 6 km pattern be seen?
columns 8 km	the same in 8 × 8 km columns	and an 8 km one?

Two numbers score each result, both on the lateral part, on the nodes the rays sample, at 1–4 km where the earthquakes are:

• the false structure: the rms of the lateral part of the found model, for the null truth (whose lateral part is zero, so every bit of it is invented);
• the recovery $R^2 = 1 - \sum (\text{found} - \text{true})^2 / \sum \text{true}^2$: 1 is perfect, 0 is no better than the 1D model, and a negative value is worse than drawing nothing.

null: rms false structure (%); the others: R^2
                Vp  Vp/Vs
truth                    
null          0.16   0.11
anomaly       0.65   0.50
columns 6 km  0.35   0.09
columns 8 km  0.44   0.29

Now look at them, not only at their scores. Each row is one truth: the planted lateral Vp and what the pipeline found, on the A–A′ section and on the 2 km map.

Read the rows one at a time.

• Null. The found model is nearly blank: what the method invents is a fraction of a per cent, far below the real model's few per cent in B.3.
• Anomaly. The NW–SE gradient is recovered: slow in the north-west, fast in the south-east. The shallow patches under the stations come out as the strongest found anomalies, smeared down to about 1 km. The reservoir at 1–4 km, where the earthquakes are, appears only faintly: rays that start inside a body cross little of it, and a slow body also moves the earthquakes in it.
• Columns. Both patterns are found in the right places where the earthquakes are, but weaker; the 8 km columns more completely than the 6 km ones. Outside the cloud of earthquakes, few rays cross and the model stays at the background.

The same for Vp/Vs, which is harder: it comes from the difference between S and P times.

Vp/Vs is recovered less well than Vp at every scale: the reservoir's low Vp/Vs appears, weak and broad; the 6 km columns hardly at all. In words a newcomer can carry to any tomography paper:

• lateral Vp patterns of 8 km can be seen here, 6 km ones only partly; Vp/Vs patterns only from about 8 km, and only where the earthquakes are; nothing smaller was tested;
• small anomalies come back weaker than they are: a found anomaly of 3 % may be a true one of 5 % (B.6 shows that broad ones lose less);
• anything near the size of the false structure in the null test means nothing.

Station terms or 3D structure?

B.3 left a question: the 3D model and the 1D model with station terms fit the picks equally well. What if the 3D fit also had station terms? Then each station's delay could go into its term instead of into the model. The twin tells us what that does when the truth is known. These runs used a damping scale of 0.182 (three times stronger than the chosen one); both rows of each pair used the same, so the comparison is fair.

                           null false Vp (%)  null false Vp/Vs (%)  anomaly R2 Vp  anomaly R2 Vp/Vs
3D, no station terms                    0.15                  0.10           0.59              0.47
3D + fitted station terms               0.12                  0.08           0.31              0.38

In the twin, fitted terms cut the recovery of the anomaly: the truth's shallow structure, which a station term can mimic, goes into the terms, and the model loses it. They also lower the false structure a little. In the twin, then, the choice costs recovery and buys little.

Now the real data, both ways.

Vp: with terms, the lateral pattern correlates 0.61 with the one without, at 22 % of its amplitude
Vp/Vs: with terms, the lateral pattern correlates 0.29 with the one without, at 9 % of its amplitude
real_terms3d_none: robust sd P 39 ms, S 74 ms
real_terms3d_fit: robust sd P 34 ms, S 62 ms

On the real data the terms take much more. The Vp pattern keeps its shape but only about a fifth of its size; the Vp/Vs pattern almost disappears. The fit improves, as it must with 82 more unknowns.

So the real lateral pattern, and above all the real Vp/Vs pattern, is in large part interchangeable with a delay at each station. Leaving the terms out is the better choice in the twin, where the truth is known to be structure; for the real Earth it is an assumption: that the stations' delays come from rock the rays of many earthquakes cross, not from the top few hundred metres under each station alone. Keep this in mind when B.6 reads the real model.

B.5 · How strong should the regularisation be?

Part A, A.6, showed the trade-off on a planted truth: weaker regularisation recovers more and invents more. Part A chose the strength so that the data were fitted to their noise level. Here the noise level itself is uncertain (B.3: bad picks, and delays a model may or may not hold), so we use three referees:

• the twin: on the null truth the false structure must stay small, and on the anomaly the recovery should be as high as possible;
• the real data split in two: run the whole pipeline on the even-numbered earthquakes, then on the odd-numbered ones (all earthquakes are still located; only the velocity fit is split). Two independent halves of the data should give the same model. Structure that one half sees and the other does not was fitted to noise or to a few earthquakes;
• prediction: the half of the earthquakes left out of the velocity fit is still located in the model; how well the model predicts their picks is a test on data it was not fitted to.

Four settings, each the full pipeline (1D fits, cut, 3D fit) on the twin (two noise draws) and on both real halves:

A scale 9 times smaller is a squared penalty 81 times stronger.

setting	damping scale	lateral smoothing scale
damping ×9 stronger	0.061	as chosen
as chosen	0.546	0.099 (Vp), 0.030 (Vp/Vs)
lateral smoothing ×3 stronger	0.546	0.033, 0.010
no damping	none	as chosen

                       null false Vp (%)  null false Vp/Vs (%)  anomaly R2 Vp  anomaly R2 Vp/Vs  real halves agree Vp  real halves agree Vp/Vs  held-out P (ms)  held-out S (ms)
setting                                                                                                                                                                         
damping x9 stronger                 0.11                  0.07           0.42              0.38                  0.97                     0.98            51.14            93.20
as chosen                           0.16                  0.11           0.65              0.50                  0.98                     0.97            41.75            80.94
smoothing x3 stronger               0.10                  0.06           0.54              0.32                  0.99                     0.99            47.70            91.92
no damping                          0.17                  0.11           0.67              0.52                  0.96                     0.96            41.17            79.75

Read the columns from left to right.

• Stronger damping or stronger smoothing invents less, but recovers clearly less of the anomaly and predicts the left-out earthquakes clearly worse: it holds real structure back toward the 1D model.
• No damping recovers the anomaly and predicts the left-out picks about as well as the chosen setting (within about 1 ms), invents as much, and its real halves agree a little less.
• So recovery and prediction improve as the regularisation weakens, up to the chosen setting, and then stop improving.

The rule: among the settings that predict and recover as well as the best, take the strongest regularisation. That is the chosen setting. It is Part A's rule in another form: do not let the model fit more than the data can support. It is a judgement made with numbers, not a theorem: a study that cared more about false structure than about recovery would choose a stronger one.

What the settings do to the picture is clearest on a section. The top row is the anomaly twin, the bottom row the real data (all earthquakes), for the three settings that differ.

Stronger damping lowers the amplitude; stronger lateral smoothing spreads each anomaly sideways. The real picture changes with the setting in the same way as the twin's, which is why the twin can choose the setting. And the main real features, slow in the north-west and fast in the south-east, are there in every setting.

B.6 · The real 3D model

The full pipeline's model: maps at four depths (Vp above, Vp/Vs below), then both sections. Earthquakes are drawn where the 3D fit located them.

lateral rms at 1-4 km on the sampled nodes: Vp 5.1 %, Vp/Vs 1.1 %  (null test: 0.16 %, 0.11 %)
lateral Vp at 0 km, 2nd to 98th percentile of the sampled nodes: -15 % to +12 %
lateral Vp at 1 km, 2nd to 98th percentile of the sampled nodes: -10 % to +15 %
lateral Vp at 2 km, 2nd to 98th percentile of the sampled nodes: -7 % to +12 %
lateral Vp at 3 km, 2nd to 98th percentile of the sampled nodes: -4 % to +4 %

Three features stand out:

1. Vp is low in the north-west of the field and high in the south-east at 1–3 km, by about −7 to +12 % at 2 km. The change is sharp, along a straight east–west edge near 38.80° N. That edge follows a row of the 1 km grid, and the L1 smoothing (B.3) favours such steps, so its exact position and sharpness are not resolved; only that the change happens there.
2. Vp/Vs is high in the north-west and low in the south-east.
3. Near the surface (0–1 km) there are strong blocks under single stations, and a very slow corner north-east of B′, outside the earthquakes, beyond the ±10 % of the colour scale.

The lateral rms is about thirty times the false structure of the null test for Vp and ten times for Vp/Vs, so features 1 and 2 are not inventions of the method in the sense of B.4. That is necessary, not sufficient. Two more questions: how does the method render features like these, and are they stable?

The first has a direct test. Take the real model itself as the truth of the twin, make synthetic picks from it exactly as in B.4, and run the pipeline again. The difference between this truth and what comes back is what the method does to features of this shape and size.

the real model as truth: R^2 Vp 0.91, Vp/Vs 0.78

The pattern comes back: the large NW–SE contrast in Vp almost entirely, the Vp/Vs pattern in its shape but weaker and broader. Broad features like these lose much less amplitude than the smaller ones of B.4. So the method can carry features of this size from the data into the model. It does not prove the earth has them; it says that if it does, this is what we would see.

B.7 · Is the result stable?

A feature that a few earthquakes create would change when those are left out. B.5 split the earthquakes into even and odd halves. Here are both, next to each other, at 2 km.

halves agree (correlation at 1-4 km): {'Vp': np.float64(0.98), 'Vp/Vs': np.float64(0.97)}

Even and odd are interleaved in time and space, so each half still covers the whole field. Two harder tests:

• ten random halves, each run through the whole pipeline: at every node, the spread of the ten models measures how much the answer depends on which earthquakes were used;
• north-west against south-east: only the earthquakes on one side of the field's centre shape the velocity. Each side's rays sample its own side best, so this asks whether a feature is seen from two different directions.

nodes at 1-4 km with |lateral Vp| > 1 %: 789; anomaly larger than twice the spread of the random halves at 98 % of them
random halves agree with each other: correlation 0.99 on average, 0.98 at worst
north-west vs south-east earthquakes, Vp: correlation 0.78
north-west vs south-east earthquakes, Vp/Vs: correlation 0.70

The anomaly is many times the spread almost everywhere it is large; the right map is dark only where the anomaly itself is near zero, at the edges between slow and fast. The north-west and south-east earthquakes agree less than random halves do, as they should: each side sees its own part of the field best.

Stable is not the same as correct. A station delay drawn into the model would be stable too: every half of the earthquakes is recorded at the same stations. B.4 showed that this is the one alternative the tests here cannot rule out.

One more limit: absolute depth. B.2 found the earthquakes 1.5 km deeper than the catalogue. A delay common to all S picks (for example from how S is picked) looks to the data like a slower model with deeper, earlier earthquakes (Part A, A.8; Session 2, Location & relocation). The test: shift every S pick by a fixed amount and run the pipeline again (these runs used an earlier setting: damping 0.182, one grid-correction round, one start).

                        median depth (km)  Vp at 0 km (km/s)  1D loss  3D loss
S picks earlier by (s)                                                        
-0.1                                3.542              4.030    0.178    0.198
-0.05                               3.315              4.048    0.177    0.198
0                                   3.142              4.276    0.175    0.201
0.05                                2.944              4.476    0.175    0.202
0.1                                 2.637              4.657    0.176    0.207
0.2                                 2.087              5.047    0.181    0.217

Over a shift of ±0.1 s the loss changes by only a few per cent, while the median depth moves by about a kilometre and the shallow Vp by more than half a km/s. The data do not fix these well. So the absolute depth of the earthquakes and the absolute speed near the surface are less certain than the lateral changes this session is about.

What may be read from the 2016 Geysers model

Putting B.4–B.7 together:

• Supported, with one condition: lateral Vp at 1–3 km, slow in the north-west and fast in the south-east. It is many times the false structure, comes back when used as a planted truth, and is the same from random halves of the earthquakes. The condition: that the stations' delays come from the rock the rays cross, not only from the ground right under each station; with fitted station terms it shrinks to about a fifth (B.4).
• Weaker: the Vp/Vs pattern, high in the north-west and low in the south-east. It passes the same tests, but fitted station terms remove almost all of it.
• Not reliable: the blocks at 0–1 km under single stations and the north-east corner. Only rays to one or a few stations sample them.
• Not resolved: anything below about 4 km (no earthquakes, so no rays); lateral patterns smaller than about 8 km; the exact position of sharp edges; absolute depth and shallow speed.

How does this compare with published work? Julian, Ross, Foulger & Evans (1996) imaged The Geysers with travel-time tomography and found a strong low in Vp/Vs, about −9 %, over the most intensively exploited part of the reservoir, which they note is not evident in Vp alone. They read it as low pore pressure and dry rock, partly from pore water boiling as steam is extracted. Our south-east Vp/Vs low is much weaker, and it sits on our strongest Vp high, so it is not simply the same feature. Their data are from the early 1990s and ours from 2016; their reference is a different model; and B.4 showed that our Vp/Vs pattern is the part most easily traded for station delays. Whether the two lows coincide with the most exploited part of the field is a question for a project; the project notebook has the production and injection records to start from.

Takeaways

Machine learning

• The adjoint of the eikonal solver is reverse-mode differentiation: one extra solve per station gives the gradient of the misfit for every node, which is what lets a 3D model of thousands of unknowns be fitted at all.
• When the data cannot decide every unknown, the regularisation decides the rest. Its strength is a hyperparameter, chosen by a planted twin of the real problem, by whether independent halves of the real data agree, and by how well the model predicts earthquakes it was not fitted to.
• A model that fits better is not a better model: here a 3D model and a 1D model with one delay per station fit the picks equally well, and only a twin with a known truth can say what each choice recovers.

Seismology

• Tomography of a local network sees structure only where rays cross from many directions: here, among the earthquakes, between about 1 and 4 km depth, and at sizes of about 8 km and larger.
• Every tomographic image should come with its null test, its recovery test and a stability test, and should say how its result changes when station delays are fitted instead.
• The 2016 Geysers picks show slower rock in the north-west of the field and faster rock in the south-east at 1–3 km depth, with high and low Vp/Vs respectively; the Vp pattern is the firmer of the two.