The water-level parameter nobody tunes

There is one number in the receiver-function pipeline that sets the resolution of every Moho pick you will ever make from a station, and almost nobody tunes it. It is the water level: a single scalar, usually inherited verbatim from a script written by someone who inherited it from someone else. Its default is 0.01, or 0.001, or whatever the last paper used. It is a bias-variance dial, and across most of the community it is set with a shrug.

What the water level actually does

Spectral water-level deconvolution recovers the crustal impulse response by dividing the radial-component spectrum by the vertical-component spectrum. The vertical component is your best estimate of the incident P wavelet plus source and path effects; the radial component is that same wavelet reverberating through the crust beneath the station. The division looks like this:

$$\hat{H}(\omega) = \frac{R(\omega)\, V^{*}(\omega)}{\max\!\left(\lvert V(\omega)\rvert^{2},\; c \cdot \max_{\omega}\lvert V(\omega)\rvert^{2}\right)}$$

FormulaLaTeX

Copy LaTeX

Wherever the vertical spectrum has a null, a frequency where the incident wavelet carried no energy, that division blows up. The fix, introduced for teleseismic body waves by Clayton and Wiggins [1] and made standard by Langston [2], is to clamp the denominator so it never drops below a fraction $c$ of its peak value. That fraction is the water level. It is not a physical quantity. It is a regularization strength, and like every regularization strength it lives on a tradeoff curve.

Small c: a sharp pulse on a bed of ringing

Set $c$ very small and you barely regularize at all. The Ps conversion comes back sharp and full-amplitude, which would be ideal if the pulse were the whole trace. It is not. The noisy spectral nulls now divide through almost unclamped, and each one inverse-transforms into an oscillation smeared across the entire receiver function. This is variance in the textbook sense: run the same deconvolution on ten events from the same source region and the ringing is different every time, because it is driven by noise, not structure.

Stacking helps less than you would hope. Random ringing partly averages out, but only partly, and what survives can align into a coherent side-lobe that reads as a mid-crustal discontinuity that is not there. A picker, human or automated, will happily pick it.

Large c: a quiet trace with the wrong amplitude

Turn $c$ up and the ringing vanishes. The trace looks clean, professional, publishable. You have paid for it in bias. A large water level over-regularizes the whole spectrum, broadening the Ps pulse and biasing its amplitude low. Receiver-function amplitude is not decorative. It is the observable that constrains the velocity contrast across the interface, as Ammon showed when he pinned down the absolute amplitude content of P waveforms [3], and in H-k stacking it feeds the weighted sum over the converted phase and its multiples [6]. Bias the amplitude and you bias the impedance contrast, the composition inference, and the confidence you attach to the pick. The trace is quiet because you have thrown away signal along with the noise.

The curve has a minimum, and it moves

Total error is the sum of variance and squared bias, and that sum has a minimum somewhere between the ringing regime and the over-smoothed regime. The instrument below makes the tradeoff draggable. The left panel plots variance, bias, and their sum against $\log_{10} c$; the right panel shows the recovered Ps pulse morphing from sharp-and-ringing to broad-and-quiet as you sweep the dial. The hero number tracks how much of the true Ps amplitude survives.

The second lever is the point of this piece. Drag the station noise up and the amber optimum slides roughly a decade along the c axis, while the grey "script default" marker stays exactly where the script put it. The sweet spot is not universal. It depends on the station's noise level, the source spectra of the events in your catalog, and the frequency band of the target conversion. A single hard-coded $c$ is close to right for one noise regime and wrong everywhere else, and a continental deployment spans every noise regime you can name.

What this does to a network-scale pipeline

The failure is quiet because both failure modes look plausible. The noisy-station traces ring, the stack inherits coherent side-lobes, and the map gains structure. The quiet-station traces are over-smoothed, the Ps amplitudes drop, and the inferred contrasts flatten. Neither failure announces itself in QC the way a dead channel does. A pipeline that hard-codes one water level for a thousand stations is optimizing for the average station and quietly corrupting the tails on both ends.

