Bottom Shear Modulus Profiler

This invention relates to a method and apparatus for making measurements of a seabed and converting the measurements to a shear modulus profile of a seabed. More particularly, the invention relates to a method and apparatus for making use of the gravity waves found in the sea to measure the seabed shear modulus profile in a passive environment.

The determination of the shear modulus profile within the seabed is an important and challenging problem in geophysics, underwater acoustics, hydrodynamics and geotechnical engineering. Such information is needed for modeling the propagation of seismic waves, acoustic waves and gravity water waves in the ocean. Such information is also important for the design of foundations for offshore structures.

Existing direct methods using penetrometers and boreholes are costly and time consuming thereby being impractical for surveying a large area of seafloor. On the other hand, the bottom shear properties are very difficult to extract by indirect methods using acoustic waves because the shear waves are difficult to excite acoustically. Active seismic methods using interface waves are limited in that the necessary forward and inverse analyses are quite involved numerically.

It is an object of the present invention to provide a passive means for measuring the shear modulus profile of an ocean bed.

It is a further object of the present invention to provide a method and apparatus for measuring the shear modulus profile of a seabed without the need to disturb the seabed by extracting bores or inducing shear waves on the seabed.

To this end, an apparatus has been created to measure the shear modulus imposed on a seabed by naturally occurring surface waves.

In effect, the invention is comprised of a device housing three seismometers and a pressure transducer arranged to provide the seabed admittance (seabed displacement divided by wave pressure amplitude) produced by the surface or gravity waves that naturally occur on the surface of the ocean. In the preferred embodiment, the coordinates of the seabed are measured by the seismometers that are arranged in three dimensional coordinates. The displacement produced in the seabed by the surface waves are measured by the seismometer in units of velocity or acceleration of the seafloor over a discrete time inverval. A pressure transducer is provided to measure the pressure of the water induced by surface waves in the local area of the seabed being measured. The measured units are transmitted through a cable to an amplifier and ultimately to a computer, which is programmed to convert the velocity or acceleration measurements of the bed and pressure of the water into a shear modulus of the seabed.

FIG. 1 is a schematic drawing of the composite shear modulus profiler apparatus;

FIG. 2 is an elevational drawing of the bottom shear modulus profiler sensor taken through line 2--2 of FIG. 1;

FIG. 3 is a sectional elevational view taken in plan through line 3--3 of FIG. 2;

FIG. 4 is a sectional plan view taken through line 4--4 of FIG. 2;

FIG. 5 is a schematic view of the process of the subject invention;

FIG. 6 is a portion of a histogram of 3-components (X, Y, Z) of ground motion and wave pressure measured by the BSMP sensor;

FIG. 7 is a set of power spectra of X, Y, Z and P calculated from the measured histograms shown in FIG. 6;

FIG. 8 is a set of vertical and horizontal admittance spectra calculated from the power spectra in FIG. 7; and

FIG. 9 is a typical bottom shear modulus profile as measured by the apparatus of the present invention and calculated from the histogram of FIG. 6.

FIG. 10 is a definition of shear modulus profile (MODEL) and the incremental change for calculation of the kernel C.

FIG. 11 is a definition of the bed structure for numerical integration, and the two exact shear modulus profiles A' and C' used for the numerical tests.

FIG. 12 is a test for linearity between the changes in the bed admittances .xi. and .eta. and the incremental change in the shear modulus .DELTA.G at point B(z=-11 m) for the exact shear modulus profile A' and rigid basement.

FIG. 13 is synthetic data of bed admittance versus wave period used for numerical tests. Rigid basement for cases 1 and 2 and elastic half space basement for test case 3.

FIG. 14 is the numerical test results case 1; the initial shear modulus profile, the inverted shear modulus profile and the exact shear modulus profile. The number on the inverted profile indicates the number of iteration.

FIG. 15 is the numerical test results case 2; the definition is the same as FIG. 14.

FIG. 16 is the numerical test results case 3; the definition is the same as FIG. 14.

