Nonlinear analysis of drainage systems to examine surface deformation : an example from Potwar Plateau ( Northern Pakistan )

We devise a procedure in order to characterize the relative vulnerability of the Earth’s surface to tectonic deformation using the geometrical characteristics of drainage systems. The present study focuses on the nonlinear analysis of drainage networks extracted from Digital Elevation Models in order to localize areas strongly influenced by tectonics. We test this approach on the Potwar Plateau in northern Pakistan. This area is regularly affected by damaging earthquakes. Conventional studies cannot pinpoint the zones at risk, as the whole region is characterized by a sparse and diffuse seismicity. Our approach is based on the fact that rivers tend to linearize under tectonic forcing. Thus, the low fractal dimensions of the Swan, Indus and Jehlum Rivers are attributed to neotectonic activity. A detailed textural analysis is carried out to investigate the linearization, heterogeneity and connectivity of the drainage patterns. These textural aspects are quantified using the fractal dimension, as well as lacunarity and succolarity analysis. These three methods are complimentary in nature, i.e. objects with similar fractal dimensions can be distinguished further with lacunarity and/or succolarity analysis. We generate maps of fractal dimensions, lacunarity and succolarity values using a sliding window of 2.5 arc minutes by 2.5 arc minutes (2 .5×2.5). These maps are then interpreted in terms of land surface vulnerability to tectonics. This approach allowed us to localize several zones where the drainage system is highly structurally controlled on the Potwar Plateau. The region located between Muree and Muzaffarabad is found to be prone to destructive events whereas the area westward from the Indus seems relatively unaffected. We conclude that a nonlinear analysis of the drainage system is an efficient additional tool to locate areas likely to be affected by massive destructing events affecting the Earth’s surface and therefore threaten human activities. Correspondence to: F. Shahzad (geoquaidian@gmail.com)


