Peng 2025 W-phase plus IDBP joint inversion

This directory is an independent implementation of the algorithm described by Peng et al. (2025), "Slip Modeling of Large Earthquakes by Joint Inversion of W-Phase and Back-Projected Images". It does not modify ../formal_ff_release/ or its result directories. The copied peng2025_ff package contains the data/QC/geometry helpers needed by this program.

Reproducible event workflow

For the organized us7000srb1 workflow, start with final_workflow_20260721/README.md. The result inventory is in RESULTS_INVENTORY_20260721.md; the deletion audit is in DELETE_CANDIDATES_20260721.md.

Use final_workflow_20260721/events/us7000srb1/results/ for main-model figures and final_workflow_20260721/events/us7000srb1/sensitivity/ for explicitly labeled TauP/MCCC, NP1/NP2, NT, and normalization tests.

The long command history below is retained as methodological notes. Some older commands intentionally refer to deleted development outputs and are not reproducible entry points. Use the numbered scripts in final_workflow_20260721/scripts/ for current runs.

What is implemented

The program follows the five steps in the paper.

1. It reads the selected Grond point-source W-phase CMT, using the effective centroid location, centroid time, full six-component moment tensor and source duration.
2. It reuses the formal Z/R observation preparation and station-level QC.
3. It reads an independent regional body-wave BP image and aggregates it onto the finite-fault grid and MTW time windows.
4. Its paper-faithful mode builds the approximate IDBP array-response kernel G_I from cumulatively stacked mainshock BP energy in the first 6 seconds. The paper does not search 4/6/8/10 snapshots. Kernel sensitivity is represented separately by comparing the mainshock early image, an aftershock image and a foreshock image, when those files are available. taup_replay is retained only as a non-paper diagnostic implementation.
5. The separate run_idbp.py module implements the EPSL 2016 image-level reference-event deconvolution d_I ~= G_I A_I with NNLS. It outputs radiators, source duration, propagation trend and reference-image resolution. This standalone A_I is a diagnostic/source-image product; it is not treated as slip.
6. build_physical_idbp_kernel.py constructs a station-level unit-source operator. G_I_TauP uses iasp91 P travel times as a physical baseline and replays a positive pulse through the BP mean/semblance stack. G_I_MCCC is an optional sensitivity product using bounded, station-wide group-QC shifts; it never introduces a free shift for each subfault or time window.
7. It solves the joint W-phase plus IDBP plus smoothing system with scipy.optimize.nnls, which is the Lawson-Hanson non-negative least-squares algorithm used in the paper.
8. It converts the solved coefficients into normalized moment density A^J and one scalar moment M0_j per time window.
9. It performs FFWI-M by holding A^J fixed and solving six moment-tensor components per time window. The tensor is constant over grid points and may vary between time windows, as in Eq. (6) of the paper.
10. It repeats the joint solve and FFWI-M update until the mean Kagan-angle change is below 10 degrees and total scalar-moment change is below 2%, or until --iterations is reached.

Equations implemented

For the Peng rapid mode, the first-six-second field is the approximation to the 3-D array response G_I; the observed mainshock BP image remains d_I. For the standalone EPSL mode, the reference-event image is used to form shifted columns of G_I, and NNLS solves for non-negative A_I. After the image and W-phase densities are normalized, the practical Eq. (1) factor is reported as X=1. A_I is never assigned a slip unit.

The joint matrix is assembled as three separately normalized blocks:

G_W / tr(G_W.T G_W)
gamma * G_I / tr(G_I.T G_I)
lambda * F / tr(F.T F)

The right-hand side uses the corresponding three denominators. The smoothing operator is exactly Eq. (5): the diagonal is -1, and each existing left, right, up or down neighbor has coefficient +1/4. It is not scaled by the physical subfault dimensions and is not boundary-renormalized.

The paper also states that G^I and d^I are normalized separately before this block scaling, but does not specify the norm in the main text or the supporting information. The solver therefore exposes two independent switches:

• --block-normalization trace|frobenius selects the denominator used for the three Eq. (4) blocks. trace is the literal tr(G.T G) form printed in the paper; frobenius uses its square root as a diagnostic alternative.
• --idbp-pre-normalization none|frobenius|sum_abs controls the separately normalized IDBP operator/data pair before the Eq. (4) denominator. none reproduces the earlier raw-operator path, frobenius uses the matrix Frobenius norm and data L2 norm, and sum_abs is a transparent sensitivity test. The exact paper norm remains undocumented, so no option is silently declared to be the unique paper baseline.

For the current NP2 run, raw trace scaling gives W-phase and IDBP traces of 1.306e6 and 5.844; with gamma=4.5, their effective Frobenius norm ratio is about 2.13e3. Consequently, gamma=4.5 cannot be interpreted as a literal 4.5-to-1 numerical block balance on the raw operator. Every output records the selected modes and both preprocessing scale factors.

The audit command is:

/opt/anaconda3/envs/grond/bin/python audit_joint_operator.py \
  --joint-dir final_workflow_20260721/events/us7000srb1/sensitivity/03_normalization/trace_operator_audit_NP2 \
  --gi-solution final_workflow_20260721/events/us7000srb1/sensitivity/01_kernel/physical_taup_NP2_NT6_STF6/physical_gi_solution.npz \
  --output-dir final_workflow_20260721/events/us7000srb1/sensitivity/03_normalization/trace_operator_audit_NP2

It verifies that the saved predictions reproduce the design matrices, that the BP observation flattening is exact, and that the active columns of the 90-column physical G_I are selected as subfault-major/time-window-minor. The current NP2 audit gives machine-precision errors for all these checks; the active W-phase source columns are present only for windows 1--4.

For each time window j, the non-negative joint coefficients are separated as

M0_j = m_unit * sum_i x_ij
A^J_ij = x_ij / sum_i x_ij
S^J_ij = A^J_ij * M0_j / (mu * S_q)

FFWI-M uses the fixed A^J weights to form the aggregate W-phase matrix and solves six unrestricted moment-tensor components per time window. The iteration diagnostics are:

Misfit = sum((d - syn)^2) / sum(d^2)
FD Index = (W1 - W2) / W1

Here W1 is the joint-inversion W-phase Misfit and W2 is the FFWI-M W-phase Misfit. FD Index is a frequency-dependence diagnostic, not a general quality score.

Configuration

The default peng2025_joint_config.json is a Philippines event template. The current template uses one broad time window for the first direct test. For a Peng-style MTW run, set the explicit MTW fields to the desired event-specific values, for example:

"nt": 6,
"window_duration_s": 20.0,
"window_offset_s": 10.0,
"stf_duration_s": 20.0

nt is the number of time windows. window_duration_s controls the BP image bin and triangular STF duration is specified independently by stf_duration_s. window_offset_s is the shift between adjacent windows and allows overlapping MTW windows. The legacy th_s field remains as a label and fallback for old one-window configurations; explicit paper-mode fields should always be provided for a multi-window run.

The CMT path can be set in JSON or overridden with --grond-cmt. The BP image must be an NPZ containing image, along_km, down_km, and source_time_s. The BP station CSV is required by the physical station-level kernel and the legacy taup_replay kernel.

Physical G_I validation before joint inversion

Run the theoretical baseline first. This example uses the NP2 hypocenter/fault-start geometry and origin time; the 18.27 s Grond CMT centroid shift is not inserted into the finite-fault time axis:

/opt/anaconda3/envs/grond/bin/python build_physical_idbp_kernel.py \
  --config peng2025_joint_config.json --plane-name NP2 \
  --source-lat 5.5918 --source-lon 125.0469 --source-depth-km 55.193 \
  --bp-image ../bodywave_backprojection/results/np2_bp_all_peng2025_waveform_aligned/bodywave_bp_image.npz \
  --formal-subfaults final_workflow_20260721/events/us7000srb1/inputs/formal_subfaults_depth40km_NP2.csv \
  --station-csv ../bodywave_backprojection/results/np2_bp_all_peng2025_waveform_aligned/prepared_bodywave_stations.csv \
  --output-dir final_workflow_20260721/events/us7000srb1/sensitivity/01_kernel/physical_taup_NP2_NT6_STF6

