diff --git a/literature/figs/demo_10_deviations.png b/literature/figs/demo_10_deviations.png index ee29663..1907e8e 100644 Binary files a/literature/figs/demo_10_deviations.png and b/literature/figs/demo_10_deviations.png differ diff --git a/literature/figs/demo_11_multiverse.png b/literature/figs/demo_11_multiverse.png index efb8324..8638dc2 100644 Binary files a/literature/figs/demo_11_multiverse.png and b/literature/figs/demo_11_multiverse.png differ diff --git a/literature/figs/demo_12_bayes.png b/literature/figs/demo_12_bayes.png index 786eb55..0b34186 100644 Binary files a/literature/figs/demo_12_bayes.png and b/literature/figs/demo_12_bayes.png differ diff --git a/literature/figs/demo_13_branch_asymmetry.png b/literature/figs/demo_13_branch_asymmetry.png index bb7779d..12dd942 100644 Binary files a/literature/figs/demo_13_branch_asymmetry.png and b/literature/figs/demo_13_branch_asymmetry.png differ diff --git a/literature/figs/demo_14_network_pairs.png b/literature/figs/demo_14_network_pairs.png index 89c0256..86aec58 100644 Binary files a/literature/figs/demo_14_network_pairs.png and b/literature/figs/demo_14_network_pairs.png differ diff --git a/literature/figs/demo_15_window_envelope.png b/literature/figs/demo_15_window_envelope.png index b84a961..cb7987a 100644 Binary files a/literature/figs/demo_15_window_envelope.png and b/literature/figs/demo_15_window_envelope.png differ diff --git a/literature/figs/demo_1_methods.png b/literature/figs/demo_1_methods.png index 90ee1a3..2eb5a2c 100644 Binary files a/literature/figs/demo_1_methods.png and b/literature/figs/demo_1_methods.png differ diff --git a/literature/figs/demo_2_aggregation.png b/literature/figs/demo_2_aggregation.png index 35bb134..2105dbb 100644 Binary files a/literature/figs/demo_2_aggregation.png and b/literature/figs/demo_2_aggregation.png differ diff --git a/literature/figs/demo_3_uncertainty.png b/literature/figs/demo_3_uncertainty.png index 8ea193c..bd05e7b 100644 Binary files a/literature/figs/demo_3_uncertainty.png and b/literature/figs/demo_3_uncertainty.png differ diff --git a/literature/figs/demo_4_frequency_depth.png b/literature/figs/demo_4_frequency_depth.png index ce8f421..ec3b1c1 100644 Binary files a/literature/figs/demo_4_frequency_depth.png and b/literature/figs/demo_4_frequency_depth.png differ diff --git a/literature/figs/demo_5_window_band.png b/literature/figs/demo_5_window_band.png index d90fb58..58e37f0 100644 Binary files a/literature/figs/demo_5_window_band.png and b/literature/figs/demo_5_window_band.png differ diff --git a/literature/figs/demo_6_stacking.png b/literature/figs/demo_6_stacking.png index 9672151..cefa721 100644 Binary files a/literature/figs/demo_6_stacking.png and b/literature/figs/demo_6_stacking.png differ diff --git a/literature/figs/demo_7_reference.png b/literature/figs/demo_7_reference.png index 5884bfb..7fc7638 100644 Binary files a/literature/figs/demo_7_reference.png and b/literature/figs/demo_7_reference.png differ diff --git a/literature/figs/demo_8_artifacts.png b/literature/figs/demo_8_artifacts.png index a47f949..a3b23ed 100644 Binary files a/literature/figs/demo_8_artifacts.png and b/literature/figs/demo_8_artifacts.png differ diff --git a/literature/figs/demo_9_multiverse.png b/literature/figs/demo_9_multiverse.png index a8d0d01..2726374 100644 Binary files a/literature/figs/demo_9_multiverse.png and b/literature/figs/demo_9_multiverse.png differ diff --git a/literature/figs/realdata_1_validation.png b/literature/figs/realdata_1_validation.png new file mode 100644 index 0000000..cd2c753 Binary files /dev/null and b/literature/figs/realdata_1_validation.png differ diff --git a/literature/figs/realdata_2_interferograms.png b/literature/figs/realdata_2_interferograms.png new file mode 100644 index 0000000..53c13ad Binary files /dev/null and b/literature/figs/realdata_2_interferograms.png differ diff --git a/literature/figs/realdata_3_warmup.png b/literature/figs/realdata_3_warmup.png new file mode 100644 index 0000000..500b7f6 Binary files /dev/null and b/literature/figs/realdata_3_warmup.png differ diff --git a/paper/manuscript_marine.qmd b/paper/manuscript_marine.qmd index c223e58..5a49842 100644 --- a/paper/manuscript_marine.qmd +++ b/paper/manuscript_marine.qmd @@ -53,12 +53,12 @@ Time-series analysis; Interferometry; Seismic noise; Coda waves; Inverse theory; # Introduction {#sec:intro} -Changes in subsurface properties occur due geodynamics, which drive earthquake damage -and volcanic eruption, and hydrodynamics, which controls fluid exchange between the atomsphere -and the solid Earth. These processes influence the mechanical property of Earth materials, -which directly affect the speed at which seismic waves propagates. Changes in seismic -velocity, often measured and referred to as \dvv\, can be tracked by measuring changes in arrival -times of seismic waves, especially scatterd waves such as coda waves, +Changes in subsurface properties occur due to geodynamics, which drive earthquake damage +and volcanic eruption, and hydrodynamics, which controls fluid exchange between the atmosphere +and the solid Earth. These processes influence the mechanical properties of Earth materials, +which directly affect the speed at which seismic waves propagate. Changes in seismic +velocity, often measured and referred to as \dvv$\,$can be tracked by measuring changes in arrival +times of seismic waves, especially scattered waves such as coda waves, provided that the source and receivers are at the same location. ```{=latex} @@ -73,13 +73,13 @@ the current coda; the exact relation is \] The first-order shortcut $\dvv \approx -\varepsilon$ is accurate to order $\varepsilon^2$ --- negligible below $1\,\%$ but a $0.17\,\%$ absolute bias at -the $4\,\%$ changes seen on landslides. +a true $\dvv$ of $4\,\%$, the changes seen on landslides. \end{minipage}} \end{center} ``` Due to the sensitivity of coda waves to small perturbations in the material properties, -\dvv\ was discovered as an effective method to monitor chanes during volcanic unrest since its +\dvv\ was discovered as an effective method to monitor changes during volcanic unrest since its discovery at Piton de la Fournaise [@Brenguier2008] and is now calculated in continuous along side of more conventional seismic monitoring methods in the same volcano observatory [@Duputel2009] and was a determining early warning parameter used by the Icelandic Meteorological @@ -93,39 +93,39 @@ along side of more conventional seismic monitoring methods in the same volcano o Each of these deployments rests on the same fragile assumption: that the \dvv\ curve an operator acts on is a property of the subsurface, not of the analyst's processing choices. -The elevated sensitivity comes at the price a long series of processing choices +The elevated sensitivity comes at the price of a long series of processing choices, and at almost every step the analyst makes a choice. Among these are choices of - estimators between windowed phase measurements or strething [@Mikesell2015, @Mao2020, @Yuan2021], + estimators between windowed phase measurements or stretching [@Mikesell2015, @Mao2020, @Yuan2021], frequency band and coda window (which together set the sampled depth; @Obermann2013 [@Obermann2016]), reference window [@Brenguier2014, @Ermert2023, @Okubo2024], how much to substack to increase coherence among windows (Hadzianou Celine), and how to aggregate and weight the many cross-component and station-pair measurements that make up a single reported \dvv\ time series (e.g., @Hobiger2012). - These choices are typically made by habit, justified briefly if at all , and rarely reported in + These choices are made by habit, justified briefly if at all, and rarely reported in enough detail to reproduce. The community has long flagged *individual* pitfalls (spurious changes from non-stationary noise, @Zhan2013; measurement-error formulae, @Clarke2011 [@Weaver2011]), but the *cumulative* effect of the full choice set on both the value and its stated uncertainty has not been quantified -comprehensively in the literature. +in the literature. This is a reproducibility problem of exactly the "garden of forking paths" type identified in the statistical sciences [@Gelman2013; @Steegen2016]: many individually reasonable analyses of the same data give different answers, and -without full reporting parameter choices we cannot fully interpret \dvv\ as robust values. -. Here we make the problem concrete for \dvv\ monitoring. We use purely synthetic correlation +without full reporting of parameter choices we cannot interpret a reported +\dvv\ as a robust value. Here we make the problem concrete for \dvv\ monitoring. We use purely synthetic correlation time series to test the methods (estimators) and parameter choices that the community -makes to estimate \dvv\, which we report over 103 studies Appendix~\ref{app:survey}. +makes to estimate \dvv$\,$which we report over 103 studies in Appendix~\ref{app:survey}. We do not aim to report the "best" pipeline, which is most often the one reported in -scientific papers, but instead document the various parameter impacts (Section~\ref{sec:results}). -We then propose a new measurement error that incoporates these effects into a data +scientific papers, but instead document the parameter impacts (Section~\ref{sec:results}). +We then propose a new measurement error that incorporates these effects into a data covariance matrix $C_d$ (Section~\ref{sec:bayes}). -One example of propagating such error in downstream scientific insights is the +One example of propagating such error into downstream science is the migration of the surface \dvv\ measurement to depth profiles of perturbations in shear wave velocity $\Delta V_S(z)/V_S(z)$, which often depends on the wavefield constituting the coda waves, such as surface waves or body waves, and that depend on the source-receiver pair geometry. We illustrate the propagation of errors to a -depth profile (Section~\ref{sec:depth}). We only utilize synthetic examples for +depth profile (Section~\ref{sec:depth}). We only use synthetic examples for ground truthing on the signal processing parameters, since the concepts behind the observations of phase lags in scattered waves is well established [@Obermann2013]. We package this new methodology in a Python software, ``codameter``, which we also @@ -170,7 +170,7 @@ reproduces the observation that high frequencies are retained only at short lag times, so a fixed late window samples different depths at different bands. A homogeneous velocity change is imposed exactly by stretching the lapse-time -axis, $u_{\rm cur}(t) = u_{\rm ref}\!\left(t/(1+\dvv)\right)$, and a repeated time +axis, $u_{\rm cur}(t) = u_{\rm ref}\!\left(t\,(1+\dvv)\right)$, and a repeated time series is produced by generating this stretched coda with a prescribed ground-truth $\dvv(t)$ and additive band-limited noise at a controlled signal-to-noise ratio. The concept has been demonstrated using full waveform @@ -179,13 +179,13 @@ we do not repeat that full-waveform modeling here. Because the imposed $\dvv(t)$ a recovered series from it is an artefact of the processing, not of the data. -We utilize seven \dvv\ estimators that were implemented in ``noisepy`` in the +We use seven \dvv\ estimators that were implemented in ``noisepy`` in the `monitoring_methods` module [@Jiang2020]: trace stretching (TS, @lobkis01), windowed cross-correlation (WCC, @poupinet84), dynamic time warping (DTW, @Mikesell2015), -te moving-window cross-spectrum (MWCS; @Clarke2011), and three wavelet-domain methods, the +the moving-window cross-spectrum (MWCS; @Clarke2011), and three wavelet-domain methods, the wavelet cross-spectrum (WCS, @Mao2020) and the two wavelet stretching (WTS) and wavelet DTW (WTDTW) -introduced and benchmarked numerically by @Yuan2021. The framework, the figures below, and a - implementation are released in the open ``codameter`` package +introduced and benchmarked numerically by @Yuan2021. The framework, the figures below, and an +implementation are released in the open ``codameter`` package (Section~\ref{sec:discussion}). Formal definitions of all seven methods, of the two aggregation pathways, and of the uncertainty conventions are described in Appendices~\ref{app:estimators} and~\ref{app:aggregation}. @@ -260,7 +260,7 @@ impacts the recovered \dvv\ and its error. The deliverable of this section is a measurement covariance, which every subsequent step in the inference chain (Sections~\ref{sec:depth}--\ref{sec:stress}) consumes. -Each of the parametric components is implemented in a single comprehensive package ``codameter`` +Each of the parametric components is implemented in a single package, ``codameter``, that borrows from ``msnoise`` [@Lecocq2014] and ``noisepy`` [@Jiang2020] to extract the measurement of \dvv\ from the ambient noise monitoring workflow and generalize it to any coda wave from repeated source-receiver paths. @@ -283,8 +283,8 @@ its own dedicated synthetic scenario (see the cross-referenced subsection).} \toprule \textbf{Axis} & \textbf{Best-case RMS} & \textbf{Worst-case RMS} & \textbf{Section} \\ \midrule -Estimator (family split) & $<\!0.01\,\%$ (TS/WTS, up to 5\,\% true \dvv) & -cycle-skip $>\!0.5\,\%$ past $\sim\!1.4\,\%$ true \dvv\ (MWCS) & +Estimator (family split) & $<\!0.01\,\%$ (TS, up to 5\,\% true \dvv) & +cycle-skip $>\!0.5\,\%$ past $\sim\!1.5\,\%$ true \dvv\ (MWCS) & Section~\ref{sec:methods-fig} \\ \hline Cross-component aggregation & $\sim\!0.03\,\%$ (Approach B, averaged images) & $\sim\!0.31\,\%$ (Approach A, unweighted) & Section~\ref{sec:aggregation} \\ \hline @@ -293,8 +293,8 @@ $\sim\!0.05\,\%$ (individual-pair range) & Section~\ref{sec:uncertainty} \\ \hli Frequency band & $\sim\!0.003$--$0.03\,\%$ (matched to depth) & $\sim\!0.10\,\%$ (mismatched) & Section~\ref{sec:params} \\ \hline Coda window & $\sim\!0.01\,\%$ (adapted to band) & -$\sim\!3.8\,\%$ (fixed, wrong band) & Section~\ref{sec:window} \\ \hline -Reference scheme & $\sim\!0.03$--$0.04\,\%$ (fixed stack / joint inversion) & +$\sim\!3.9\,\%$ (fixed, wrong band) & Section~\ref{sec:window} \\ \hline +Reference scheme & $\sim\!0.03\,\%$ (fixed stack / joint inversion) & $\sim\!0.16\,\%$ (60-day moving) & Section~\ref{sec:params} \\ \hline Stack length & $\sim\!0.018\,\%$ (10-day) & $\sim\!0.044\,\%$ (1-day, noisy) & Section~\ref{sec:params} \\ @@ -313,7 +313,7 @@ every estimator to this convention in both signs, and checks the full `run_pipeline` call path end to end, so the convention cannot silently drift back. -On small, clean \dvv\, all seven estimators agree, demonstrating robustness of the methods +On small, clean \dvv$\,$all seven estimators agree (Fig.~\ref{fig:methods}a). Sweeping the same clean recovery out to $\pm 5\,\%$ (Fig.~\ref{fig:methods}b) shows exactly where and how each family first departs from the 1:1 line, and the estimator choice becomes consequential @@ -322,8 +322,8 @@ is split according to their phase measurement approaches. The stretching family coda and remain accurate for high SNR coda waves; the phase methods (MWCS) read a wrapped phase and may cycle-skips, while the *same* cross-wavelet phase (WCS), once unwrapped in 2-D, recovers the change. The warping methods (DTW, WTDTW) track but under-shoot the largest -strains. No estimator is simply "right" and choices in the estimators -can alters the \dvv\ measurements for larger strain changes. +strains. No estimator is simply "right"; the choice of estimator materially changes +the \dvv\ measurement at larger strain. \begin{figure} @@ -331,33 +331,42 @@ can alters the \dvv\ measurements for larger strain changes. \includegraphics[width=\textwidth]{demo_1_methods.png} \caption{Estimator choice across the seven NoisePy methods. (a) Clean, small \dvv: all agree. (b) The same clean recovery swept over $\pm 5\,\%$ - true \dvv: MWCS cycle-skips almost immediately past $\pm 1$--$2\,\%$, DTW/WTDTW - break past $\pm 3$--$4\,\%$, WCS degrades smoothly, while TS/WTS/WCC track the - 1:1 line throughout. (c) Large, noisy \dvv\ (a pre-failure landslide signal): + true \dvv: MWCS cycle-skips past $\sim\!1.5\,\%$ on either branch; TS and WTS + track the 1:1 line throughout; WCC tracks just as tightly on the negative + branch but breaks sharply near the positive edge; DTW and WTDTW break + asymmetrically, WTDTW near $+1\,\%$ but only near $-3\,\%$ on the other + branch; WCS degrades smoothly, crossing $1\,\%$ error beyond $\pm 4\,\%$. + (c) Large, noisy \dvv\ (a pre-failure landslide signal): MWCS cycle-skips, 2-D-unwrapped WCS and the stretching family track, the warping methods under-shoot.} \label{fig:methods} \end{figure} -On a clean, noiseless sweep of true \dvv\ from 0 to 5\,\% (Fig.~\ref{fig:methods}b), -the phase-wrapped MWCS estimator is the first to break: its error stays below -$0.1\,\%$ up to a true \dvv\ of $\sim\!1.2\,\%$, then jumps discontinuously past -$0.5\,\%$ error by $\sim\!1.4\,\%$ true \dvv\ --- the cycle-skip. The stretching -family (TS, WTS) stays below $0.01\,\%$ error out to the full 5\,\% tested, and -WCC stays below $0.1\,\%$; the warping methods (DTW, WTDTW) are accurate at small -\dvv\ but develop intermittent, large ($>\!1\,\%$) errors from $\sim\!3\,\%$ -true \dvv\ onward as the warp path becomes ill-conditioned, while WCS degrades -smoothly rather than catastrophically, crossing $1\,\%$ error only beyond -$\sim\!4\,\%$ true \dvv. +On a clean, noiseless sweep of true \dvv\ from $-5$ to $5\,\%$ (Fig.~\ref{fig:methods}b), +the phase-wrapped MWCS estimator is the first to break on either branch: its +error stays below $0.1\,\%$ out to $\sim\!1.3\,\%$ true \dvv, then exceeds +$1\,\%$ error by $\sim\!1.5\,\%$ --- the cycle-skip, essentially symmetric in +sign. The stretching family (TS, WTS) stays below $0.1\,\%$ error out to the +full $5\,\%$ tested, on both branches. WCC is just as accurate on the negative +branch (error stays below $0.1\,\%$ throughout) but breaks sharply on the +positive branch, crossing both $0.1\,\%$ and $1\,\%$ error abruptly at the +edge of the tested range ($\sim\!4.75\,\%$) --- a sign asymmetry from the +physical convention itself (Section~\ref{sec:intro}), not a processing +artefact. The warping methods develop large ($>\!1\,\%$) errors asymmetrically +as the warp path becomes ill-conditioned: WTDTW crosses $1\,\%$ error already +at $\sim\!1\,\%$ true \dvv\ on the positive branch but only at $\sim\!3\,\%$ on +the negative branch, and DTW crosses at $\sim\!2.5\,\%$ versus $\sim\!3\,\%$. +WCS degrades smoothly rather than catastrophically, crossing $1\,\%$ error +beyond $\pm4\,\%$ true \dvv\ on either branch. ## Aggregating cross-component results {#sec:aggregation} -Each three-component seismic station (e.g., Z, N, E) carries 6 cross-componet correlations -(ZZ, NN, EE, ZE, ZN, NE), whether they are calculated at the single station or -a inter-station pair. Each carry a signature of the changes in velocity, components -may be dominated by Love or Rayleigh waves [@Lin2008; @Stehly2006], but there scattering and non-straight -ray path induce cross-component leakage between modes [@Hennino2001; @Margerin2019], and thus it is often -assumed in practice that coda waves of cross-components with multi-scattering characteristixcs +Each three-component seismic station (e.g., Z, N, E) carries 6 cross-component correlations +(ZZ, NN, EE, ZE, ZN, NE), whether they are calculated at a single station or +an inter-station pair. Each carries a signature of the changes in velocity; components +may be dominated by Love or Rayleigh waves [@Lin2008; @Stehly2006], but scattering and non-straight +ray paths induce cross-component leakage between modes [@Hennino2001; @Margerin2019], and thus it is often +assumed in practice that coda waves of cross-components with multi-scattering characteristics (e.g., no clearly separated phases) are composed of "surface waves" with strong S-wave sensitivity. Combining them together requires parameter choices, such as averaging them directly [@Liu2014], or weighted (e.g., using coherence-based weighting @Hobiger2012, @DePlaen2016). Combining is another workflow @@ -446,7 +455,7 @@ The remaining choices are no less consequential. The **frequency band** sets the sampled depth: in a two-layer medium, high frequencies recover shallow, often seasonal signals and low frequencies recover a deeper, maybe more tectonic, signal (Fig.~\ref{fig:params}a). On our synthetic two-layer groundwater scenario, a band -matched to each layer (1.5--6\,Hz for the shallow, seasonal signal; 0.2--0.8\,Hz +matched to each layer (1.5--6$\,$Hz for the shallow, seasonal signal; 0.2--0.8$\,$Hz for the deep, drought-trend signal) recovers each with RMS error $\sim\!0.003$--$0.03\,\%$; using the wrong band for a given depth inflates the error to $\sim\!0.10\,\%$ for both --- a $\sim\!4$--30$\times$ degradation. @@ -455,8 +464,8 @@ The **reference** defines what survives: a moving reference re-baselines continuously and erases slow trends that a fixed reference or a joint inversion [@Brenguier2014] preserve (Fig.~\ref{fig:params}c). On the volcano synthetic, a fixed total-stack reference gives RMS $\sim\!0.03\,\%$ and a Brenguier-style -joint inversion $\sim\!0.04\,\%$ (both preserve the pre-eruptive trend), while a -60-day moving reference gives RMS $\sim\!0.16\,\%$ --- roughly 4--5$\times$ worse +joint inversion a comparable $\sim\!0.03\,\%$ (both preserve the pre-eruptive trend), while a +60-day moving reference gives RMS $\sim\!0.16\,\%$ --- roughly $5\times$ worse --- because it re-baselines away the very trend being measured. The **stacking length** trades noise against temporal resolution, rounding off @@ -489,13 +498,14 @@ The **reference** choice also has a construction axis beyond fixed-versus-moving how much of the record a "fixed" reference actually stacks. On the volcano synthetic, a reference built from the whole pre-eruptive record gives RMS $\sim\!0.03\,\%$, one built from only the earliest 15\,\% of that record (an -"early stack," noisier for having fewer days) gives $\sim\!0.04\,\%$ -($\sim\!1.5\times$ worse), a 60-day moving stack gives $\sim\!0.16\,\%$ -(erasing the trend, as above), and the Brenguier-style joint inversion -[@Brenguier2014] gives $\sim\!0.04\,\%$ while additionally preserving the trend ---- so among the reference-construction choices, whole-period and joint -inversion are comparably good, early-stacking is a modest ($\sim\!1.5\times$) -noise penalty, and only the moving reference actively destroys signal. +"early stack") gives a comparable $\sim\!0.03\,\%$ --- a shorter, earlier +reference is not automatically worse: it avoids averaging across the +developing pre-eruptive trend the way a longer stack does. A 60-day moving +stack gives $\sim\!0.16\,\%$ (erasing the trend, as above), and the +Brenguier-style joint inversion [@Brenguier2014] gives $\sim\!0.03\,\%$ while +additionally preserving the trend --- so among the reference-construction +choices, whole-period, early, and joint inversion are all comparably good, +and only the moving reference actively destroys signal. ## Coda window and its covariation with frequency band {#sec:window} @@ -506,9 +516,9 @@ and scattering attenuation both grow with frequency, high-frequency coda energy falls into the noise floor much sooner than low-frequency coda: a coda window that is well past the direct arrival for a $\sim\!0.3$--$0.8\,$Hz band is, at $\sim\!3$--$6\,$Hz, sampling almost pure noise (Fig.~\ref{fig:params}b). On our -synthetic, a fixed 20--40\,s window at the high band gives RMS error -$\sim\!3.8\,\%$ (the noise floor, not the signal), while a window hand-adapted -to the band (3--12\,s) recovers the truth at $\sim\!0.01\,\%$ --- nearly a +synthetic, a fixed 20--40$\,$s window at the high band gives RMS error +$\sim\!3.9\,\%$ (the noise floor, not the signal), while a window hand-adapted +to the band (3--12$\,$s) recovers the truth at $\sim\!0.01\,\%$ --- nearly a 400$\times$ difference from this one choice alone. This is itself an instance of a literature-documented tension: band-matched windowing is recommended, but Table~\ref{tab:survey} shows many surveyed studies instead reuse one fixed @@ -530,11 +540,11 @@ correctly shrinking at the high band, though the low-versus-mid ordering is not perfectly monotonic on this synthetic (an artifact of how the fixed additive noise floor interacts with each band's filter, not a claim that the detector is exact). Recovering \dvv\ with each band's own detected window instead of one -universal fixed (10--30\,s) window (Fig.~\ref{fig:window-envelope}b) gives RMS -$\sim\!0.032\,\%$ vs.\ $\sim\!0.036\,\%$ at the low band (a modest, -$\sim\!1.1\times$ gain), $\sim\!0.017\,\%$ vs.\ $\sim\!0.035\,\%$ at the mid -band ($\sim\!2\times$), and $\sim\!0.019\,\%$ vs.\ $\sim\!1.8\,\%$ at the high -band ($\sim\!93\times$) --- the fixed window is adequate at low frequency and +universal fixed (10--30$\,$s) window (Fig.~\ref{fig:window-envelope}b) gives RMS +$\sim\!0.030\,\%$ vs.\ $\sim\!0.036\,\%$ at the low band (a modest, +$\sim\!1.2\times$ gain), $\sim\!0.017\,\%$ vs.\ $\sim\!0.035\,\%$ at the mid +band ($\sim\!2\times$), and $\sim\!0.020\,\%$ vs.\ $\sim\!2.0\,\%$ at the high +band ($\sim\!100\times$) --- the fixed window is adequate at low frequency and catastrophic at high frequency, while the envelope-derived window is close to the best achievable at every band without ever being told what band it is measuring. @@ -582,7 +592,7 @@ recorded at the two stations, the correlated wavefield in the coda may differ [@ While the interpretation of such coda in terms of Earth's structure effect is difficult [@Snieder2002], the stability of the wavefield excited in the coda is the main requirement for stable \dvv\ measurements [@Hadziioannou2009]. Given the challenge in interpreting both sides independently, -researchers typically would measure \dvv\ on each lag and then report its average [@Kidiwela2026]. +researchers typically measure \dvv\ on each lag and then report its average [@Kidiwela2026]. The causal (positive-lag) and acausal (negative-lag) branches sample opposite-direction paths with different source-side illumination, and in a 3D medium their sensitivity kernels sample partly different volumes, so the two can report @@ -591,11 +601,11 @@ choice (Fig.~\ref{fig:branches}). When a change is *localized* to the volume one branch samples, symmetrizing or averaging the branches --- the common default --- dilutes it toward zero, while the branch that carries the change recovers it (Fig.~\ref{fig:branches}a); here preferring the branch of greatest change is -an researchers' judgement. +a researcher's judgement. When instead both branches share the *same* change, their difference is measurement noise, and selecting the branch of greatest change over-reports it --- a max-of-two-estimators selection bias that grows as -SNR falls (Fig.~\ref{fig:branches}b). When both side exhibit the same change +SNR falls (Fig.~\ref{fig:branches}b). When both sides exhibit the same change (sign, coherence) but with different magnitude, it is reasonable to use the \dvv\ of greatest change given the already low sensitivity in coda waves [@Obermann2013]. An acceptable workflow is to measure both branches, evaluate on that consistency, @@ -626,7 +636,7 @@ the error-bar change it induces. The scenario is a single representative station pair monitoring a shallow volcanic edifice, in the style of the permanent broadband deployments used at effusive/dome volcanoes such as Piton de la Fournaise [@Brenguier2008]: a coda band matched to -shallow depths (0.4--1.0\,Hz), a coda window past the direct arrival (10--30\,s), +shallow depths (0.4--1.0$\,$Hz), a coda window past the direct arrival (10--30$\,$s), and daily correlations sampled every 3 days over 2.5 years at a per-day correlation-coefficient SNR of 7, typical of a continuously operating station. The synthetic ground truth combines an annual seasonal \dvv\ cycle (as from @@ -636,19 +646,21 @@ scenario isolates the measurement-step choices from the network-aggregation choi already covered in Sections~\ref{sec:aggregation}--\ref{sec:uncertainty}. The previous sections only identified single choices, but the overall research workflow -involves them all. We now estimate the combined effects of these parametic choices. +involves them all. We now estimate the combined effects of these parametric choices. Starting from a single **best-practice baseline** (trace stretching, a band matched to the target depth, a coda window well past the direct arrival, a 10-day stack, a long stable reference, coherence gating; the cross-cutting rules of @Brenguier2014 [@Weaver2011; @Clarke2011] as distilled in our survey), we change one parameter at a time to a deviation from best practice documented in the literature and measure the resulting error against the known truth (Fig.~\ref{fig:deviations}). The ranking is unambiguous: relative -to a best-practice RMS error of $\sim\!0.02\,\%$, an unwrapped phase estimator that -cycle-skips is catastrophic, a moving reference and a wrapped-phase MWCS each inflate the -error roughly an order of magnitude, an over-long stack and a late low-SNR window -distort the co-eruptive drop, while the frequency band (which mostly sets -precision here) and the gating are second-order. The joint inversion is -nearly as good as the fixed reference and preserves the trend. +to a best-practice RMS error of $\sim\!0.03\,\%$, the two-dimensionally-unwrapped +WCS estimator is catastrophic here ($\sim\!50\times$ worse), a wrapped-phase MWCS +inflates the error by roughly an order of magnitude ($\sim\!12\times$), and DTW by +$\sim\!6\times$; among the non-estimator choices a moving reference is worst +($\sim\!5\times$, and most distorts the recovered drop), while stack length, coda +window, and frequency band deviations are each more modest +($1.5$--$2.5\times$). The joint-inversion reference stays closest to the +fixed-reference baseline among the reference-scheme deviations. \begin{figure} \centering @@ -672,10 +684,12 @@ widest exactly at the sharp co-eruptive drop; the RMS error against the known truth spans over two orders of magnitude across pipelines ($\sim\!0.02$--$2.8\,\%$). Attributing the variance of the outcome to each axis with a first-order (main-effect) sensitivity index -(Fig.~\ref{fig:multiverse}b) shows that, for this dataset, the **coda window and -the stack length** control the RMS error most, the **stack length** dominates the -recovered drop amplitude, and the estimator is third; the first-order indices sum -to well under one, so a large part of the spread is *interaction* between choices +(Fig.~\ref{fig:multiverse}b) shows that, for this dataset, the **coda window** +controls the RMS error most, with the estimator second and the stack length +third; for the recovered drop amplitude the order changes to **stack length** +first, window second, and estimator third. The first-order indices sum to +well under one in both cases ($\sim\!0.66$ for RMS, $\sim\!0.63$ for the +drop), so a large part of the spread is *interaction* between choices compounding. The appropriate object is therefore not a single curve but a distribution of \dvv\ over the processing choices. @@ -687,8 +701,9 @@ but a distribution of \dvv\ over the processing choices. each coloured by its RMS error against the truth on a colourblind-safe scale (bright accurate, dark biased; the worst run off the clipped axis); the grey band is the 10--90\% inter-pipeline spread, widest at the velocity drop. (b) First-order variance attribution: which - choice controls the RMS error and the recovered drop. Window and stack length - dominate; the sub-unity sum signals strong interactions.} + choice controls the RMS error and the recovered drop. Window, estimator, and + stack length dominate, in that order for RMS and reordered for the drop; + the sub-unity sum signals strong interactions.} \label{fig:multiverse} \end{figure} @@ -828,9 +843,76 @@ separate reimplementation. The retrospective run is validated against the already-published @Clements2023 CI.LJR result before any new claim is drawn from it. - +After the correction, single-station \dvv\ (NoisePy correlations, a +codameter 5-member ensemble, 2--4$\,$Hz, 2018--2019) validates against the +published @Clements2023 product. The comparison is not trivial to get right: +the CD2022 90-day-comp product is a *trailing* 90-day stack, so it lags a +centered-smoothed daily series by about 45 days, and comparing without +matching that smoothing caps the correlation near 0.7 even on a real annual +cycle (Fig.~\ref{fig:realdata-validation}). We match by applying the same +trailing 90-day mean to the daily series, compare demeaned --- the two +products reference different epochs, and a constant offset is bookkeeping, +not error --- and exclude the first 150 days of each station's series as +reference burn-in. Under this matched comparison, CI.LJR reaches $r=0.990$ +(681 overlapping days); CI.ARV reaches $r=0.66$--$0.92$ depending on the join +method, reflecting real 2018 data gaps; CI.RXH reaches $r=0.68$, the weakest +of the three, on a station whose recovered \dvv\ is nearly flat and whose +signal is low to begin with. + +\begin{figure} + \centering + \includegraphics[width=\textwidth]{realdata_1_validation.png} + \caption{Single-station \dvv\ at three CI stations, 2018--2019, against the + published Clements \& Denolle (2022) product (dashed, reference-shifted). + Daily \dvv\ (points) with the ensemble spread (shaded, epistemic) and + measurement error bars, and a 45-day-smoothed curve. The annotated $r$ is + computed on this 45-day-smoothed comparison, not the smoothing-matched + comparison reported in the text; matching CD2022's own trailing 90-day + window instead of a centered 45-day one raises $r$ at every station (to + 0.990/0.66--0.92/0.68 for LJR/ARV/RXH), which is why the two numbers + differ.} + \label{fig:realdata-validation} +\end{figure} + +Why validation quality differs so much by station is visible in the +correlation function itself (Fig.~\ref{fig:realdata-interferograms}). CI.LJR +shows a stable, narrow coda near zero lag through the year. CI.RXH shows +multipath, several coherent arrivals spread across the full $\pm 8\,$s of +lag shown, and a visible shift in that pattern in April--May 2019, +consistent with a site change rather than a processing artefact. CI.ARV's +coherent energy is compact and concentrated near zero lag but comparatively +sparse, consistent with its higher scatter and join-method sensitivity. + +\begin{figure} + \centering + \includegraphics[width=\textwidth]{realdata_2_interferograms.png} + \caption{Daily north-south/east-west cross-component correlation + functions, 2--4$\,$Hz, $\pm 8\,$s lag, 2019. CI.LJR's coda is stable and + narrow. CI.RXH shows multipath (several coherent bands across the full lag + range) and a shift in that pattern in April--May 2019. CI.ARV's coherent + energy is compact and near zero lag but comparatively sparse. The visual + difference is the reason validation quality differs by station in + Fig.~\ref{fig:realdata-validation}.} + \label{fig:realdata-interferograms} +\end{figure} + +As an optional supplement, Fig.~\ref{fig:realdata-warmup} shows the +ensemble's warm-up behaviour on a separate 90-day smoke run at CI.LJR: \dvv$\,$is +undefined until enough history has accumulated for the moving-reference +member to compute a trailing reference, and each band's per-epoch stretching +correlation coefficient is reported alongside the recovered series. + +\begin{figure} + \centering + \includegraphics[width=\textwidth]{realdata_3_warmup.png} + \caption{CI.LJR single-station \dvv, a separate 90-day smoke run (January--April + 2023), four frequency bands. The ensemble spread (shaded) and measurement + error bars are reported once at least one member is defined; the annotated + member count grows as the moving-reference member's warm-up period + elapses, reaching all 5 members later in the window. Bottom panel: + per-band stretching correlation coefficient.} + \label{fig:realdata-warmup} +\end{figure} # Depth: propagating the measurement covariance through sensitivity kernels {#sec:depth} @@ -843,7 +925,7 @@ yet report. The executable stages are in development in the open codameter packa The field has converged on one physical rule for depth --- *depth is set by the frequency band and the coda lapse time, it is not assumed* [@Obermann2013; @Obermann2016] --- and, increasingly, on multi-band measurement as the means to -resolve it [@Takano2017; @Feng2020; @Mao2022, @Mao2025]. The step from a set of per-band +resolve it [@Takano2017; @Feng2020; @Mao2022; @Mao2025]. The step from a set of per-band \dvv\ series to a depth profile is where reporting is least consistent: many studies read a single band as a single depth, only a few invert several bands against surface-wave sensitivity kernels, and the propagation of the measurement @@ -1064,8 +1146,10 @@ package in which the depth stage is being built. Let $r(t)$ and $c(t)$ be the reference and current cross-correlations, band-passed to $[f_1,f_2]$ (central frequency $f_c$, bandwidth $B$) and read over -the coda window $W=[t_1,t_2]$ on one or both branches. We use the convention -$\delta t/t = -\,\dvv$, and write $\varepsilon$ for the recovered \dvv. +the coda window $W=[t_1,t_2]$ on one or both branches. Each method below +estimates a trial stretch factor $\varepsilon$, the fractional dilation of +the current coda relative to the reference; physical \dvv\ follows via the +convention boxed in Section~\ref{sec:intro}, $\dvv=-\varepsilon/(1+\varepsilon)$. **Trace stretching (TS).** For a trial stretch $s$, define $r_s(t)=r\!\left(t/(1+s)\right)$ and the windowed correlation coefficient @@ -1112,16 +1196,19 @@ per-scale estimates are pooled with cross-wavelet-power weights. # Aggregation and uncertainty conventions {#app:aggregation} For component $k$ of a station pair, stretching yields the image -$\mathrm{CC}_k(\varepsilon,t)$. - -**Approach A (average the \dvv).** -$\varepsilon_k(t)=\arg\max_{\varepsilon} \mathrm{CC}_k(\varepsilon,t)$ and the pair -estimate is the (possibly weighted) mean $x_p(t)=\sum_k w_k\,\varepsilon_k/\sum_k -w_k$, with $w_k=\max_{\varepsilon} \mathrm{CC}_k$ (coherence-weighted) or $w_k=1$ -(unweighted). - -**Approach B (average the images).** -$x_p(t)=\arg\max_{\varepsilon}\big[\tfrac1{N_c}\sum_k \mathrm{CC}_k(\varepsilon,t)\big]$. +$\mathrm{CC}_k(\varepsilon,t)$; each per-component stretch factor converts to +physical $\dvv_k$ via the boxed relation. + +**Approach A (average the per-component \dvv).** +$\dvv_k(t)$ follows from $\varepsilon_k(t)=\arg\max_{\varepsilon} +\mathrm{CC}_k(\varepsilon,t)$, and the pair estimate is the (possibly +weighted) mean $x_p(t)=\sum_k w_k\,\dvv_k/\sum_k w_k$, with $w_k=\max_{\varepsilon} +\mathrm{CC}_k$ (coherence-weighted) or $w_k=1$ (unweighted). + +**Approach B (average the images).** The pair stretch factor +$\varepsilon_p(t)=\arg\max_{\varepsilon}\big[\tfrac1{N_c}\sum_k +\mathrm{CC}_k(\varepsilon,t)\big]$ converts to $x_p(t)=\dvv_p(t)$ via the same +relation. **Network over pairs.** With pair weights $W_p$ (e.g.\ mean coherence), \begin{equation} diff --git a/paper/manuscript_marine.tex b/paper/manuscript_marine.tex index 0626455..b32f377 100644 --- a/paper/manuscript_marine.tex +++ b/paper/manuscript_marine.tex @@ -246,7 +246,7 @@ \title{The reproducibility cost of ad-hoc processing choices in ambient-noise seismic velocity-change monitoring} \author{M. A. Denolle} -\date{2026-08-08} +\date{2026-08-11} \begin{document} \maketitle \begin{abstract} @@ -284,20 +284,35 @@ \section{Introduction}\label{sec:intro} -Changes in subsurface properties occur due geodynamics, which drive +Changes in subsurface properties occur due to geodynamics, which drive earthquake damage and volcanic eruption, and hydrodynamics, which -controls fluid exchange between the atomsphere and the solid Earth. -These processes influence the mechanical property of Earth materials, -which directly affect the speed at which seismic waves propagates. -Changes in seismic velocity, often measured and referred to as \dvv, can -be tracked by measuring changes in arrival times of seismic waves, -especially scatterd waves such as coda waves, provided that the source -and receivers are at the same location. +controls fluid exchange between the atmosphere and the solid Earth. +These processes influence the mechanical properties of Earth materials, +which directly affect the speed at which seismic waves propagate. +Changes in seismic velocity, often measured and referred to as +\dvv\(\,\)can be tracked by measuring changes in arrival times of +seismic waves, especially scattered waves such as coda waves, provided +that the source and receivers are at the same location. + +\begin{center} +\fbox{\begin{minipage}{0.95\linewidth} +\textbf{Sign convention.} \dvv\ is the fractional seismic velocity change, +positive for a velocity \emph{increase}. The stretching family of estimators +measures the stretch factor $\varepsilon$ that maps the reference coda onto +the current coda; the exact relation is +\[ +\dvv = -\frac{\varepsilon}{1+\varepsilon}. +\] +The first-order shortcut $\dvv \approx -\varepsilon$ is accurate to order +$\varepsilon^2$ --- negligible below $1\,\%$ but a $0.17\,\%$ absolute bias at +a true $\dvv$ of $4\,\%$, the changes seen on landslides. +\end{minipage}} +\end{center} Due to the sensitivity of coda waves to small perturbations in the material properties, \dvv~was discovered as an effective method to -monitor chanes during volcanic unrest since its discovery at Piton de la -Fournaise \citep{Brenguier2008} and is now calculated in continuous +monitor changes during volcanic unrest since its discovery at Piton de +la Fournaise \citep{Brenguier2008} and is now calculated in continuous along side of more conventional seismic monitoring methods in the same volcano observatory \citep{Duputel2009} and was a determining early warning parameter used by the Icelandic Meteorological Office used @@ -314,50 +329,48 @@ \section{Introduction}\label{sec:intro} acts on is a property of the subsurface, not of the analyst's processing choices. -The elevated sensitivity comes at the price a long series of processing -choices and at almost every step the analyst makes a choice. Among these -are choices of estimators between windowed phase measurements or -strething \citep[\citet{Mao2020}, \citet{Yuan2021}]{Mikesell2015}, -frequency band and coda window (which together set the sampled depth; -\citeauthor{Obermann2013} +The elevated sensitivity comes at the price of a long series of +processing choices, and at almost every step the analyst makes a choice. +Among these are choices of estimators between windowed phase +measurements or stretching \citep[\citet{Mao2020}, +\citet{Yuan2021}]{Mikesell2015}, frequency band and coda window (which +together set the sampled depth; \citeauthor{Obermann2013} \citetext{\citeyear{Obermann2013}; \citealp{Obermann2016}}), reference window \citep[\citet{Ermert2023}, \citet{Okubo2024}]{Brenguier2014}, how much to substack to increase coherence among windows (Hadzianou Celine), and how to aggregate and weight the many cross-component and station-pair measurements that make up a single reported \dvv~time -series (e.g., \citet{Hobiger2012}). These choices are typically made by -habit, justified briefly if at all , and rarely reported in enough -detail to reproduce. The community has long flagged \emph{individual} -pitfalls (spurious changes from non-stationary noise, \citet{Zhan2013}; +series (e.g., \citet{Hobiger2012}). These choices are made by habit, +justified briefly if at all, and rarely reported in enough detail to +reproduce. The community has long flagged \emph{individual} pitfalls +(spurious changes from non-stationary noise, \citet{Zhan2013}; measurement-error formulae, \citeauthor{Clarke2011} \citetext{\citeyear{Clarke2011}; \citealp{Weaver2011}}), but the \emph{cumulative} effect of the full choice set on both the value and -its stated uncertainty has not been quantified comprehensively in the -literature. +its stated uncertainty has not been quantified in the literature. This is a reproducibility problem of exactly the ``garden of forking paths'' type identified in the statistical sciences \citep{Gelman2013, Steegen2016}: many individually reasonable analyses -of the same data give different answers, and without full reporting -parameter choices we cannot fully interpret \dvv~as robust values. . +of the same data give different answers, and without full reporting of +parameter choices we cannot interpret a reported \dvv~as a robust value. Here we make the problem concrete for \dvv~monitoring. We use purely synthetic correlation time series to test the methods (estimators) and -parameter choices that the community makes to estimate \dvv, which we -report over 103 studies Appendix\textasciitilde{}\ref{app:survey}. We do -not aim to report the ``best'' pipeline, which is most often the one -reported in scientific papers, but instead document the various -parameter impacts (Section\textasciitilde{}\ref{sec:results}). We then -propose a new measurement error that incoporates these effects into a -data covariance matrix \(C_d\) -(Section\textasciitilde{}\ref{sec:bayes}). - -One example of propagating such error in downstream scientific insights -is the migration of the surface \dvv~measurement to depth profiles of +parameter choices that the community makes to estimate \dvv\(\,\)which +we report over 103 studies in Appendix\textasciitilde{}\ref{app:survey}. +We do not aim to report the ``best'' pipeline, which is most often the +one reported in scientific papers, but instead document the parameter +impacts (Section\textasciitilde{}\ref{sec:results}). We then propose a +new measurement error that incorporates these effects into a data +covariance matrix \(C_d\) (Section\textasciitilde{}\ref{sec:bayes}). + +One example of propagating such error into downstream science is the +migration of the surface \dvv~measurement to depth profiles of perturbations in shear wave velocity \(\Delta V_S(z)/V_S(z)\), which often depends on the wavefield constituting the coda waves, such as surface waves or body waves, and that depend on the source-receiver pair geometry. We illustrate the propagation of errors to a depth profile -(Section\textasciitilde{}\ref{sec:depth}). We only utilize synthetic +(Section\textasciitilde{}\ref{sec:depth}). We only use synthetic examples for ground truthing on the signal processing parameters, since the concepts behind the observations of phase lags in scattered waves is well established \citep{Obermann2013}. We package this new methodology @@ -409,7 +422,7 @@ \section{Synthetic framework}\label{sec:methods} A homogeneous velocity change is imposed exactly by stretching the lapse-time axis, -\(u_{\rm cur}(t) = u_{\rm ref}\!\left(t/(1+\dvv)\right)\), and a +\(u_{\rm cur}(t) = u_{\rm ref}\!\left(t\,(1+\dvv)\right)\), and a repeated time series is produced by generating this stretched coda with a prescribed ground-truth \(\dvv(t)\) and additive band-limited noise at a controlled signal-to-noise ratio. The concept has been demonstrated @@ -419,16 +432,16 @@ \section{Synthetic framework}\label{sec:methods} every departure of a recovered series from it is an artefact of the processing, not of the data. -We utilize seven \dvv~estimators that were implemented in -\texttt{noisepy} in the \texttt{monitoring\_methods} module -\citep{Jiang2020}: trace stretching (TS, \citet{lobkis01}), windowed -cross-correlation (WCC, \citet{poupinet84}), dynamic time warping (DTW, -\citet{Mikesell2015}), te moving-window cross-spectrum (MWCS; -\citet{Clarke2011}), and three wavelet-domain methods, the wavelet -cross-spectrum (WCS, \citet{Mao2020}) and the two wavelet stretching -(WTS) and wavelet DTW (WTDTW) introduced and benchmarked numerically by -\citet{Yuan2021}. The framework, the figures below, and a implementation -are released in the open \texttt{codameter} package +We use seven \dvv~estimators that were implemented in \texttt{noisepy} +in the \texttt{monitoring\_methods} module \citep{Jiang2020}: trace +stretching (TS, \citet{lobkis01}), windowed cross-correlation (WCC, +\citet{poupinet84}), dynamic time warping (DTW, \citet{Mikesell2015}), +the moving-window cross-spectrum (MWCS; \citet{Clarke2011}), and three +wavelet-domain methods, the wavelet cross-spectrum (WCS, +\citet{Mao2020}) and the two wavelet stretching (WTS) and wavelet DTW +(WTDTW) introduced and benchmarked numerically by \citet{Yuan2021}. The +framework, the figures below, and an implementation are released in the +open \texttt{codameter} package (Section\textasciitilde{}\ref{sec:discussion}). Formal definitions of all seven methods, of the two aggregation pathways, and of the uncertainty conventions are described in @@ -507,12 +520,11 @@ \section{\texorpdfstring{parameter-dependent \dvv~and its subsequent step in the inference chain (Sections\textasciitilde{}\ref{sec:depth}--\ref{sec:stress}) consumes. -Each of the parametric components is implemented in a single -comprehensive package \texttt{codameter} that borrows from -\texttt{msnoise} \citep{Lecocq2014} and \texttt{noisepy} -\citep{Jiang2020} to extract the measurement of \dvv~from the ambient -noise monitoring workflow and generalize it to any coda wave from -repeated source-receiver paths. +Each of the parametric components is implemented in a single package, +\texttt{codameter}, that borrows from \texttt{msnoise} +\citep{Lecocq2014} and \texttt{noisepy} \citep{Jiang2020} to extract the +measurement of \dvv~from the ambient noise monitoring workflow and +generalize it to any coda wave from repeated source-receiver paths. Table\textasciitilde{}\ref{tab:results-synthesis} previews the RMS error against the known synthetic ground truth for the best- and worst-case @@ -533,8 +545,8 @@ \section{\texorpdfstring{parameter-dependent \dvv~and its \toprule \textbf{Axis} & \textbf{Best-case RMS} & \textbf{Worst-case RMS} & \textbf{Section} \\ \midrule -Estimator (family split) & $<\!0.01\,\%$ (TS/WTS, up to 5\,\% true \dvv) & -cycle-skip $>\!0.5\,\%$ past $\sim\!1.4\,\%$ true \dvv\ (MWCS) & +Estimator (family split) & $<\!0.01\,\%$ (TS, up to 5\,\% true \dvv) & +cycle-skip $>\!0.5\,\%$ past $\sim\!1.5\,\%$ true \dvv\ (MWCS) & Section~\ref{sec:methods-fig} \\ \hline Cross-component aggregation & $\sim\!0.03\,\%$ (Approach B, averaged images) & $\sim\!0.31\,\%$ (Approach A, unweighted) & Section~\ref{sec:aggregation} \\ \hline @@ -543,8 +555,8 @@ \section{\texorpdfstring{parameter-dependent \dvv~and its Frequency band & $\sim\!0.003$--$0.03\,\%$ (matched to depth) & $\sim\!0.10\,\%$ (mismatched) & Section~\ref{sec:params} \\ \hline Coda window & $\sim\!0.01\,\%$ (adapted to band) & -$\sim\!3.8\,\%$ (fixed, wrong band) & Section~\ref{sec:window} \\ \hline -Reference scheme & $\sim\!0.03$--$0.04\,\%$ (fixed stack / joint inversion) & +$\sim\!3.9\,\%$ (fixed, wrong band) & Section~\ref{sec:window} \\ \hline +Reference scheme & $\sim\!0.03\,\%$ (fixed stack / joint inversion) & $\sim\!0.16\,\%$ (60-day moving) & Section~\ref{sec:params} \\ \hline Stack length & $\sim\!0.018\,\%$ (10-day) & $\sim\!0.044\,\%$ (1-day, noisy) & Section~\ref{sec:params} \\ @@ -554,12 +566,19 @@ \section{\texorpdfstring{parameter-dependent \dvv~and its \subsection{Estimator family}\label{sec:methods-fig} -On small, clean \dvv, all seven estimators agree, demonstrating -robustness of the methods (Fig.\textasciitilde{}\ref{fig:methods}a). -Sweeping the same clean recovery out to \(\pm 5\,\%\) -(Fig.\textasciitilde{}\ref{fig:methods}b) shows exactly where and how -each family first departs from the 1:1 line, and the estimator choice -becomes consequential at large, noisy +As of codameter v0.4.0, all seven estimators return physical \dvv~under +the sign convention above rather than the raw stretch factor +\(\varepsilon\); the synthetic generator imposes changes in the same +convention, so a positive imposed \dvv~recovers as positive. +\texttt{tests/test\_sign\_convention.py} holds every estimator to this +convention in both signs, and checks the full \texttt{run\_pipeline} +call path end to end, so the convention cannot silently drift back. + +On small, clean \dvv\(\,\)all seven estimators agree +(Fig.\textasciitilde{}\ref{fig:methods}a). Sweeping the same clean +recovery out to \(\pm 5\,\%\) (Fig.\textasciitilde{}\ref{fig:methods}b) +shows exactly where and how each family first departs from the 1:1 line, +and the estimator choice becomes consequential at large, noisy \dvv~(Fig.\textasciitilde{}\ref{fig:methods}c), where the effect of the methods is split according to their phase measurement approaches. The stretching family (TS, WTS) and WCC match the whole dilated coda and @@ -567,46 +586,56 @@ \subsection{Estimator family}\label{sec:methods-fig} wrapped phase and may cycle-skips, while the \emph{same} cross-wavelet phase (WCS), once unwrapped in 2-D, recovers the change. The warping methods (DTW, WTDTW) track but under-shoot the largest strains. No -estimator is simply ``right'' and choices in the estimators can alters -the \dvv~measurements for larger strain changes. +estimator is simply ``right''; the choice of estimator materially +changes the \dvv~measurement at larger strain. \begin{figure} \centering \includegraphics[width=\textwidth]{demo_1_methods.png} \caption{Estimator choice across the seven NoisePy methods. (a) Clean, small \dvv: all agree. (b) The same clean recovery swept over $\pm 5\,\%$ - true \dvv: MWCS cycle-skips almost immediately past $\pm 1$--$2\,\%$, DTW/WTDTW - break past $\pm 3$--$4\,\%$, WCS degrades smoothly, while TS/WTS/WCC track the - 1:1 line throughout. (c) Large, noisy \dvv\ (a pre-failure landslide signal): + true \dvv: MWCS cycle-skips past $\sim\!1.5\,\%$ on either branch; TS and WTS + track the 1:1 line throughout; WCC tracks just as tightly on the negative + branch but breaks sharply near the positive edge; DTW and WTDTW break + asymmetrically, WTDTW near $+1\,\%$ but only near $-3\,\%$ on the other + branch; WCS degrades smoothly, crossing $1\,\%$ error beyond $\pm 4\,\%$. + (c) Large, noisy \dvv\ (a pre-failure landslide signal): MWCS cycle-skips, 2-D-unwrapped WCS and the stretching family track, the warping methods under-shoot.} \label{fig:methods} \end{figure} -On a clean, noiseless sweep of true \dvv~from 0 to 5,\% +On a clean, noiseless sweep of true \dvv~from \(-5\) to \(5\,\%\) (Fig.\textasciitilde{}\ref{fig:methods}b), the phase-wrapped MWCS -estimator is the first to break: its error stays below \(0.1\,\%\) up to -a true \dvv~of \(\sim\!1.2\,\%\), then jumps discontinuously past -\(0.5\,\%\) error by \(\sim\!1.4\,\%\) true \dvv~--- the cycle-skip. The -stretching family (TS, WTS) stays below \(0.01\,\%\) error out to the -full 5,\% tested, and WCC stays below \(0.1\,\%\); the warping methods -(DTW, WTDTW) are accurate at small \dvv~but develop intermittent, large -(\(>\!1\,\%\)) errors from \(\sim\!3\,\%\) true \dvv~onward as the warp -path becomes ill-conditioned, while WCS degrades smoothly rather than -catastrophically, crossing \(1\,\%\) error only beyond \(\sim\!4\,\%\) -true \dvv. +estimator is the first to break on either branch: its error stays below +\(0.1\,\%\) out to \(\sim\!1.3\,\%\) true \dvv, then exceeds \(1\,\%\) +error by \(\sim\!1.5\,\%\) --- the cycle-skip, essentially symmetric in +sign. The stretching family (TS, WTS) stays below \(0.1\,\%\) error out +to the full \(5\,\%\) tested, on both branches. WCC is just as accurate +on the negative branch (error stays below \(0.1\,\%\) throughout) but +breaks sharply on the positive branch, crossing both \(0.1\,\%\) and +\(1\,\%\) error abruptly at the edge of the tested range +(\(\sim\!4.75\,\%\)) --- a sign asymmetry from the physical convention +itself (Section\textasciitilde{}\ref{sec:intro}), not a processing +artefact. The warping methods develop large (\(>\!1\,\%\)) errors +asymmetrically as the warp path becomes ill-conditioned: WTDTW crosses +\(1\,\%\) error already at \(\sim\!1\,\%\) true \dvv~on the positive +branch but only at \(\sim\!3\,\%\) on the negative branch, and DTW +crosses at \(\sim\!2.5\,\%\) versus \(\sim\!3\,\%\). WCS degrades +smoothly rather than catastrophically, crossing \(1\,\%\) error beyond +\(\pm4\,\%\) true \dvv~on either branch. \subsection{Aggregating cross-component results}\label{sec:aggregation} Each three-component seismic station (e.g., Z, N, E) carries 6 -cross-componet correlations (ZZ, NN, EE, ZE, ZN, NE), whether they are -calculated at the single station or a inter-station pair. Each carry a -signature of the changes in velocity, components may be dominated by -Love or Rayleigh waves \citep{Lin2008, Stehly2006}, but there scattering -and non-straight ray path induce cross-component leakage between modes +cross-component correlations (ZZ, NN, EE, ZE, ZN, NE), whether they are +calculated at a single station or an inter-station pair. Each carries a +signature of the changes in velocity; components may be dominated by +Love or Rayleigh waves \citep{Lin2008, Stehly2006}, but scattering and +non-straight ray paths induce cross-component leakage between modes \citep{Hennino2001, Margerin2019}, and thus it is often assumed in practice that coda waves of cross-components with multi-scattering -characteristixcs (e.g., no clearly separated phases) are composed of +characteristics (e.g., no clearly separated phases) are composed of ``surface waves'' with strong S-wave sensitivity. Combining them together requires parameter choices, such as averaging them directly \citep{Liu2014}, or weighted (e.g., using coherence-based weighting @@ -707,20 +736,21 @@ \subsection{Frequency band, reference and stacking}\label{sec:params} recover shallow, often seasonal signals and low frequencies recover a deeper, maybe more tectonic, signal (Fig.\textasciitilde{}\ref{fig:params}a). On our synthetic two-layer -groundwater scenario, a band matched to each layer (1.5--6,Hz for the -shallow, seasonal signal; 0.2--0.8,Hz for the deep, drought-trend -signal) recovers each with RMS error \(\sim\!0.003\)--\(0.03\,\%\); -using the wrong band for a given depth inflates the error to -\(\sim\!0.10\,\%\) for both --- a \(\sim\!4\)--30\(\times\) degradation. +groundwater scenario, a band matched to each layer (1.5--6\(\,\)Hz for +the shallow, seasonal signal; 0.2--0.8\(\,\)Hz for the deep, +drought-trend signal) recovers each with RMS error +\(\sim\!0.003\)--\(0.03\,\%\); using the wrong band for a given depth +inflates the error to \(\sim\!0.10\,\%\) for both --- a +\(\sim\!4\)--30\(\times\) degradation. The \textbf{reference} defines what survives: a moving reference re-baselines continuously and erases slow trends that a fixed reference or a joint inversion \citep{Brenguier2014} preserve (Fig.\textasciitilde{}\ref{fig:params}c). On the volcano synthetic, a fixed total-stack reference gives RMS \(\sim\!0.03\,\%\) and a -Brenguier-style joint inversion \(\sim\!0.04\,\%\) (both preserve the -pre-eruptive trend), while a 60-day moving reference gives RMS -\(\sim\!0.16\,\%\) --- roughly 4--5\(\times\) worse --- because it +Brenguier-style joint inversion a comparable \(\sim\!0.03\,\%\) (both +preserve the pre-eruptive trend), while a 60-day moving reference gives +RMS \(\sim\!0.16\,\%\) --- roughly \(5\times\) worse --- because it re-baselines away the very trend being measured. The \textbf{stacking length} trades noise against temporal resolution, @@ -755,14 +785,15 @@ \subsection{Frequency band, reference and stacking}\label{sec:params} fixed-versus-moving: how much of the record a ``fixed'' reference actually stacks. On the volcano synthetic, a reference built from the whole pre-eruptive record gives RMS \(\sim\!0.03\,\%\), one built from -only the earliest 15,\% of that record (an ``early stack,'' noisier for -having fewer days) gives \(\sim\!0.04\,\%\) (\(\sim\!1.5\times\) worse), -a 60-day moving stack gives \(\sim\!0.16\,\%\) (erasing the trend, as -above), and the Brenguier-style joint inversion \citep{Brenguier2014} -gives \(\sim\!0.04\,\%\) while additionally preserving the trend --- so -among the reference-construction choices, whole-period and joint -inversion are comparably good, early-stacking is a modest -(\(\sim\!1.5\times\)) noise penalty, and only the moving reference +only the earliest 15,\% of that record (an ``early stack'') gives a +comparable \(\sim\!0.03\,\%\) --- a shorter, earlier reference is not +automatically worse: it avoids averaging across the developing +pre-eruptive trend the way a longer stack does. A 60-day moving stack +gives \(\sim\!0.16\,\%\) (erasing the trend, as above), and the +Brenguier-style joint inversion \citep{Brenguier2014} gives +\(\sim\!0.03\,\%\) while additionally preserving the trend --- so among +the reference-construction choices, whole-period, early, and joint +inversion are all comparably good, and only the moving reference actively destroys signal. \subsection{Coda window and its covariation with frequency @@ -777,10 +808,10 @@ \subsection{Coda window and its covariation with frequency past the direct arrival for a \(\sim\!0.3\)--\(0.8\,\)Hz band is, at \(\sim\!3\)--\(6\,\)Hz, sampling almost pure noise (Fig.\textasciitilde{}\ref{fig:params}b). On our synthetic, a fixed -20--40,s window at the high band gives RMS error \(\sim\!3.8\,\%\) (the -noise floor, not the signal), while a window hand-adapted to the band -(3--12,s) recovers the truth at \(\sim\!0.01\,\%\) --- nearly a -400\(\times\) difference from this one choice alone. This is itself an +20--40\(\,\)s window at the high band gives RMS error \(\sim\!3.9\,\%\) +(the noise floor, not the signal), while a window hand-adapted to the +band (3--12\(\,\)s) recovers the truth at \(\sim\!0.01\,\%\) --- nearly +a 400\(\times\) difference from this one choice alone. This is itself an instance of a literature-documented tension: band-matched windowing is recommended, but Table\textasciitilde{}\ref{tab:survey} shows many surveyed studies instead reuse one fixed window across bands. @@ -803,12 +834,12 @@ \subsection{Coda window and its covariation with frequency this synthetic (an artifact of how the fixed additive noise floor interacts with each band's filter, not a claim that the detector is exact). Recovering \dvv~with each band's own detected window instead of -one universal fixed (10--30,s) window +one universal fixed (10--30\(\,\)s) window (Fig.\textasciitilde{}\ref{fig:window-envelope}b) gives RMS -\(\sim\!0.032\,\%\) vs.~\(\sim\!0.036\,\%\) at the low band (a modest, -\(\sim\!1.1\times\) gain), \(\sim\!0.017\,\%\) vs.~\(\sim\!0.035\,\%\) -at the mid band (\(\sim\!2\times\)), and \(\sim\!0.019\,\%\) -vs.~\(\sim\!1.8\,\%\) at the high band (\(\sim\!93\times\)) --- the +\(\sim\!0.030\,\%\) vs.~\(\sim\!0.036\,\%\) at the low band (a modest, +\(\sim\!1.2\times\) gain), \(\sim\!0.017\,\%\) vs.~\(\sim\!0.035\,\%\) +at the mid band (\(\sim\!2\times\)), and \(\sim\!0.020\,\%\) +vs.~\(\sim\!2.0\,\%\) at the high band (\(\sim\!100\times\)) --- the fixed window is adequate at low frequency and catastrophic at high frequency, while the envelope-derived window is close to the best achievable at every band without ever being told what band it is @@ -862,9 +893,9 @@ \subsection{Causal and acausal branches}\label{sec:branches} \citep{Snieder2002}, the stability of the wavefield excited in the coda is the main requirement for stable \dvv~measurements \citep{Hadziioannou2009}. Given the challenge in interpreting both sides -independently, researchers typically would measure \dvv~on each lag and -then report its average \citep{Kidiwela2026}. The causal (positive-lag) -and acausal (negative-lag) branches sample opposite-direction paths with +independently, researchers typically measure \dvv~on each lag and then +report its average \citep{Kidiwela2026}. The causal (positive-lag) and +acausal (negative-lag) branches sample opposite-direction paths with different source-side illumination, and in a 3D medium their sensitivity kernels sample partly different volumes, so the two can report genuinely different \dvv~without either being wrong. Two regimes bound the choice @@ -873,11 +904,11 @@ \subsection{Causal and acausal branches}\label{sec:branches} averaging the branches --- the common default --- dilutes it toward zero, while the branch that carries the change recovers it (Fig.\textasciitilde{}\ref{fig:branches}a); here preferring the branch -of greatest change is an researchers' judgement. When instead both +of greatest change is a researcher's judgement. When instead both branches share the \emph{same} change, their difference is measurement noise, and selecting the branch of greatest change over-reports it --- a max-of-two-estimators selection bias that grows as SNR falls -(Fig.\textasciitilde{}\ref{fig:branches}b). When both side exhibit the +(Fig.\textasciitilde{}\ref{fig:branches}b). When both sides exhibit the same change (sign, coherence) but with different magnitude, it is reasonable to use the \dvv~of greatest change given the already low sensitivity in coda waves \citep{Obermann2013}. An acceptable workflow @@ -912,20 +943,20 @@ \section{The multiverse of processing choices}\label{sec:multiverse} shallow volcanic edifice, in the style of the permanent broadband deployments used at effusive/dome volcanoes such as Piton de la Fournaise \citep{Brenguier2008}: a coda band matched to shallow depths -(0.4--1.0,Hz), a coda window past the direct arrival (10--30,s), and -daily correlations sampled every 3 days over 2.5 years at a per-day -correlation-coefficient SNR of 7, typical of a continuously operating -station. The synthetic ground truth combines an annual seasonal -\dvv~cycle (as from near-surface thermoelastic/hydrologic effects), a -slow pre-eruptive inflation ramp, and a sharp co-eruptive velocity drop -with partial recovery. This single-pair scenario isolates the -measurement-step choices from the network-aggregation choices already -covered in +(0.4--1.0\(\,\)Hz), a coda window past the direct arrival +(10--30\(\,\)s), and daily correlations sampled every 3 days over 2.5 +years at a per-day correlation-coefficient SNR of 7, typical of a +continuously operating station. The synthetic ground truth combines an +annual seasonal \dvv~cycle (as from near-surface +thermoelastic/hydrologic effects), a slow pre-eruptive inflation ramp, +and a sharp co-eruptive velocity drop with partial recovery. This +single-pair scenario isolates the measurement-step choices from the +network-aggregation choices already covered in Sections\textasciitilde{}\ref{sec:aggregation}--\ref{sec:uncertainty}. The previous sections only identified single choices, but the overall research workflow involves them all. We now estimate the combined -effects of these parametic choices. Starting from a single +effects of these parametric choices. Starting from a single \textbf{best-practice baseline} (trace stretching, a band matched to the target depth, a coda window well past the direct arrival, a 10-day stack, a long stable reference, coherence gating; the cross-cutting @@ -935,13 +966,16 @@ \section{The multiverse of processing choices}\label{sec:multiverse} deviation from best practice documented in the literature and measure the resulting error against the known truth (Fig.\textasciitilde{}\ref{fig:deviations}). The ranking is unambiguous: -relative to a best-practice RMS error of \(\sim\!0.02\,\%\), an -unwrapped phase estimator that cycle-skips is catastrophic, a moving -reference and a wrapped-phase MWCS each inflate the error roughly an -order of magnitude, an over-long stack and a late low-SNR window distort -the co-eruptive drop, while the frequency band (which mostly sets -precision here) and the gating are second-order. The joint inversion is -nearly as good as the fixed reference and preserves the trend. +relative to a best-practice RMS error of \(\sim\!0.03\,\%\), the +two-dimensionally-unwrapped WCS estimator is catastrophic here +(\(\sim\!50\times\) worse), a wrapped-phase MWCS inflates the error by +roughly an order of magnitude (\(\sim\!12\times\)), and DTW by +\(\sim\!6\times\); among the non-estimator choices a moving reference is +worst (\(\sim\!5\times\), and most distorts the recovered drop), while +stack length, coda window, and frequency band deviations are each more +modest (\(1.5\)--\(2.5\times\)). The joint-inversion reference stays +closest to the fixed-reference baseline among the reference-scheme +deviations. \begin{figure} \centering @@ -970,12 +1004,14 @@ \section{The multiverse of processing choices}\label{sec:multiverse} Attributing the variance of the outcome to each axis with a first-order (main-effect) sensitivity index (Fig.\textasciitilde{}\ref{fig:multiverse}b) shows that, for this -dataset, the \textbf{coda window and the stack length} control the RMS -error most, the \textbf{stack length} dominates the recovered drop -amplitude, and the estimator is third; the first-order indices sum to -well under one, so a large part of the spread is \emph{interaction} -between choices compounding. The appropriate object is therefore not a -single curve but a distribution of \dvv~over the processing choices. +dataset, the \textbf{coda window} controls the RMS error most, with the +estimator second and the stack length third; for the recovered drop +amplitude the order changes to \textbf{stack length} first, window +second, and estimator third. The first-order indices sum to well under +one in both cases (\(\sim\!0.66\) for RMS, \(\sim\!0.63\) for the drop), +so a large part of the spread is \emph{interaction} between choices +compounding. The appropriate object is therefore not a single curve but +a distribution of \dvv~over the processing choices. \begin{figure} \centering @@ -984,8 +1020,9 @@ \section{The multiverse of processing choices}\label{sec:multiverse} each coloured by its RMS error against the truth on a colourblind-safe scale (bright accurate, dark biased; the worst run off the clipped axis); the grey band is the 10--90\% inter-pipeline spread, widest at the velocity drop. (b) First-order variance attribution: which - choice controls the RMS error and the recovered drop. Window and stack length - dominate; the sub-unity sum signals strong interactions.} + choice controls the RMS error and the recovered drop. Window, estimator, and + stack length dominate, in that order for RMS and reordered for the drop; + the sub-unity sum signals strong interactions.} \label{fig:multiverse} \end{figure} @@ -1132,6 +1169,83 @@ \section{Toward deployment: a real-data retrospective the already-published \citet{Clements2023} CI.LJR result before any new claim is drawn from it. +After the correction, single-station \dvv~(NoisePy correlations, a +codameter 5-member ensemble, 2--4\(\,\)Hz, 2018--2019) validates against +the published \citet{Clements2023} product. The comparison is not +trivial to get right: the CD2022 90-day-comp product is a +\emph{trailing} 90-day stack, so it lags a centered-smoothed daily +series by about 45 days, and comparing without matching that smoothing +caps the correlation near 0.7 even on a real annual cycle +(Fig.\textasciitilde{}\ref{fig:realdata-validation}). We match by +applying the same trailing 90-day mean to the daily series, compare +demeaned --- the two products reference different epochs, and a constant +offset is bookkeeping, not error --- and exclude the first 150 days of +each station's series as reference burn-in. Under this matched +comparison, CI.LJR reaches \(r=0.990\) (681 overlapping days); CI.ARV +reaches \(r=0.66\)--\(0.92\) depending on the join method, reflecting +real 2018 data gaps; CI.RXH reaches \(r=0.68\), the weakest of the +three, on a station whose recovered \dvv~is nearly flat and whose signal +is low to begin with. + +\begin{figure} + \centering + \includegraphics[width=\textwidth]{realdata_1_validation.png} + \caption{Single-station \dvv\ at three CI stations, 2018--2019, against the + published Clements \& Denolle (2022) product (dashed, reference-shifted). + Daily \dvv\ (points) with the ensemble spread (shaded, epistemic) and + measurement error bars, and a 45-day-smoothed curve. The annotated $r$ is + computed on this 45-day-smoothed comparison, not the smoothing-matched + comparison reported in the text; matching CD2022's own trailing 90-day + window instead of a centered 45-day one raises $r$ at every station (to + 0.990/0.66--0.92/0.68 for LJR/ARV/RXH), which is why the two numbers + differ.} + \label{fig:realdata-validation} +\end{figure} + +Why validation quality differs so much by station is visible in the +correlation function itself +(Fig.\textasciitilde{}\ref{fig:realdata-interferograms}). CI.LJR shows a +stable, narrow coda near zero lag through the year. CI.RXH shows +multipath, several coherent arrivals spread across the full \(\pm 8\,\)s +of lag shown, and a visible shift in that pattern in April--May 2019, +consistent with a site change rather than a processing artefact. +CI.ARV's coherent energy is compact and concentrated near zero lag but +comparatively sparse, consistent with its higher scatter and join-method +sensitivity. + +\begin{figure} + \centering + \includegraphics[width=\textwidth]{realdata_2_interferograms.png} + \caption{Daily north-south/east-west cross-component correlation + functions, 2--4$\,$Hz, $\pm 8\,$s lag, 2019. CI.LJR's coda is stable and + narrow. CI.RXH shows multipath (several coherent bands across the full lag + range) and a shift in that pattern in April--May 2019. CI.ARV's coherent + energy is compact and near zero lag but comparatively sparse. The visual + difference is the reason validation quality differs by station in + Fig.~\ref{fig:realdata-validation}.} + \label{fig:realdata-interferograms} +\end{figure} + +As an optional supplement, +Fig.\textasciitilde{}\ref{fig:realdata-warmup} shows the ensemble's +warm-up behaviour on a separate 90-day smoke run at CI.LJR: \dvv\(\,\)is +undefined until enough history has accumulated for the moving-reference +member to compute a trailing reference, and each band's per-epoch +stretching correlation coefficient is reported alongside the recovered +series. + +\begin{figure} + \centering + \includegraphics[width=\textwidth]{realdata_3_warmup.png} + \caption{CI.LJR single-station \dvv, a separate 90-day smoke run (January--April + 2023), four frequency bands. The ensemble spread (shaded) and measurement + error bars are reported once at least one member is defined; the annotated + member count grows as the moving-reference member's warm-up period + elapses, reaching all 5 members later in the window. Bottom panel: + per-band stretching correlation coefficient.} + \label{fig:realdata-warmup} +\end{figure} + \section{Depth: propagating the measurement covariance through sensitivity kernels}\label{sec:depth} @@ -1146,13 +1260,13 @@ \section{Depth: propagating the measurement covariance through is set by the frequency band and the coda lapse time, it is not assumed} \citep{Obermann2013, Obermann2016} --- and, increasingly, on multi-band measurement as the means to resolve it -\citep{Takano2017, Feng2020, Mao2022}. The step from a set of per-band -\dvv~series to a depth profile is where reporting is least consistent: -many studies read a single band as a single depth, only a few invert -several bands against surface-wave sensitivity kernels, and the -propagation of the measurement error into the depth estimate is seldom -shown. The depth assignment is frequently the scientific claim itself ---- whether the change lies in the aquifer or the overlying soil, +\citep{Takano2017, Feng2020, Mao2022, Mao2025}. The step from a set of +per-band \dvv~series to a depth profile is where reporting is least +consistent: many studies read a single band as a single depth, only a +few invert several bands against surface-wave sensitivity kernels, and +the propagation of the measurement error into the depth estimate is +seldom shown. The depth assignment is frequently the scientific claim +itself --- whether the change lies in the aquifer or the overlying soil, whether coseismic softening is a shallow site response or slip on the fault at depth \citep{Rubinstein2005} --- so a depth reported without its uncertainty cannot support that claim. @@ -1322,6 +1436,30 @@ \section{Discussion}\label{sec:discussion} and lets a reader re-run a study's pipeline on the truth-known synthetic to see its bias before trusting it on data. +\textbf{Validate against something you did not generate.} Every +synthetic test in this paper passed before the discovery below, and that +is exactly the danger: a synthetic built under the same sign convention +as the estimator reading it will always agree, whether the convention is +physically correct or not. Testing codameter's estimators against a real +cross-network deployment (Section\textasciitilde{}\ref{sec:deployment}), +recovered \dvv~anticorrelated with the published Clements \& Denolle +(2022) product and with seasonal hydrology at three stations +(\(r=-0.69,-0.45,-0.40\)). Ground-truthing through the exact call path +made the cause obvious: imposing a \(+0.5\,\%\) velocity change returned +\(-0.50\,\%\). The synthetic generator and all seven estimators had +consistently used the stretch factor \(\varepsilon\) (positive for a +coda dilation, i.e.~a slowdown), not physical \dvv~(positive for a +speedup) --- internally coherent, so every synthetic-recovery test in +the sections above passed, but opposite to the sign convention the field +expects and to the published product it was compared against. The fix is +the convention boxed in Section\textasciitilde{}\ref{sec:intro}, shipped +as codameter v0.4.0. Internal consistency is not correctness: a pipeline +that only checks itself will confirm whatever convention it started +with, and only the comparison against an independently-produced result +caught this one. Full details and the ground-truthing procedure are in +codameter PR \#36 and its commit history on branch +\texttt{fix/dvv-sign-convention}. + \textbf{Propagate the covariance, do not truncate it.} The measurement covariance is not the end of the analysis but its first input. Section\textasciitilde{}\ref{sec:depth} sets out the next step of the @@ -1363,8 +1501,11 @@ \section{Estimator definitions}\label{app:estimators} Let \(r(t)\) and \(c(t)\) be the reference and current cross-correlations, band-passed to \([f_1,f_2]\) (central frequency \(f_c\), bandwidth \(B\)) and read over the coda window \(W=[t_1,t_2]\) -on one or both branches. We use the convention \(\delta t/t = -\,\dvv\), -and write \(\varepsilon\) for the recovered \dvv. +on one or both branches. Each method below estimates a trial stretch +factor \(\varepsilon\), the fractional dilation of the current coda +relative to the reference; physical \dvv~follows via the convention +boxed in Section\textasciitilde{}\ref{sec:intro}, +\(\dvv=-\varepsilon/(1+\varepsilon)\). \textbf{Trace stretching (TS).} For a trial stretch \(s\), define \(r_s(t)=r\!\left(t/(1+s)\right)\) and the windowed correlation @@ -1418,17 +1559,20 @@ \section{Estimator definitions}\label{app:estimators} \section{Aggregation and uncertainty conventions}\label{app:aggregation} For component \(k\) of a station pair, stretching yields the image -\(\mathrm{CC}_k(\varepsilon,t)\). - -\textbf{Approach A (average the \dvv).} -\(\varepsilon_k(t)=\arg\max_{\varepsilon} \mathrm{CC}_k(\varepsilon,t)\) -and the pair estimate is the (possibly weighted) mean -\(x_p(t)=\sum_k w_k\,\varepsilon_k/\sum_k -w_k\), with \(w_k=\max_{\varepsilon} \mathrm{CC}_k\) -(coherence-weighted) or \(w_k=1\) (unweighted). - -\textbf{Approach B (average the images).} -\(x_p(t)=\arg\max_{\varepsilon}\big[\tfrac1{N_c}\sum_k \mathrm{CC}_k(\varepsilon,t)\big]\). +\(\mathrm{CC}_k(\varepsilon,t)\); each per-component stretch factor +converts to physical \(\dvv_k\) via the boxed relation. + +\textbf{Approach A (average the per-component \dvv).} \(\dvv_k(t)\) +follows from \(\varepsilon_k(t)=\arg\max_{\varepsilon} +\mathrm{CC}_k(\varepsilon,t)\), and the pair estimate is the (possibly +weighted) mean \(x_p(t)=\sum_k w_k\,\dvv_k/\sum_k w_k\), with +\(w_k=\max_{\varepsilon} +\mathrm{CC}_k\) (coherence-weighted) or \(w_k=1\) (unweighted). + +\textbf{Approach B (average the images).} The pair stretch factor +\(\varepsilon_p(t)=\arg\max_{\varepsilon}\big[\tfrac1{N_c}\sum_k +\mathrm{CC}_k(\varepsilon,t)\big]\) converts to \(x_p(t)=\dvv_p(t)\) via +the same relation. \textbf{Network over pairs.} With pair weights \(W_p\) (e.g.~mean coherence), \begin{equation} diff --git a/src/codameter/deviations.py b/src/codameter/deviations.py index c2a184e..b42069d 100644 --- a/src/codameter/deviations.py +++ b/src/codameter/deviations.py @@ -431,6 +431,7 @@ def fig_deviation_ranking(rows=None): def fig_multiverse_full(mv=None): """The ultimate multiverse: every pipeline + the variance attribution.""" import matplotlib.pyplot as plt + from matplotlib.colors import Normalize if mv is None: mv = multiverse() @@ -444,7 +445,7 @@ def fig_multiverse_full(mv=None): # (a) fan of pipelines, coloured by RMS error with a colourblind-safe, # perceptually uniform sequential map (dark = accurate, bright = biased). order = np.argsort(-np.nan_to_num(rms)) - norm = plt.Normalize(np.nanpercentile(rms, 5), np.nanpercentile(rms, 95)) + norm = Normalize(np.nanpercentile(rms, 5), np.nanpercentile(rms, 95)) cmap = plt.cm.viridis_r for i in order: ax[0].plot(yrs, curves[i] * PCT, color=cmap(norm(rms[i])), lw=0.3, alpha=0.16) @@ -463,13 +464,10 @@ def fig_multiverse_full(mv=None): ) ax[0].plot(yrs, truth * PCT, color=C["truth"], lw=2.6, label="ground truth") ax[0].axvline(2.0, color="0.6", ls="--", lw=1) - # Clip tightly to the truth scale; the cycle-skipping pipelines run off-axis - # (that is the point — the colourbar flags them) but would otherwise swamp the + # Fixed, symmetric range: the cycle-skipping pipelines run off-axis (that is + # the point -- the colourbar flags them) but would otherwise swamp the # signal and make the panel unreadable. - span = (np.nanmax(truth) - np.nanmin(truth)) * PCT - ax[0].set_ylim( - (np.nanmin(truth) * PCT - 0.35 * span, np.nanmax(truth) * PCT + 0.35 * span) - ) + ax[0].set_ylim((-0.8, 0.8)) ax[0].set( xlabel="time (years)", ylabel="dv/v (%)", diff --git a/src/codameter/uq_bayes.py b/src/codameter/uq_bayes.py index d202852..e9a1025 100644 --- a/src/codameter/uq_bayes.py +++ b/src/codameter/uq_bayes.py @@ -409,11 +409,11 @@ def _fig_bayes(res, run): truth = run.truth sd_cd = np.sqrt(np.diag(res.Cd)) sd_post = np.sqrt(np.diag(res.mu_cov)) - fig = plt.figure(figsize=(7.2, 3.0), layout="constrained") - gs = fig.add_gridspec(1, 3, width_ratios=[1.5, 1.0, 1.1]) + fig = plt.figure(figsize=(7.2, 5.6), layout="constrained") + gs = fig.add_gridspec(2, 2, height_ratios=[1.0, 1.0]) # (a) ensemble + posterior + the two bands. - ax0 = fig.add_subplot(gs[0]) + ax0 = fig.add_subplot(gs[0, 0]) for k in range(run.members.shape[0]): ax0.plot(yrs, run.members[k] * 100, lw=0.5, color="0.7", alpha=0.6) if truth is not None: @@ -448,7 +448,7 @@ def _fig_bayes(res, run): ) # (b) the data covariance matrix. - ax1 = fig.add_subplot(gs[1]) + ax1 = fig.add_subplot(gs[0, 1]) vmax = float(np.percentile(np.diag(res.Cd), 85)) # robust to the warm-up spike im = ax1.imshow( res.Cd, @@ -472,7 +472,9 @@ def _fig_bayes(res, run): fig.colorbar(im, ax=ax1, fraction=0.046) # (c) time-dependent sigma_d(t) and the effective-sample-size collapse. - ax2 = fig.add_subplot(gs[2]) + # Full width on its own row: it carries four legend entries and was too + # horizontally squeezed sharing a row with (a) and (b). + ax2 = fig.add_subplot(gs[1, :]) ax2.plot( yrs, sd_cd * 100, color=C["volcano"], lw=1.6, label=r"$\sigma_d(t)$ (total)" )