Introduction
Quantification of spatial patterns and their rigorous analysis using nonlinear analyses is increasing in geology, information sciences, forestry, life sciences and peripheric disciplines (Dombradi et al., 2007;Gloaguen et al., 2007;Dong, 2009;Melo et al., 2006;Marwan et al., 2007;Dougherty and Henebry, 2001;Dong, 2000a;Martínez et al., 2007;Feagina et al., 2007;Turcotte, 1992;Watterson, 1986).The main goal of these approaches is to identify heterogeneities in natural patterns and study their causes.In tectonic geomorphological analysis, the heterogeneities take the form of breaks, linearity and asymmetries in drainage systems that are due to the ongoing surface and subsurface geological processes.These processes are episodic in nature and the drainage patterns can preserve their cumulative effects (Shahzad et al., 2009;Bull, 2007).
Previous studies have proposed that the linear (Wobus et al., 2006;Shahzad et al., 2009) and nonlinear analysis (Dombradi et al., 2007;Gloaguen et al., 2007;Guillermo et al., 2004;Shahzad and Gloaguen, 2010a) of individual streams as well as the whole drainage network may be used to investigate surface deformation.The linear techniques mainly focus on the secondary parameters e.g.contributing area, channel lengths, elevation etc. but ignore the fractal nature of drainage systems.The usual linear analysis e.g.stream profile analysis (Wobus et al., 2006;Shahzad et al., 2009) uses the slope-area relationship which can give similar results for different contributing effects.
The non-uniqueness of usual approaches can be reduced if we take the spatial distribution of the drainage network into account.The nonlinear analysis of the drainage system is very important because patterns with different space filling properties can discriminate areas which would have had similar signatures using linear or similar analyses.We propose that the way the drainage system fills the space is a strong indicator of the area vulnerability to surface deformation.We have shown that drainage systems under tectonic forcing tend to be reorganized and linearized (Gloaguen et al., 2008).This reorganization has a strong effect on the space filling of the F. Shahzad et al.: Nonlinear analysis of drainage systems drainage system.Therefore we propose to apply a nonlinear analysis including three fractal measures i.e. fractal dimension, lacunarity analysis and succolarity analysis in order to characterize the disequilibrium of drainage systems and to quantify the changes from a dendritic system into a structurally controlled one.These measures allow the examination of textural properties of the drainage system and thus should be a strong support for the evaluation of surface deformation localization and intensity.
In this paper, we show that the analysis of the river network complexity allows the localization of zones prone to tectonic induced damages even in areas where seismic deformation occurs.In other words, by examining the vulnerability to surface deformation using the spatial variability of drainage patterns, we may get insights how surface processes operate.Fractal measures were applied to the drainage network of the Potwar Plateau and its surrounding areas (Fig. 1).The drainage pattern in the study area is the result of spatially variable tectonic and erosional processes (Shahzad et al., 2009) and hence represent areas with variable vulnerability to surface deformation (Burbank and Beck, 1991).The disconnected, linearized and anomalous drainage patterns in a low seismic context give motivations for this study.We first extracted the drainage network from Shuttle Radar Topographic Mission's digital elevation data (SRTM DEM -90-m resolution) using D8 algorithm (O'Callaghan and Mark, 1984).The Box Counting method (Foroutan-pour et al., 1999;Dombradi et al., 2007;Guillermo et al., 2004) was applied to calculate the fractal dimension of the Indus, Swan and Jehlum rivers in order to identify their anomalous drainage behavior.The fractal dimension (D) distribution map was prepared using the extracted drainage network.
We selected eight regions with similar D values in order to evaluate their behaviors using the lacunarity and the succolarity approaches.This shows that the distribution of translation invariance in these regions is somehow similar.Nonetheless, drainage locations with similar D values represent very different drainage textures, which should be differentiated further.The lacunarity analysis (∧) was performed using a gliding box method (Melo et al., 2006;Zhang et al., 2007;Dombradi et al., 2007) that allowed us to prepare lacunarity maps.The succolarity analysis (σ ) was used to measure the rotation symmetry of the drainage patterns.It is an attempt to investigate the rotation of the area with similar translation behavior.Physically, it actually measures the percolation or draining ability of a textural system (Melo and Conci, 2008) in the horizontal (left to right and right to left) and vertical directions (top to bottom and bottom to top).Like for fractal dimension and lacunarity analysis, we prepared similar spatial distribution maps.An average succolarity value was used to characterize the relative heterogeneity of four selected sites where the drainage pattern shows similar translation invariance and cannot be distinguished by lacunarity and fractal analysis.
Three parameters of drainage patterns are associated with the tectonic deformation i.e. distribution, translation invariance and orientation.The distribution of drainage patterns is studied using fractal analysis and a low fractal dimension value i.e.D=1 is expected in the region of intense deformation.The lacunarity analysis is used to study the translation invariance in the drainage patterns, and in a relative analysis, a low value of lacunarity represents a less deformed region.The third and last parameter i.e. orientation is very important, as the drainage patterns are also affected by the orientation of the tectonic structures.In a relative analysis, this orientation can be characterized using succolarity analysis, where high mean succolarity values represent regions of intense deformation.

Study area
Analysis and mapping of surface deformation on the Potwar Plateau is of significant importance both geologically and geographically.The study area gathers some of the major settlements in Pakistan including the capital city Islamabad (Fig. 1).During the recent earthquake of 8 October 2005, major deformations were observed in different parts of this region.It therefore becomes highly important to investigate the vulnerability to surface deformation and provide reconnaissance observations about the seismic hazard.The plateau is highest to the north, where altitudes reach up to 2688 m a.s.l.The Swan River bisects and drains most of the water of the Potwar Plateau.Other major rivers draining the study area are the Indus and Jehlum rivers on the western and the eastern sides respectively.The relief of this area is rugged (Shahzad et al., 2009).
Tectonically, the Potwar Plateau is part of the Himalayan fold and thrust belt, a prominent feature along the India-Eurasia collision zone.From north to south, this collision zone consists of the Main Karakoram Thrust (MKT), the Main Mantle Thrust (MMT), the Main Central Thrust (MCT), the Main Boundary Thrust (MBT) and the Salt Range Thrust (SRT) as major tectonic boundaries.As a result of strain locking the deformation started to offset towards the Main Central Thrust (MCT) and then to the Main Boundary Thrust (MBT).The north-western part of this mountain belt is referred to as north-western Himalayan fold and thrust belt (NWHFTB) (Shahzad et al., 2009;Kamp et al., 2008;Jadoon and Frisch, 1997;Kazmi and Jan, 1997;Monalisa et al., 2007).The Potwar Plateau extends from the MBT in the north to the SRT in the south, while the Jehlum and Kalabagh Faults mark the eastern and western boundaries (Fig. 1).The study area is characterized by variable deformation patterns, as the Swan Syncline divides this region into two parts i.e. the Northern Potwar Deformed Zone (NPDZ) and the Southern Potwar Platform Zone (SPPZ) (Shahzad et al., 2009;Jadoon and Frisch, 1997).(Kamp et al., 2008).(b) Tectonic setting of the India Eurasia collision zone (Kamp et al., 2008).(c) Seismotectonic map of the study area (Kazmi and Jan, 1997;Kamp et al., 2008;Shahzad et al., 2009;Monalisa et al., 2007).

