## Abstract

One of the central problems of a collapse hazard evaluation is the answer to the question: at what point on a collapse-prone site will a future collapse sink be located? If we have information about underground cavities, which can be the reason for a collapse, this problem may be solved fairly easily; otherwise, prediction is more difficult. This paper describes an approach to define the most probable location of a collapse at the stage of preparation for this event. The situation when we assume the presence of underground cavities, but have no exact information about them, is considered. This is typical for karst regions or for piping events. The approach includes the mathematical analysis of data from cone penetration tests or dynamic probing tests, carried out during standard site investigation. The basic principle is the recognized correspondence between the spatial distribution of these test data and the stress distribution in soils above an underground cavity. The results of practical application of this approach are given for some examples in Russia. In particular, in the city of Dzerzhinsk (near Nizhniy Novgorod) the forecast was completely confirmed, when inside a hazardous zone delineated in 1994 a collapse sink of karst origin was formed in 1996.

The existence of weakened zones in a soil around a collapse sink is a well-known phenomenon, which is clearly shown by the results of cone penetration tests or dynamic probing tests (see Fig. 1). These zones have been described by various workers (e.g. Lebegue 1973). There are two tendencies in the spatial distribution of static or dynamic point resistance of soils near collapse sinks: if the distance between a test point and a collapse sink is small, the resistance is relatively small at a particular depth and the gradient of increase of resistance during cone penetration is also small.

It would be possible to explain these characteristics by unloading of a collapse sink's slopes, but the same phenomenon has been observed near very small collapse sinks. For example, in Figure 2, cross-section I–I, a very clear indication of this effect in unsaturated sands can be seen, although the collapse sink is too small for a marked slope unloading. In cross-section II–II, the spatial distribution of static point resistance is very similar; however, it is connected not with a collapse sink but with a subsurface cavity found during borehole drilling. The same result has been obtained by Rey & Suderlau (1984), who carried out dynamic probing tests in loess overlapping a cavity.

Thus, it is possible to conclude that this phenomenon is not explained by unloading in the sloping region of a collapse sink, but is the result of a change in stress distribution in soil overlapping a cavity during its enlargement. This change is the initial stage of a collapse. As has been shown by Anikeev *et al*. (1991), the appearance of an underground cavity is accompanied by a decrease of normal stresses in overlapping soil not only on the cavity's contour, but also near the surface, where it may be easily found by cone penetration tests or dynamic probing tests. It is possible to assume that the pattern of spatial stress distribution in the soil above a cavity forms long before a collapse; therefore the occurrence of a collapse sink does not cause radical changes in this pattern.

This hypothesis has been confirmed by the mathematical modelling of a collapse during the initial stages of its development (Khomenko & Kolomensky2000). The modelled picture (Fig. 3), in its upper part, is identical to the contour pattern of static or dynamic point resistance of soil near a collapse sink and above a cavity (see Figs 1 and 2). These results led the author to propose that by analysing the structure of a stress field, indirectly expressed through the spatial distribution of static or dynamic point resistance, it is possible to forecast a collapse location.

The field of static or dynamic point resistance is 3D, but for this study it is possible to assess horizontal sections, treating them as 2D fields. However, there are difficulties in this operation because the change of cone resistance with depth shows considerable fluctuation. This problem can be avoided by using an integrated parameter, the average cone resistance for a penetration point (*R**). The field of this parameter can be simply interpreted as a 2D surface; for example, as a plane with depressions.

## Principles

About 100 years ago Protodyakonov (1912) noted an interesting fact: the form of the horizontal section of a collapsed soil body is rounder in plan, the higher the position of the section. As a result of this, a newly formed collapse sink always looks like a circle or an oval in bird's eye view, irrespective of the configuration of the cavity triggering the collapse. This idea is well illustrated by Figure 4.

As the form of a collapsed soil body is determined by the configuration of the stress field in soil overlapping a cavity, it is logical to assume that in a homogeneous soil mass this 3D field can be approximated by a family of surfaces of revolution with a common vertical axis of symmetry. Obviously, the bottom of the axis passes through the orthocentre of the greatest circle that can be inscribed in the horizontal section of a cavity in contact with overlapping soil. The top of the axis passes through the centre of the future collapse sink.

