Evaluating imaging uncertainty associated with the near surface and added value of vertical arrays using Bayesian seismic refraction tomography

Traditional tomographic methods do not consider the uncertainties associated with near-surface velocities and static corrections and provide a deterministic solution to the estimation problem. However, these uncertainties significantly affect structural mapping and interpretation of seismic imaging results. On the other hand, Bayesian first-arrival tomography provides multiple near-surface models that fit observed traveltimes equally well and enable the study of potential solution distributions. We demonstrate this approach on a complex synthetic near-surface model, representative of arid environments, to quantify associated velocity and statics uncertainties. We evaluate two different parameterizations for subsurface velocities in the context of near-surface Bayesian tomography: Voronoi tessellation with natural neighbor interpolation and the more conventional Delaunay triangulation with linear interpolation. Our analysis shows that the Voronoi cell parameterization with natural neighbor interpolation is more appropriate for this problem. Finally, the new approach is applied to compare two alternative acquisition geometries comprising conventional surface receivers and surface receivers augmented with vertical receiver arrays. The results demonstrate that adding vertical receiver arrays to conventional surface receivers can significantly reduce the near-surface velocity uncertainty and thus increases the accuracy of the seismic imaging results. Furthermore, the study shows that Bayesian tomography can be used as a tool for evaluating different source and receiver geometries during the acquisition design stage