FIG. 17 is the spectrum of the bed admittance versus wave period measured at the Mississippi River Delta used in the inversion and the calculated spectra for the two inverted shear modulus profiles cases 1 and 2 obtained by inversion.

FIG. 18 is inversion results for the Mississippi River Delta, case 1; starting with initial profile 1. The number on the inverted profile indicates the number of iteration.

FIG. 19 is inversion results for the Mississippi River Delta, case 2; starting with initial profile 2. The number on the inverted profile indicates the number of iteration.

FIG. 20 is the shear modulus profile measured at boring site 5 about 5 km. east of the test site and inverted shear modulus profiles cases 1 and 2.

The embodiment of the present invention is described herein in terms of a bottom shear modulus profiler in which seismometers are used to provide bed velocity or acceleration data and a pressure transducer to measure the wave induced water pressure at the locality at which the velocity or acceleration data is being measured.

As best seen in FIG. 1, the apparatus of the present invention is comprised essentially of a Bottom Shear Modulus Profile sensor assembly 2, a cable 4 to deliver the sensor assembly signals to a computer 6 and an amplifier assembly 8 to amplify the signals obtained by the sensor assembly 2.

As best seen in FIG. 2, the bottom shear modulus profiler sensor assembly 2 is comprised of a waterproof cylindrical member 10 and a cap 12 in which seismometers 14, 16, 18 and a pressure transducer 20 are housed. The cap 12 is screwed into the threaded opening 13 in the cylindrical member 10 and provided with O-ring seals 15 to insure the waterproof condition.

In the preferred embodiment, the seismometers are arranged at 90.degree. angles to each other to form x-y-z coordinates. The seismometer 14 is secured to the bottom 22 of the cylindrical member 10 and extends perpendicularly upwardly from the bottom 22. Seismometers 16 and 18 are secured to the wall 24 of the cylindrical member 10 normal to each other and parallel to the bottom 22. In practice, it has been found that Teledyne Geotech Model S-500 Short Period seismometers with a modification in electrical circuits are suitable for the application.

Each of the seismometers 14, 16 and 18 are connected electronically through cables 50, 52, 54 respectively to a plug 26 in the cap 12. The exterior of the plug 26 connects to the common cable 4 (as shown in FIG. 4) that delivers the signals to the amplifier assembly 8.

The pressure transducer 20 is threaded into the cap 12 with one end 28 exposed to the sea water. The other end of the pressure transducer 20 is a plug 30 from which a cable 32 extends to connect to the plug 26. A Kistler Instrument Corp. Quartz Pressure Transducer Model 206 is used in the application.

The amplifier assembly 8 is comprised of three discrete amplifiers 34, 36, 38 to amplify the signals from the seismometers 14, 16, 18 and an amplifier 40 to amplify the signals from the pressure transducer 20. In practice it has been found that Analog Module Model LNVA/0.01 with 0.1 to 1.0 Hz bandwidth 40 & 80 dB BNC outputs and + and -15 volts power supply amplifiers are suitable to amplify the seismometer signals and an A/C voltage amplifier serves well to amplify the pressure transducer signals.

The process of the present invention proceeds by dropping the bottom shear modulus profiler 2 to the seabed and allowing it to rest on the seabed surface. The movement of the seafloor is detected by the seismometers 14, 16 and 18 which provide velocity readings in x-y-z coordinates, regardless of the orientation of the shear modulus profiler sensor assembly 2 on the seabed. The outputs of the seismometers 14, 16 and 18 are in volts and are converted to either velocity or acceleration.

The voltage signals are passed through the cable 4 to amplifiers 34, 36 and 38 where the signals are amplified and delivered to the computer 6 for analysis. Typically, the voltage output from the seismometers 14, 16 and 18 will be amplified in the range of 26 dB to 64 dB. Amplification is necessary to obtain a high enough voltage level to enable processing of the data.

The pressure transducer 20 measures the pressure of the sea water at the location at which the seismometer signals are being generated. The output of the pressure transducer 20 is also in volts and is converted at the computer 6 to Pascals or Newtons per square meter. The amplifier 40 typically amplifies the pressure transducer signal 20 db to 64 dB.