This writes physical_gi_solution.npz, physical_gi_summary.json, checkerboard_recovery.csv, leave_group_out_summary.csv and physical_gi_validation.png/pdf. The checkerboard is a noiseless self-consistency test of the operator and NNLS implementation; it is not an independent resolution proof. The leave-group-out test removes one of eight azimuth groups at a time and measures operator sensitivity. It is a station-geometry test, not a held-out observed-waveform prediction, because the archived BP NPZ does not contain station-wise images.

The script also solves the actual aggregated BP observation d_I ~= G_I A_I. In the current NP2 test, the TauP kernel gives an image Misfit of about 0.0095 and correlation about 0.972; this shows that the physical BP operator can reproduce the archived image after solving for relative A_I. It does not imply that the same kernel is already compatible with the W-phase data. That compatibility must be checked by the joint solve and is a separate test.

The constrained MCCC sensitivity run is separate:

/opt/anaconda3/envs/grond/bin/python build_physical_idbp_kernel.py \
  --config peng2025_joint_config.json --plane-name NP2 \
  --source-lat 5.5918 --source-lon 125.0469 --source-depth-km 55.193 \
  --bp-image ../bodywave_backprojection/results/np2_bp_all_peng2025_waveform_aligned/bodywave_bp_image.npz \
  --formal-subfaults final_workflow_20260721/events/us7000srb1/inputs/formal_subfaults_depth40km_NP2.csv \
  --station-csv ../bodywave_backprojection/results/np2_bp_all_peng2025_waveform_aligned/prepared_bodywave_stations.csv \
  --station-qc-csv ../bodywave_backprojection/results/np2_bp_all_peng2025_waveform_aligned_smooth/bodywave_bp_station_qc.csv \
  --static-correction-mode mccc --output-dir final_workflow_20260721/events/us7000srb1/sensitivity/01_kernel/mccc_static_correction_NP2_NT6_STF6

The formal baseline remains G_I_TauP unless the MCCC version improves held-out station predictions without degrading conditioning and cross-band/group stability. For the current waveform_linear BP operator, the TauP kernel is full rank but has condition number about 5.5e4; this is already a warning that the 15-subfault/6-window parameterization is ill-conditioned. The constrained MCCC sensitivity run reduced the condition number to about 4.9e3 and changed the aggregated BP-image Misfit from about 0.0094 to 0.0092. This is encouraging but not sufficient to promote MCCC to the primary kernel because its shifts were estimated from the same event; it remains a sensitivity product until independent held-out waveform tests are available.

After validation, either physical kernel can be passed to the joint solver:

/opt/anaconda3/envs/grond/bin/python run_peng2025_joint_inversion.py \
  --config peng2025_joint_config.json \
  --bp-image ../bodywave_backprojection/results/np2_bp_all_peng2025_waveform_aligned/bodywave_bp_image.npz \
  --kernels physical_taup \
  --idbp-solution final_workflow_20260721/events/us7000srb1/sensitivity/01_kernel/physical_taup_NP2_NT6_STF6/physical_gi_solution.npz \
  --gamma-bp 10 --lambda-smooth 0.3 --iterations 10 \
  --output-dir final_workflow_20260721/events/us7000srb1/results/03_final_model_NP2_depth40_gamma10_lambda0p3

Gamma-lambda search and current result

The first controlled search used the same NP2 hypocenter/fault-start geometry, origin time, six 9-s MTW windows, W-phase band 0.00167-0.01 Hz, and the same waveform_linear BP stack as the physical-kernel construction. It evaluated all 48 combinations of

gamma = 0, 0.1, 0.5, 1, 2, 4.5, 10, 20
lambda = 0.01, 0.1, 0.5, 1, 3, 10

with --iterations 1, separately for TauP and MCCC kernels. The exact search command is:

/opt/anaconda3/envs/grond/bin/python search_peng2025_parameters.py \
  --config peng2025_joint_config.json \
  --bp-image ../bodywave_backprojection/results/np2_bp_all_peng2025_waveform_aligned/bodywave_bp_image.npz \
  --kernels physical_taup \
  --idbp-solution results_physical_gi_tauP_np2_nt6/physical_gi_solution.npz \
  --bp-station-csv ../bodywave_backprojection/results/np2_bp_all_peng2025_waveform_aligned/prepared_bodywave_stations.csv \
  --plane-name NP2 --source-depth-km 55.193 \
  --anchor-mode hypocenter_fault_start --time-reference origin \
  --wphase-fmin 0.00167 --wphase-fmax 0.01 \
  --nt 6 --window-duration 9 --window-offset 9 --stf-duration 9 \
  --gamma-values 0,0.1,0.5,1,2,4.5,10,20 \
  --lambda-values 0.01,0.1,0.5,1,3,10 \
  --iterations 1 --output-dir results_gamma_lambda_search_physical_taup

For the normalization audit, keep the physical kernel, geometry, frequency band and MTW settings fixed and repeat a small grid with only the IDBP preprocessing changed:

for mode in none frobenius; do
  /opt/anaconda3/envs/grond/bin/python search_peng2025_parameters.py \
    --config peng2025_joint_config.json \
    --bp-image ../bodywave_backprojection/results/np2_bp_all_peng2025_waveform_aligned/bodywave_bp_image.npz \
    --bp-station-csv ../bodywave_backprojection/results/np2_bp_all_peng2025_waveform_aligned/prepared_bodywave_stations.csv \
    --kernels physical_taup \
    --idbp-solution results_physical_gi_tauP_np2_nt6/physical_gi_solution.npz \
    --plane-name NP2 --source-depth-km 55.193 \
    --anchor-mode hypocenter_fault_start --time-reference origin \
    --wphase-fmin 0.00167 --wphase-fmax 0.01 \
    --nt 6 --window-duration 9 --window-offset 9 --stf-duration 9 \
    --gamma-values 0.5,1,2,4.5,8 --lambda-values 0.1,1,10 \
    --block-normalization trace --idbp-pre-normalization "$mode" \
    --iterations 1 --output-dir "results_gamma_lambda_norm_${mode}"
done

The two grids are a normalization sensitivity experiment, not two independent physical models. Compare W-phase Misfit, BP-image Misfit/correlation, the recorded preprocessing scales, and the A_I/A^J spatial diagnostic before selecting a formal convention.

For the current Philippines NP2 test, the event-specific refined search used the Frobenius block diagnostic only to expose the Pareto trade-off; it did not assume the Nepal value gamma=4.5, lambda=1. With raw G^I and Frobenius block normalization, the coarse Pareto knee was near gamma=2.5, lambda=0.05. The refined front ranged from approximately (W Misfit, BP CC)=(0.414, 0.01) at gamma=1.5 to (0.922, 0.467) at gamma=6. The gamma=2.5, lambda=0.05 candidate was then tested with five FFWI-M iterations. Its final W-phase Misfit was 0.170, but the Kagan-angle and total-moment changes had not both met the paper's convergence thresholds. A gamma=4, lambda=0.05 comparison ended with W-phase Misfit 0.171 and BP correlation 0.434, also without convergence after five iterations. Repeating gamma=4 with lambda=0.1 and 0.3 gave final W-phase Misfits of 0.171 and 0.173, and BP correlations of 0.434 and 0.440; neither converged. This small change in lambda therefore does not explain the FFWI-M instability.

