www.nonlin-processes-geophys.net/16/393/2009/ © Author(s) 2009. This work is distributed under the Creative Commons Attribution 3.0 License. Nonlinear Processes

Water level fluctuations in deep bore wells in the vicinity of seismically active Koyna region in western India provides an opportunity to understand the causative mech- anism underlying reservoir-triggered earthquakes. As the crustal porous rocks behave nonlinearly, their characteristics can be obtained by analysing water level fluctuations, which reflect an integrated response of the medium. A Fractal di- mension is one such measure of nonlinear characteristics of porous rock as observed in water level data from the Koyna region. It is inferred in our study that a low nonlinear dynam- ical system with three variables can predict the water level fluctuations in bore wells.


Introduction
Water table fluctuations in the upper crust represent an integrated response of the medium in a region, induced due to prevailing stresses.Thus, any variation in the water table before, during, or after an earthquake offers important clues in understanding the genesis of earthquakes.The piezometric head in a well contains signals of boundary stresses active around the zone of influence of the well, body forces active within this zone besides the hydrological properties of the affective porous media, such as porosity and permeability, physical and chemical properties of matrix and fluids, etc.The well influence zone, thus, needs to be treated as a complex multiphase and multi-component evolving media.It would be highly simplifying to assume such a media as linear Correspondence to: D. V. Ramana (dvr@ngri.res.in)poroelastic, a frequent assumption underlying many recent studies on interpretation of water levels in wells (Grecksch et al., 1999;Chadha et al., 2005).The water level fluctuations due to Spitak Earthquake (Gavrilenko et al., 2000), and the hydrological response due to earthquakes in central Japan (Matsumoto et al., 2003a, b) have been earlier studied.In this paper, we made an attempt to study the nature of the media in the Koyna region by applying non-linear analysis of the water level variations observed due to a local earthquake in the region.We performed a simple test to know whether the dimension of this water level time series is fractal and what is the nature of the underlying dynamical system.

Seismological and geological set up of the Koyna region
The Koyna region in western India is located within the interior of the Indian continental plate.The Koyna-Warna area (Fig. 1) is 225 km south of Mumbai in western India and lies 50 km east of the Arabian Sea.The Koyna and Warna dam reservoirs lie east of an elevated N-S escarpment parallel to the west coast of India.The mean elevation varies between 600 m on the escarpment to 100-200 m in the Konkan plains towards the west coast.Earthquakes in the region were first observed with the construction of Koyna dam and simultaneous impoundment during early 1960s (Gupta, 2002).Since then, earthquakes continue to occur in the region and show very good correlation with the annual cycle of loading and unloading of Koyna and Warna reservoirs (Pandey and Chadha, 2003), though a comprehensive model is yet to be developed.The case of triggered earthquakes in Koyna is unique as these have continued for the last 45 years, whereas other such cases reported worldwide do not show Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.such activity over such a long period.This would mean that, other than reservoir loading and unloading, there are other factors influencing the triggering phenomenon.These factors could be categorised in terms of deeper sources of inducing stresses or changes in the constitution of the media.This region lies in the Deccan Traps rocks of 67.4 Ma (Duncan and Pyle, 1988) comprising of several sequences of lava flows.
The thickness of the Deccan Traps varies from 1.5 km in the Koyna region to less than 100 m in the peninsular shield in central India (Kailasam et al., 1976).The massive, compact lava sequences have low permeability but there is significant migration of water through faults, fractures, columnar jointing or through vesicles.Several lineaments have been mapped from Landsat satellite imageries in the region; two NNE-SSW faults are delineated from Koyna in the north to Warna in the south.Cretaceous sediments with a basement complex of Proterozoic times underlie the Deccan trap rocks.
The constitutive relationship for such inhomogeneous and fractured materials would thus be quite complex.Water level fluctuations will have signatures of this complexity.
In our work, we analysed the water level variations in a bore well drilled in the region, as this indicator of tectonics will respond to porous structures of the region.Thus, one would expect its value to lie within 2 and 3.This has been borne out by the work reported in sequel.