The computer 6 is programmed to convert the signals from the pressure transducer 20 and the seismometers 14, 16 and 18 into a shear modulus profile in accordance with the linear inverse theory combined with the Yamamoto theory of wave-seabed interaction described in detail in "Seabed Shear Modulus Profile Inversion Using Surface Gravity (Water) Wave-Induced Bottom Motion", Geophys Jour. Royal Astr. Soc. (1986) pp. 630-1 to 630-19 published May 1986 and incorporated herein by reference. Also incorporated herein by reference is Determination of Shear Modulus Profile By Inverse Analysis of Wave-Induced Bottom Motion by Tsuyosh Torii. A thesis submitted to the Faculty of the Univerity of Miami--Coral Gables, Fla., June, 1985.

As described in the above-cited Yamamoto article, "Seabed Shear Modulus Profile Inversion Using Surface Gravity (Water) Wave-Induced Bottom Motion", there are two types of hydrodynamic pressure generated by gravity water waves which can excite the seabed. The first type is the linear component which has the same frequency as the surface water waves and attenuate exponentially with depth. This component is usually referred to as the single frequency component. The second type is the non-linear component resulting from the non-linear interaction between waves of opposing directions of propagation. The type usually referred to as the double frequency component does not attenuate with depth and has twice the frequency of the surface water waves. This non-linear component is responsible for generating micorseisms as confirmed by observation.

In the above article, only the situation where the single frequency component of the wave-induced pressure dominates at the sea bottom is considered. This situation usually exists on the continental shelves and in the shallow seas. Recent measurements of the bottom pressure and the bottom motion in the deep ocean indicated that even the motion of the deep sea bottom at frequency lower than about 0.03 Hz is clearly due to the single frequency component of the bottom pressure which is induced by long-period gravity water waves.

In the theory of wave-seabed interaction, the forward problem is to calculate the motion of a seabed with prescribed layered structures of physical properties which is forced by gravity water waves of a given period propagating in a constant water depth. The problem is complicated in that the sediments are permeable, that the stress-strain relation of a sediment is quite non-linear and that the motion of the seabed affects the dispersion relationship of water waves.

Because the permeability of sediments is very small and because the frame structure of sediments is soft, the solid frame and the pore liquid within the sediments are completely coupled in their motion under gravity water waves. In addition, because the stiffness of pore water is much larger compared with the rigidity of sediment frame, the wave-induced motion of the water-saturated sediments is entirely due to shear deformation. The wave-induced shear strain in the sediments is usually quite small so that the sediments behave quite linearly except in the extreme conditions of very soft mud beds under large storm waves. Furthermore, because water waves have a very low propagation velocity compared with that of shear waves through most of the sediments, inertia forces are negligibly small compared with elastic restoring forces in the equations of motion. In other words, under usual circumstances the sediment behaves as a massless incompressible solid under water waves. These theoretical findings have been substantiated by laboratory experiments and field experiments. Therefore, the wave-induced motion of a sediment bed is independent of the density of the sediment but only dependent on the shear modulus of the sediment. This simple behaviour of a sediment bed under water waves is the basis for the bottom shear modulus profiler proposed.

The dispersion relation of the gravity water waves propagating over a sediment bed is given by Yamamoto (1982) as ##EQU1## Here, k=wavenumber; h=depth; .omega.=2.pi./T=angular frequency; T=wave period; g=gravitational acceleration; W=amplitude of vertical displacement at the bed surface; W'=amplitude of vertical seepage displacement of pore water at the bed surface and a=amplitude of surface water waves. Seepage displacement W' is usually negligible. As the vertical bed admittance (.eta.=W/a) is of the order of a maximum of 1 percent or so even for soft clay beds, the dispersion relation is reduced to the familar one given by

This means that the forcing of a seabed by water waves is completely prescribed if the amplitude spectrum of waves is measured.

The wave-induced motion of an arbitrarily layered bed may be represented by the coupling matrix. Let the displacement-stress vector r=(r.sub.1, r.sub.2, r.sub.3, r.sub.4) represent the depth functions of the horizontal displacement u, the vertical displacement w, the shear stress .tau..sub.xz, and the vertical normal stress .tau..sub.zz as

