Study of the factors affecting the karst volume assessment in the Dead Sea sinkhole problem using microgravity field analysis and 3-D modeling

Thousands of sinkholes have appeared in the Dead Sea (DS) coastal area in Israel and Jordan during two last decades. The sinkhole development is recently associated with the buried evaporation karst at the depth of 25– 50 m from earth’s surface caused by the drop of the DS level at the rate of 0.8–1.0 m/yr. Drop in the Dead Sea level has changed hydrogeological conditions in the subsurface and caused surface to collapse. The pre-existing cavern was detected using microgravity mapping in the Nahal Hever South site where seven sinkholes of 1–2 m diameter had been opened. About 5000 gravity stations were observed in the area of 200 ×200 m2 by the use of Scintrex CG3M AutoGrav gravimeter. Besides the conventional set of corrections applied in microgravity investigations, a correction for a strong gravity horizontal gradient (DS Transform Zone negative gravity anomaly influence) was inserted. As a result, residual gravity anomaly of –(0.08 ÷0.14) mGal was revealed. The gravity field analysis was supported by resistivity measurements. We applied the Emigma 7.8 gravity software to create the 3-D physical-geological models of the sinkholes development area. The modeling was confirmed by application of theGSFCprogram developed especially for 3-D combined gravity-magnetic modeling in complicated environments. Computed numerous gravity models verified an effective applicability of the microgravity technology for detection of karst cavities and estimation of their physical-geological parameters. A volume of the karst was Correspondence to: L. V. Eppelbaum (levap@post.tau.ac.il) approximately estimated as 35 000 m 3. The visual analysis f large sinkhole clusters have been forming at the microgravity anomaly site, confirmed the results of microgravity mapping and 3-D modeling.


Introduction
The Dead Sea (Fig. 1) is situated within a unique region, which has interest for many scientific disciplines -above all for Earth sciences.At about 415 m below the mean sea level, the Dead Sea is the lowest point on the face of the Earth.It is located at the centre of the Dead Sea Transform Zone (DSTZ), a large feature that extends over 1000 km from the southern tip of the Sinai Peninsula to the Taurus Mountain in Turkey, and is the most prominent tectonic feature in the Middle East (Ben-Avraham, 2001).
Sinkholes developed during two last decades along the Dead Sea western and eastern shorelines in Israel and Jordan is the main concern of the region.The sinkhole development is recently associated with the buried evaporation karst at the depth range of 25-50 m from surface caused by the drop in the Dead Sea level at the rate of 0.8-1.0m/yr.Drop in the Dead Sea level has changed hydrogeological conditions in the subsurface and caused surface to collapse (Yechieli et al., 2006).Ezersky et al. (2005) proposed to develop an integrated approach for prediction of natural hazards caused by the development of sinkholes in the Dead Sea region of Israel and Jordan through the joint application of geophysical and hydrogeological studies.One of the problems resolved in the framework of this investigation is an improvement of the microgravity method.In this article we consider the possibilities of 3-D microgravity modeling aimed to estimate qualitative and quantitative parameters of buried karst cavities, detect the buried salt dissolution caverns and verify our calculations with in-situ geodetic measurements.Our studies demonstrate the possibilities of gravity field analysis on models and allow to estimate the real parameters of sinkholes formed in the Nahal Hever South area (western Dead Sea coast).We performed also geodetic mapping of the sinkhole cluster (formed at the place of detected microgravity anomaly) to verify the results of microgravity mapping and 3-D gravity modeling.

Site geology: a brief description
A group of sinkholes about 20-50 m apart appeared in the Nahal Hever South area in July 1998 (Fig. 2).These sinkholes (0.5-2.0 m deep and up to 3 m in diameter) did not change their appearance over approximately 8-9 months, but then their sizes began to vary over some time (Fig. 3).
The geological section of this area is composed of alluvial fan sediments down to 18 m deep, a marl layer of 5 m thick and a salt layer of 11 m thick (a generalized geological section is presented in Fig. 4).A 5 m clay and gravel layer underlies the salt layer.The elevation of the area is approximately -393 m.The top of the salt layer is located at the depth of 24 m.In one of boreholes, a cavity filled by dense mud was detected at the depth of 23-29 m, which is assumed to be the result of the salt layer dissolution.The visible section of sinkholes consists of sand-gravel intercalating claysand layers; iron oxide mineralization is also visible (Fig. 4).(after Yechieli, 2002).