Nature of the experimental set-up
Figure 1 shows twenty-one bore wells drilled in the Koyna-Warna region to study the variation of in-situ pore pressure in and around the seismic source volume (Gupta et al., 2000).The well at Ukalu (UKA) near the Warna reservoir has a depth of 150 m with its lower part left uncased and is in contact with the rock (pore fluids).The distance between the dam and the well is nearly 1 km down stream.Battery powered pressure transducers (OTT-Hydrometrie, 1992) were installed in these bore wells.These transducers read the hydraulic pressure of the water column above the sensor element with a resolution of ±1 mm at a sample rate of 15 min.As changes in rock stresses occur, these are partly transferred into pore pressure.Changes in pore pressure are transferred into changes in the depth of the ground water level as seen in the observation well.The UKA well data show very good response to earth tides.Signal-to-noise ratio of the largest tidal constituent M2 are above 50 for most recordings, with peakto-peak tidal amplitudes of up to 24 cm.Presence of tidal signals in UKA wells indicate that it is connected to a confined, that is, hydraulically locked and fully water-saturated aquifer.Accordingly, the well-aquifer system is capable of detecting weak pore pressure anomalies caused by crustal strain of the order of nanostrain (Gupta et al., 2000).

Time series analysis methodology
The correlation dimension describes the complexity of the system by providing the information about the dimension of the attractor.It may be related to the minimum number of variables describing the dynamic system (Sugihara and May, 1990).Several methods have been given in the literature to compute the correlation dimension of the time series data (Grassberger and Procaccia, 1983;Theiler, 1987).The most commonly used algorithm of Grassberger and Procaccia (1983) is used in our study to compute the correlation dimension of water table variations.The algorithm is based on the phase space reconstruction of the time series.According to the Takens embedding theorem (Takens, 1981), an m-dimensional phase space can be reconstructed from the single available time series x i of water table variations as, where τ and m are the time delay and embedding dimension, respectively.The above delay coordinate embedding provides a one-to-one image of the original system for sufficiently chosen m (Takens, 1981).For an m-dimensional phase space, the correlation dimension is obtained by the correlation integral C(r) denoting the fraction of the pairs of vectors with the distance smaller than a chosen radius r as, (2) where H (x) is the Heviside function defined as H (x)=1 for x>0 and H (x)=0 for x≤0, vectors X i and X j can be obtained using Eq. ( 1), ||.|| denotes the Euclidean distance between the two vectors, n is the number of data points and r is the radius of the sphere centered on X i .It has been shown that for sufficiently small values of r and in the limit of infinite amount of data, C(r) follows a power law as C(r)∼r D , where D is the correlation exponent.D is generally obtained from the slope of the graph of log C(r) against log r over a sufficient range of small inter distances and over the particular scaling region of r.The presence of chaos in the system can be inferred by plotting the correlation exponent D against the respective m.If it saturates to some finite value, the system generally considers being governed by deterministic dynamics.The saturation value of the correlation exponent is the correlation dimension of the attractor of the time series.If the saturation occurs at the small value of correlation exponent, then it is considered that the process generating the time series is governed by low-dimensional dynamics.The increasing correlation exponent without any saturation indicates that the process generating the time series is stochastic.
The computation of correlation exponent is influenced by the selection of the time delay τ .It should not be too large as the system may lose memory of its initial state.For small values of τ , the resulting phase space would be highly dependent and the new information about the evolution of the system cannot be produced (Frazer and Swinney, 1986).Hence, choice of τ should provide low correlation between adjacent points in the embedded phase space.The computation of mutual information between x i and x i+τ is suggested to estimate the reasonable τ (Frazer and Swinney, 1986), which can be given as, (3) where p(x i ) and p(x i+τ ) are the probabilities to find a time series value in the ith and i+τ th interval, respectively, and p(x i , x i+τ ) is the joint probability that an observation x i falls into the ith and x i+τ into the i+τ th interval.These probabilities can be obtained by plotting the histogram of the data.
After obtaining the mutual information function, the usual practice is to plot this function against varying τ .The reasonable value of τ can be selected by locating the time at which the first minimum in the mutual information function occurs.The details can be found in Frazer and Swinney (1986).