The residual audit also found no edge-localized residual: edge cells contain about 77% of the residual energy, but they are 12 of the 15 cells and have lower mean residual energy than the three interior cells. In the pre-fix diagnostic run, the late-window residual increase was accompanied by a structural source-column mask: only source windows 1--4 were generated for the W-phase matrix, while the BP vector still contained six windows. This was not evidence that the W-phase waveform ended early. The MTW builder has now been corrected so that rupture-velocity bounds set onset delays but do not delete later independent MTW windows. The corrected NT=6 and a reduced NT=4 run must therefore be compared before drawing a physical conclusion about late-window residuals. These are diagnostic candidates, not yet a final published model.

The MCCC run uses the same command with physical_mccc and results_physical_gi_mccc_np2_nt6/physical_gi_solution.npz. The outputs are in results_gamma_lambda_search_physical_taup/, results_gamma_lambda_search_physical_mccc/, and the comparison plots are in results_gamma_lambda_comparison/.

The current search does not show a parameter combination that fits both data types well:

kernel	selection	W-phase Misfit	BP correlation
TauP	maximum BP correlation with W Misfit <= 1 (gamma=4.5, lambda=1)	0.905	0.494
TauP	W-phase-only minimum (gamma=0, lambda=0.01)	0.203	-0.624
MCCC	maximum BP correlation (gamma=0.1, lambda=0.01)	21.567	0.509
MCCC	W-phase-only minimum (gamma=0, lambda=0.01)	0.203	-0.634

These values are a diagnostic result, not a failed numerical solve. They show that, with the current BP image, physical operator, geometry, MTW parameterization, and block normalization, the two data sets impose incompatible spatial constraints. In particular, the independent fit d_I ~= G_I A_I with correlation about 0.973 cannot be used to predict that the joint coefficient vector A^J will fit the BP image equally well. The independent A_I fit has more freedom and is a different inverse problem.

The Pareto plots and CSV files are generated with:

/opt/anaconda3/envs/grond/bin/python analyze_gamma_lambda_search.py \
  --search-dirs results_gamma_lambda_search_physical_taup results_gamma_lambda_search_physical_mccc \
  --output-dir results_gamma_lambda_comparison

For the short-listed TauP point and the MCCC maximum-BP point, ten-step FFWI-M runs were also completed. FFWI-M converged in both cases by the paper's stopping tests and reduced the W-phase-only FFWI-M Misfit to about 0.187; however, the joint-solve BP correlations remained about 0.494 and 0.509, respectively. This means that FFWI-M updates the moment tensors and improves the W-phase explanation, but does not remove the underlying W-phase/BP spatial conflict. The detailed outputs are in results_physical_taup_ffwi10_g4p5_l1/ and results_physical_mccc_ffwi10_g0p1_l0p01/.

The next diagnostic should therefore compare the independent A_I field with the joint A^J field and audit the column normalization, trace scaling, time-window mapping, and BP observation/operator convention. More extreme gamma-lambda values should not be interpreted as a physical improvement until that comparison is understood.

NP1 versus NP2 geometry diagnostic

The NP1 control experiment was regenerated from the same 134 stations, the same hypocenter/fault-start anchor, source depth 55.193 km, origin time, NT=6, 9-s windows, and gamma=4.5, lambda=1. NP1 was not evaluated by reusing the NP2 BP image: its BP projection, subfault coordinates, and TauP physical kernel were rebuilt independently.

The independent physical-kernel fits were similar:

             G_I condition       d_I ~= G_I A_I Misfit    correlation
NP1             4.80e3                   0.0319              0.9713
NP2             5.52e4                   0.0094              0.9726

However, the joint solutions did not improve for NP1:

             W-phase Misfit       joint BP Misfit       joint BP correlation
NP1              1.105                 0.573                    0.473
NP2              1.112                 0.344                    0.494

The normalized difference between the independent A_I and joint A^J fields was 0.559 for NP1 and 0.280 for NP2. Under this controlled test, NP2 is more compatible with the joint BP constraint even though its physical kernel is more ill-conditioned. This is not proof that NP2 is the true fault plane; it only means that changing to NP1 does not remove the present W-phase/BP incompatibility.

The complete diagnostic products are in results_joint_spatial_diagnostic_np1_np2/:

• joint_spatial_diagnostic.png: observed d_I, independent G_I A_I, joint G_I A^J, and both residual maps;
• A_I_A_J_comparison.png: independent and joint density fields and their difference for all six MTW windows;
• joint_window_diagnostics.png and joint_window_diagnostics.csv: BP residual/correlation by data window and the change in W-phase Misfit when a source window is removed;
• joint_spatial_diagnostic.json: machine-readable summary.

The diagnostic can be regenerated with:

/opt/anaconda3/envs/grond/bin/python diagnose_joint_spatial_consistency.py \
  --np1-joint results_physical_taup_joint_np1_nt6_diag \
  --np1-gi results_physical_gi_tauP_np1_nt6/physical_gi_solution.npz \
  --np2-joint results_physical_taup_joint_np2_nt6_diag \
  --np2-gi results_physical_gi_tauP_np2_nt6/physical_gi_solution.npz \
  --output-dir results_joint_spatial_diagnostic_np1_np2

The window test is diagnostic only. In particular, windows 5-6 have no active joint columns in this run and therefore cannot be interpreted as evidence that those physical time intervals contain no source process.

Run commands

Always run from this directory with the Grond environment:

cd /Users/ansel/Documents/New\\ project/grond/Philippines_eq/us7000srb1_wphase/formal_ff_peng2025_joint
export PYTHONPATH=.

/opt/anaconda3/envs/grond/bin/python run_peng2025_joint_inversion.py \
  --config peng2025_joint_config.json \
  --bp-image ../bodywave_backprojection/results/np2_bp_all_first_pass/bodywave_bp_image.npz \
  --dry-run

One joint solve without FFWI-M:

/opt/anaconda3/envs/grond/bin/python run_peng2025_joint_inversion.py \
  --config peng2025_joint_config.json \
  --bp-image ../bodywave_backprojection/results/np2_bp_all_first_pass/bodywave_bp_image.npz \
  --gamma-bp 4.5 --lambda-smooth 1.0 --iterations 1 \
  --output-dir results_np2_joint_iter1

Full iterative run using cumulative BP energy from the first six seconds:

/opt/anaconda3/envs/grond/bin/python run_peng2025_joint_inversion.py \
  --config peng2025_joint_config.json \
  --bp-image ../bodywave_backprojection/results/np2_bp_all_first_pass/bodywave_bp_image.npz \
  --idbp-kernel first6_empirical --idbp-kernel-duration-s 6 \
  --gamma-bp 4.5 --lambda-smooth 1.0 --iterations 5 \
  --output-dir results_np2_joint_peng2025

For an independent image to define the IDBP kernel, use --idbp-kernel-image /path/to/kernel_image.npz; the image passed to --bp-image remains the observed d_I.

Independent EPSL 2016 IDBP stage

The reference event must be small enough to approximate a point source and must be processed with the same stations, filter, source grid, travel-time model and time convention as the mainshock. The program solves

d_I ~= G_I A_I

with NNLS and writes idbp_solution.npz, idbp_radiators.csv, idbp_summary.json and idbp_diagnostics.png:

/opt/anaconda3/envs/grond/bin/python run_idbp.py \
  --main-image /absolute/path/mainshock/bodywave_bp_image.npz \
  --reference-image /absolute/path/reference/bodywave_bp_image.npz \
  --time-step-s 2 \
  --potential-threshold 0.8 \
  --output-dir results_idbp_np2

Use the reference-image-derived kernel in Peng Eq. 4 with:

/opt/anaconda3/envs/grond/bin/python run_peng2025_joint_inversion.py \
  --config peng2025_joint_config.json \
  --bp-image /absolute/path/mainshock/bodywave_bp_image.npz \
  --idbp-solution results_idbp_np2/idbp_solution.npz \
  --idbp-kernel standalone_idbp \
  --gamma-bp 4.5 --lambda-smooth 1 --iterations 5 \
  --output-dir results_np2_joint_with_reference_idbp

This uses the reference-image-derived G_I while retaining the mainshock BP image as d_I.

Parameter search

Use the search driver for a coarse scan. It writes one complete result directory per candidate and can resume completed candidates:

/opt/anaconda3/envs/grond/bin/python search_peng2025_parameters.py \
  --config peng2025_joint_config.json \
  --bp-image ../bodywave_backprojection/results/np2_bp_all_first_pass/bodywave_bp_image.npz \
  --gamma-values 0,0.5,1,2,3,4.5,6,8,10 \
  --lambda-values 0,0.05,0.3,1,3 \
  --kernels first6_empirical --kernel-images 6 \
  --iterations 1 --resume \
  --output-dir results_peng2025_parameter_search

The command above searches only gamma and lambda while using the fixed six-second mainshock kernel. To compare the paper's three kernel sources, provide independent BP image files for the aftershock and foreshock:

  --kernel-source-map \
    "aftershock=/absolute/path/aftershock_bp_image.npz,foreshock=/absolute/path/foreshock_bp_image.npz"

The --bp-image argument remains the mainshock observation d_I; the source map changes only the image used to construct the IDBP kernel. Do not interpret the optional --idbp-kernel-images override as a paper parameter search.

The default ranking minimizes W-phase Misfit. This is the paper-compatible primary criterion for the waveform. To exclude candidates that do not reproduce the BP image at all, add for example --bp-correlation-min 0.2. Alternatively, --metric joint_data_misfit ranks the trace-normalized W plus gamma-weighted BP data residual. The search does not hide the trade-off: parameter_search_results.csv retains W-phase Misfit, BP Misfit, combined data Misfit, BP correlation, VR, FD Index and convergence status for every case. For a paper-faithful search, omit --bp-correlation-min and use the reported W-phase, BP-image and combined Misfit curves to select the trade-off. The correlation option is only an optional diagnostic filter and is not a criterion from Peng et al. After selecting the best coarse region, run a narrower grid with smaller parameter spacing and then rerun the best candidate with --iterations 5 or more. This is computational staging only: every finalist is still evaluated with the same in-loop FFWI-M convergence test, and kernel source is reported as a separate sensitivity comparison rather than mixed with the early-image count.

Outputs

The output directory contains:

• joint_summary.json: complete parameter provenance, per-iteration metrics, Kagan angles, scalar moments, Misfit and FD Index.
• joint_solution.npz: raw NNLS coefficients x, normalized A^J, per-window M0_j, moment tensors and subfault coordinates.
• joint_subfaults.csv: fault-start coordinates, A^J, per-window moments, total moment and slip from S^J_ij.
• joint_predictions.npz: observed and synthetic W-phase/BP vectors plus their Misfit values. idbp_design_matrix and bp_data_raw preserve the physical input; idbp_design_matrix_solver and bp_data are the exact preprocessed arrays used by NNLS.
• joint_iteration_history.npz: the (A^J) density and moment-tensor fields saved after every joint solve, for testing spatial patch movement across FFWI-M iterations.
• bp_observation.npz: normalized finite-grid d_I.
• bp_array_response_kernel.npz: raw G_I used as the physical input. In first6_empirical mode this is Peng's six-second response approximation; in standalone_idbp mode it is projected from the EPSL reference-event image. The response amplitudes are preserved in this file; the solver-side global preprocessing and Eq. 4 block normalization are recorded in joint_summary.json. They are not normalized independently column by column.
• idbp_solution.npz, idbp_radiators.csv, idbp_summary.json and idbp_diagnostics.png: products of the independent EPSL IDBP stage.
• wphase_qc_all.csv and wphase_qc_kept.csv: the exact W-phase QC audit.
• component_weight_summary.csv: Z/R balancing information.

The residual concentration audit is run with:

/opt/anaconda3/envs/grond/bin/python audit_joint_residuals.py \
  --joint-dir results_np2_event_candidate_g4_l0p05_ffwi5 \
  --output-dir results_np2_event_candidate_g4_l0p05_ffwi5/residual_audit \
  --n-along 5 --n-down 3

It reports residual energy per cell and per MTW window, together with the scalar-moment contribution of every MTW window from joint_solution.npz. Edge-cell residual energy must be compared with the number of edge cells; an edge fraction near the edge-cell fraction is not evidence of an edge-localized error.

MTW column audit and corrected NT=4/6 comparison

The original MTW builder incorrectly removed a source-window column when the window center fell later than the slow-rupture-time estimate. This was not a W-phase data-duration test. It created a dimensionally inconsistent situation: the BP vector retained all requested windows, while the W-phase operator had only early-window columns. The builder now keeps every requested independent MTW window. The rupture-velocity setting supplies the onset delay; it does not delete later MTW basis functions.

For the NP2, 30x20 km, 9-s-window diagnostic, the pre-fix run had 45 W-phase columns (15 subfaults by only 3--4 available windows) and zero fitted moment in windows 5--6. After the correction, NT=6 has 90 columns, 15 per window. A matched NT=4 run has 60 columns, also 15 per window. The corrected first-pass results with gamma=4, lambda=0.1, Frobenius block normalization, and five FFWI-M iterations are:

run	W Misfit	BP Misfit	BP CC	final FFWI-M W Misfit	converged
NT=4	0.518	0.029	0.939	0.165	no
NT=6 corrected	0.604	0.043	0.873	0.164	no

The NT=4 result is not automatically preferred: it has fewer temporal degrees of freedom. The corrected NT=6 result distributes nonzero moment over all six windows (approximately 14.7%, 1.7%, 19.2%, 29.8%, 14.5%, and 20.1%). These results establish that the previous zero late-window moment was an implementation artifact, but they do not yet establish FFWI-M convergence or the physical uniqueness of six 9-s windows.

A corrected nine-point local search over gamma={2,4,6} and lambda={0.05,0.1,0.3} was then run with fixed NP2 geometry and one mechanism iteration. The W-phase-priority end was gamma=2, lambda=0.05 (W Misfit 0.514, BP CC 0.570). The BP-compatible branch was gamma=4, lambda=0.05--0.3 (W Misfit about 0.923--0.926, BP CC about 0.879--0.880). These values define an event-specific trade-off and are not the Nepal parameters copied into this event.

Five FFWI-M iterations were run for two representative candidates. The gamma=2, lambda=0.05 model reached a final W-phase Misfit of 0.169 but its BP CC fell to 0.530 and the Kagan changes remained large. The gamma=4, lambda=0.1 model reached a final W-phase Misfit of 0.164 while retaining BP CC 0.873; it also did not meet both the 10-degree Kagan and 2-percent scalar-moment convergence criteria within five iterations. At this stage, gamma=4, lambda=0.1 is the better diagnostic compromise, not a proven unique optimum. More iterations should only be run after checking whether the mechanism sequence is oscillatory; the spatial and temporal residual audits should remain part of model selection.

The optional analyze_peng2025_stability.py reports active-area overlap and centroid/peak scatter across a selected candidate set. These are supplementary sensitivity descriptions only. They are not used as formal acceptance criteria. The paper's formal iteration criteria remain mean Kagan-angle change below 10 degrees and total scalar-moment change below 2 percent.

FFWI-M damping diagnostic

FFWI-M now accepts --ffwi-update-damping alpha. Before this test, the initial Grond CMT tensors (stored in N m) and FFWI-M coefficients (stored relative to solver.m_unit_nm) were compared in different units, which could make the first moment-change value appear as 100%. The internal tensors are now converted to a common dimensionless solver unit before Kagan and moment diagnostics; reported scalar moments remain in N m.