Consequently, the 2D field of *R** above a cavity should become a symmetrical surface of revolution with a vertical axis. This circumstance further simplifies interpretation of the mathematical model. It is possible to estimate the configuration of this surface from common sense: at the centre (above a cavity) *R** should have the minimal value corresponding to maximal unloading of stresses and at the periphery its value should asymptotically approach a certain maximal value, which corresponds to the relict field of *R**. This exists before the occurrence of an underground cavity as an object unloading the overlapping soil.

Proceeding from these preconditions, it is possible to come to the conclusion that the field of *R** above a cavity should be similar to the surface formed by rotation around an ordinate axis of an upside-down Gaussian curve, which is expressed by the function $$mathtex$$\[R{\ast}{\,}{=}{\,}R{\ast}_{f}{\,}{-}{\,}{\Delta}R{\ast}_{c}exp[{\,}{-}{\,}({\rho}/{\rho}_{k})^{2}]\]$$mathtex$$(1) where *R** is the average of static or dynamic point resistance of soil (MPa), *R**_{f} is the relict maximal value of *R** (MPa), Δ*R**_{c} is the decrease in the *R**_{f} value above the centre of a cavity (MPa), ρ is the distance measured in plan between a penetration point and the centre of the cavity (m), and ρ_{k} is the empirical coefficient (m). The Δ*R**_{c} and ρ_{k} values are determined by the cavity's size, the thickness of the covering soils, and their properties.

The model expressed by equation (1) repeatedly proved to be true, in particular at sites in the Nizhniy Novgorod region. For example, the model was confirmed by the data from dynamic probing tests made near a collapse sink formed in fine alluvial Middle Quaternary unsaturated sands overlapping a cavity in Permian carbonate and sulphate rocks (Fig. 1). Figure 5a demonstrates very good correlation between *R** and ρ at this site using equation (1).

Thus, mathematical analysis is necessary in the search for cavities ready to collapse, and hence for the forecast of future collapse sink locations. The analysis of a 2D field of *R** allows anomalies to be found corresponding to the standard criterion expressed by equation (1). As this field is represented in the geometrical expression of the surface of certain configurations, the forecast requires a search for depressions on this surface, to determine their exact locations and to estimate the characteristics of these depressions from the point of view of their conformity to certain requirements, details of which will be given below.

It is possible to carry this out as follows. For each point of cone penetration tests or dynamic probing tests made during certain times on a certain site, the size of *R** must be calculated as an average. The depth from which averaging begins should be a constant for the site. This initial depth must be chosen to exclude a first layer below the surface of variable thickness, if this layer consists of soils with excessively high or excessively low static or dynamic point resistance (made ground, organic matter, weak material, etc.). If this is not done, the structure of the field of *R** can be very complicated. Then a contour map of *R** must be made as a field picture represented as a 3D surface. On this map using the ‘watershed’ principle, the boundaries of depressions on this surface should be delineated. The above-mentioned anomalies caused by cavities can correspond to separate sites of depression.

The location and identification of anomalies is carried out separately for different groups of points of cone penetration or dynamic probing tests. Each group is situated within the boundary of each delineated site of depression. Such groups should include no fewer than five test points, as discussed in the next section. For the mathematical processing of the test data for the points of any group, two approaches can be used. Each of the models on which they are based possesses certain limitations, but despite this the forecasts constructed with their help have proved to be successful.

## Approach

### Location and identification of the central part of the anomaly

The basis of this variation of the approach is the precondition that the form of the central part of the anomaly can be expressed by a circular paraboloid (Fig. 5b) or by a circular cone (Fig. 5c). For the realization of the forecast that is based on this precondition, the contour map of *R** can use a Cartesian coordinate system, which must be inserted optionally. Then, within one group for each test point a serial number is allocated. As a result, in each group, each point with number *i* will have three initial parameters: *R*** _{i}*,

*x*(abscissa), and

_{i}*y*(ordinate).