Ring structures identification and classification
Ring structures (RS) phenomenon is widely presented in the Earth's environments (Fig. 4).The RS may be generally classified as Terrestrial, Extraterrestrial (natural RS) and Archaeological (artificial RS).In the developed block-scheme (Fig. 5) natural RS are classified using Khain (1995) ar-rangement and the author's developments (Eppelbaum et al., 1998).Taking into account the majority of natural and artificial RS, sometimes springing up a problem of not only RS identification, but also of it correct classification.The karst (sinkhole) phenomenon in the RS scheme (Fig. 5) is classified as a component of erosional RS.

A general methodology of RS delineation
The first problem arising by RS geophysical studying is their revealing (delineation) against the background noise.Very frequently it is difficult to single out RS in complex geological-geophysical environment, especially taking into account the commonly employed rectangular network of geophysical observations.Apparently, Aitken (1961) firstly proposed the idea of circular graticules set application in detailed magnetic survey for RS delineation.A computing of directional derivatives is realized (in a few various modifications) in the well-know software WinSurfer (Golden Software).
Obviously, the most clearly approach of summing up horizontal gradients (we will designate horizontal gradient as U x ) using a circular "apparent" graticule was presented in Khesin et al., 1983, 1996 (Fig. 6).The apparent graticule radii drawn with the interval of 45 • determine horizontal gradients (Fig. 6C).When summing up the gradients in various horizontal directions, the presence of circular features should be intensified, whereas other signals are leveled.Here the correlation of the sum of gradients (or the average gradient) for a circle with the radius R n and a ring external to this circle limited by R n and R n+1 radii makes it possible to determine whether the circular feature revealed reflects the centric Fig. 5. Classification of ring structures in the Earth's environments (after Eppelbaum et al., 1998;Eppelbaum, 2007a).
or ring structure (Fig. 6D).The value of n i=1 U x →0 inside a circle in the absence of a centric structure (Khesin et al., 1996).Application of this method is explained on a model of the inclined circular cylinder magnetized along its dipping (Fig. 6A, B).
The abovementioned procedure may be significantly modified on the basis of so-called "multimodel approach" (Eppelbaum, 2005;2008) by the way of performing geophysical measurements at different levels over the earth's surface.Undoubtedly, computing difference between the horizontal gradients of geophysical fields along z-axis (we will designate this parameter as U xz ) will provide additional important information about the studied geological section.
It should be noted that the present thriving developing of Remote Operated Vehicle (ROV) industry (which allowing cheaply and operatively carry out geophysical measurements at low altitudes over the earth's surface) makes it possible to apply this methodology for rapid and economic studying different geological phenomena (Eppelbaum, 2008) including investigation of karst terranes.

Field investigations and preliminary data processing
To reveal potential of the microgravity method, the feasibility study was carried out throughout the Nahal Hever South sinkhole development area (Rybakov et al., 2001).About 5000 gravity stations were observed in this area by the use of Scintrex CG-3M AutoGrav gravity meter.Microgravity mapping of 200×200 m 2 area (with the observation step of 3 m) was carried out in 1999 after the first sinkholes appeared.Location of the gravity observations points in the Nahal Hever South area is displayed in Fig. 7.All necessary corrections (including complete terrain correction) were inserted (Rybakov et al., 2001).The residual microgravity anomaly of -(0.8÷0.14)mGal, extending in the north-south direction, was revealed in the area.The anomaly consists of three sections with sizes from about 20×20 m 2 to 50×50 m 2 .
These investigations were carrying out in the zone of strong horizontal gradient ( g B(RHG) ) (up to 10 mGal/km) caused by the influence of the regional negative Dead Sea gravity minimum (Fig. 8).The values of the g B(RHG) were polynomial computed and subtracted from the observed g B map (Rybakov, 2007).The residual scheme of the horizontal g B gradient for the Nahal Hever South area is presented in Fig. 9.As shown in this figure, g B(RHG) reaches of high values for such a small area.
The residual gravity map is shown in Fig. 10B.The resistivity map (Fig. 10A) illustrates the resistivity distribution at a depth between 10 and 15 m (constructed on the basis of resistivity field modeling).Two anomalies can be clearly seen in the both maps.The first (No. 1) is a low resistivity anomaly at the center of the resistivity map (Fig. 10A) and corresponding negative gravity anomaly in the center of gravity map (Fig. 10B).The second similar pair of anomalies (marked as No. 2 in Fig. 10A) is legibly seen in both these maps.Most of the early sinkholes (before 1999) were formed within the area of anomaly No. 1.
Figure 11A, B repeats the maps presented in Fig. 10A, B, but with location of the dissolute salt contour.The previous surveys enable us to reach a complete understanding of the physical-geological model of the area and create the necessary precursors to mapping the salt layer (Ezersky et al., 2006).

Gravity field transformations
The recent studies have shown that the major limitation of the microgravity method is the disturbing effect caused by  Khesin et al., 1996).Isolines of a model field (B) and the sum of its gradients (D): (1) positive, (3) zero, (4) negative; cylinder edge projection: (5) upper, (6) lower; (7) contour of the portion treated on the (B).
bodies underlying and surrounding the studied targets.It is especially important in zones of the salt layer existence.For instance, gravity effect from the salt layer distorts the observed gravity anomaly because salt density exceeds the surrounding medium density.We propose to apply different procedures of gravity field transformations allowing underlining the desired effects and suppressing the background noise.