For the NP2 gamma=4, lambda=0.1 diagnostic, the undamped five-iteration control had final BP CC 0.873 and FFWI-M W Misfit 0.164. A ten-iteration run with alpha=0.5 had final BP CC 0.856 and FFWI-M W Misfit 0.167. Its mean Kagan angle reached 17.4 degrees, but total moment changes still oscillated between roughly 30% and 38% in later iterations. Damping therefore reduces some directional jumps but does not establish convergence for this event; it is retained as a sensitivity option, not used in the formal model.

Scale and moment-tensor decomposition audit

The joint solver currently separates each MTW coefficient vector into a normalized A^J and a window scalar moment M0_j; the sum of A^J in every window is one by construction. The new audit records these sums over the iteration history, so the proposed inverse A^J-versus-M0 scale seesaw can be tested rather than assumed. For the present damping run, all six sums stay at one to numerical precision. Thus the observed moment oscillation is not caused by an unconstrained overall A^J scale in this implementation.

The audit also decomposes every FFWI-M tensor into ISO, DC, CLVD and deviatoric fractions using Pyrocko's standard decomposition. In the current full six-component run, several windows show ISO fractions of roughly 0.3--0.65 and CLVD fractions up to roughly 0.45, and these fractions vary between iterations. This is evidence that the unconstrained six-component FFWI-M update is using non-double-couple degrees of freedom; it is not yet proof that those components are physical. The outputs are moment_tensor_decomposition.png/pdf and the corresponding JSON fields in audit_joint_residuals.py.

The next model test should therefore be an explicitly labelled zero-trace/deviatoric FFWI-M sensitivity run. It must be compared with the full-tensor run using W-phase Misfit, BP CC, Kagan/moment stability, and the spatial slip pattern. The deviatoric constraint should not become the formal model unless it improves stability without materially changing the robust source features; it is not part of the current baseline.

That matched test has now been run for the NP2 gamma=4, lambda=0.1 case (six 9-s windows, origin time, 55.193 km fault-start anchor, TauP G^I). The full-tensor control gives W-phase Misfit 0.604, BP CC 0.873 and FFWI-M W Misfit 0.164 after five iterations. The zero-trace test gives W-phase Misfit 0.710, BP CC 0.896 and FFWI-M W Misfit 0.156. Its BP fit improves slightly, but the W-phase joint fit worsens by about 17%, and the Kagan sequence (56.4, 60.1, 37.5, 32.0 degrees) still does not converge; total moment changes also remain as large as 25.8%. Therefore zero trace removes the ISO leakage by construction but does not solve the FFWI-M instability and is not accepted as the formal model. It remains a documented sensitivity result. The output is results_np2_event_nt6_g4_l0p1_deviatoric_iter5/, with the decomposition and residual audit under residual_audit/.

Use the optional mode explicitly when reproducing this test:

--ffwi-tensor-space deviatoric

The default is --ffwi-tensor-space full, which preserves the six-component FFWI-M parameterization used by the formal Peng-style implementation.

FFWI-M SVD and inactive-window guard

The ten-iteration full-tensor baseline was used to test whether the moment oscillation came from a generally rank-deficient FFWI-M matrix. Five of the six MTW windows remained full rank with local condition numbers of about 8--11. In one iteration, however, the second MTW window had sum(A^J)=0; its six-column FFWI-M block consequently had rank zero and infinite condition number. This is an identifiability failure caused by the joint NNLS solution assigning zero density to that window, not a failure of the station Green functions themselves. The other windows remained full rank.

The solver now records global and per-window singular values, condition numbers, effective ranks and active-window flags in joint_iteration_history.npz and joint_summary.json. The residual audit also writes ffwi_svd_audit.png/pdf. When a window is exactly inactive, its FFWI-M update is skipped and its previous tensor is retained; this prevents an undefined zero-sensitivity least-squares solution from contaminating the next iteration. This is an identifiability guard, not a manually imposed physical time-window freeze. The guard does not make the present model converged: the ten-iteration NP2 run still has substantial Kagan/moment variation, so it is not accepted as evidence of a unique time-dependent mechanism.

An earlier NT=5 versus NT=6 comparison used 9-s BP aggregation windows but left stf_duration_s=32 s. Those outputs are retained only as a historical diagnostic and must not be used to decide whether a 9-s triangular W-phase source-time function is appropriate.

The valid Peng-style comparison uses the same 9-s duration for the BP window, the window offset, and the triangular W-phase STF (window_duration=9, window_offset=9, stf_duration=9). The independently rebuilt NT=5 kernel is 75 x 75, full rank, and fits the independent image with correlation 0.976 (NT=6: 90 x 90 and 0.973). After ten iterations, NT=5 gives W-phase Misfit 0.736 and BP CC 0.853; NT=6 gives 0.767 and 0.851. Both runs make the second 9-s window inactive after the fifth iteration, while the other windows remain active. Their outputs are results_np2_event_nt5_g4_l0p1_stf9_full_iter10_guard/ and results_np2_event_nt6_g4_l0p1_stf9_full_iter10_guard/.

This result does not support the claim that the sixth window alone is too many: reducing NT from six to five does not remove the same second-window active-set transition, and the final integrated spatial patterns are similar. It does show that the current data/model combination cannot stably identify that 9-s time window. The likely next test is therefore a coarser, matched MTW parameterization (for example 15-s or 20-s duration and offset), followed by the same SVD, active-window, Kagan-angle and scalar-moment audits. A 9-s parameterization should not be selected as the formal time resolution merely because it is used in Peng et al. (2025); it is retained as a Peng-consistent sensitivity test until the coarser-window comparison is made.

No-rupture-speed control

Peng et al. state that the main Nepal joint inversion does not impose a rupture-speed constraint; the 2.7 km/s result is an additional sensitivity test. The previous runs in this section used a legacy distance/vr_max onset delay and therefore were not an exact reproduction of that main-mode timing. The code now exposes rupture_time_mode:

• none (default): source time is j * window_offset_s for every subfault, matching the unconstrained MTW interpretation;
• fastest: applies the legacy distance / vr_max delay and is retained only for an explicit rupture-speed sensitivity test.

With rupture_time_mode=none, matched 9-s STF tests at the current NP2 geometry, TauP kernel, origin-time anchor, gamma=4.0, lambda=0.1, and Frobenius block normalization give:

model	joint W Misfit	joint BP CC	FFWI-M W Misfit	FFWI-M convergence
NT=5	0.448	0.869	0.138	yes, iteration 4
NT=6	0.687	0.807	0.133	yes, iteration 4

All five or six windows remain active in these no-speed runs. This reverses the earlier interpretation: the zero second window and non-convergence were caused mainly by the extra rupture-time delay, not by a demonstrated failure of 9-s MTW itself. The numbers above are still a diagnostic comparison, not the final Philippine parameter selection, because gamma, lambda, and the block-normalization convention must be searched together under the same no-speed timing mode.

The 450% weighting reported by Peng et al. corresponds to gamma=4.5 only under the paper's own trace-normalization convention. It is a Nepal-case choice, not a universal Philippine value; it should be tested rather than adopted as a prior. In particular, a gamma value from the Frobenius-normalized diagnostic runs is not numerically interchangeable with a trace-normalized 450% value.

Trace-normalized Philippine parameter search

The first trace-normalized, no-rupture-speed NT=6 search tested 56 (gamma, lambda) combinations. All candidates completed successfully. After the one-iteration screening, BP correlation reached a plateau near gamma=4.5; increasing gamma to 6--10 improved the correlation by less than 0.001. A five-iteration FFWI-M check of representative candidates found that all tested models converged and had nearly identical FFWI-M W-phase Misfits (about 0.132--0.134), while the BP correlation was already saturated for gamma >= 4.5.

