{"id":"7f42c508-48a6-4389-b849-b13eac5282ae","arxiv_id":"2501.13532","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"A discrete adjoint state method, built on a second-order fast marching solver, provides self-consistent gradients for velocity and source location in eikonal traveltime tomography and enables deterministic and probabilistic inversion.","lead":"This paper derives a discrete adjoint method for eikonal-equation traveltime tomography that computes gradients for both subsurface velocity and earthquake source locations in one coherent framework. It demonstrates the method on synthetic 2D and 3D tests with L-BFGS optimization and Hamiltonian Monte Carlo, showing that probabilistic uncertainty quantification is feasible for this class of inverse problems.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The discrete adjoint is derived for a fixed stencil configuration, but the actual FMM forward map switches stencils and is nondifferentiable at ties, so the claim of exactly consistent gradients is unverified; without a finite-difference check the central claim lacks numerical support.","rationale":"The reader's weakest assumption identifies exactly the concern I consider load-bearing: the discrete adjoint is derived for a fixed stencil configuration, while the FMM forward solver is nondifferentiable at stencil switches and ties, and no numerical verification compares the adjoint gradient to the derivative of the actual forward map. This is not an attack on the authors; it is a precise gap between the theory as written and the implemented forward model. The paper is otherwise coherent: the derivation is algebraically reasonable for a fixed stencil, the triangular-structure argument exploiting FMM ordering is plausible, and the code is publicly available, which is real evidence. The empirical inversions in Section 4 are suggestive, but they rely on the very gradients that are unvalidated. I considered whether the more serious concern might be the inverse-crime nature of the synthetic tests, where the same forward solver generates data and is used for inversion, so forward-modeling error is not tested. That is real and worth noting, but it is secondary: even under inverse-crime conditions, incorrect gradients would degrade the inversions, so the absence of gradient validation is more fundamental. I also noted the text/figure iteration-count inconsistencies (Fig. 4 caption says 150 iterations, text says 80; text for the NUTS 2D example says 50000 iterations and 150000 saved models, while the caption says 46000 iterations; the 3D NUTS text says 10000 iterations, caption says 1454). These inconsistencies undermine reproducibility but do not affect the correctness of the central mathematical claim, so I treat them as minor. My recommended verdict is therefore unchanged from the reader's CONDITIONAL: the paper is likely sound but needs a gradient-validation experiment to support the exact-consistency claim. The concrete test above would settle the concern directly.","tokens_in":21842,"tokens_out":4238,"duration_ms":43392,"concrete_test":"Run a finite-difference verification in the released EikonalSolvers.jl code: on a small 2D grid (e.g., 20×20) with a random smooth velocity model and an off-grid source, compute the adjoint gradient from Section 3 (both dψ/dv and dψ/ds_r) and compare to central finite differences of the actual forward FMM traveltime code, perturbing each velocity parameter and source coordinate by ±ε with the same code path, including receiver interpolation. Report relative L2 error and maximum pointwise error as ε→0. Then repeat on models deliberately chosen to trigger stencil switches or ties (e.g., a constant-velocity starting model or geometry with equal traveltime neighbors) to check whether errors spike at those points.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim of exact discrete consistency rests on the adjoint derivation in Section 3.1.1, which differentiates the implicit equation (9) and, via eq. (20), treats the stencil matrices D^x, D^y, D^z as fixed. In the actual forward solver, those stencils are selected at runtime by the FMM upwind rule (eq. 5), including the max(·,0) terms and the fallback between second-order (eq. 3) and first-order (eq. 4) stencils when upwind values are unavailable. The stencil selection is a piecewise-constant function of the traveltime field, so the implemented forward map is only piecewise differentiable; eq. (20) gives a one-sided, stencil-conditional derivative, not the derivative of the full algorithm. At ties or stencil switches, one-sided derivatives disagree, and the derivative in eq. (22) does not correspond to any unique linearization of the forward code. The paper never validates the adjoint gradient against finite differences of the actual FMM implementation, so the size and frequency of any mismatch are completely unquantified. This is load-bearing because the paper's stated advantage over continuous adjoint methods is exactly that the discrete adjoint is consistent with the discretized forward model, and that consistency fails precisely at nondifferentiable points. For L-BFGS, an inaccurate gradient may still behave like a generalized subgradient, but convergence guarantees are lost. For HMC/NUTS, a biased gradient biases the stationary distribution, undermining the claimed uncertainty quantification. The empirical inversions in Section 4 could be robust in spite of this, but the paper provides no evidence to establish that.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops a discrete adjoint method for traveltime tomography based on the eikonal equation, with second-order fast marching used as the forward solver. Gradients are derived with respect to both velocity structure and source location, including a treatment of grid refinement around the source, and are then embedded in deterministic (L-BFGS) and probabilistic (HMC/NUTS) inversion frameworks. Synthetic 2D and 3D examples are presented for both approaches. The central claim is that the discrete adjoint derivation yields gradients that are exactly consistent with the discretized FMM forward model, enabling efficient joint inversion for velocity and source position.","tokens_in":22058,"tokens_out":2592,"duration_ms":24123,"significance":"If the consistency claim holds, the paper would be a useful contribution: it unifies velocity and source-location gradients in one discrete adjoint framework, allows arbitrary source/receiver positions, handles local grid refinement, and demonstrates both deterministic and probabilistic inversion, with open-source code provided. The derivations are self-contained and the numerical experiments illustrate the intended use. However, the exact-consistency claim is not numerically verified against the actual piecewise-differentiable FMM forward map, and this gap is load-bearing for both the L-BFGS and HMC applications. The paper is therefore promising but requires additional validation and careful qualification before the central claim can be accepted.","major_comments":[{"comment":"The adjoint derivation differentiates the implicit equation (9) while treating the stencil matrices D^x, D^y, D^z as fixed. In the implemented FMM (Section 2.1, Eq. (5)), the stencils are selected at runtime by the max(·,0) upwind rule, with fallback between second-order (Eq. (3)) and first-order (Eq. (4)) stencils depending on the availability of upwind values. The stencil configuration is a piecewise-constant function of the traveltime field, so Eq. (20) gives a branch-local derivative, not the derivative of the full implemented forward map. At stencil switches or ties, the gradient in Eq. (22) does not correspond to a unique linearization of the forward code. This is load-bearing because the paper's stated advantage over continuous adjoint methods is exactly discrete consistency, and the mismatch is never quantified.","section":"Section 3.1.1, Eqs. (9)-(22)"},{"comment":"No finite-difference validation of the adjoint gradients is presented in any of the deterministic or probabilistic experiments. Given the nondifferentiable stencil selection in the forward solver, a direct comparison of the adjoint gradient against finite differences of the actual FMM implementation is essential to support the claim of exact consistency stated in the Introduction, Section 5, and the abstract. Without such a test, the magnitude and frequency of any gradient error remain unknown, which also affects the HMC/NUTS results because a biased gradient biases the stationary distribution.","section":"Section 4 (numerical experiments)"},{"comment":"The gradient with respect to source location differentiates Eq. (2) while assuming a fixed set S of nodes surrounding the source. As the source moves across a grid-cell boundary, the set of enclosing nodes changes discontinuously, making the source-location forward map piecewise differentiable. The derivation does not address this, and no numerical test checks the behavior of dψ/ds_r at such boundaries. Since joint inversion for source location is a central contribution, this branch-dependence should be discussed and ideally validated numerically.","section":"Section 3.2, Eqs. (59)-(60)"}],"minor_comments":[{"comment":"The phrases “exactly the counterparts of the forward model” and “machine precision” overstate the result given the piecewise-differentiable stencil selection; the text should be qualified to refer to differentiability within fixed stencil configurations.","section":"Introduction, Section 5"},{"comment":"The velocity v_s in Eq. (2) is not uniquely defined when the four (or eight) surrounding grid points have different velocities; the paper should specify how v_s is chosen.","section":"Section 2.1, Eq. (2)"},{"comment":"The caption states that the inversion ran for 150 iterations, while the text reports 80 iterations; this discrepancy should be corrected.","section":"Section 4.1, Figure 4 caption"},{"comment":"The text reports 50000 NUTS iterations but then states that 150000 models were saved after burn-in; the relationship between these numbers should be clarified.","section":"Section 4.2, first example"},{"comment":"The caption refers to a random model “after 1454 iterations,” while the text says 10000 iterations were run; the reported iteration counts should be reconciled.","section":"Section 4.2, Figure 7 caption"},{"comment":"The right-hand side written as ∂u/∂τ is ambiguous; it should be written as ∂û_h/∂τ_p to match the notation introduced in Eq. (42).","section":"Section 3.1.4, Eq. (44)"}],"recommendation":"major_revision","confidential_remarks":"The manuscript is technically interesting and the code availability is a strong point. The main risk is the unquantified gap between the discrete adjoint gradient and the actual piecewise-differentiable FMM forward map; a finite-difference check is essential and is likely to be expected by readers. The iteration-count inconsistencies in Section 4 should also be fixed. The paper fits the scope of physics.geo-ph and could become acceptable after these revisions."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The paper gives a discrete adjoint method for eikonal traveltime tomography that covers velocity and source location in one framework, with off-grid sources/receivers and local grid refinement. It also shows HMC-based uncertainty quantification. That last part is the genuine addition over Li et al. (2013) and Treister & Haber (2016): the source-location gradient is derived through the same discrete adjoint machinery, and the probabilistic examples are more than a toy. The derivation is clear and self-contained, the code is public, and the synthetic examples are easy to follow.\n\nThe soft spot is the one the stress test flags: the adjoint is derived by differentiating the implicit equation with the stencil matrices treated as fixed, but the actual FMM selects stencils through max() and fallback logic, so the forward map is only piecewise differentiable. At ties or stencil switches, the derivative in eq. (20) is a one-sided conditional derivative, not the derivative of the code. The paper never checks the gradient against finite differences, so the size of any mismatch is unknown. This matters because the stated advantage over continuous adjoints is exactly consistency with the discrete forward model. For L-BFGS it might still work as a subgradient, but for HMC a biased gradient biases the target distribution. So the central claim is plausible but unverified.\n\nA second, smaller complaint: the tests are inverse-crime synthetic experiments—data generated with the same solver on the same grid—so they demonstrate self-consistency, not real-world accuracy. That is acceptable for a methods paper, but it should be stated. There are also small inconsistencies in iteration counts (text says 80, caption says 150; text says 50, caption says 40; the HMC text mentions 50000 iterations and then \"150000 were saved\"), which should be cleaned up.\n\nOverall, I think the method is worth engaging with. The missing gradient check is fixable and should be required. I would send this to a serious reviewer, mainly for the source-location gradient and the HMC application, and ask for a finite-difference validation plus discussion of the nondifferentiability.","headline":"Clean derivation, public code, and a real extension to source-location gradients and HMC, but the discrete adjoint is only proved for fixed FMM stencils and the missing finite-difference check leaves the central consistency claim unverified.","tokens_in":22669,"tokens_out":2780,"would_cite":true,"duration_ms":24383,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper derives gradient formulas for velocity and source location directly from the discretized fast-marching eikonal solver, giving a single adjoint framework that works for both L-BFGS inversion and Hamiltonian Monte Carlo sampling.","keywords":["traveltime tomography","eikonal equation","discrete adjoint method","fast marching method","velocity inversion","source location inversion","Hamiltonian Monte Carlo","uncertainty quantification"],"falsifier":"Run central-difference checks of the discrete-adjoint gradients against the same FMM code: perturb each velocity voxel and each source coordinate by a small epsilon, recompute traveltimes, and compare the ratio (psi(theta+epsilon)-psi(theta-epsilon))/(2 epsilon) with the adjoint gradient at models near stencil switches and ties in the max operation. If pointwise relative errors are far above round-off where the solver's stencil selection changes, the claim that the gradient is exactly consistent with the implemented forward map is refuted; if the errors stay near machine precision everywhere, the claim is supported.","tokens_in":21570,"feed_emoji":"🌊","tokens_out":11507,"duration_ms":93901,"temperature":0.7,"pith_summary":"The paper's central proposition is that traveltime tomography based on the eikonal equation can be turned into a coherent inverse problem by deriving gradients directly from the discretized forward solver, rather than from a separately discretized continuous adjoint equation. Working from a second-order fast-marching method, the authors write the discrete eikonal equation as an implicit system, differentiate it to obtain an adjoint equation, and use the same adjoint variable for both the velocity gradient and the source-location gradient. The resulting gradients are, by construction, consistent with the discrete forward model, and the fast-marching update order turns the adjoint system into a triangular one that is cheap to solve. This makes it practical to run joint inversions for velocity and source position, and because the gradients are cheap and accurate, they can feed both deterministic optimization (L-BFGS) and gradient-based Bayesian sampling (Hamiltonian Monte Carlo). The paper demonstrates the approach on synthetic 2D and 3D experiments, including one where Bayesian sampling reveals a bimodal posterior that a deterministic point estimate would miss.","feed_headline":"One adjoint equation drives both velocity and source tomography","feed_subtitle":"Fast-marching traveltime tomography now has a single gradient for both deterministic and Bayesian inversion.","key_machinery":"The load-bearing object is the discrete implicit form of the eikonal equation, $f_i = (\\sum_j D^x_{ij} u_j)^2 + (\\sum_j D^y_{ij} u_j)^2 + (\\sum_j D^z_{ij} u_j)^2 - 1/v_i^2 = 0$, where the $D$ matrices carry the second-order upwind finite-difference stencils selected by the fast-marching solver. The adjoint state $\\lambda$ is defined by $(\\partial f / \\partial u)^T \\lambda = -(\\partial \\psi / \\partial u)^T$, and the same $\\lambda$ is reused in the velocity-gradient and source-location-gradient chain rules. The structural trick is that ordering grid points by the fast-marching wavefront arrival makes $\\partial f / \\partial u$ lower triangular, so the transposed adjoint system is upper triangular and solvable in a single pass; the stencils are sparse, with only two or three nonzeros per row. For the refined source-region grid, a second adjoint system is solved on the fine grid with the coarse-grid traveltime acting as the adjoint source, and the chain-rule factor $N_{qa}$ maps the velocity interpolation back to the coarse grid.","core_discovery":"The contribution the paper is trying to establish is that all the pieces of eikonal traveltime tomography—velocity update, source-location update, local grid refinement around the source, and arbitrary source and receiver placement—can be handled by one discrete adjoint formalism built on the same second-order fast-marching forward model. The forward solver is written as an implicit discrete equation; differentiating it gives the adjoint equation whose solution is a single field lambda. The same lambda is then folded into the chain rules for both velocity parameters and source coordinates, with the case of a refined source-region grid handled by solving a second, analogous adjoint problem on the fine grid. Because the fast-marching method solves grid points in increasing traveltime order, the matrix in the adjoint system is triangular, so the gradient computation costs roughly two forward solves and is independent of the number of receivers. The paper claims this yields gradients that correspond exactly to the discretized forward model and demonstrates that these gradients drive successful deterministic L-BFGS reconstructions and Hamiltonian Monte Carlo posteriors in 2D and 3D.","pith_inferences":["The same discrete-adjoint construction should extend to other grid-based eikonal solvers, such as factored eikonal forms or higher-order stencils, by replacing the derivative matrices in the forward relation; the chain-rule structure of the derivation is solver-agnostic.","A practical validation step would be a finite-difference check of the gradients near points where the forward solver switches stencils or ties in the max selection; that would quantify how often the fixed-stencil differentiability assumption matters in practice.","For field applications, the probabilistic branch could be used to produce posterior covariance maps for both velocity and hypocenter locations, which would feed directly into seismic hazard or event-relocation studies.","The triangular structure of the adjoint system suggests the gradient cost scales like the forward solve, so large 3D regional models with many sources and receivers are the natural next target; the paper's 3D examples use relatively small grids."],"forward_implications":["Gradients produced by the discrete adjoint are consistent with the FMM forward model by construction, removing the forward/adjoint discretization mismatch that can arise when the continuous adjoint equation is discretized separately.","The same single adjoint field supplies both velocity and source-location gradients, so joint inversion becomes a natural extension of the single-parameter case rather than a separate two-step procedure.","Ordering grid points along the fast-marching wavefront makes the adjoint matrix triangular, so gradient computation is efficient and the cost is nearly independent of the number of receivers.","The formalism accommodates arbitrary source and receiver positions plus a refined grid around the source, which reduces the dominant traveltime errors near the source without dropping those terms from the gradient.","Because gradients are cheap and accurate, the same framework feeds both L-BFGS deterministic inversion and HMC probabilistic sampling, yielding uncertainty quantification that can expose multimodality in the posterior."],"supporting_citations":[{"why":"Introduces the fast marching method; the forward solver the adjoint is built on.","marker":"Sethian, 1996"},{"why":"Applies FMM to seismic traveltime computation; source of second-order upwind finite differences and the entropy-satisfying update.","marker":"Sethian and Popovici, 1999"},{"why":"Proposes a refined grid around the source in FMM; basis for the refined-grid treatment in the gradient derivation.","marker":"Rawlinson and Sambridge, 2004b"},{"why":"Provides the discrete upwind eikonal formulation based on max operations that eq. (5) and the implicit forward model follow.","marker":"Treister and Haber, 2016"},{"why":"Review of the adjoint-state method; supplies the general derivation used to write the discrete adjoint equation.","marker":"Plessix, 2006"},{"why":"Continuous adjoint eikonal tomography; representative earlier approach that the discrete adjoint method is contrasted with.","marker":"Leung and Qian, 2006"},{"why":"Recent continuous-adjoint traveltime tomography workflow for velocity and source location that the paper extends with a unified discrete formalism.","marker":"Tong, 2021a"},{"why":"Provides the Hamiltonian Monte Carlo software used to run the probabilistic experiments.","marker":"Zunino et al., 2023"},{"why":"Supplies the No-U-turn sampler used in the HMC sampling examples.","marker":"Hoffman and Gelman, 2014"}],"fun_headline_variants":["Discrete adjoint unifies velocity and source tomography","One adjoint field drives eikonal traveltime inversion","Adjoint gradients for joint velocity-source inversion","Eikonal inversion with a single discrete adjoint","Triangular adjoint accelerates traveltime tomography"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The whole gradient-consistency claim rests on treating the fast-marching solver's selection of neighboring grid points as a fixed smooth configuration when differentiating, although the solver actually switches stencils through max operations and order fallbacks.","fun_headline_variants_meta":{"raw":{"variants":["Discrete adjoint unifies velocity and source tomography","One adjoint field drives eikonal traveltime inversion","Adjoint gradients for joint velocity-source inversion","Eikonal inversion with a single discrete adjoint","Triangular adjoint accelerates traveltime tomography"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000257,"raw_usage":{"total_tokens":1589,"prompt_tokens":969,"completion_tokens":620,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":585,"completion_tokens_details":{"reasoning_tokens":547}},"tokens_in":585,"tokens_out":620,"duration_ms":6046,"temperature":1.0,"reasoning_tokens":547,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-10T15:51:11.350855+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run central-difference checks of the discrete-adjoint gradients against the same FMM code: perturb each velocity voxel and each source coordinate by a small epsilon, recompute traveltimes, and compare the ratio (psi(theta+epsilon)-psi(theta-epsilon))/(2 epsilon) with the adjoint gradient at models near stencil switches and ties in the max operation. If pointwise relative errors are far above round-off where the solver's stencil selection changes, the claim that the gradient is exactly consistent with the implemented forward map is refuted; if the errors stay near machine precision everywhere, the claim is supported.","supporting_citations":[{"cited_title":"A fast marching algorithm for the factored eikonal equation","cited_arxiv_id":null,"evidence_quote":"Provides the discrete upwind eikonal formulation based on max operations that eq. (5) and the implicit forward model follow."},{"cited_title":"An adjoint state method for three-dimensional transmission traveltime tomography using first-arrivals","cited_arxiv_id":null,"evidence_quote":"Continuous adjoint eikonal tomography; representative earlier approach that the discrete adjoint method is contrasted with."},{"cited_title":"HMCLab : A framework for solving diverse geophysical inverse problems using the Hamiltonian Monte Carlo method","cited_arxiv_id":null,"evidence_quote":"Provides the Hamiltonian Monte Carlo software used to run the probabilistic experiments."}],"review_version":1}