Then, the wave-induced motion of layered massless incompressible elastic bed is given as ##EQU2## in which G(z) is shear modulus which may be an arbitrary function of depth z. z is the vertical distance from the mudline and is positive in the upward direction. Note that the coupling matrix is made up of wavenumber k and shear modulus G(z) alone. This simple structure of the coupling matrix together with the explicit dispersion relation given by equation (2.2) provides stability and precision for inversion.

Various methods are available for solving the eigenvalue problem for the displacement-stress vector involving the coupled linear differential equations such as equation (2.4). For example, the coupled linear differential equation (2.4) may be integrated directly by numerical methods as an initial value problem. A fourth-order RungeKutta method is adopted as our numerical integration technique, because the method is simple and well suited for modelling a continuously layered seabed. Two independent sets of initial values may be defined at the basement (z=-d).

If the basement at z=-d is a rigid boundary, the two initial vectors may be given as

On the other hand, if the basement at z=-d is an elastic half-space, the two initial vectors may be given as

The displacement-stress vector in the sediment bed is given as a linear combination of the two particular vectors as

The proportionality constants A and B can be determined from the boundary conditions at the seabed interface. At the interface (z=0) shear stress vanishes;

and vertical stress must be balanced by wave-induced bottom pressure P(0) as

in which P(z) is the wave-induced pressure amplitude at z=z as

in which .rho.=water density. In the experiments, bottom pressure may be measured directly or deduced from measured wave amplitude through equation (2.10).

The linear inverse problem most often encountered in geophysics can be written in the discrete form, ##EQU3## or in matrix notation

where e.sub.i is one of N observations (DATA), C is a kernel derived from theory, and m is the unknown function (MODEL) defined on the I interval.

The singular decomposition of matric C is used to solve equation (3.3) following Lanczos (1961).

U.sub.i, the columns U are the N-dimensional eigenvectors of a coupled eigenvector-eigenvalue problem which spans the data space. V.sub.i, the columns V are the M-dimensional eigenvectors of a coupled eigenvector-eigenvalue problem which spans the solution (MODEL) space. is the diagonal matrix with the corresponding eigenvectors .lambda..sub.i which are sorted into descending order. Zero valued eigenvectors are not included, In particular

and the U.sub.i and V.sub.i are orthonormal such that

In terms of U, V and , the solution to equation (3.3) is ##EQU4##

For the data with noise, the following is considered for the construction of the solution using the truncated singular expansion, ##EQU5##

There are good reasons for using the truncated expansion. With the assumption of Gaussian white noise one can show that the uncertainty of each coefficient multiplying V.sub.i in equation (3.7) is proportional to .lambda..sub.i.sup.-1/2. Thus the expansion equation (3.7) is in terms of vectors whose coefficients increase in uncertainty. It is also usually the case that the vectors V.sub.i become more oscillatory as i increases (exceptions to this rule are very rare). Thus, the truncated expansion equation (3.8) simultaneously removes those parts of the solution which are sensitive to noise an gives smoothed estimate of the solution. This point is demonstrated in this paper using synthetic data without noise and measured data with noise.

The resolution matrix has been adopted to determine the quality of the above solution using the truncated singular value expansion ##EQU6## In the ideal case, in which C.sup.-1 C results in the identity matrix, the resolution is perfect. When C.sup.-1 C has near unity values on the diagonal and small values elsewhere, the computed solution represents a smoothed version of the true solution.

All of the above analyses apply to the linear inverse problem. For the non-linear problem such as the one under consideration, an iterative linearization scheme has been used by which successive approximations to the model can computed. Starting with the non-linear relation between data e and model m, and recalling equation (3.3)

In order to minimize .phi.=.vertline.e-f(m.sub.o).vertline..sup.2, it is necessary to assume an initial value m.sub.o and then linearize. Expanding equation (3.10) in a Taylor series form: ##EQU7## The derivatives .delta.f/.delta.m are approximated by finite differences in the last equation. As will be shown in the following section, the finite difference approximations to the derivatives are quite good for our problem.