The current formal candidate is therefore:

rupture_time_mode = none
block_normalization = trace
gamma = 4.5
lambda = 1.0
NT = 6
window_duration = window_offset = STF_duration = 9 s

For this candidate, BP CC is approximately 0.972, all six windows remain active, and the FFWI-M W-phase Misfit is approximately 0.133. The local lambda=0.1 sensitivity at the same gamma gives BP CC approximately 0.973 and FFWI-M W-phase Misfit approximately 0.135. This small difference does not justify treating lambda=0.1 as a uniquely better value; both should be retained in the stability comparison. The one-iteration search CSV is in results_np2_trace_novr_gamma_lambda_nt6_iter1/, and the FFWI-M screening results are in results_np2_trace_novr_gamma_lambda_nt6_iter5/.

Stable first-six-second IDBP path

When no suitable foreshock or aftershock reference image is available, the current Philippines workflow uses the station-subset/frequency-tested waveform-linear BP products as an auditable mainshock ensemble. Each product is normalized by its total image energy, and the median is used as the robust mainshock image. The cumulative first six seconds of this image is then used as Peng et al.'s rapid array-response approximation.

Construct and validate the formal-grid (G^I) and (A^I) with:

/opt/anaconda3/envs/grond/bin/python run_stable_idbp.py \
  --stability-root ../bodywave_backprojection/results/np2_bp_stability_waveform_linear \
  --formal-subfaults results_peng2025_bp05_2hz_hypocenter_nt6/joint_subfaults.csv \
  --nt 6 --window-duration-s 9 --window-offset-s 9 \
  --kernel-duration-s 6 \
  --output-dir results_stable_idbp_np2_nt6

The program solves (d_I \simeq G_I A_I) with NNLS, reports the image Misfit/correlation, matrix rank and condition number, and performs a leave-one-BP-product-out validation. A_I is relative radiated-energy density, not slip. The output idbp_solution.npz contains the raw formal-grid G_i, d_i and a_i arrays. A high same-sample fit is insufficient; the cross-validation metrics must also be reported.

The validated matrix can be passed explicitly to the joint solver:

/opt/anaconda3/envs/grond/bin/python run_peng2025_joint_inversion.py \
  --config peng2025_joint_config.json \
  --bp-image results_stable_idbp_np2_nt6/stable_bp_image.npz \
  --idbp-solution results_stable_idbp_np2_nt6/idbp_solution.npz \
  --idbp-kernel stable_first6 \
  --plane-name NP2 --anchor-mode hypocenter_fault_start \
  --source-depth-km 55.193 \
  --wphase-fmin 0.00167 --wphase-fmax 0.01 \
  --nt 6 --window-duration 9 --window-offset 9 --stf-duration 9 \
  --gamma-bp 4.5 --lambda-smooth 1 --iterations 1 \
  --output-dir results_stable_first6_joint_np2_nt6_iter1

The stable_first6 mode uses the stored (G^I) exactly and checks that its dimensions match the finite-fault grid and MTW count. It is intentionally separate from first6_empirical, which rebuilds a kernel directly from one input BP image. The current validation found a full-rank (90\times90) matrix, so the joint result must still be assessed with held-out image fits, W-phase Misfit, FFWI-M Kagan/moment convergence and geometry checks. A boundary peak in the six-second response is recorded as a warning rather than treated as proof of a physical rupture endpoint.

The current NP2 diagnostic output is results_stable_first6_joint_np2_nt6_full_v2/. With gamma=4.5, lambda=1, six 9-s windows and the 0.00167--0.01 Hz W-phase band, the first joint solution has W-phase Misfit 0.894, BP Misfit 0.240 and BP correlation 0.237. FFWI-M reduces the W-phase Misfit to 0.262 while the BP correlation remains 0.237; the update converges with an average Kagan-angle change of approximately 0.026 degrees and total scalar moment change of approximately 0.0006%. This is a valid diagnostic result, but not yet a jointly consistent final model: FFWI-M improves the W-phase mechanism fit without resolving the BP/W-phase spatial conflict.

Interpretation limits

The first-six-second field is the rapid array-response approximation stated in Peng et al. It is not a direct conversion of BP energy into slip. The BP term constrains a spatial-temporal image, while the W-phase term constrains long-period waveform moment release. A low joint W-phase VR or a high Misfit must therefore be investigated together with gamma, lambda, the selected IDBP image/kernel, MTW timing, and the separate W-phase/BP residuals.

The paper's published gamma and lambda values are not automatically portable to another event. Because every data/operator block is normalized by its own trace denominator here, they must be searched for the present event and station set.

Depth search

search_peng2025_depth.py recomputes the matching BP image, six-snapshot mainshock Kernel and joint inversion at every reference depth. It keeps the station list, plane, grid, gamma and lambda fixed and writes separate W, BP, joint and FFWI-M metrics. The CSV also reports the shallowest edge of the full fault plane; a positive subfault-center depth alone does not guarantee that the whole plane remains below the surface.

/opt/anaconda3/envs/grond/bin/python search_peng2025_depth.py \
  --joint-config peng2025_joint_config.json \
  --bp-config ../bodywave_backprojection/bodywave_backprojection_config_all.json \
  --plane-name NP1 --depths 20,25,30,35,40,45,50,55 \
  --station-csv ../bodywave_backprojection/results/mainshock_np1_depth21p6_common128/prepared_bodywave_stations.csv \
  --gamma-bp 4.5 --lambda-smooth 1 --iterations 1 \
  --skip-invalid-geometry --output-dir results_depth_search_np1_coarse

Use --iterations 1 for a coarse scan and rerun Pareto candidates with --iterations 5. The one-iteration table uses the Joint W Misfit as a proxy for the unavailable FFWI-M final Misfit, and marks this in ffwi_m_available.

Independent depth diagnostics

diagnose_slab2_geometry.py compares subfault centers with the local Slab2 depth, dip, and strike grids. The files in this project are comma-delimited lon,lat,value grids without headers; NaN cells are ignored. For example:

/opt/anaconda3/envs/grond/bin/python diagnose_slab2_geometry.py \
  --joint-config peng2025_joint_config.json \
  --plane-name NP2 --source-depth-km 35.5 \
  --slab-depth-file cot_slab2_dep_02.24.18.xyz \
  --slab-dip-file cot_slab2_dip_02.24.18.xyz \
  --slab-strike-file cot_slab2_str_02.24.18.xyz \
  --output-dir results_slab2_np2_d35p5

Its nearest 3-D distance is only a proxy unless the input is a surface mesh with normals. The Slab2 interface depth must not be silently substituted for the earthquake hypocentral or W-phase centroid depth.

phase_depth_scan.py uses TauP P, pP and sP predictions and reports station- median local energy/SNR over candidate depths. Phase windows overlapping the direct P window are excluded. It is a screening diagnostic, not a full depth-phase inversion; the output must be reviewed together with waveform plots and phase-pick quality.

Kinematic search

After choosing a depth/plane candidate, use the small kinematic search. The first pass should keep 40 s MTW windows and compare NT=1,2,3 plus moderate rupture-velocity bounds:

/opt/anaconda3/envs/grond/bin/python search_peng2025_kinematics.py \
  --config peng2025_joint_config.json \
  --bp-image /absolute/path/to/matching/bodywave_bp_image.npz \
  --bp-station-csv /absolute/path/to/prepared_bodywave_stations.csv \
  --plane-name NP1 --source-depth-km 25 \
  --anchor-mode centroid_centered --time-reference origin \
  --nt-values 1,2,3 \
  --vr-windows 1.5:3.5,2:3.5,2:4.5,2:5 \
  --gamma-bp 4.5 --lambda-smooth 1 \
  --window-duration 40 --window-offset 40 --stf-duration 80 \
  --iterations 1 --output-dir results_kinematic_search_np1_depth25