Computing of informational characteristics
Sometimes, for removing the noise of different origin, presented as quasi-white noise, it is useful to calculate parameter called as "informational" one.Khesin (1974) suggested using the following informational parameter: and where J i is the informational parameter, P j is the relative frequency of the j -th interval of the i-th indicator on the histogram of its distribution, U i and U i are the amplitude and the error of this indicator determination, respectively.Expressions (1) and (2) were presented also in the book of Khesin et al. (1996).
Eppelbaum (Eppelbaum, 2001;Eppelbaum et al., 2003) proposed to apply the following expression: where U i is the geophysical observation (in our case -gravity) at ith point at the area under study and n is the total number of observations, K is the some coefficient.
An experimental testing indicates that application of expression (3) is more effective, especially in the complex geological media.

Computing of gravity fField derivatives
It is well known (e.g., Veselov and Sagitov, 1966;Blakely, 1995) that computing of vertical and horizontal derivatives of gravity field g is a powerful tool for delineation of some peculiarities of studied targets.The gravity field modeling was performed using Emigma 7.8 (PetRos EiKOn Inc) and GSFC (Khesin et al., 1996) software.
Let's consider the following example (Fig. 12).Two sinkholes of typical form occur closely at the depth of 7 m (upper edge).Upper part of these sinkholes has the density contrast (-1850 kg/m 3 ) and lower part -(-2000 kg/m 3 ).Figure 12A displays the results of model computation of gravity fields at two levels: 0.3 m (red line) and 1.5 m (blue line).Amplitude of these anomalies is (-0.22÷0.24mGal (220÷240 microGal)) and only very small positive anomaly (about 0.012 mGal) in the center of the negative anomaly testifies to the presence of some positive mass (between two sinkholes).Meanwhile for solving some applied geological problems it is high important to recognize, do we have one anomalous target or two closely occurred objects.
An absolute ratio between the main (negative) and small (positive) gravity anomalies (Fig. 12A) is about 20 (i.e., against the noise background this small positive anomaly could be not detected).At the same time the g z distribution (Fig. 12B) clearly indicates to the presence of some additional object (in our case -allows to delineate two different anomalies from two narrowly located objects).A ratio between the main and secondary anomalies in the graph of g z (Fig. 12B) is 3.4.Thus, probability of this anomaly detection in the graph g z increases in almost 6 times comparing with g curve.
The case of computation of gravity effects from the nearsurface closely locating karstic bodies (their upper edge is located at the depth of 0.8 m) is shown in Fig. 13.Upper part of these sinkholes has the density contrast (-1800 kg/m 3 ) and lower part -(-1900 kg/m 3 ).Figure 13A also presents the results of model computation of gravity fields at two levels: 0.3 m (red line) and 1.5 m (blue line).The corresponding ratio between the main and small gravity anomalies (Fig. 13A) is 4.52.The calculated ratio between the main and secondary anomalies of g z (Fig. 13B) consists of 1.54.Thus, probability of this anomaly detection for g z increases in 2.94 times comparing with g anomaly.
Figure 14 shows computation of gravity effects from the same geological section, but second level of "observation" is not 1.5, but 10 m (Fig. 14A) (here is proposed a future ROV exploration).Interestingly, that the changing of the second "observation" level for the value of 10 m does not allow to increase essentially the value of vertical derivative amplitude comparing with the field observed at level of 0.3 m.
The geological model presented in Figs.13C and 14C and the gravity effect from this section (computed for the level of 0.3 m over the earth's surface) was used for computing the values of g x and g xx (Fig. 15B and C, respectively).Analysis of the graph g x allows to determine projection of left boundary of the left occurring target and right boundary of the right occurring target.A positive anomaly arising between two extremums testify to a presence of some additional target (in our case it is a surrounding medium, i.e. we could made a suggestion about the geometrical parameters of these targets).Behaviour of the graph g xx clearly reflects location of vertical boundaries of two closely occurred objects with some small interval (surrounding medium) between these bodies.

Estimation of gravity effect due to salt layer occurring the coastal plain of the Dead Sea
Many investigators (e.g., Frumkin and Raz, 2001;Yechieli, 2000Yechieli, , 2002;;Yechieli et al., 2006) note the presence of salt layer in the western coastal plain of the Dead Sea.We compiled some averaged parameters of the salt body and surrounding medium (e.g., see Fig. 4) for construction of geological section presented in Fig. 16B.The gravity effects from this salt body (the main peculiarity is that this salt layer with density of 2150 kg/m 3 occurs in the surrounding medium with density of 1900 kg/m 3 and underlying surrounding layer has density of 2100 kg/m 3 ).Thus, this salt body produces a positive gravity effect.The values of the gravity effects were estimated depend upon a position of the right end (butt-end).The employed parameters along the strike were (-80 m, +80 m).The computations were performed for 5 variants of the right end location of this salt body (Fig. 16A): (1) infinite, (2) at 300 m (within the geological section presented in Fig. 16C), (3) at 350 m, (4) at 400 m and (5) at 700 m.It is clear that the curves (1) and (2) display a different behavior in the right part of the geological section.For a first view, the curves (1) and ( 5) are very similar (they have only different computation levels).However, calculation of difference curve ( g (infinity)g (700 m) ) shows (Fig. 16B) that we have no difference in the left and right parts of the profile, but in the center of the section the dif- ference anomaly reaches the values of more than -0.12 mGal (-120 microGal) over the left butt-end of the salt body.This fact once more indicates that only 3-D computing of gravity anomalies will allow solving the desired geological problems.

Application of Fisher-Lustikh expression
Fisher-Lustikh expression (Veselov and Sagitov, 1968) is half-forgotten, but effective procedure (especially when the contrast density σ is a priori known).Fisher-Lustikh expression is written as where H 1 and H 2 are the minimal and maximal distances to the earth's surface of some structure with contrast density, G is the universal gravity constant, g x(max) is the maximal value of the horizontal derivative of gravity field and σ is the density contrast.
Example of g x computing is presented in Fig. 17 (the gravity effect computed for the variant (3) in Fig. 16 (right end of the anomalous body is at 300 m) was selected as a basis curve).The first advantage of this computation is that the max value of g x testifies to the projection of the left end (butt-end) of the salt body to the earth surface (this parameter from the visual analysis of g graph might be not exactly determined).
Taking into account that g x =11.85• 10 −9 1/s 2 , and σ =250 kg/m 3 , we can calculated the value H 2 (depth to the lower edge of the salt body), H 2 ∼ =27 m.It should be noted that the real depth H 2 is 28 m (hence we have a difference of  1 m).Interestingly, that this difference is decreasing by analysis of curves with the right butt-end location at 350, 400 and 700 m and reducing to zero for the case of infinite location of the right end.This fact is explained that the expression (4) was developed for the case of infinite location of the anomalous body right end.

Downward continuation of microgravity data
Analytical continuation of gravity field in the lower semispace should underline anomalous effects from small geological objects.However, it is necessary to take into account the possible arising of disturbances of two kinds: (1) noise component in the observed gravity field may be transformed to significant fictitious anomalies, (3) possible erroneous realization of the downward continuation below the upper edge of anomalous body may initiate the effect of "field destruction" (Strakhov, 1976).Among the majority of downward continuation methods developed for potential geophysical fields, we will shortly consider two: (1) method based on the Gauss' theorem of the averaged value of harmonic function, and (3) Strakhov's (1976) method.
(1) Procedure of downward continuation of gravity field, based on the known Gauss theorem of the averaged value of harmonic function (e.g., Blakely, 1995): where S is the sphere with the radius R, point 0 is the sphere's center and M is the point located at sphere.
If we will select six points at the sphere with radius R=h, we can write the following equation for the downward continuation of gravity field at some levelh: g (0, 0, −h) = 6 g (0, 0, 0) − g (0, h, 0) − g (h, 0, 0) where g(0,0,h) is the observation at the level h over the earth's surface (we may receive this value, for instance, by direct measurements or by the way of performing upward continuation of gravity data).
The equality is more exacting with decreasing of the sphere radius (thus, this procedure is preferable for detailed gravity investigations).Example of such downward continuation to the levels of 0.7 and 1.0 m below the earth's surface is presented in Fig. 18.It is obvious that intensity of gravity field is increasing for 30 and 50%, respectively comparing with the surface (0.3 m over the earth's surface) computation.(3) Strakhov's scheme of downward continuation (2-D) consists of applying the following expression: where h is the depth of the downward continuation (and simultaneously length of horizontal intervals where g values are selecting).Preference of Strakhov's method is that employment of this procedure does not demand any additional data besides the gravity observations along a studied profile.However, it is a 2-D procedure with its known imperfections. Figure 19 shows comparison of Strakhov's and "averaged value" methods computed for more deeply occurring target (comparing with the previous example).We can see that the amplitudes of the computed curves (Fig. 19B) differ for a significant value.At the same time it necessary to underline that application of these procedures permits to perform only some qualitative estimation of karst structures and its classification.