Our problem is to extract the shear modulus profile, G(z), within the seabed from the measured bed admittance spectra. Let us define the horizontal bed admittance, .xi., and the vertical bed admittance, .eta., by

Here U and W are amplitudes of horizontal and vertical displacement of the bed at the mudline z=0 and a is the amplitude of the water wave. For our problem, equation (3.12) must be expressed in terms of the bed admittances .xi., .eta. and the bed shear modulus G. Using the finite difference approximations to derivatives .delta..xi./.delta.G and .delta..eta./.delta.G, based on small discrete intervals .DELTA.G=.alpha.G, equation (3.12) can be expressed as, ##EQU8## or equation (3.14) can be written as: ##EQU9## In equations (3.14) and (3.15), .xi..sub.j and .eta..sub.j are the data of the horizontal and vertical bed admittances for the jth wave. .xi..sub.jo and .eta..sub.jo are the calculated horizontal and vertical bed admittances of the seabed with the initial shear modulus profile (G.sub.1, G.sub.2, . . . , G.sub.i, . . . , G.sub.m) for the jth wave. .xi..sub.ji and .eta..sub.ji are the calculated horizontal and vertical bed admittances of the seabed with the shear modulus profile with an incremental change .DELTA.G.sub.i =.alpha.G.sub.i at the jth grid point as shown in FIG. 10 for the jth wave, G.sub.i is the initial shear modulus at the jth grid point. .alpha. is a small real number. .alpha.=0.01 is used to calculate the kernel given by the matrix in equation (3.14). (G.sub.1, . . . , G.sub.2, . . . , G.sub.i, . . . , G.sub. m) is the unknown shear modulus profile to be found by inversion. (S.sub.i, S.sub.2, . . . , S.sub.i, . . . , S.sub.M) is the function and expressed as: ##EQU10## then

M is the number of unknowns in the shear modulus profile in the bed and N(=2n) is the number of admittance data induced by water waves at n(=N/2) different wave periods.

For the numerical tests and the inverse analysis of the Mississippi River Delta data shown in the following sections, the matrix equation (3.15) consists of the data space N=20, the model space M=12, and the eigenvalues .lambda..sub.i and the eigenvectors V.sub.i (i=1,2, . . . , 12). Note that this is an over-determined system (N>M). The under-determined systems (M<N) and the completely determined systems (M=N) can be also solved by this method without modification.

Equation (3.14) can be shown in matrix notation as .DELTA.e=C m.(3.17)

The linear system expressed by equation (3.17) is solved for .DELTA.m using equations (3.7) and (3.8), and then repeating the process after replacing m.sub.o by m.sub.o +.DELTA.m. This process is repeated until .DELTA.m becomes acceptably small. In practice the method does not always converge.

The numerical tests were conducted where the consistency and the uniqueness of the inversion procedures described in the previous section are ascertained by numerical examples using the synthetic data without noise which are generated by the forward theory for two prescribed shear modulus profiles and basement boundary conditions. To begin with, we investigate the linearity of the forward theory relative to its acceptability for the linear inverse theory. The bed structure is defined in FIG. 11 about the layered system, for the grid system used for numerical integration, the bed thickness and the basement system. This bed structure is used throughout the numerical tests and the example applications presented in this paper. The model used consists of a wave with T=6.0 s, a=1.0 m and d=19.2 m and the shear modulus profile type A' on a rigid basement defined in FIG. 11. The 50 m thick bed is represented by 51 grid points at a constant interval of 1 m in the numerical integration of equation (2.4).

The method used to test the linearity of the forward theory is the same as the method to generate the kernel C given by equation (3.14). The increment defined in FIG. 10 is varied from 0.01 to 8.0 to the positive side and -0.01 to 0.96 to the negative side. As an example, the relation of .vertline..DELTA..xi..vertline., .vertline..DELTA..eta..vertline. versus .vertline..alpha..vertline.=.vertline..DELTA.G/G.vertline. evaluated at point B(z=-11 m) is in FIG. 12. Here, .DELTA..xi. and .DELTA..eta. are the changes in the horizontal and vertical admittances respectively at the mudline (z=0) due to the incremental change in the shear modulus at point B. We have found the similar linearity relations as FIG. 12 exist for all the other grid points.