Introduction
Complex near-surface geology is one of the biggest challenges for land seismic data processing and interpretation (Robinson and Al-Husseini, 1982 ;Keho & Kelamis 2012 ).The accuracy of seismic imaging results strongly depends on the near-surface velocity models (Nosjean et al. 2017 ), especially when mapping low-relief structures with typical vertical closures of < 60 m (or ∼30 ms).In such cases, reliable structural imaging of target horizons requires very accurate near-surface velocities and corresponding long-wavelength static corrections.One of the main instruments for nearsurface velocity reconstruction is a first-arrival traveltime tomography that relies on arrival times of diving or refracted waves (Zhang & Toksöz 1998 ).The obtained tomography model can be used to compute static corrections or directly embedded into a depth velocity model for seismic migration.
In industry practice, seismic tomography is typically conducted deterministically, resulting in a single velocity model highly dependent on the initial model and various userdefined parameters such as regularization or smoothing.As a result, the uncertainty associated with the estimated solution is not evaluated.At the same time, seismic uncertainty is a critical factor in quantifying different exploration risks (Osypov et al. 2013 ;Bakulin et al. 2017 ;Nosjean et al. 2017 ).
Estimating uncertainty in near-surface tomography is essential in accurately mapping subsurface structures and reducing exploration risks.One approach is to conduct multiple tomographic runs with different parameters, providing preliminary insights into solution variability (Wang 2003 ;Ley et al. 2013 ;Silvestrov et al. 2023 ).However, a more rigorous and comprehensive approach is through a Bayesian framework, resulting in a vast collection of models that fit observed traveltimes equally wel l and al low for calculating various solution statistics (Rawlinson et al. 2014 ).Such a method enables a more accurate estimation of uncertainties and helps decision-making in exploration activities.A variety of methods have been used previously to obtain probabilistic solutions for the tomographic problem.These methods include Hamiltonian Monte Carlo (Fichtner et al. 2019 ), variational inference techniques (Zhang & Curtis 2020 ) and reversiblejump Markov Chain Monte Carlo (rj-MCMC) (Bodin et al. 2012 ;Egorov et al. 2020 ).Among these, the rj-MCMC method has been extensively employed to analyze different field datasets in the context of global seismology (Bodin & Sambridge 2009 ;Galetti et al. 2015 ;Ryberg & Haberland 2018 ).This method involves defining velocities on a set of nodes where their positions are also treated as unknowns.As a result, the tomographic engine conducts inversion simultaneously for the values of velocities and node locations.To make the algorithm computationally feasible, the number of nodes in this method is typically much smaller than in conventional ray-based tomography with regular grids.As a result, velocity values must be estimated using interpolation and extrapolation methods based on a chosen parameterization away from these nodes.
The original formulation of the Bayesian tomography with rj-MCMC uses Voronoi tessellation as a parameterization method and implies a nearest-neighbor interpolation.This approach results in a velocity field represented as a piecewise-constant function, often failing to provide geologically feasible models.Several other parameterizations were suggested for Bayesian tomography to more accurately describe subsurface velocities.For example, Belhadj et al. ( 2018 ) invoke Johnson-Mehl tessellation, which also results in piecewise-constant velocities but with more complex shapes of cell boundaries.The authors also introduced a Gaussian kernel-based parameterization, which yields smoother models.However, both parameterizations require additional per-node parameters, such as cell implan-tation times for the former and Gaussian kernels widths for the latter, leading to a more complex inversion process.Hawkins et al. ( 2019 ) examined two parameterizations that use Delaunay triangulation and either linear or Clough-Tocher interpolation.Both methods yielded smooth velocity models and reduced parameterization-induced errors when the true velocity was smooth.Egorov et al. ( 2021 ) suggested an alternative parameterization based on Voronoi cells but utilizing the natural neighbor interpolation (Sibson 1981 ), and demonstrated its benefits through a simplified nearsurface model example.
One of the possible applications of the Bayesian tomography framework is to investigate how different source/receiver configurations affect the accuracy and uncertainty of tomography results (Egorov et al. 2022 ;Watts et al. 2023 ).To reduce the near-surface uncertainty in complex areas, a viable option is to employ a combined acquisition geometry that consists of surface sensors and vertical receiver arrays.Distributed acoustic sensing (DAS) can be used to easily implement such a survey with a single continuous fiber-optic cable connecting multiple vertical arrays or upholes.This survey method is called smart DAS uphole acquisition and has been studied by Bakulin et al. ( 2017Bakulin et al. ( , 2018 ) ).The method allows for the simultaneous recording of surface data and vertical data from the same shots.As a result, such a dataset contains additional transmitted/refracted raypaths and delivers traveltimes at multiple buried receivers inside the near surface at many different angles.This results in rich i l lumination coverage for depth imaging and accurate near-surface characterization.The smart DAS uphole acquisition technique also improves velocity estimation accuracy through deterministic first-arrival tomography (Alshuhail et al. 2019 ).Egorov et al. ( 2022 ) utilized Bayesian tomography to conduct an initial assessment of the precision and certainty enhancements made possible by DAS vertical receiver arrays.
Building on the prior research by Egorov et al. ( 2020Egorov et al. ( , 2021Egorov et al. ( , 2022 ) ), we perform a more rigorous investigation of the application of Bayesian first-arrival tomography for nearsurface characterization in the exploration context.First, we apply the method to a synthetic dataset generated from a realistic SEAM Ar id model (Or istaglio 2012 ; Bakulin & Silvestrov 2021 ), representing near-surface challenges in a desert environment.We evaluate the performance of a conventional linear interpolation method using Delaunay triangles for the velocity model and compare it with an alternative natural neighbor interpolation technique based on Voronoi cells.The latter approach yields better results than the former.Subsequently, we apply the Bayesian approach to investigate the effect of vertical arrays on reducing velocity estimation uncertainty.Finally, we incorporate this uncertainty into the seismic processing workflow and evaluate the accuracy of longwavelength static corrections, thereby analyzing their impact on imaging and mapping deep low-relief structures.

Method
We utilize a 2D hierarchical reversible-jump Markov Chain Monte Carlo (rj-MCMC) technique proposed by Bodin et al. ( 2012 ).The method is based on Bayes' theorem, which links p ( m | d ) , the posterior probability of the model m given the data d , with prior knowledge p ( m ) and the likelihood function p ( d | m ) (Wang 2016 ): (1) The likelihood defines how well the model fits the data.We use the Gaussian likelihood (Tarantola 2005 ;Bodin et al. 2012 ): ) } . (2) For the described problem, d is the vector of measured traveltimes of length M , g( m ) are the traveltimes computed for the model m and C e is the covariancematrix of the data noise.
If the noise is assumed to be Gaussian and uncorrelated, this matrix can be replaced by a diagonal matrix, C e =  2 n I , where  n is the standard deviation of arrival picking error and I is the identity matrix.
We parameterize the subsurface by N unstructured nodes.Each of the nodes is defined by the coordinates x i and z i and the velocity value v ( x i , z i ) , all of which are the parameters updated by the algorithm.The velocity model in the whole domain is computed from this set of nodes by a chosen parameterization algorithm.Traditionally, Voronoi parameterization has been used for this purpose, where the velocities between the nodes are interpolated by the nearest-neighbor algorithm.An alternative parameterization involves computing a Delaunay triangulation from the nodes' coordinates and interpolating the speed linearly within the triangles, which is effective for smooth velocity models and refraction geometry (Ryberg & Haberland 2018 ).As an alternative, we propose parameterizing the considered problem using a discrete version of the natural neighbor algorithm.The natural neighbor is an interpolation technique proposed by Sibson ( 1981 ), which produces smooth interpolant functions without requiring additional user-defined parameters.The natural neighbor method computes the interpolant velocity function v at an arbitrary location ( x, z ) as follows: Here, w i ( x, z ) are the weights computed from the areas u i .Such areas define the extent to which one Voronoi cell (corresponding to x in the tessellation { ( x i , z i ) , ( x, z ) } N i = 1 ) borr ows from the Voronoi cells of the original { ( x i , z i ) } N i = 1 tessellation.Note that we do not compute these weights geometrically but use a more effective scatter approach proposed by Park et al. ( 2006 ).The scatter approach is a numerically effective way of calculating the interpolant velocity in equation ( 3), which increases a computation speedup by looping over the tessellation cells { ( x i , z i ) } instead of output grid locations { ( x, z ) } .Models provided by natural neighbor interpolation are smooth and do not require additional input parameters.They contain no overshoots or undershoots, meaning that the velocity values in the model are always between the minimum and maximum velocities at the nodes, controlled by the priors.Finally, in contrast to linear interpolation based on the Delaunay parameterization (Ryberg & Haberland 2018 ), the natural neighbor parameterization seamlessly extrapolate velocity values away from the nodes and does not require extra conditioning or addition of the nodes on the model boundaries.
In addition to the velocities at the nodes and node coordinates, we also need to estimate the number of nodes in the model N (a small number compared to the number of regular grid cells) and the standard deviation of picking noise  n .To improve the convergence of the Markov chains, we use the coordinate scaling as suggested by Zhang et al. ( 2018 ) for all the compared parameterizations.We apply a scaling factor of four to the vertical coordinate when obtaining the velocity grid from the nodes, which makes the grid cells along the vertical axis artificially 'further away' from each other for the interpolation algorithm.This scaling promotes horizontal structures in the stochastic models.
We use multiple noninteracting Markov chains.Each chain gets an initial model from the priors.We use wide uniform prior distributions for velocity,  n and N .We allow the velocities to change from 700 to 8000 m s − 1 ,  n from 0.1 to 100 ms, and N from 11 to 901.At each step of the chain, either a new position or velocity is proposed for one of the nodes, or the birth/death of one of the nodes happens, or  n changes.The new model is then accepted or rejected by comparing its probability to the probability of the current model.A version of the Metropolis-Hastings acceptance cr iter ion modified for the reversible-jump method (Green 1995 ) is used for comparison.The cr iter ion sometimes may accept a model with the larger misfit, which results in a chain's random-wal k-li ke behavior and allows the chain to continue sampling from the posterior model distribution instead of converging to a single model.It is important to note that we use log-likelihoods and log-probabilities for practical calculations and comparisons to avoid the issues of numerical overflow.After collecting enough chains, we analyze the accepted models and estimate the velocity uncertainty by computing the mean, various moments (e.g. standard deviation), or by directly looking at velocity distributions for each of the locations in the subsurface.Finally, we compute the traveltimes with an eikonal solver.We do a complete eikonal computation at every step of the algorithm, in contrast to the approximation by Bodin et al. ( 2012 ), who used inner and outer loops, computed rays in the outer loop and updated the model with fixed rays in the inner loop.Accurate recomputation takes extra time but allows us to avoid tuning the inner loop length.Galetti et al. ( 2015 ) applied a similar algorithm modification.They demonstrated that recomputing the rays for each model produces different estimated uncertainties compared to the approximate method.
Fundamental restrictions of traveltime tomography also need to be considered when obtaining velocity estimates.For example, it is well-known that refraction first-arrival tomography with surface sources and receivers struggles in the presence of internal low-velocity layers when the so-called 'shingling' phenomenon occurs.In that case, early arrivals become discontinuous, forming several branches referred to as 'shingles' (Knox 1967 ;Shen et al. 2012 ).First, 'shingling' makes the arrivals harder to pick.Second, traditional eikonal solvers cannot model the correct traveltimes after the first discon- In this study, we use an eikonal solver as its fast computation of traveltimes enables Bayesian tomography in the first place.Several strategies can be used when picking shingling.The most obvious option is to disregard the algorithm shortcomings and invert all the picked traveltimes, including shingled branches.This approach may be acceptable when the shingling is not severe.Another option is to invert only the first 'shingle' and disregard the picked traveltimes after the first discontinuity, thus building a velocity model only down to the first velocity inversion.Alternatively, one may invert only the last 'shingle' , effectively replacing the upper part of the section with a different structure.Finally, other approaches exist to picking in severe shingling, e.g.picking a consistent velocity gradient that merges with a deep refractor (Diggins & Hampson 2021 ).Here, we follow the simplest strategy and keep all the picks in the tomography: this allows for a simplified automatic picking workflow without manual editing, at least for the synthetic dataset used.While it is theoretically possible to design an algorithm aware of such modeling errors (Meles et al. 2022 ) or change the problem to waveform inversion, such approaches remain too costly for practical applications.Therefore, we focus on a more straightforward approach based on an eikonal solver and demonstrate that it could be helpful in practical applications in models with velocity inversions.
After Bayesian tomography, we compute static corrections from each near-surface velocity model sampled by the rj-MCMC, analyze their variation and estimate associated structural uncertainty for a selected low-relief structure of interest.buried receivers (figure 2 c).We use 2000 m as a maximum lateral offset for picking on surface data and 1000 m maximum lateral offset for the buried array data.These traveltimes picked on the synthetic data are input traveltimes to the tomography algorithm.

Velocity uncertainty
First, we run two tomography experiments for the surface seismic acquisition to compare the linear interpolation based on Delaunay triangulation and natural neighbor interpolations based on Voronoi cells.We compare the mean of velocity models, their standard deviations (figure 3 ), traveltimes in the mean model (figure 4 ) and the stochastic velocity models (figure 5 ).One can observe that the standard deviations of the velocities identified by the linear interpolation are higher than the natural neighbor interpolation show-ing that the solution is less certain.Such behavior is likely due to the treatment of the model edges at the surface or because the chains have not fully converged in this case.The mean velocity model for the natural neighbor interpolation has a lower root-mean-square traveltime misfit (7 ms for natural neighbor vs. 14 ms for linear interpolation), showing that such a mean model is more accurate.This can also be observed by comparing the computed and picked traveltimes in Figure 4 .The histograms of the traveltime errors in the mean models shown in Figure 6 demonstrate overall lower errors of the natural neighbor parameterization.The histogram for the linear case also displays skewed behavior in the direction of negative errors, which suggests that this mean model provides overall larger traveltimes (slower velocities) compared to the ground truth.All of this indicates that the natural neighbor scheme possesses a lower 'parameterization error' (Hawkins et al. 2019 ). Figure 5 shows the stochastic models for the tested parameterizations.The boundaries between the triangles in the linear parameterization can be observed in the stochastic models.Also, while the stochastic models in figure 5 osci l late strongly (particularly at depth, where ray coverage is diminishing), they all provide a similar fit and possess identical traveltimes.The traveltime misfit plots for all the accepted samples in the chains are shown in figure 7 for both parameterization types.These plots suggest quicker and more stable convergence for the case of natural neighbor parametrization.
It is clear from the rapidly increasing velocity's standard deviations with depth that velocity inversions and lowvelocity anomalies complicate the inverse problem solution.Therefore, we add the first arrivals picked on the vertical receiver array data to the surface dataset to improve the velocity estimation at the deeper part of this near-surface model.Then, based on previous results showing the super ior ity of  the natural neighbor interpolation, we apply it to all remaining tomographic experiments.The mean velocity models and standard deviations corresponding to the tomography results with the added upholes are displayed in figure 8 .The mean velocity models show a more pronounced low-velocity anomaly at X = 5000 m (red arrow) corresponding to a karstified field and a low-velocity layer that starts at ∼200 m (blue arrow).With the addition of receivers at depth, the standard deviation at depth decreases.
These results demonstrate that the vertical arrays provide higher resolution and improve the accuracy of the estimated velocity models: the mean velocities are closer to their true values from the actual model, whereas the standard devia-tions are lower.Another way to confirm this is to assess the posterior probability distribution functions (PDF) estimated from the tomography between three types of acquisition and the true velocities (figures 9 and 10 ).PDFs for surface data show a rapid uncertainty increase below the low-velocity layers and anomalies (figure 9 ).The addition of vertical arrays improves the velocity estimate at depth and decreases the uncertainty, particularly for the velocity inversion occurring at a 200 m depth for X = 3000 and X = 4700 m (figure 10 ).The higher amount of detail provided by the tomography with added vertical arrays manifests itself in the increased number of nodes in the models estimated by tomography with the vertical arrays (figure 11 ).Note that even in the part of the model estimated with high accuracy according to the predicted uncertainty (very near surface), the tomography sti l l recovers only a smooth approximation of the actual velocity.Due to strong vertical heterogeneity, some portions of the true velocity profiles fall outside the estimated uncertainty ranges (figures 9 and 10 ).This effect is not unique to our results and can be observed in other works on transdimensional Bayesian tomography (e.g.Huang et al. 2021 ).This is caused by the nature of the considered tomographic problem capable of inverting relatively smooth velocity models without thin layers.Sti l l, whi le the uncertainty estimates do not assign the true structure a high probability, they are sti l l instrumental in assessing the results.In particular, near-surface tomography is often used to estimate static corrections.The multiple models provided by the algorithm, which all fit the observed traveltimes equally well, can be used to analyze these corrections and their impact on the mapping of deep structures.We provide an example of such an analysis in the next section.Considering a hypothetical low-relief structure with a vertical closure of 30 ms, the 170 ms uncertainty in the longwavelength statics solution obtained from only surface data makes such an anticline a very risky prospect.By contrast, the inversion of data augmented with vertical arrays possesses decreased uncertainty of 20 ms (less than the expected closure) and provides a significantly improved structure's location and shape.This demonstrates that the combined acquisition of surface seismic and vertical receiver arrays can decrease the risks when exploring for low-relief structures under a complex near surface.

Discussion
The primary limitation of Bayesian seismic tomography is its high computational cost.In the cases presented, a single chain requires ∼6 hours to run on one Intel Xeon 2.6 GHz processor core for surface acquisition in the aforementioned model.Overall, the current 2D tomography code, with optimized parameters, takes approximately eight hours to run on a 32-core machine.To simplify the prediction of full-scale 3D Bayesian tomography runtime, we compared the computation times of the 2D eikonal solver with the corresponding 3D solver.We generated 2D and 3D near-surface models with typical exploration seismic dimensions of 10 km lateral extent and a maximum depth of 1 km for both models.One shot point computation in 2D takes ∼0.007 seconds, on average, based on 100 computations.However, the corresponding computation time for the 3D version is ∼100 times longer, at 0.73 seconds.3D tomography requires modeling more shots (at least an order of magnitude higher comparing to 2D) and more complex models.The latter wi l l result in more unknowns required for the transdimensional method, leading to an additional increase in computational time, which may be the most significant factor.Estimating this factor in advance is challenging.However, assuming it is at least another order of magnitude, the computational complexity of 3D tomography can be considered at least four orders of magnitude higher than the 2D example presented.Although a 3D application was demonstrated by Zhang et al. ( 2018 ), the feasibility of using rj-MCMC tomography for exploration-scale 3D seismic studies is yet to be established.Fichtner et al. ( 2019 ) recently demonstrated that an alternative MCMC method, Hamiltonian Monte Carlo, may improve the inference speed.
This study is focused on estimating static corrections that are commonly used in time processing and imaging.Although widely used in practice, statics is a simplified approximation.A depth imaging study that aims to estimate structural uncertainty would require conducting full-scale migrations of a synthetic dataset with the numerous velocity models generated by the Bayesian tomography and analyzing the corresponding structural uncertainty.
The use of the assumed Gaussian distribution for uncorrelated picking errors may be considered simplistic, and further studies are necessary to determine the extent of the model's applicability.Likewise, incorporating first-break picking errors and distribution (manual and auto-picking algorithms) could enable more realistic uncertainty estimation.The emergence of machine-learning-based solutions for first-break picking (Duan & Zhang 2020 ) may facilitate further investigations in this area.
The parameterization and algorithms used in Bayesian inference for geophysical problems can significantly affect the results, as shown in studies by Zhang & Curtis ( 2022 ) and Zhao et al. ( 2022 ).To avoid subjective analysis of the results, Arnold & Curtis ( 2018 ) suggested interrogating the Bayesian solutions to answer general questions about subsurface structures.In this study, we compare and contrast the static correction distributions obtained from different parameterizations of the Bayesian tomography results, which can be considered a form of interrogation.

Conclusions
We investigate the efficacy of Bayesian first-arrival tomography for estimating near-surface velocity models in arid environments.The tomography yields a range of models that accurately match the observed traveltimes and identifies the velocity uncertainty at each location on the grid.The estimation process does not require an initial model except for broad prior intervals utilized to limit the unknown veloci-ties.We demonstrate the benefits of combining the Voronoi cell parameterization with the natural neighbor parameterization in modeling.Compared to the conventional Delaunay triangulation-based linear interpolation, our method generates stochastic models that are naturally smooth, with no overshoots or undershoots, and requires no additional parameters.Both parameterizations produce velocity estimates that fit the data.However, in the realistic SEAM Arid model example we demonstrated, the natural neighbor parameterization provides lower velocity uncertainty than linear interpolation.Furthermore, the mean velocity model in this parameterization provides a lower traveltime misfit showing that the results are more accurate.
The Bayesian tomography algorithm is applied to demonstrate the reduction in uncertainty resulting from the simultaneous acquisition of surface seismic data and vertical receiver arrays.Results indicate a significant decrease in velocity uncertainty and, consequently, the uncertainty of static corrections and imaging when such arrays are used together with surface data.The observed vertical uncertainty in two-way travel times is reduced from 170 to 20 ms, which is critical for delineating low-relief structures.Additional ly, unli ke conventional upholes, vertical receiver arrays record offset data from all sources in the survey, reducing imaging uncertainty beneath and away from the array locations.The results show that Bayesian tomography can assist in acquisition design by identifying the optimal placement of vertical arrays for maximum uncertainty reduction.

Figure 1 .
Figure 1.Vertical sections of true (a) and smoothed (b) 2D SEAM Arid velocity model used in the study.Part (a) also displays the datum (288 m) used for statics computation and the locations of 300 m deep vertical arrays used in the inversion.

Figure 2 .
Figure 2. Examples of the first-arrival picks (red) on common-shot seismic gathers (CSG) from the surface and buried data: (a) surface CSG without shingling, (b) surface CSG with shingling and (c) buried array CSG.

Figure 3 .
Figure 3. Mean velocity models (a, b) and standard deviations (c, d) for linear interpolation based on Delaunay triangulation (a, c) and natural neighbor (b, d) interpolation based on Voronoi tessellation.Only surface data are used for inversion.

Figure 4 .
Figure 4.A comparison of picked and computed traveltimes in the mean velocity model for linear (a) and natural neighbor (b) interpolations.
Using the section from the SEAM Arid velocity model (figure 1 a)(Oristaglio 2012 ; Bakulin & Silvestrov 2021 ), we create the 2D acoustic synthetic dataset to evaluate the Bayesian tomography algorithm.The model contains complexities typical of arid environments with multiple velocity inver-sions.The model includes alternating high and low-velocity layers (carbonates vs. sands/shales) and an additional lowvelocity zone in the central part of the model corresponding to a karst field.First, we model the seismic gathers with 250 m shot spacing and 25 m surface receiver spacing.Then, to estimate the impact of vertical receiver arrays, we add them to the acquisition with 500 m lateral spacing and 6.25 m receiver spacing along the vertical arrays.The maximum depth of the arrays is 300 m, and their geometry is shown in figure1a.Fi gure 1 b shows a smoothed version of this velocity model, provided for visual comparison with tomography results obtained next.The smoothing took place in the slowness domain, and the standard deviation of the Gaussian used for smoothing was chosen so that this model is visually comparable to the velocity models reconstructed by tomography.We pick first arrivals using a threshold-based algorithm on common-shot gathers with the surface (figure 2 a and b) and

Figure 5 .
Figure 5.A comparison of several stochastic velocity models from the linear (left) and natural neighbor (right) parameterizations.

Figure 6 .
Figure 6.Histograms of traveltime errors computed using mean models for linear (red) and natural neighbor (blue) parameterizations.

Figure 7 .
Figure 7. Comparing the traveltime misfit progression for the linear (red) and natural neighbor (blue) parameterizations.The increasing sample numbers represent the iterations.Observe smoother behavior and lower misfit achieved by the natural neighbor approach.The plot only displays the accepted elements (models) of the chains.

Figure 8 .
Figure 8.The mean velocity model (a) and standard deviations (b) resulting from the Bayesian tomography that jointly use surface data and 300-m deep vertical receiver arrays as shown in figure 1 a.

Figure 9 .
Figure9.A comparison between the true velocities versus depth (red) and the mean velocity obtained from tomography (solid green).The velocity errors are estimated as three standard deviations and shown as green fil l between dashed green lines.These are overlaid over the posterior probability density functions (PDFs, shown in grayscale colors) of velocity estimated by tomography at three lateral locations in the SEAM Arid model, using surface receivers only.

Figure 10 .
Figure 10.Same as figure 9 but for joint inversion of surface and vertical array (300 m) data.

Figure 12
Figure12compares the true static corrections for the chosen datum of 288 m (green line) with the static corrections derived from the stochastic velocity models obtained by the tomography.Static corrections corresponding to the mean velocity models are shown as red lines.The PDFs of the statics are plotted as the background image.We estimate the uncertainty as 99% symmetric confidence error bars taken from the PDFs and shown as dashed red lines.The considered datum of 288 m is below several low-velocity layers and anomalies, so the static corrections for this datum have high uncertainty in the case of surface acquisition (170 ms maximum uncertainty at X ≈ 1750 m in figure12 a).However, the true statics stay within the error bars, suggesting that the uncertainty estimates are meaningful.Note that the mean velocity model estimated by tomography sti l l provides reasonably accurate statics despite high uncertainty (i.e. the statics computed from the mean model shown as the red line are close to the actual statics in green).

Figure 11 .
Figure 11.Posterior probability density functions (PDFs) shown for the number of nodes in the model obtained from surface receivers only and from the combined acquisition with both surface receivers and vertical arrays of 300 m deep.The higher number of nodes in the latter case suggests that the subsurface models estimated by combined acquisition are better resolved with more detail.

Figure 12 .
Figure 12.The uncertainty of static corrections estimated from the stochastic velocity models obtained through Bayesian tomography is presented for different acquisition scenarios: (a) surface acquisition and (b) the simultaneous acquisition of surface seismic and vertical array data with a depth of 300 m.The datum used for these analyses is at 288 m.