Calculation of some quantitative parameters of karst cavity on the basis of inverse problem solution
Gravity field intensity F is expressed as where W is the gravity potential.
For anomalous magnetic field U a we can write (when magnetic susceptibility ≤0.1 SI unit) (Khesin et al., 1996): where V represents the magnetic potential.Let's consider analytical expressions of some typical models employed in magnetic and gravity fields (Table 1).
Here Z v is the vertical magnetic field component at vertical magnetization, I is the magnetization, b is the horizontal semi-thickness of TB, m is the elementary magnetic mass, z is the depth to a center of body (for HCC and sphere) and depth to the upper edge of TB and rod (point source), and M is the mass of sphere.
Taking into account all above mentioned, we can apply for the gravity field analysis the advanced interpreting methodologies developed in magnetic prospecting for complicated environments (Khesin et al., 1996).
We can also calculate the gravity moment of the object with contrast density (Khesin et al., 1996): where M g is the gravity moment, g a is the amplitude of gravity anomaly (in mGal) and h is depth of HCC occurrence (in meters).

Field Analytical expression
Magnetic Thin bed (TB) Point source (rod)

Gravity
Horizontal Circular Cylinder (HCC) Sphere It should be noted that by observations at an inclined relief we would receive some fictitious value of M g .To obtain the real value, we should apply the following expression where ω 0 is the terrain relief inclination angle along the observation profile (ω 0 >0 where the inclination is towards the positive the direction of the x-axis).
Figure 20 testifies application of improved versions of tangents and characteristic point methods for quantitative examination of gravity anomaly computed from a model of sinkhole.This disturbing body was approximated by a HCC model (center of the body was determined with required accuracy).Parameter The interpretation methods were tested on the results of microgravity survey performed at the Medford Cave site (Florida, USA).This survey was carried out with the aim to delineate the shallow subsurface air-filled cavities and tunnels (Butler, 1984).A density of 1900 kg/m 3 was used for the Bouguer and terrain corrections.One of profiles presented in Butler (1984) is shown in Fig. 21.The most significant anomaly I (almost 70 microGals) was interpreted using the abovementioned methodology.Results of interpretation have a good agreement with the geological data.The parameters M g was calculated as -0.238 [mGal•m].