From FIG. 12, it is possible to conclude that the forward theory is appropriate for the linear inversion if .alpha.>-1.0, where linearity exists in .vertline..alpha..vertline.<0.1 and weak non-linearity occurs in -1.0<.alpha.<-0.1 and 0.1<.alpha.<8, approximately. Note that the kernel C approximated by finite differences given by equation (3.14), which is calculated using .alpha.=0.01, should be quite a good approximation to derivatives. The three cases (numerical test cases 1, 2 and 3) are examined to determine the uniqueness and the consistency of the solution by the inversion. The two exact shear modulus profiles are shown in FIG. 11. For the two exact shear modulus profile-basement combinations (type A'/rigid bed and type C'/elastic half-space), the horizontal and vertical bed admittances at the bed surface (z=0) are calculated for the following 10 waves;

The synthetic data without noise of the bed admittances so calculated are plotted in FIG. 13. It is important to note that these bed admittance values without noise are on the smooth curves against wave period T.

To start the linearization-iteration procedure of the inversion, an initial shear modulus profile is needed. A good initial profile may be determined from an inverse relation between the vertical bed admittance and the shear modulus for the homogeneous poro-elastic half-space:

where z=1/k=L/2.pi. and e=2.718.

The inverted shear modulus is truncated at the Lth mode, if ##EQU11## (a) the inverted shear modulus G.sub.L, that is, where G.sub.o is the initial shear modulus and G.sub.L is the inverted shear modulus;

(b) the eigenvalue .lambda..sub.L varies largely, that is

(c) G.sub.L profile starts to oscillate.

When .DELTA.G/G.sub.o <0.1 (or .DELTA.G/(0.01G.sub.o)<0.1.times.10.sup.2), however, the model advances to the higher mode L'(>L). The iterative linear inversion procedure outlined in the previous section is implemented for the numerical test case 1 and also illustrated in FIG. 5. In the inversion, the values of shear modulus at 12 elevations z=-1, -3, -5, -8, -11, -14, -18, -22, -26, -30, -40 and -50 m are varied as independent variables. The shear modulus values at the grid points in between these points are linearly interpolated. The initial shear modulus profile shown in FIG. 14 is a good initial shear modulus profile which is obtained by equation (4.1). The entire process of how the iterative linear inversion procedure corrects the initial profile through 11 successive interactions and finally converges to the exact profile is illustrated in FIG. 14 for numerical test case 1.

In order to ascertain the uniqueness and the consistency of the inversion procedure, another and a very poor initial shear modulus profile is used in the inversion procedure (numerical test case 2) as shown in FIG. 15. As clearly illstrated in FIG. 15, the inversion procedure finds the exact profile again after 10 successive iterative inversions indicating that the proposed inversion procedure is not only unique but also usually consistent. This unusual behaviour of the inversion procedure is due to the extremely simple and linear relation between the bed motion and the bed shear modulus structure represented by the forward model given by equation (2.4).

We have found that the uniqueness and consistency of the inversion procedure hold for the case (numerical case test 3) of the basement of an elastic half-space as shown in FIG. 16.

It can be concluded that if perfect data .xi. and .eta. (i.e. .xi. and .eta. without error and on the smooth curve against T) exist, the shear modulus profile can be determined consistently and uniquely.

Therefore, as a next step, the real data .xi. and .eta. from Mississippi River Delta (i.e. .xi. and .eta. with error) are examined. Continuous simultaneous measurements of bottom oscillations and wave characteristics were made at the East Bay of the Mississippi River Delta. Wave staffs, pressure sensors and 3-component bottom accelerometers were placed in about 19 m of water at an oil platform in an area having a fine-grained clay bottom. The accelerometer package was placed about 0.3 below the mudline.