_{i}The statistical processing of the initial parameters consists in the determination of the polynomial factors of a stochastic function such as $$mathtex$$\[R{\ast}{\,}{=}{\,}K_{1}{\,}{+}{\,}K_{2}x{\,}{+}{\,}K_{3}y{\,}{+}{\,}K_{4}(x^{2}{\,}{+}{\,}y^{2})\]$$mathtex$$(2) where *K*_{1} is the first polynomial factor (MPa), *K*_{2} and *K*_{3} are the second and third polynomial factors (MN m^{−1}), *K*_{4} is the fourth polynomial factor (MN), *x* is the abscissa of the test point (m), and *y* is the ordinate of the test point (m).

Then, the coordinates of the centre of the anomaly are inserted into a Cartesian coordinate system calculated according to $$mathtex$$\[x_{c}{\,}{=}{\,}K_{2}/2K_{4}\]$$mathtex$$(3) $$mathtex$$\[y_{c}{\,}{=}{\,}K_{3}/2K_{4}\]$$mathtex$$(4) where *x*_{c} is the abscissa of the anomaly's centre (m) and *y*_{c} is the ordinate of the anomaly's centre (m).

Deviations of actual values of *R** from the axial-symmetric trend surface expressed by equation (2) cause the existence in a horizontal plane of a circular zone inside which the centre of a future collapse sink will be situated. The radius of this zone will be defined by a confidence interval and a trend in concentration. The coordinates of the zone's centre coincide with the coordinates of the centre of an anomaly on the trend's surface and so can be determined using equations (3) and (4).

For an estimation of the sizes of the circular zone inside which there should be a collapse sink, the distance on a horizontal plane from the zone's centre to each test point of any group (ρ* _{i}*) must be defined. For this numerical set the linear correlation between

*R** and ρ should be investigated. From the results of forecast verifications, this is clearly an empirical basis for the identification of an anomaly, and so this value can be accepted using a 2σ confidence interval and 0.6 as the limit of correlation coefficient. Thus, if the linear-correlation coefficient between

*R** and ρ is <0.6 it may be concluded that a collapse is not possible. If it was not confirmed, the standard error of the estimate of ρ must be calculated at its regression on

*R**. The radius of the circular zone inside which there will be a centre of a future collapse sink is equal to twice the size of the standard error.

It is impossible to know the precise location of the border of the central part of a depression on a trend surface of *R**, which is approximated by a circular paraboloid (and, especially, a circular cone). For this reason, variants of the forecast calculations should be carried out, for the various combinations of test points belonging to a single depression, if there are more than five points. All the variants of calculations must be considered equivalent, therefore the forecast zone of collapse hazard can be the point on a complex configuration representing a combination of several merged circles with various radii.

### Location and identification of whole anomaly

The second variant of realization of the approach is based on an assumption that the investigated depression can be described completely by the function expressed by equation (1), with *R**_{f} being equal to or close to the maximal value of *R*** _{i}* in a given group of test points. After the transformation of equation (1) and taking the natural logarithm, the result is $$mathtex$$\[{\,}(R{\ast}){^\prime}{\,}{=}{\,}A{\,}{+}{\,}B{\rho}^{2}\]$$mathtex$$(5)

$$mathtex$$\[{\,}(R{\ast}){^\prime}{\,}{=}{\,}ln(R{\ast}_{f}{\,}{-}{\,}R{\ast})\]$$mathtex$$(6)

$$mathtex$$\[A{\,}{=}{\,}ln({\Delta}R{\ast}_{c})\]$$mathtex$$(7)

$$mathtex$$\[B{\,}{=}{\,}{\,}{-}{\,}1/{\rho}_{k}.\]$$mathtex$$(8)

The further statistical processing of the initial data can be carried out in the same way as for the first variant of the approach, but with the following differences: (1) for each point of the group, (*R**)′ must be calculated to replace *R**; (2) linear correlation between (*R**)′ and ρ must be investigated, and should be accepted using a 1σ confidence interval and 0.8 as the correlation coefficient limit (the standard of identification); (3) the radius of the circular zone inside which there will be the centre of a future collapse sink is equal to the square root of the value of the standard error.

Because the value of (*R**)′ cannot be infinite, for the groups with five test points, the limiting condition is proposed: $$mathtex$$\[R{\ast}_{f}{\,}{=}{\,}(R{\ast}_{i})_{max}{\,}{+}{\,}0.1{\ }MPa\]$$mathtex$$(9) where (*R*** _{i}*)

_{max}is the maximal value of

*R** in a given group of test points (MPa).