Methodology
The first step of this study is to extract the drainage network image (binary).Digital Elevation Models from the Shuttle Radar Topographic Mission (SRTM DEM -3 -resolution) were used in this study.Initially, the DEM is processed in order to remove pits/holes and then the drainage network was extracted using a D8 flow grid algorithm (O'Callaghan and Mark, 1984).This algorithm calculates possible flow directions at each cell towards the 8 neighboring cells.The least cost algorithm is used to connect the stream flow directions in order to generate a continuous network of streams.TecDEM, a MATLAB based software was used to prepare the drainage networks (Shahzad and Gloaguen, 2010a).The extracted drainage network is then filtered to represent the streams with a contributing area greater than 1 km 2 in order to keep only the parts where effective flow occurs and not hill slope diffusion (Tarboton et al., 1991).We use the streams of Strahler order higher than or equal to 3 in order to reduce the effects of DEM noise.The drainage network is finally converted into a binary image where the drainage network has a pixel value of 1 (Melo et al., 2006;Plotnick et al., 1993).
The estimation of fractal dimension is a widely used approach.In the present case, we expect to detect linearized streams controlled by active tectonics.Unfortunately the fractal dimension approach is equivocal, as it only gives an indication of the intensity of the network complexity, but not a pattern description.To overcome the non-uniqueness of the approach we perform a lacunarity and succolarity analysis (Mandelbrot, 1983).These methods are explained in the following sections.The simplified flow chart of the processes is shown in Fig. 2.

Fractal analysis
The use of fractal analysis for the quantitative description of geomorphological features has increased in the last decades (Gloaguen et al., 2007(Gloaguen et al., , 2008;;Mandelbrot, 1983;Dombradi et al., 2007;Ben-Zion and Sammis, 2003).The drainage patterns are characterized by irregularities with self-similar characteristics.In this investigation, the fractal dimension for the selected rivers (Guillermo et al., 2004) and whole drainage network (Dombradi et al., 2007) is calculated using a box counting method.The box counting method uses an approach of placing boxes of variable sizes on an image and then counting the number of drainage pixels within the box (Foroutan-pour et al., 1999;Dombradi et al., 2007;Guillermo et al., 2004).Box sizes s in each grid and respective number of boxes N (s) are counted.The fractal dimension value is calculated by using the following equation.
where s is the length of the box side and N (s) is the number of boxes.The slope of the bes fitted line for the double logarithmic plot of N (s) and 1/s is equal to D. This method is used to examine the D value of three major rivers in this region.
In order to investigate the deformational processes, a grid based approach was used to investigate the spatial distribution of fractal dimensions.A moving window/grid of 2.5 arc minutes by 2.5 arc minutes (2.5 ×2.5 ) on the binary image of the drainage network was used to calculate D values for each window using Eq.(1).A spatially distributed D value map was prepared with color coding so that the regions with smaller D values can be identified easily.
The fractal dimension distribution map represents an overall picture of the linearity of the drainage pattern.The low fractal dimension i.e. when D→1, corresponds to highly linearized drainage patterns (Gloaguen et al., 2008) and thus corresponds to the high vulnerability to surface deformation.Similarly, the regions with high fractal dimension i.e. when D→2, represents the dendritic drainage behavior, less affected by surface deformation.In a relative analysis of the spatial distribution of fractal dimensions, D can be classified into three categories.The values of D less than 1.3 are associated with regions of high vulnerability to surface deformation, the values between 1.3 to 1.4 represent medium vulnerability and values higher than 1.4 are classified as almost invariant regions.

Lacunarity analysis
The lacunarity analysis is used to improve the textural description of the river patterns by studying not only the distribution and sizes of empty sets (Dong, 2000b) but also informs about the deviation of fractal objects from translation invariance (Plotnick et al., 1996).This approach is very useful in distinguishing between the different textural patterns with similar fractal dimension.The lacunarity analysis has been in use for image processing and analysis applications in geology, forestry, medicine, biology, and related fields (Melo et al., 2006;Plotnick et al., 1993;Zhang et al., 2007;Marwan et al., 2007;Dougherty and Henebry, 2001).There are plenty of methods (Melo et al., 2006) to calculate lacunarity of binary images e.g."Gliding Box method, Differential Box Counting", "three-term local quadrant variance" (3TLQV).For its simplicity, we used the Gliding Box method where a square box of side r is glided along all possible direction of the texture.The mass s , i.e., the total number of drainage pixels counted during this process, is calculated.This procedure is repeated with sliding boxes of gradually increasing size i.e. r+i.Plotnick et al. (1993) suggested that the gliding box should be of size r=1 to some fraction of image length (M).As a result, a frequency distribution of mass s with variable box size r is obtained.This frequency distribution is then converted into a probability distribution P (s,r) by normalizing with the number of total boxes N(r) of size r.The dimensionless lacunarity (r) is calculated using the first and second moments of this distribution, as shown in the following equations (Plotnick et al., 1996(Plotnick et al., , 1993) The lacunarity curves are analyzed by using log-log plots between box sizes (r) and lacunarity values ( (r)).
A lacunarity distribution map with a moving window of 2.5 arc minutes by 2.5 arc minutes (2.5 ×2.5 ) is prepared on the binary image of the drainage network.In each moving window, the underlying image of 2.5 ×2.5 (i.e.51×51 pixels for 90 m DEM) was taken as sub image and the box size of r=1 to 20 (less than half of the moving window length as suggested by Plotnick et al. (1993)) was used to measure the lacunarity values.This process provides us with a distribution of box sizes vs. lacunarity values.Now, instead of the whole distribution, we need a single lacunarity value which can be used as representation of the underlying area.This single value of lacunarity should be carefully selected, because, as it can observed on the log-log plots (Fig. 4d), it varies with the resolution of the data set, i.e. high values of box sizes (e.g.r=20, 25) will give low lacunarity and low values of r (e.g. 1, 2, 3) will give high values of lacunarity (Plotnick et al., 1993).In order to have less mean error, we decided to select a central value of box size at log(r)=1 and a corresponding value of lacunarity was obtained.This value can vary from case to case depending upon data resolution and moving window size.We obtained a single lacunarity value at each pixel position.After repeating this method for all moving windows, a spatial distribution of lacunarity values was obtain.This lacunarity distribution map of the study area is shown in Fig. 4b.This method was also repeated for all the regions of interest where similar fractal dimension values have been observed in Fig. 4a.A combined plot of the lacunarities of all the 8 identified regions was prepared.This plot was used to identify the relative vulnerability of the sites using high and low lacunarity values.
The drainage textures can display high deviations from the translation invariance of the systems and therefore high lacunarity values.This high lacunarity is attributed to the heterogeneous drainage patterns i.e. highly deformed regions.The lowest value of lacunarity represents the area with low neotectonic activity in the selected regions.In a relative sense, low values of lacunarity represent the stable regions or areas with low neotectonic activity.The spatial distribution map of the lacunarity values can be divided into three categories i.e. low, medium and high lacunarity values.The values of (r) less than 1.5 represent homogeneous drainage patterns and are associated with regions of low vulnerability of surface deformation.The (r) values between 1.5 to 2.5 represent medium vulnerability and values higher than 2.5 represent heterogeneous drainage patterns and are classified as highly vulnerable regions.

Succolarity analysis
Succolarity analysis is a complimentary approach to fractal and lacunarity analysis.It measures the percolation or draining capacity of the image, i.e., how much a given fluid can flow though this image (Mandelbrot, 1983).Melo and Conci (2008) suggested that this type of analysis is very useful when the input texture has direction and/or flow information www.nonlin-processes-geophys.net/17/137/2010/ Nonlin.Processes Geophys., 17, 137-147, 2010 associated with it e.g.drainage patterns.In relation to surface deformation, succolarity analysis measure the rotational effect of the geological structures.Succolarating fractals include the filaments that would have allowed percolation, i.e., the amount of interconnected pixels in drainage textures.These textures consist of two types of pixels, i.e. empty spaces and impenetrable mass i.e. drainage.The first approach to calculate this fractal measure was provided by Melo and Conci (2008) using a box counting approach.He suggested that it can be measured along any possible flow direction (0 • ∼360 • ).In order to measure the mean succolarity of the pattern, we selected four possible directions, i.e. along 0 • , 90 • , 180 • , and 270 • .First, the rotated image is prepared at the desired angle and the succolarity is calculated using the following procedure.
1.The image is flooded by checking all its top pixels.If the pixel represents any empty space on the image, it means that the fluid can flood this pixel and its neighbors (4 neighbors for each pixels: top; bottom; left and right).This process is applied recursively on all the pixels until impenetrable mass is encountered.The flooded image provides information about the length of filaments in the image along that flow direction.
2. Each flooded image is analyzed using a box counting approach similar to the box counting method for fractal dimension estimation.This method uses an approach of placing boxes of variable size k (k=1 to n, where n is the number of possible integer divisions Melo and Conci, 2008) on the flooded image and then counting the number of flooding pixels (NP(k)) within the box.
3. A pressure matrix (PR) is created which is equal to the size of the flooded image and consists of linearly increasing weight from top to bottom.The contents of this matrix depends upon the box size (k) because it is applied on the centroid of the box and indicates the amount of pressure over it.The PR at any location pc for box size of k i.e.PR(pc,k) represents the centroid of the box on the scale of pressure.
4. Finally, the succolarity for a given direction (dir) is calculated using the following equation.
The dir represents the required flow/percolation direction of interest.This equation provides the succolarity values such that it is a dimensionless value like D. The results of σ can be analyzed in a variety of ways, i.e., along each direction or using the mean value.If we have large amount of data for comparison, then we can use the mean value to distinguish between their draining properties.The average value can be calculated using the arithmetic average of succolarity along available directions.
where σ represents mean succolarity and is calculated by taking the mean of the succolarities along the four directions (Melo and Conci, 2008).The value of this measure ranges from 0 to 1 and represents the percentage of percolation.We prepared the succolarity map of the area using the minimum of mean succolarity ( σmin ) with a moving window of 2.5 arc minutes by 2.5 arc minutes (2.5 ×2.5 ).The succolarity plot of the selected regions was also prepared using the mean succolarity value.Again the succolarity values were divided into three categories i.e. low (less than 0.5, medium (0.5 to 0.75) and high succolarity (greater than 0.75).The highly succolarating textures represent the farthest location of impenetrable mass over a filament i.e. the streams are at distant location.In the neotectonics scenario, the high values of succolarity represent an increased vulnerability to surface deformation because of the high percentage of empty space and the existence of long filaments.These filaments are probably caused by the tectonic lineaments present in the tectonically active regions.The regions where the succolarity values are low, represent the places where the geological conditions did not allow the growth of lineaments and thus they consist of high drainage density.

Results and discussion
Assessing drainage systems in the context of surface deformation is getting very important because they represent the witnessing of erosive and tectonic processes which may have disconnected, linearized and changed the dendritic behavior of the drainage system.The morphology of the drainage network of the Potwar Plateau is clearly influenced by geological forces at different stages (Fig. 1).The study area features mainly two types of spatially distributed drainage patterns, parallel and dendritic.The transition from one type to another is characterized by climatic and geological processes.The parallel drainage patterns characterize the flow along steepened regions and are mainly controlled by local geological settings.The steepened regions are produced as a result of local uplift, thus identifying the linearized drainage network reveals uplifted regions and vice versa.The undisturbed drainage networks follow basic fractal geometry and thus their existence and differentiation is studied using various fractal characteristics, i.e., fractal dimension, heterogeneity and connectivity (succolarity).
Fractal dimension analysis was used to study the effects of neotectonic activities.The idea here is to quantify the influence of neotectonics on the drainage system by measuring a reduction in complexity as the deformation intensity increases.This means that for a highly vulnerable site, lower fractal dimension values should be obtained.In other words, it is expected that the drainage patterns change their original complexity (dendritic pattern) and linearize (parallel pattern) as topographical changes triggered by faults occur and thus modify the flow network geometry.Estimated fractal dimensions for three major rivers in the study area are shown in Fig. 3.They all display values close to the Euclidian dimension 1, except for the Indus River where it is slightly higher i.e. 1.14.These low values indicate the existence of controlling processes i.e. erosion and uplift on the tectonic evolution (Guillermo et al., 2004).These values suggest that due to deformational processes they have experienced a change in their natural meandering pattern, which made them more linearized.In a relative analysis based upon D values, it is concluded that the Swan River has continuously evolved on a relatively flat surface.
The fractal dimension distribution map (Fig. 4a) was prepared using the box counting method, which shows the spatial anomaly in the drainage patterns.This map shows that the boundaries of the Potwar Plateau are characterized by low fractal dimension values i.e. highly deformed regions.Low fractal dimensions are found in the interior of the plateau.It is probably due to the two stage deformation of the plateau i.e. the interior of the plateau has been deformed only at a late stage.Very low values of fractal dimensions are observed in the SPPZ near Domeli Diljabba Triangle zone, which is an independently known active region (Aamir and Siddiqui, 2006).High fractal values are observed near the main channel of the Indus and Swan Rivers.This shows that here the patterns are more dendritic and are controlled by erosion processes and thus show low vulnerability.For a detailed relative analysis, the lowest values of fractal dimension (approaching to 1) were identified in eight regions of the study area.These regions are labeled and shown in Fig. 4a and their fractal values are shown in Table 1.All these regions have potentially high vulnerability of surface deformation as the fractal dimension values are less than 1.30.As these regions are characterized by diverse neotectonics en- vironments, these values imply a likely correlation between the complex scaling properties of the evolved river network and the vulnerability of surface deformation.The drainage patterns of these sites approach linearity by diversifying the filling spaces.The low values of these regions also indicate the recent steepening of the landscape.Along region 1, near the Attock gorge, the Indus River as well as its tributaries are highly linearized.In region 2, the Jehlum River takes a sharp bend and changes its flow direction from East-West to North-South.This change in drainage pattern causes the diversity in fractal dimension value.In region 7, the Indus River is again linearized and produces very low fractal dimension values.
Although the quantification of the space filling property of drainage patterns can reveal the vulnerability to deformation, it may still be difficult to differentiate between the patterns with different empty spaces distribution.It means that different drainage patterns may have similar fractal dimension values, but their empty spaces are distributed differently.Thus, the regions can be further characterized using their deviation from translation invariance.The lacunarity analysis is very suitable for active orogens where the drainage system has been deviated from translational invariance.High lacunarity of the drainage patterns indicates that empty spaces in an image have high diversity of sizes i.e. the image is variant in terms of its empty space distribution.The drainage networks are controlled by numerous mechanisms of space filling which depend upon the active processes operating.These processes are spatially variable not only within local regions but also in different tectonogeomorphic environments.Thus, lacunarity analysis can be used to identify the steepened regions that received uplift and deformed homogeneously.The spatial distribution of lacunarity in the Potwar area was prepared using the method described in the methodology section.The high values of lacunarity represent the regions with high deviation from translation invariance.In general, medium to high lacunarity was observed in the areas with high fractal dimension.The deviation from translation invariance of different sites indicates that they have been differently influenced by geological processes.If they are controlled by tectonic processes, then they display relatively high deviations as compared to area mainly controlled by erosion processes.The locations with low fractal dimension and high lacunarity are not treated as highly deformed regions because they miss the basic assumption of vulnerability to surface deformation i.e. high fractal dimension.All eight highlighted regions with similar fractal dimension values are studied using lacunarity analysis (Fig. 4b).The lacunarity values for these regions are shown in Table 1.This table shows that the regions with high fractal dimension generally display medium to high lacunarity values.The regions 1, 2, 5 and 6 can be identified easily.The low fractal dimension values and high relative lacunarity shows that regions 2 and 6 are intensely deformed, while regions 1 and 5 are relatively less deformed as compared to all eight regions.Region 2 is located to the west of Muree and is at a junction point between the Jehlum, Kishanganga and Kunhar Rivers.This area is highly deformed and has received intense deformation during the 8 October 2005 earthquake (Shahzad et al., 2007).Geologically, this is the region of Hazara Kashmir Syntaxis, an active portion of the India-Eurasia collision.Region 6 is the area between Talagang and Chakwal, small but densely populated towns in the Potwar region.The surface of this region is developed at the interaction between the Domeli Diljabba triangle zone, the Salt Range and the Swan syncline.It is progressively deformed due to the joint interaction and activation of these tectonic features (Aamir and Siddiqui, 2006).Region 5 is a relatively less deformed zone but has been strongly influenced by the activities along the Domeli Diljabba triangle zone.The drainage linearization in this region may be attributed to a two phase deformation of the Himalayan orogeny i.e. initial disconnection and then rough steepening.This analysis shows that regions 4, 3, 7 and 8 have similar fractal dimensions and lacunarity values.The lacunarity plot (Fig. 4d) shows that it is difficult to differentiate between these four regions because the lacunarity curves are overlapping or intersecting each other.This is because of the fact that the drainage textures in these four regions display similar space filling behavior.Thus the connectivity i.e. succolarity analysis proves to be a useful tool for this purpose, as it is a complimentary approach to fractal and lacunarity analysis.The confluence of two streams at high angles is usually attributed to the localized tectonic control.If the river systems have similar patterns of occupying space and display translation invariance, then a succolarity analysis helps to quantify their orientation and contribute to a further characterization.
In this study, this analysis was applied on those four regions where fractal dimension and lacunarity analysis cannot differentiate the relative vulnerability.In this case, the texture is differentiated not on the basis of occupation space and translation invariance, but on the basis of the ordination.This analysis was performed using the minimum of mean succolarity values along four possible directions, i.e., horizontally and vertically (Fig. 4c).The spatial distribution of succolarity values divides the Potwar Plateau in three zones i.e. eastern, central and western zone.These three zones have different structural orientation and thus the extreme low values of succolarity can be used to mark their rough boundaries.This is because the drainage orientation abruptly changes along the faults.Recent studies (Shahzad et al., 2009;Aamir and Siddiqui, 2006) also suggest that the structural orientation of Potwar Plateau varies from east to west as is clear from the seismotectonic map of the study area.The percolation property of the drainage networks of the regions 3, 4, 7, and 8 (Fig. 4d) suggest that the relative vulnerability of these regions is descending as follow: 3, 8, 4 and 7 as shown in the Table 1.Region 3, which lies near Khair-i-murat Fault in the NPDZ is more deformed.As succolarity changes with factor of division or box size k, the gradual decrease in succolarity indicates that the draining capability of this region is very high at a low factor of division.The drainage heterogeneity is attributed to the presence of anticlines and synclines in this region.The major cause of the drainage linearization is the lithological control and drainage in a plain area.In Fig. 4d, the succolarity values of regions 4 and 7 cut each other at a factor of division 3.This means that only after this threshold value we can differentiate between the two patterns.Region 7 is linearized due to the effects of the Indus River erosion and tectonics of the Kalabagh Fault.Thus the eight regions with high vulnerability for surface deformation have been graded with decreasing vulnerability from 2, 6, 3, 4, 8, 7, 5, and 1.In other words, a prediction can be made about the relative surface deformation from non-linear analysis of drainage patterns.These analyses can also inform us about the steepness, homogeneity and orientation of relative surface deformation.

Conclusions
Nonlinear analyses are useful methods for the quantitative description of drainage system evolution.The variation in fractal dimension of the Indus, Swan and Jehlum Rivers suggests a variability in the structural control in the region.The fractal dimension distribution map shows the spatial anomaly in the drainage patterns near eight selected regions identified as vulnerable to surface deformation.The lacunarity analysis was applied to differentiate between these eight regions as it is a complimentary method to fractal dimension analysis.Since the lacunarity analysis could only distinguish the vulnerability between four regions, i.e., regions 1, 2, 5, www.nonlin-processes-geophys.net/17/137/2010/ Nonlin.Processes Geophys., 17, 137-147, 2010

Fig. 3 .
Fig. 3. Fractal dimension estimation of Swan, Indus and Jehlum Rivers with insets showing the orientation of the rivers.

Fig. 4 .
Fig. 4. (a) Fractal dimension distribution map of the study area.The low values of fractal dimension represent highly deformed regions.Eight sites with low values of fractal dimension are highlighted in all maps.(b) Lacunarity distribution map of the study area where high lacunarity values represent the regions with high vulnerability of surface deformation.(c) Succolarity distribution map of the study area.The drainage network is also shown in the background.The low values of succolarity are in regions of high drainage density.(d) Lacunarity and succolarity curves of the selected regions were used to differentiate the relative vulnerability of surface deformation.

Table 1 .
Summary of the nonlinear measures i.e. fractal dimension, lacunarity and succolarity analyses of eight selected regions.