Results and discussions
Figure 2 shows the water level fluctuations during January-March 2005 at UKA well after removing the tidal effects.This is a complex time series which has been analysed by following the non-linear time series methodology.Before carrying out all the computations, the data is standardized by subtracting the average and dividing the output of that with the standard deviation.It is argued that the larger data size is required to compute the correlation exponents (Sivakumar, 2000).A small data set produces a smaller scaling region in calculating the slope of the straight line fitted to logC(r) against log(r), which leads to underestimation of the correlation exponent.Hence, it is aimed at keeping the data size larger (∼8640).The computation of correlation exponent is also affected by other factors such as the presence of noise.
The noise in the data affects the scaling behaviour.The other effects of noise on the computation of correlation dimension are discussed elsewhere (Sivakumar, 2000).The real time series is always expected to be contaminated by some level of noise.It is, however, shown by Sivakumar et al. (1999c) that the presence of noise significantly affects the prediction accuracy and the correlation dimension estimates are not influenced even for high noise levels.The authors argue that the correlation dimension estimation could be attempted for preliminary investigation of character of chaos in the hydrological data.To compute the correlation dimension, the delay time was computed using the mutual information function for different lags (Fig. 3) and the lag time, where the first minimum in the mutual information function occurs, can be considered as the choice of optimal τ .It is observed that the first minimum occurred at lag 24.The correlation integral is computed against different r and m values, as discussed earlier.Fig. 4a shows the relationship between some values of log(r) and logC(r) for various m ranging from 1 to 10.The slopes are computed and then plotted in Fig. 4b as correlation exponent against various m.It is observed that correlation exponent increases with m and at a certain m (∼6) it gets saturated.The saturation value of correlation exponent, called correlation dimension, is approximately 2.7.To describe the dynamics of the water table variations, the minimum number of required variables is 3.The number of variables sufficient to embed the phase space of the attractor, in case of water table variations, is 6.The low and non-integer dimension shows the presence of low dimensional chaos in water table variations.Costain and Bollinger (1996) have inferred fractal dimension of water levels at a depth by downward continuation of stream flow data, using pore pressure diffusion equation with an assumed value of a diffusion constant.Their estimates of fractal dimension range between 1.2-2.4.The value reported here is larger than these values.This difference can happen when, in extrapolation by diffusion equation it is possible that high frequency behaviour of stream flow would be attenuated.The value of fractal dimension being 2.7 implies that the underlying dynamics of the water level fluctuation is non-linear.These values would give clues for construction of non-linear dynamical systems to explain and predict water table fluctuations.Thus we need a low dimensional nonlinear dynamical system with three variables to describe and predict the water level fluctuations.

Conclusions
The water levels in the wells in porous crustal rocks of Koyna region respond to changes in the prevailing stress/strain regime.The present analysis of the water level data showed that the fractal dimension of the underlying nonlinear system in 2.7 indicating that a minimum number of variables describing the system would be 3.The tectonic system is shown to behave as a low dimensional non-linear dynamical system to construct a prediction model for water level variation in bore wells in this region.

Fig. 1 .
Fig. 1.Map showing Koyna-Warna region in western India.The 21 bore wells drilled in the region are shown in bold circles.Small dots are the epicenters of the earthquakes recorded for one year in the seismically active Koyna region.

Fig. 2 .
Fig. 2. Water level fluctuations with time (during January-March 2005 at Ukalu after removing the tidal effects).

Fig. 3 .
Fig. 3. Mutual information function for water level fluctuations time series.

Fig. 4 .
Fig. 4. (a) Correlation integral C vs. r for the different values of m.(b) Correlation dimension of the water level fluctuation time series.