A stricter, but also more difficult modification of the second method of the approach is possible if the field of *R** can be considered as a combination of the relict field, which existed before the occurrence of a cavity, and of local depression-like anomalies, which were imposed on it later. To a relict component of a field of *R**, expressed as a 2D surface, to which the test points located on ‘watersheds’ belong, are the points with greatest values of this parameter. The relict field can be described by a series of polynomials. This operation must be carried out step by step; during each step the area of ‘watersheds’ expands, and the area of depressions narrows.

This method can be carried out by gradually involving in statistical processing the increasing quantity of test points with maximal values of *R**. For the approximation of a relict field by polynomials the quantity of such points should be not fewer than four for a polynomial of first degree, seven for a polynomial of second degree, and eleven for a polynomial of third degree. Approximation by polynomials of higher degrees is inexpedient for many practical and theoretical reasons. From this point of view, the number of test points within an investigated site should be not fewer than nine; otherwise, the realization of the approach by the modification of its second method will not be possible.

The description of a relict field of *R** by polynomial expressions allows, in turn, the investigation of its structure by statistical analysis. On certain occasions, such research can represent a particular interest that is not connected with a collapse. Only trend surfaces with a high enough correlation of trend are representative, when the correlation ratio of the analytic trend expression *R**_{f} = f(*x*,*y*) is >0.8.

## Results

In recent years dynamic probing tests and especially cone penetration tests have often been used for the evaluation of collapse hazard, particularly in karst areas. Examples of such research include successful investigations of subsurface cavities (Rey & Suderlau 1984), of buried sinkholes (Chang & Basnett 1997), and of depressions on the surface of soluble rocks under overburden (Kaufmann & Quinif 2001). However, the approach presented in this paper not only represents post-formation investigation, but also provides a method of prediction.

There are many examples of the use of this approach in some regions of European Russia, mainly in karst areas (Table 1 and Fig. 6). Most of the investigated sites are situated in Moscow and Dzerzhinsk (Nizhniy Novgorod region), and the geological conditions of these cities are presented in Figures 7 and 8. In all the sites listed in Table 1 the origin of the underground cavity triggering a possible collapse is karst or piping. The mechanisms of the corresponding collapse-sink formation are illustrated by schemes given in Figures 9 and 10. The approach described here makes it possible to recognize stages 1 and 2 of this process, which are otherwise hard to identify.

Also, verifications of forecasts that used the approach were made. For instance, one verification had an inverse character (Khomenko 1992): the results of dynamic probing tests made in 1982 during the site investigation for construction of two houses in Moscow allowed a retrospective forecast of a hazardous zone in part of the site, and this was then compared with the location where in 1987 a collapse sink of karst origin appeared (Fig. 11a). Another example is an indirect verification carried out unexpectedly in the situation shown in Figure 2. As a result of the forecast made in December 1992 in the city of Dzerzhinsk, a hazardous zone was delineated (Fig. 11b). Borehole drilling executed in January 1993 inside this zone found at a depth of between 8.8 and 10.0 m an underground cavity in unsaturated sands. In the third case, also in Dzerzhinsk (Fig. 11c), verification showed a clear direct correlation of the centre of a collapse sink that occurred in December 1996 (shown in Fig. 12) with a predicted hazardous zone determined 2 years earlier.

## Conclusions

Attempts to forecast collapse using data from dynamic probing tests or cone penetration tests in soils above underground cavities (overburden) were undertaken, with a focus on karst terrains. The novelty of the proposed approach was demonstrated by the use of statistical trend models on the basis of certain principles of collapse development observed *in situ* and interpreted in terms of soil mechanics. Moreover, the studies on which this approach is based have significantly enhanced our understanding of collapse preparation with depth. The approach can be used successfully in karst areas and on piping-prone sites.

The approach proposed in this paper has proved to be effective and allows forecasts to be made. Its use does not demand any information other than that which can be obtained during standard engineering site investigation. The immediate prospects of the development of this approach are connected with a transition from clear spatial to spatial–temporal forecasting.

## Acknowledgements

The author would like to thank the enterprise ‘Antikarst and River Bank Protection’ (Dzerzhinsk, Nizhniy Novgorod region) for its co-operation and support.

- © 2008 The Geological Society of London