The results of the experiments indicated that the bottom undergoes an elastic response to wave-induced bottom pressure in a quite linear manner. The amplitudes of the wave-induced pressure at z=13.2 m, the bottom horizontal acceleration and the bottom vertical acceleration under 10 waves of nearly sinusoidal shapes were depicted by Dawson, Suhayda & Coleman (1981). Using this set of data, the horizontal bed admittance (.xi.) and the vertical bed admittance (.eta.) are calculated for the 10 waves and are plotted in FIG. 17. It should be noted that the signal-to-noise ratio of the data is of the order of 3-4. These very noisy values of the bottom admittance are used as input data for the inversion analysis.

The iterative linear procedure outlined in the previous section is implemented for the bed admittance data and the results are illustrated in FIG. 18. The grid system of the upper 50 m of the bed is represented at a constant integration step of 1 m as before. The basement below -50 m is assumed to be an elastic half-space whose shear modulus is equal to that of the grid point at z=-50 m. The numerical integration procedure and the inversion procedure are the same as for the numerical tests.

FIG. 18 shows how the initial shear modulus profile corrects itself through eight successive iterative inversions and finally converges to the final profile. Because of the large noise, only the first two, three or four modes are included in the solution.

In order to evaluate the uniqueness and the consistency of the inverted shear modulus profile from noisy data, another independent iterative inversion procedure is conducted using a completely different initial shear modulus profile as shown in FIG. 19. Again only the first three, four or five modes are included in the solution because of the large noise. Comparison between FIGS. 18 and 19 indicates that the overall feature of the two final profiles are very similar except near the peaks of about z=-18 m. Both profiles clearly reveal the existence of a very soft sublayer below about z=-25 m. This is an indication of the consistency of the inversion. The fact that the two final profiles differ from each other weakens the uniqueness of the inverted solution to some degree. This is largely due to the noise or the scatter in the measureed data as may be observed in Fig. 17. The calculated bottom admittance spectra for the two inverted final shear modulus profiles are compared with the measured data in FIG. 17. It may be concluded that the iterative inversion procedure tends to smooth the noisy data and converges to one of many possible averaged lines through data. This indicates some degree of uncertainty introduced by the noise in the real data which is absent in the inversion using the synthetic data without noise.

The geological characteristics of the core taken in the boring from this area was determined. FIG. 20 shows the shear strength profile in the bottom taken at boring 5 which is about 5 km east of the experimental site. Note that a approximately linear strength profile exists above a thick sublayer of soft mud below about z=-25 m. The two final shear modulus profiles obtained from the inverse analysis are compared with the shear profile obtained from the boring in FIG. 20. A conversion factor C=200 is used to convert shear strength into shear modulus for the comparisons. It must be noted that the boring site and the test site are separated by about 5 km so that the comparisons are only very approximate. Nonetheless, good overall agreement between the inverted profiles and the profile by direct measurements is obvious. In particular, the very soft sublayer below about z=-25 m revealed by the inverse analysis is most remarkable.

Having no control over the excitation mechanism (natural water waves), it may be useful to assess the depth of penetration and the resolution one can obtain for the shear modulus profile using this method. As may be seen in euquation (2.4), the depth of penetration is proportional to the wavelength of water waves that excite the seabed. The resolution of the shear modulus profile depends on the bandwidth and the resolution of the bottom admittance spectrum. Therefore the penetration depth and the resolution of the shear modulus profile depend on the water depth and surface wave spectrum at the time and place of measurement. We can only assess these qualities based on the actual data. In the case of Mississippi River Delta experiments discussed, the bandwidth of wave period covered is from 5.2 to 8.0 s which corresponds to the bandwidth from 42 to 80 m for a water depth of 19.2 m. As shown in FIG. 20, the penetration depth is 50 m and may be extended a bit more. Therefore one may say that the penetration depth of this method is approximately one wavelength. According to FIG. 20, the resolution is the best (say 1 m) near the mudline and decreases with depth. Therefore, one may say that the resolution of the shear modulus profile is a few percent of the wavelength.