Use --iterations 5 only for the short list. Increasing NT or reducing the window below the W-phase time resolution is a model-complexity test, not an automatic improvement; inspect convergence, total moment and the BP residual alongside W-phase Misfit.

For the catalog-hypocenter experiment, use the matching hypocenter-anchored BP image and make the convention explicit:

/opt/anaconda3/envs/grond/bin/python search_peng2025_kinematics.py \
  --config peng2025_joint_config.json \
  --bp-image ../bodywave_backprojection/results/np2_hypocenter_anchor/bodywave_bp_image.npz \
  --bp-station-csv ../bodywave_backprojection/results/np2_hypocenter_anchor/prepared_bodywave_stations.csv \
  --plane-name NP2 --source-depth-km 55.193 \
  --anchor-mode hypocenter_fault_start --time-reference origin \
  --nt-values 2,3 --vr-windows 2:3.5,2:4.5 \
  --window-duration 30 --window-offset 15 --stf-duration 60 \
  --iterations 1 --output-dir results_kinematic_hypocenter_coarse

Source-anchor modes

The code supports two explicit reference systems:

• centroid_centered with time_reference=origin is the current recommended formal baseline. The grid remains centered on the effective Grond CMT source, but finite-fault W-phase Green functions and BP windows share PDE origin time.
• The finite-fault program no longer exposes time_reference=centroid. The previous centroid-time result is obsolete for the Peng finite-fault workflow and is not retained as a production option.
• hypocenter_fault_start with time_reference=origin uses the catalog event latitude, longitude, and depth as the first fault edge and uses PDE origin time as the finite-fault time zero. The grid coordinates then run from 0..fault_length_km and 0..fault_width_km.

The second mode is a hypothesis to be tested, not a silent replacement for the CMT-centered model. Use the matching BP anchor when generating its image:

/opt/anaconda3/envs/grond/bin/python \
  ../bodywave_backprojection/run_us7000srb1_bodywave_backprojection.py \
  --config ../bodywave_backprojection/bodywave_backprojection_config_all.json \
  --anchor-mode hypocenter_fault_start \
  --output-dir ../bodywave_backprojection/results_np2_hypocenter_anchor

/opt/anaconda3/envs/grond/bin/python run_peng2025_joint_inversion.py \
  --config peng2025_joint_config.json \
  --bp-image ../bodywave_backprojection/results_np2_hypocenter_anchor/bodywave_bp_image.npz \
  --plane-name NP2 --source-depth-km 55.193 \
  --anchor-mode hypocenter_fault_start --time-reference origin \
  --gamma-bp 4.5 --lambda-smooth 1 --iterations 5 \
  --output-dir results_peng2025_np2_hypocenter_anchor

The BP image and joint inversion must use the same anchor mode, source depth, and time reference. Translating signed centered coordinates to 0..L/0..W for plotting is not sufficient to change the physical anchor.

Matched intermediate-anchor test

An intermediate latitude/longitude is useful as a diagnostic of spatial registration, but it must not be described in the paper as a newly measured hypocenter. The following driver regenerates both the BP image and the finite- fault W-phase Green matrix for every candidate:

/opt/anaconda3/envs/grond/bin/python search_peng2025_anchor.py \
  --joint-config peng2025_joint_config.json \
  --bp-config ../bodywave_backprojection/bodywave_backprojection_config_all.json \
  --plane-name NP2 --source-depth-km 21.56465 \
  --anchors 5.727169:125.105374,5.693327:125.090756,5.659484:125.076137,5.625642:125.061519,5.591800:125.046900 \
  --station-csv ../bodywave_backprojection/results/np2_bp_all_first_pass/prepared_bodywave_stations.csv \
  --gamma-bp 4.5 --lambda-smooth 1 --iterations 1 \
  --output-dir results_anchor_search_np2

The first pass is deliberately one iteration. Rerun the best two or three anchors with --iterations 5; compare W-phase Misfit and BP correlation separately. A candidate is informative only if the source coordinates, depth, BP image and W-phase Green functions are all regenerated together.

Hypocenter-anchor gamma-lambda search

For the current Hypocenter fault-start experiment, the parameter search must carry the anchor and origin-time flags explicitly. Otherwise a search started from the default JSON can silently evaluate the wrong spatial reference:

/opt/anaconda3/envs/grond/bin/python search_peng2025_parameters.py \
  --config peng2025_joint_config.json \
  --bp-image ../bodywave_backprojection/results/np2_hypocenter_anchor/bodywave_bp_image.npz \
  --bp-station-csv ../bodywave_backprojection/results/np2_hypocenter_anchor/prepared_bodywave_stations.csv \
  --plane-name NP2 --source-depth-km 55.193 \
  --anchor-mode hypocenter_fault_start --time-reference origin \
  --gamma-values 1,2,3,4.5,6,8 \
  --lambda-values 0.1,0.3,1,3,10 \
  --iterations 1 --output-dir results_hypocenter_gamma_lambda_coarse

Use the resulting W-phase Misfit, BP correlation, and frequency-dependence index to select a short list. Then rerun only those candidates with --iterations 5. Do not rank candidates by BP correlation alone: the current Hypocenter result already demonstrates that a high BP correlation can coexist with poor W-phase fitting.

Paper-style diagnostic figures

Before the final fine search, generate the paper-style diagnostic set from one complete FFWI-M candidate. This is a plotting step only; it does not modify the solution or select a final parameter set:

MPLCONFIGDIR=/tmp/codex-mplconfig \
/opt/anaconda3/envs/grond/bin/python plot_peng2025_paper_figures.py \
  --joint-dir results_np2_trace_novr_gamma_lambda_nt6_iter5/g4p5_l1_physical_taup_mainshock \
  --search-csv results_np2_trace_novr_gamma_lambda_nt6_iter1/parameter_search_results.csv \
  --output-dir results_np2_trace_novr_gamma_lambda_nt6_iter5/g4p5_l1_physical_taup_mainshock/paper_figures

The script writes PNG and PDF versions of four figures:

1. figure1_pareto_tradeoff: W-phase Misfit versus BP-image correlation for the coarse gamma-lambda search. The star marks the current diagnostic candidate, not a result that is valid independently of the normalization and search grid.
2. figure2_source_and_bp_overlay: a geographic projection of the discrete source cells and a fault-plane view of the time-integrated normalized A^J field with relative observed BP-image contours. The map uses a grayscale shaded-relief background so that the warm source color scale is visually dominant. The BP contours stay on the native 3 x 5 grid; no interpolation is used in the formal figure.
3. figure3_temporal_moment_and_mechanism: scalar-moment fraction and the FFWI-M mechanism for each MTW window. The vertical placement of the beachballs is a layout device, not a second physical coordinate.
4. figure4_wphase_and_bp_fit: selected W-phase waveform fits plus observed, predicted, and residual integrated BP images. Non-negative BP images use a sequential scale; only the signed residual uses a diverging scale centered at zero. Waveform traces are station-normalized because the saved inversion data have already been weighted by observed RMS.

For the current candidate, the post-FFWI-M summary is approximately W-phase Misfit 0.13 and BP correlation 0.97. These values should be copied from the candidate's joint_summary.json and reported with the exact definition used; the pre-FFWI joint Misfit and the fixed-A^J FFWI-M Misfit are different diagnostics. In particular, a visually smooth BP panel must not be presented as evidence of fine spatial resolution when the integrated BP image itself has weak dynamic range.

IDBP kernel validation figure

Generate the kernel-validation figure separately from the final joint-model figures:

MPLCONFIGDIR=/tmp/codex-mplconfig \
/opt/anaconda3/envs/grond/bin/python plot_peng2025_kernel_validation.py \
  --taup-dir results_physical_gi_tauP_np2_nt6 \
  --mccc-dir results_physical_gi_mccc_np2_nt6 \
  --output-dir results_physical_gi_kernel_validation \
  --image-case taup

