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.
Results of dynamic probing tests near a collapse sink of karst origin that appeared in 1979 west of Dzerzhinsk, Nizhniy Novgorod region. 1, Initial (a) and smoothed (b) curves of dynamic point resistance (p) v. depth (z). 2, Fine alluvial Middle Quaternary unsaturated sands. 3, Contour of p (MPa). 4, Point of dynamic probing test made in 1979 (with its number), shown in plan view (a) and cross-section (b). 5, Cross-section line. 6, Collapse sink.
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.
Results of cone penetration tests made at a site in the industrial zone of Dzerzhinsk, Nizhniy Novgorod region. 1, Initial (a) and smoothed (b) curves of static point resistance (qc) v. depth (z). 2, Contour of qc (MPa). 3, Point of cone penetration test carried out in 1992 (with its number), shown in plan view (a) and cross-section (b). 4, Collapse sink of karst origin, which appeared in 1992 (a) and orthogonal projection of its axis on a cross-section (b). 5, Borehole which was drilled in 1993 and found a cavity in unsaturated sands. 6, Predicted configuration of the cavity in plan view (a) and cross-section (b). The others symbols are as in Figure 1.
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.
Results of computer modelling of soil unloading above a cavity triggered by its appearance: the spatial distribution of soil density in a vertical section at different times (t).
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.
A newly formed collapse sink of mining origin triggered by an abandoned underground gallery near Podolsk, Moscow region (photograph taken in 1984): it is possible to see the change with depth of the form of the sink's horizontal section, which is rectangular at the bottom and circular at the top.
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).
Quadratic exponential (a), quadratic binomial (b), and linear (c) correlation of average of dynamic point resistance (R*) and distance between a test point and the centre of a collapse sink (ρ), which appeared in 1979 near Dzerzhinsk (Fig. 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, xi (abscissa), and yi (ordinate).
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 K1 is the first polynomial factor (MPa), K2 and K3 are the second and third polynomial factors (MN m−1), K4 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 xc is the abscissa of the anomaly's centre (m) and yc 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.
Examples of use of the approach in European Russia
Location of sites investigated by the approach. 1, Single site and its number (see Table 1). 2, Group of sites in a city. 3, Area underlain by limestone and dolomite. 4, Area underlain by limestone, dolomite, and gypsum. 5, Area underlain by chalk and marl.
Schematic cross-section through the central part of Moscow. Adapted from Anon (1997). 1, Carboniferous limestone, dolomite, and clay. 2, Jurassic sand. 3, Jurassic clay. 4, Cretaceous sand. 5, Quaternary loam. 6, Quaternary sand. 7, Groundwater level.
Schematic cross-section through the western part of Dzerzhinsk. Adapted from Tolmachyov et al. (1986). 1, Lower Permian gypsum. 2, Upper Permian weatheredlimestone. 3, Upper Permian clay. 4, Quaternary sand. 5, Groundwater level.
Mechanism of typical collapse-sink formation of karst origin in investigated conditions: stage 1, cavity formation; stage 2, subsidence formation; stage 3, collapse-sink formation. 1, Soluble rock. 2, Clay. 3, Sand. 4, Level of unconfined groundwater. 5, Level of confined groundwater.
Mechanism of typical collapse-sink formation of karst or piping origin in investigated conditions: stage 1, cavity formation; stage 2, subsurface collapse; stage 3, collapse-sink formation. 1, Soluble rock or permeable saturated soil. 2, Impermeable or unsaturated permeable soil.
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.
Some examples of forecasts constructed by the approach: (a) retrospective forecast verification in Moscow; (b) forecast confirmed by discovery of cavity in Dzerzhinsk; (c) forecast confirmed by collapse formation in Dzerzhinsk. 1, Point of dynamic probing test (a) and of cone penetration test (b) with their numbers. 2, Hazardous zones of the collapse, where the centre of a future collapse sink will be situated, which are forecast by identification of the central part of the anomaly in the field of R* (a) and by identification of the whole of this anomaly (b). 3, Collapse sink with year of its formation. 4, Old sinkhole. 5, Borehole that has located an underground cavity (Fig. 2). 6, Buildings and structures.
Collapse sink of karst origin, which appeared in Dzerzhinsk in 1996 within the hazardous zone forecast in 1994 by the proposed approach (Fig. 11c).
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
