###### Abstract  
In many regions, particularly coastal areas, population growth, overuse of water, and climate change have put quality and availability of water under threat. While accurate, predictive groundwater flow models are essential for effective water resource management, the precision of these models relies on the ability to determine the geometries of geological structures and hydrogeologic systems accurately. In complex geological settings or with deep aquifers, the drilling of observation wells becomes too costly and shallow seismic surveys become the method of choice. Common-Reflection-Surface stacking is being used by major oil companies for hydrocarbon exploration but can serve also as an advanced imaging method for near-surface seismic data. Its spatial stacking aperture covers a whole group of neighboring common midpoint gathers and, as such, a multitude of traces contribute to every single stacking process. Since the level of noise suppression is proportional to the number of contributing traces, Common-Reflection-Surface stacking generates a large increase in signal-to-noise ratio. In addition, the data-driven velocity analysis is a statistical process and is, as such, the more stable the more input traces it has. At the beginning, however, the spatial operator was only used for stacking, not for velocity analysis, since limiting computational demand was mandatory to obtain results within a reasonable time frame. Today's computing facilities are thousands of times faster and even large efficiency gains do not justify the loss of effectiveness anymore that comes with a truncated velocity analysis. We show that this is particularly true for near-surface data with low signal-to-noise ratio and modest common midpoint fold. For the spatial velocity analysis, we present two options: (1) as reference, a global search of all three parameters of the Common-Reflection-Surface operator, and (2) as a quicker solution, a strategy that uses the two-parameter Common-Diffraction-Surface operator to obtain initial values for a local three-parameter optimization. For shallow P-wave data from a hydrogeological survey, we show that the computational cost of option (2) is one order of magnitude smaller than the cost of option (1), while the stack and corresponding normal-moveout velocities are very similar. Comparing the results of the spatial velocity analysis to those of preceding, computationally lighter, strategies, we find a significant improvement, both in stack section resolution and stacking parameter accuracy. seismic methods; groundwaters; near-surface; CRS stack; velocity analysis; multi-parameter optimization +
Footnote †: journal: Applied Sciences  
structure, which extends vertically below the acquisition line. Since data is recorded shot-by-shot, a dataset initially consists of a sequence of so-called shot gathers. However, it is common practice in seismic imaging to reorganize the traces into so-called common midpoint (CMP) gathers, each composed of records from different experiments which have source and receiver centered around the same midpoint, but a variable source-receiver distance, called offset. In the very basic scenario of having a planar measurement surface and a parallel planar interface separating two homogeneous subsurface layers, we would record the same trace for every midpoint location. The reflection traveltime would depend only on the offset. If we represent all measured traces in a midpoint, offset, traveltime space, all reflection amplitudes would lie on a plane with hyperbolic curvature in offset direction and zero curvature in midpoint direction. The shortest reflection traveltimes would be found in the smallest-offset plane. If we had a reflector dip relative to the measurement surface, we would still find the same traveltime surface but tilted with respect to the midpoint axis by an angle proportional to the reflector dip.