Estimation of gravity effects from various karst cavity volumes by the use of set of computed curves
It is known that 3-D modeling in microgravity is a powerful interpretation tool (e.g., Branston and Styles, 2006;Debeglia et al., 2006;Abad et al., 2007).3-D gravity modeling over four models of sinkholes (with different size and depth of the sinkholes occurrence) is presented in Fig. 22. Upper part (major) of these sinkholes occurs in the surrounding medium of 2150 kg/m 3 and lower part (minor) -in the surrounding medium of 2250 kg/m 3 .The geometrical parameters along strike were selected in such a way that to be quasi-isometric to dimension shown in the geological section (for instance, for the model 1 parameters along the strike were selected as (-10 m, +10 m), and for the model 4 -(-3.5 m, +3.5 m).The computed gravity anomalies vary from -0.53 mGal (for the model 1) to -0.03 mGal (for the model 4).It allows to predict preliminary the expected gravity effects from desired buried sinkholes.
Another method of the karst volume estimation is computing of the same model by x and z coordinates, but with different parameters along the strike (y coordinate) (Fig. 23).The gravity effect from sinkhole with negative contrast density of (-1900 kg/m 3 ) was computed in eight variants: (1) -10 m, +10 m, (3) -20 m, +20 m, (4) -30 m, +30 m, (5) -40 m, +40 m, (6) -60 m, +60 m, (7) -80 m, +80 m, (8) -100 m, +100 m, and (9) -200 m, +200 m.For simplicity we selected here symmetric variant of the sinkholes butt-ends locations along the y-axis.An analysis of computed curves indicates that the gravity effects computed from models from (-10 m, +10 m) up to (-60 m, +60 m) have a clear distinction between themselves.Gravity effects between (-60 m, +60 m) and (-80 m, +80 m) has a small difference (about 15 microGal).Then, effects from (-80 m, +80 m) and (-100 m, +100 m) have very close values (discrepancy is about 4 microGal).Finally, difference between (-100 m, +100 m) and (-200 m, +200 m) consists of about 15 microGal.Thus, we can conclude that estimation of karst volume is much better realizing at small parameters of y than by large values of butt-ends along the strike.It is clear that gravity effect difference between the computed parameters [(-20 m, +20 m) -(-10 m, +10 m)] in more than one order exceeds the difference between the parameters [(-200 m, +200 m) -(-100 m, +100 m)].Thus, the microgravity estimation using a single profile analysis is more effective by small geometrical parameters along the strike.Obviously, if we should detect The final model (Fig. 24) demonstrates the karst volume estimation from the same model, but with asymmetric location of the butt-ends along the strike.Here we present 10 variants of computation.First variant is a symmetric model from the previous figure (-20 m, +20 m).Subsequent variants with a step of 5 m (the length along the strike has a constant value=40 m) are moving from the symmetric variant and beginning from 6th model (+5 m, +45 m) the models are outside the geological section.At the same time, the gravity effects of these bodies are not so small (for example, >200 microGal for 6th model).Thus, only 3-D advance gravity computing of all effects (even not available in the studied geological section) (Eppelbaum, 2006) will facilitate to development of real exact karst models.3.5 Some results of gravity field analysis and 3-D modeling in the nahal hever south area Some procedures described in Sect.3.4 were tested in the present investigation of Nahal Never South area.A map the informational parameter J i (computed by the use of Eq. 3) is shown in Fig. 25.Examination of this map indicates that location of J i negative anomalies have a good agreement with the contours of karst development for the period of 1999 (when the gravity survey was carried out).
It is necessary to underline that the quality of final gravity maps is strongly depending on the correct computation of the regional gravity gradient (Dead Sea Transform Zone (DSTZ) influence).Therefore, besides the application of polynomial and other approximations, some physical-geological approaches should be tested (e.g., computing of 3-D gravity effects from DSTZ model and subtracting this effect from the  observed field, removing the regional background by the way of downward continuation, etc.).
We applied the Emigma 7.8 software (PetRos EiKOn Inc) to create the 3-D physical-geological models of the sinkholes area development and to estimate the volume of the buried caverns.It should be noted that application of the GSFC program, developed specially for 3-D combined gravity- magnetic modeling in complicated environments (Khesin et al., 1996), gave the similar results.
The development of the physical-geological models has been carried out on the basis of the next main stages: 1. Generalization of all geological (geological mapping and drilling), geophysical (seismics, ERT and remote sensing images) and geochemical information for compiling the initial geological section, 2. Calculation of geological bodies density using the known relationships between the seismic velocities and density (e.g., Telford et al., 1995), 3. Application of the procedures of qualitative and quantitative analysis of gravity anomalies for obtaining certain parameters of disturbing bodies and introducing these parameters to initial physical-geological section,

Integration of microgravity with other geophysical methods
The complexity of current geological/geophysical problems in the Dead Sea region, and the ambiguity of the geophysical observations interpretation call for an integration of different geophysical methods and their integration with other methods.However, an extension of a set of methods is at variance with its economical efficiency and complicated from both organizational and technical viewpoints.Besides, there is a basic limitation imposed on the number of methods.
For estimation of integration preferences we can apply some elements of the theory of information (Khesin and Eppelbaum, 1997).If a set of methods is focused on investigat-   ing some independent indicators of equal value, the anomaly detection reliability γ can be described by an error function (probability integral) as: where υ is the ratio of the anomaly square to the noise dispersion for each i-th geophysical field.Now let's assume that the anomaly is indicated by three points and that the mean square of the anomaly for each field is equal to noise dispersion.For a single method, the reliability of the detection of an anomaly of a known form and intensity by Kotelnikov's criterion (Kharkevich, 1965) is expressed by erf √ ν i 2 .Then, the value of reliability for indi-vidual methods is 0.61 and those of a set of two and three methods is 0.77 and 0.87, respectively (according to Eq. 15).Risk of an erroneous solution (q value) is calculated as This means that the q value at the integration of two or three methods decreases by the factors of 1.7 and 3, respectively.A comparison of the risk with the expenditures enables one to find an optimum set of methods.
Besides the microgravity, GPR, seismic refraction and ERT, we can recommend also employment of high-precision magnetic survey (obviously, the first micromagnetic monitoring at the Dead Sea coast was performed by Rybakov et al., 2005) -some examples of advanced detailed magnetic field analysis in Israel are given by Eppelbaum et al. (2001Eppelbaum et al. ( , 2003Eppelbaum et al. ( , 2004)), and in separate cases -seismoelectric method (Neishtatd et al., 2006).
The last achievements in the field of geophysical surveying on the basis of remote operated vehicles (ROV) (e.g., Eppelbaum, 2008) allow suggesting a sharp increasing of combined examination of the ROV and land geophysical observations.

Verification of the microgravity study
Verification of any geophysical investigation is an important stage of the study.The large sinkhole cluster has formed at the Nahal Hever South microgravity anomaly site that allows to verify directly the results of the microgravity mapping and 3-D modeling (Fig. 27).We have performed geodetic measurements along the cluster contour in March 2007.Measurements were carried out using the non-differential Garmin V GPS system with accuracy of 3-4 m.Maps of the present contours of the sinkhole cluster are shown in Figs.28 and 29.
One can see in Fig. 28 that measured contour of the sinkhole cluster presented in Fig. 3 coincides with the residual gravity anomaly marked as 1.At the place of second anomaly (2 in Fig. 28) has formed a coupled sinkhole shown in Fig. 29.A seismic refraction survey carried out in 2007 did not reveal salt unit at the site of this anomaly (Ezersky et al., 2008).Our model suggests an existing of large preexisting dissolution cavern that refilled by the fines, which washed out by recharged waters and brought out to the preexisting cavern (Legchenko et al., 2008).However, it is out of the scope of this paper.

Conclusions
The pre-existing cavern was detected using microgravity mapping in the Nahal Hever South site where several sinkholes 1-2 m in diameter had been opened (Rybakov et al., 2001).The cavern was located under thin salt layer shaped like peak-cap detected by the seismic refraction method.Residual gravity anomaly of -(0.08÷0.14)mGal was revealed through area of approximately 40×70 m 2 .
Analysis of gravity field transformations (computing of horizontal and vertical derivatives, Fisher-Lustikh expression and downward continuation) indicates that the procedures are reliable tools for localization of buried karst caverns in com- We applied the Emigma 7.8 software (PetRos EiKOn Inc) to create the 3-D physical-geological models of the sinkholes development area and to estimate the volume of the buried caverns.The modeling was verified by application of the GSFC program, developed especially for 3-D combined gravity-magnetic modeling in complicated environments (Khesin et al., 1996).Volume of the karst domain was roughly estimated as 35 000 m 3 .
Results of the performed 3-D modeling along profile C'-C clearly show a potential of the microgravity method for the large buried cavern detection at the depths of 25-50 m.The microgravity study should be combined with other geophysical methods (first of all, seismic refraction and CVES) permitting specifying the geological structure of the site (salt layer borders and its surface topography, decompaction of the uppermost part of section, etc.).
It is shown the applicability of gravity measurements at different levels both by the land survey (0.3-1.5 m over the earth's surface) and remote operated vehicles observations for future investigations of the Dead Sea sinkhole problem.

Fig. 2 .
Fig. 2. Development of sinkholes in time in the Dead Sea region.

Fig. 4 .
Fig. 4. Presumable hydrogeological model of the sinkhole formation at the Dead Sea western shore (a) and geological section of the Nahal Hever South site derived from the boreholes HS-2 and HS-3 (b) (after Yechieli, 2002).

Fig. 6 .
Fig. 6.Singling out a ring structure: (A) model field of inclined circular cylinder calculated along a profile, (B) model field complicated by random noise, (C) apparent graticule for ring structure revealing, (D) singling out a model body by summing up horizontal gradients of field within the apparent circular graticule zones (afterKhesin et al., 1996).Isolines of a model field (B) and the sum of its gradients (D): (1) positive, (3) zero, (4) negative; cylinder edge projection: (5) upper, (6) lower; (7) contour of the portion treated on the (B).

Fig. 8 .
Fig. 8. Regional gravity map in the Bouguer reduction with topographic (bathymetric) map of the studied area (after ten Brink et al., 1993, with supplements).Reproduced by permission of the American Geophysical Union.

Fig. 9 .
Fig. 9. Scheme of the residual regional horizontal gradient of g B in the studied area.

Fig. 10 .
Fig. 10.Nahal Hever south site: Comparison of gravity and resistivity maps and salt contour (A -after Ezersky et al., 2006, B -after Rybakov et al., 2001).Figure 10B reproduced by permission of the Society of Exploration Geophysicists.

Fig. 13 .
Fig. 13.Computing of vertical derivatives of gravity field for two closely disposed models of sinkholes (near-surface occurrence).Levels of computation: 0.3 and 1.5 m. (A) Computed gravity curves, (B) Computed vertical derivative of gravity field, (C) Physical-geological model.

Fig. 14 .
Fig. 14.Computing of vertical derivatives of gravity field for two closely disposed models of sinkholes (near-surface occurrence).Levels of computation: 0.3 and 10 m. (A) Computed gravity curves, (B) Calculated vertical derivative of gravity field, (C) Physical-geological model.

Fig. 15 .
Fig. 15.Computing of horizontal derivative of gravity field for the models presented in the Fig. 13.(A) Computed gravity curve (level of computation: 0.3 m), (B) Calculated first horizontal derivative of gravity field g x , (C) Calculated second horizontal derivative g xx , (D) Physical-geological model.

Fig. 16 .
Fig. 16.Computing of gravity effects from the model of salt body by different positions of the right end of the body.(A) Computed gravity curves for different positions of the salt body right end, (B) Computed difference between g B(infinity) and g B(700 m) , (C) Physical-geological model.

Fig. 17 .
Fig. 17.Computing of horizontal derivative of gravity field for determination of the left end of this body.(A) Computed gravity curve, (B): Calculated horizontal derivative of gravity field, (C) Physical-geological model.

Fig. 18 .
Fig. 18.Downward continuation of gravity field using the theorem of the averaged value of harmonic function (3-D): (A) Gravity field observed at the level of 0.3 m over the earth surface, (B) Gravity field downward continued to levels of 0.7 and 1.0 m below earth surface, (C) Physical-geological model

Fig. 19 .
Fig. 19.Downward continuation of gravity field: (A) Gravity field computed at the level of 0.3 m over the earth surface, (B) Gravity field continued downward to level of 1.7 m below earth surface using the Strakhov's method (2-D) and averaged value of harmonic function (3-D), (C) Physical-geological model.

Fig. 20 .
Fig. 20.Inverse problem solution for the model of sinkhole performed using methods developed specially for complicated environments.

Fig. 21 .
Fig. 21.Quantitative interpretation of microgravity anomaly I (underground cavity was approximated by model of HCC) at the Medford Cave site (Florida, USA).Symbol ϒ marks the determined position of the HCC center.Observed curve and geological section are taken from Butler (1984), interpretation after Eppelbaum (2007b).

Fig. 22 .
Fig. 22. 3-D gravity anomalies computed from models of sinkholes at different stages of its development.

Fig. 23 .
Fig. 23.Computation of gravity anomalies due the karstic body with different symmetric locations of the body ends along the strike.

Fig. 24 .Fig. 25 .
Fig. 24.Computation of gravity anomalies due the karstic body with different asymmetric locations of the body ends along the strike.

Fig. 26 .
Fig. 26.Results of 3-D modeling of gravity field along the fragment of profile C'-C (location of this profile is shown in Fig. 10B).

Fig. 27 .
Fig. 27.Sinkhole cluster formed at the first sinkhole locations (the photography was done in March 2006).

Fig. 28 .
Fig. 28.Residual gravity anomaly map of 1999 with contour of cluster formed in the vicinity of this anomaly.

Fig. 29 .
Fig. 29.Coupled sinkhole developed in 2004 at the place of the anomaly 2 in Fig. 28.