In concluding the referenced article, the authors stated that an amplitude inversion scheme has been developed to extract the bottom shear modulus profile from measured motion of the seabed surface at one point induced by ocean surface waves using the linear inverse theory combined with the Yamamoto theory of wave-seabed interaction. The uniqueness and the consistency of the inversion has been established by the numerical tests using synthetic data without noise. As an example application, the noisy data obtained at the Mississippi River Delta has been used in the inversion to extract the bottom shear modulus profile. A good agreement is obtained between the inverted shear modulus profile and the direct measurements.

As seen in FIG. 5, the method of the invention proceeds in essentially four basic steps.

An example of the operation of the invention is found in the following wherein: the shear modulus sensor assembly 2 was deposited in 4.11 meters of water at the Great Bahama Bank on Nov. 24, 1985:

Step 1: Signals of seabed motion velocity (X, Y and Z components), and wave pressure (P), in volts were measured by the invention BSMP sensors, amplified and recorded on a histogram (time history) as shown in FIG. 6.

Step 2: The histograms of X, Y, Z and P signals were digitized at time intervals of 0.1 second and stored in a minicomputer MINC C-11. The digital time series data of X, Y, Z and P was then used to calculate the power spectra (voltage amplitude versus period) of X, Y, Z and P as shown in FIG. 7.

Step 3: The voltage power spectra X, Y, and Z of seabed motion velocity components were converted into the power spectra of seabed motion in displacements (X, Y and Z in meters instead of volts) by multiplying conversion constants. At the same time, the voltage power spectrum (P) of wave pressure was converted into the water wave amplitude power spectrum (P in meters instead of volts). Using the amplitude spectra of Z and P, the vertical admittance spectrum and the coherence were calculated as shown in FIG. 8. Likewise, the horizontal admittance spectrum and the coherence were calculated from the amplitude spectra X, Y and P. The admittance is the ratio of the seabed displacement amplitude to the water wave amplitude, also in FIG. 8. The coherence is an indication of the quality of the admittance, the higher coherence the better quality.

Ten periods ranging from 3.0 to 5.7 seconds and corresponding values of vertical admittance and horizontal admittance were depicted from FIG. 8. This range of wave period was used because of the high coherence values. Table-I shows the ten wave periods and the corresponding admittance values:

Step 4: Thereafter a shear modulus profile, as shown in FIG. 9, was developed from the admittance data in Table-I, applying the Yamamoto theory of wave-seabed interaction and the linear inverse theory described in Reference 1 ("Shear Modulus Profile Inversion Using Surface Gravity (Water) Wave-Induced Bottom Motion", Geophys. Jour. R.A.S. (1986) pp 630-1 to 630-19). The following Yamamoto Equation was used to calculate the wave-induced seabed motion.

r.sub.1 =horizontal displacement amplitude at Z=Z

r.sub.2 =vertical displacement amplitude at Z=Z

r.sub.3 =shear stress amplitude at Z=Z

r.sub.4 =vertical normal stress amplitude at Z=Z

Z=depth below mudline

k=wave number

G=shear modulus at Z=Z

The foregoing Yamamoto equation is a shorthand equation representing the following four equations:

U=horizontal displacement of amplitude of seabed as a function of depth z.

W=vertical displacement of amplitude of seabed as a function of depth z.

.tau.=shear stress amplitude within seabed funcation of depth z.

.sigma.=vertical normal stress amplitude within seabed as a function of depth z.

z=depth of burial in seabed measured from mudline.

G=shear modulus within seabed as a function of depth z.

The Yamamoto equations relate the shear modulus profile (G vs. z) and the seabed motion U and W and the vertical stress .tau. (=-P) at the mudline which are induced by ocean waves. Therefore, Yamamoto Equations incorporated with the linear inverse theory will determine the shear modulus profile (G vs. z) from the seabed motion and wave pressure (U,W,P) which are measured by the BSMP sensors. A detailed procedure of this inverse calculation is given above. The shear modulus profile in the Bahama Bank obtained from this inversion procedure is shown in FIG. 20.