[1] Three P wave attenuation models for sedimentary rocks are given a unified theoretical treatment. Two of the models concern wave-induced flow due to heterogeneity in the elastic moduli at ''mesoscopic'' scales (scales greater than grain sizes but smaller than wavelengths). In the first model, the heterogeneity is due to lithological variations (e.g., mixtures of sands and clays) with a single fluid saturating all the pores. In the second model, a single uniform lithology is saturated in mesoscopic ''patches'' by two immiscible fluids (e.g., air and water). In the third model, the heterogeneity is at ''microscopic'' grain scales (broken grain contacts and/or microcracks in the grains), and the associated fluid response corresponds to ''squirt flow.'' The model of squirt flow derived here reduces to proper limits as any of the fluid bulk modulus, crack porosity, and/or frequency is reduced to zero. It is shown that squirt flow is incapable of explaining the measured level of loss (10 À2 < Q À1 < 10 À1 ) within the seismic band of frequencies (1-10 4 Hz); however, either of the two mesoscopic scale models easily produces enough attenuation to explain the field data.

We present a method for estimating seismic attenuation based on frequency shift data. In most natural materials, seismic attenuation increases with frequency. The high-frequency components of the seismic signal are attenuated more rapidly than the low-frequency components as waves propagate. As a result, the centroid of the signal's spectrum experiences a downshift during propagation. Under the assumption of a frequencyindependent Q model, this downshift is proportional to a path integral through the attenuation distribution and can be used as observed data to reconstruct the attenuation distribution tomographically. The frequency shift method is applicable in any seismic survey geometry where the signal bandwidth is broad enough and the attenuation is high enough to cause noticeable losses of high frequencies during propagation. In comparison to some other methods of estimating attenuation, our frequency shift method is relatively insensitive to geometric spreading, reflection and transmission effects, source and receiver coupling and radiation patterns, and instrument responses. Tests of crosswell attenuation tomography on 1-D and 2-D geological structures are presented.

[1] Accurate characterization of fractured-rock aquifer heterogeneity remains one of the most challenging and important problems in groundwater hydrology. We demonstrate a promising strategy to identify preferential flow paths in fractured rock using a combination of geophysical monitoring and conventional hydrogeologic tests. Cross-well differenceattenuation ground-penetrating radar was used to monitor saline-tracer migration in an experiment at the U.S. Geological Survey Fractured Rock Hydrology Research Site in Grafton County, New Hampshire. Radar data sets were collected every 10 min in three adjoining planes for 5 hours during each of 12 tracer tests. An innovative inversion method accounts for data acquisition times and temporal changes in attenuation during data collection. The inverse algorithm minimizes a combination of two functions. The first is the sum of weighted squared data residuals. Second is a measure of solution complexity based on an a priori space-time covariance function, subject to constraints that limit radarattenuation changes to regions of the tomograms traversed by high difference-attenuation ray paths. The time series of tomograms indicate relative tracer concentrations and tracer arrival times in the image planes; from these we infer the presence and location of a preferential flow path within a previously identified zone of transmissive fractures. These results provide new insights into solute channeling and the nature of aquifer heterogeneity at the site.

Seismic and tracer test data can be combined to estimate the spatial patterns of aquifer properties. We present an algorithm that estimates the geometry of large‐scale lithologic zones, the effective hydraulic conductivities and seismic velocities for these zones, and the effective small‐scale dispersivity. The heart of this algorithm is our split inversion method, which extracts the geometry of lithologic zones from an estimated seismic velocity field. This method determines the zonation that best matches multiple types of data, such as seismic travel times and tracer concentrations. Although the current implementation of the algorithm uses only cross‐well seismic travel times and tracer concentrations, the algorithm could incorporate other data types that are sensitive to changes in large‐scale lithology. We demonstrate the approach for two synthetic sandy aquifers, one with interbedded clay lenses and another with interbedded silt and gravel lenses. These examples differ in the uniqueness of the relationship between seismic velocities and hydraulic conductivities. For these examples our algorithm successfully maps interwell heterogeneities and accurately estimates hydraulic and seismic parameters.