Hand-tuning is not the answer either. Nobody has time to grid-search a regularization parameter by eye across a full network, re-checked every time the catalog grows. This is also why the method choice itself matters: multitaper deconvolution averages independent spectral estimates and degrades more gracefully [5], and iterative time-domain deconvolution builds the answer additively and never divides by a null at all [4]. But wherever spectral division is the workhorse, the water level is the knob that decides what your picks are worth.

The durable fix is per-trace, data-driven regularization

The water level was always a proxy for a question the data can answer directly: how much do I trust this trace? A regularization that conditions on each trace's own vertical-component spectrum and noise estimate tightens automatically on clean events and backs off on noisy ones. That is what a learned deconvolution operator does, and it is why EarthScan treats deconvolution as a trainable component rather than a script with a magic number. On clean records it behaves like a light-touch spectral division; on noisy records it supplies the damping the trace actually needs, per event, without a human sweeping a dial.

By the numbers

Regime	Ps pulse	Amplitude fidelity	Stack risk
c too small	sharp, full amplitude	high	noise ringing stacks into fake structure
c at the station optimum	near-true width	high	minimal
c too large	broadened, quiet	biased low	contrasts and H-k weights flattened
one default c, whole network	varies by station	varies	both failure modes, silently

Regime descriptions summarize the schematic tradeoff in the instrument above, not a measured survey result.

Key takeaways

1. 1The water level is a bias-variance regularization dial, not a physical constant, and it sets the effective resolution of every pick from that station.
2. 2Small c gives a sharp pulse riding on noise-driven ringing that can stack into coherent, fake mid-crustal structure.
3. 3Large c biases the Ps amplitude low, corrupting the velocity-contrast and composition inferences that depend on it.
4. 4Total error has a station-specific minimum that moves with noise, source spectrum, and target frequency, so a single default c is wrong almost everywhere.
5. 5The durable fix is a learned, per-trace regularization conditioned on each trace's own spectrum and noise estimate, not a hand-tuned scalar.

Limitations

The curves and pulse shapes in the instrument are schematic illustrations of the standard bias-variance decomposition, not error bars measured from a re-run deconvolution; the roughly one-decade shift of the optimum with noise is a property of the illustrative response surfaces, and the shift at a real station must be measured from its own catalog. The decomposition also treats noise as stationary within an event, which field noise is not. Amplitude bias interacts with stacking and with the H-k weighting in ways the single-trace view here does not resolve. And where the deconvolution method itself is the free choice, iterative time-domain and multitaper estimators change the tradeoff rather than merely moving along it, so the water-level analysis applies to spectral-division pipelines specifically.

References

[1] Clayton, R. W., and Wiggins, R. A. (1976). Source shape estimation and deconvolution of teleseismic bodywaves. Geophysical Journal of the Royal Astronomical Society, 47, 151-177. doi:10.1111/j.1365-246X.1976.tb01267.x

[2] Langston, C. A. (1979). Structure under Mount Rainier, Washington, inferred from teleseismic body waves. Journal of Geophysical Research, 84(B9), 4749-4762. doi:10.1029/JB084iB09p04749

[3] Ammon, C. J. (1991). The isolation of receiver effects from teleseismic P waveforms. Bulletin of the Seismological Society of America, 81(6), 2504-2510. doi:10.1785/BSSA0810062504

[4] Ligorria, J. P., and Ammon, C. J. (1999). Iterative deconvolution and receiver-function estimation. Bulletin of the Seismological Society of America, 89(5), 1395-1400. doi:10.1785/BSSA0890051395

[5] Park, J., and Levin, V. (2000). Receiver functions from multiple-taper spectral correlation estimates. Bulletin of the Seismological Society of America, 90(6), 1507-1520. doi:10.1785/0119990122

[6] Zhu, L., and Kanamori, H. (2000). Moho depth variation in southern California from teleseismic receiver functions. Journal of Geophysical Research, 105(B2), 2969-2980. doi:10.1029/1999JB900322