The six panels contain native-grid observed d^I, predicted G^I A^I, the image residual, an observed-predicted scatter plot, TauP/MCCC singular-value spectra, and leave-one-azimuth-group-out operator correlations. The current TauP baseline gives image correlation about 0.973 and Misfit about 0.009; MCCC gives a similar image fit while reducing the full-operator condition number. These are operator/data-mapping diagnostics, not evidence that the earthquake source is resolved at a scale finer than the 30 x 20 km cells.

The existing publication-ready tectonic-context and point-source W-phase figures remain available as independent inputs for the first main figure:

• results_wphase_ff_v18_aftershock_inset/paper_figure_18_aftershock_inset_marginal_map.png
• results_wphase_ff_v20_publication_figures/paper_figure_1_babo_cmt_posterior.png

Matched time-window search completed 2026-07-20

Before revising the paper figures, matched trace-normalized searches were run for NT=4, 5, and 6 with 9 s/9 s windows. Fixed-window-width sensitivity tests were also run for NT=6 with 6 s/6 s and 12 s/12 s windows. All 280 candidates completed successfully with the same NP2 hypocenter/fault-start anchor, origin-time reference, TauP physical G^I, and rupture_time_mode=none. The comparison files are:

• results_timewindow_search_summary/timewindow_search_comparison.csv
• results_timewindow_search_summary/timewindow_search_summary.json

The single-joint-solve results are screening values. Representative 10-iteration FFWI-M results are archived under results_timewindow_ffwi_finalists/. The current comparison is:

parameterization	final W-phase Misfit	BP correlation	FFWI-M converged
NT=4, 9/9 s, gamma=0.5, lambda=1	0.174	0.968	yes
NT=5, 9/9 s, gamma=6, lambda=0.1	0.172	0.976	yes
NT=6, 6/6 s, gamma=4.5, lambda=0.1	0.171	0.984	yes
NT=6, 9/9 s, gamma=10, lambda=0.01	0.168	0.973	yes
NT=6, 12/12 s, gamma=0.5, lambda=10	0.177	0.916	yes

These numbers do not yet define a unique final model. The 9 s candidate has the lowest W-phase Misfit, whereas the 6 s candidate has the highest BP correlation and uses a stronger, less edge-sensitive smoothing value. The 6 s/6 s TauP and MCCC sensitivity runs give final W-phase Misfits of about 0.171 and 0.173 and BP correlations of about 0.984 and 0.982, respectively. This supports kernel robustness for this candidate. Under the same 6 s/6 s parameters, NP1 gives final W-phase Misfit about 0.175 and BP correlation about 0.983; NP2 is only a marginal numerical preference, not a unique fault-plane determination.

The provisional leading candidates must still be checked against reference depth and broader geometry sensitivity before the final parameter is fixed. No publication figure should use the provisional candidate label as a final scientific conclusion.

Matched physical depth search completed 2026-07-21

The depth diagnostic was repeated with a physically matched workflow. For each NP2 reference depth, the BP image, formal subfault depths, station-level TauP G^I, and joint inversion were regenerated together. The BP estimator was fixed to the formal waveform-linear baseline: eight azimuthal groups, SNR >= 2, a 20 s correlation window, maximum lag 5 s, minimum correlation 0.2, and waveform integration. This is distinct from the older search_peng2025_depth.py, which used the empirical first-six kernel, and from the preliminary envelope-semblance diagnostic.

The matched coarse scan used NT=6, 6 s/6 s windows, gamma=4.5, lambda=0.1, trace normalization, origin time, hypocenter_fault_start, and rupture_time_mode=none:

• results_physical_depth_search_np2_nt6_stf6_waveform_linear/depth_search_results.csv
• results_physical_depth_search_np2_nt6_stf6_waveform_linear/DEPTH_SEARCH_REPORT.md
• search_peng2025_depth_physical.py

reference depth	W Misfit	BP Misfit	BP correlation	GI condition number
25 km	0.1963	0.04357	0.944	1.18e4
30 km	0.2102	0.04008	0.951	7.47e3
35 km	0.2261	0.03602	0.957	1.04e4
40 km	0.2385	0.03479	0.960	1.04e4
45 km	0.2484	0.03685	0.958	1.03e4
50 km	0.2581	0.04047	0.954	8.58e3
55 km	0.2766	0.04376	0.950	8.58e3

These are one-iteration screening values. The BP image is regenerated for each depth, so BP correlation is a matched-kernel diagnostic rather than an independent depth measurement. The coarse result favors shallower reference depths in W-phase Misfit, while the BP metric is best near 35--40 km; the two criteria do not select a unique depth.

The 25, 40, and 55 km representatives were then run through ten FFWI-M iterations under the same waveform-linear baseline:

reference depth	final W Misfit	BP correlation	mean Kagan angle	total moment change	converged
25 km	0.1328	0.9440	0.037 deg	0.015%	yes
40 km	0.1346	0.9596	0.005 deg	0.007%	yes
55 km	0.1347	0.9497	0.012 deg	0.015%	yes

The full-iteration products are under results_physical_depth_search_np2_nt6_stf6_waveform_linear_ffwi10/. The differences in final W Misfit are small, so the current data do not uniquely resolve the absolute depth. The 40 km case is the best compromise among the tested final W/BP metrics and has the smallest mechanism-update change, but it remains a provisional reference depth. Independent depth-phase or tectonic constraints remain necessary before selecting one depth for the paper.

Matched gamma-lambda search at the provisional 40 km depth

After the matched depth test, a complete 8 x 7 grid of 56 gamma-lambda combinations was run with the same physical waveform-linear BP image, TauP G^I, NT=6, 6 s/6 s windows, origin-time reference, fault-start anchor, and trace normalization. The screening table is:

• results_trace_timewindow_search_np2_depth40_nt6_stf6/parameter_search_results.csv
• results_trace_timewindow_search_np2_depth40_nt6_stf6/parameter_search_summary.json

The screening run is followed by full ten-iteration FFWI-M tests for the representative candidates in results_depth40_gamma_lambda_ffwi10/:

gamma	lambda	final W Misfit	BP correlation	Eq.4 joint Misfit	mean Kagan (deg)	moment change (%)	converged
0.5	3.0	0.1324	0.840	0.801	0.060	0.038	yes
2.0	10.0	0.1324	0.859	0.612	0.004	0.002	yes
4.5	0.3	0.1346	0.960	0.394	0.005	0.007	yes
6.0	0.3	0.1346	0.960	0.289	0.003	0.004	yes
10.0	0.3	0.1346	0.960	0.151	0.001	0.001	yes

The current provisional choice is gamma=10.0, lambda=0.3: it lies on the high-BP-correlation plateau and has the smallest Eq.4 joint Misfit among the tested candidates, while its final FFWI-M W-phase Misfit differs only slightly from the lower-gamma candidates. The joint Misfit in this table is the pre-FFWI-M Eq.4 solve metric; the post-update W-phase metric is reported separately as final W Misfit. gamma=4.5, lambda=0.3 and gamma=6.0, lambda=0.3 remain required nearby sensitivity models. This is a data-driven choice for this event, not a transfer of the Nepal values gamma=4.5, lambda=1.0.

The current model is conditional on the provisional 40 km reference depth: the matched depth tests do not uniquely resolve absolute depth. It also uses rupture_time_mode=none; no additional rupture-speed constraint has been imposed.

Canonical organized workflow

For routine reruns and paper production, use final_workflow_20260721/README.md. It contains the numbered execution scripts, clearly named result directories, copied final products, and the legacy-result index. The exploratory directories above are retained for provenance and should not be mixed into the formal result without an explicit sensitivity label.

JFWIBP