A new method is presented for identification of the permeability distribution in near‐surface aquifers. In addition to using the usual sparsely sampled pressure and permeability data, the method incorporates densely sampled seismic data, as obtained from a reflection or tomography survey, along with empirical relationships between seismic and hydraulic properties. The procedure is to first estimate by hydrologic inversion a pressure field. Then the velocity‐permeability‐pressure relationship is used to map the inverted pressure and measured seismic data to multivalued estimates of the permeability. Of those, the most probable value, based on the hydrologic inversion, is selected. In synthetic case studies the tremendous increase in coverage offered by the seismic data leads to dramatically better results in terms of both accuracy and resolution. An appealing feature is the use of relatively easy to acquire pressure data; a second is the incorporation of geophysical data which can sample an entire aquifer remotely without the need for an expensive and invasive drilling program.

The combination of differential radar tomography with conventional tracer and/or hydraulic tests facilitates high-resolution characterization of subsurface heterogeneity and enables the identification of preferential flow paths. In dynamic imaging, each tomogram is typically inverted independently, under the assumption that data sets are collected quickly relative to changes in the imaged property (e.g., attenuation or velocity); however, such "snapshot" tomograms may contain large errors if the imaged property changes significantly during data collection. Acquisition of less data over a shorter time interval could ameliorate the problem, but the resulting decrease in ray density and angular coverage could degrade model resolution. To address these problems, we propose a new sequential approach for time-lapse tomographic inversion. The method uses space-time parameterization and regularization to combine data collected at multiple times and to account for temporal variation. The inverse algorithm minimizes the sum of weighted squared residuals and a measure of solution complexity based on an a priori space-time covariance function and a spatiotemporally variable mean. We demonstrate our approach using a synthetic 2-D time-lapse (x, z, t) data set based loosely on a field experiment in which difference-attenuation radar tomography was used to monitor the migration of a saline tracer in fractured rock. We quantitatively show the benefits of space-time inversion by comparing results for snapshot and time-lapse inversion schemes. Inversion over both space and time results in superior estimation error, model resolution, and data reproduction compared to conventional snapshot inversion. Finally, we suggest strategies to improve time-lapse cross-hole inversions using ray-based inversion constraints and a modified survey design in which different sets of rays are collected in alternating time steps.

We evaluated a time-domain wave equation for modeling acoustic wave propagation in attenuating media. The wave equation was derived from Kjartansson's constant-Q constitutive stress-strain relation in combination with the mass and momentum conservation equations. Our wave equation, expressed by a second-order temporal derivative and two fractional Laplacian operators, described very nearly constant-Q attenuation and dispersion effects. The advantage of using our formulation of two fractional Laplacians over the traditional fractional time derivative approach was the avoidance of time history memory variables and thus it offered more economic computations. In numerical simulations, we formulated the first-order constitutive equations with the perfectly matched layer absorbing boundaries. The temporal derivative was calculated with a staggered-grid finite-difference approach. The fractional Laplacians are calculated in the spatial frequency domain using a Fourier pseudospectral implementation. We validated our numerical results through comparisons with theoretical constant-Q attenuation and dispersion solutions, field measurements from the Pierre Shale, and results from 2D viscoacoustic analytical modeling for the homogeneous Pierre Shale. We also evaluated different formulations to show separated amplitude loss and dispersion effects on wavefields. Furthermore, we generalized our rigorous formulation for homogeneous media to an approximate equation for viscoacoustic waves in heterogeneous media. We then investigated the accuracy of numerical modeling in attenuating media with different Q-values and its stability in largecontrast heterogeneous media. Finally, we tested the applicability of our time-domain formulation in a heterogeneous medium with high attenuation.

Reduced amplitude and distorted dispersion of seismic waves caused by attenuation, especially strong attenuation, always degrades the resolution of migrated images. To improve image resolution, we evaluated a methodology of compensating for attenuation (∼1∕Q) effects in reverse-time migration (Q-RTM). The Q-RTM approach worked by mitigating the amplitude attenuation and phase dispersion effects in source and receiver wavefields. Source and receiver wavefields were extrapolated using a previously published time-domain viscoacoustic wave equation that offered separated amplitude attenuation and phase dispersion operators. In our Q-RTM implementation, therefore, attenuation-and dispersion-compensated operators were constructed by reversing the sign of attenuation operator and leaving the sign of dispersion operator unchanged, respectively. Further, we designed a low-pass filter for attenuation and dispersion operators to stabilize the compensating procedure. Finally, we tested the Q-RTM approach on a simple layer model and the more realistic BP gas chimney model. Numerical results demonstrated that the Q-RTM approach produced higher resolution images with improved amplitude and phase compared to the noncompensated RTM, particularly beneath high-attenuation zones.