Quantitative analysis and visualization of nonplanar fault surfaces using terrestrial laser scanning (LIDAR)—The Arkitsa fault, central Greece, as a case study

Many fault surfaces commonly observed in outcrop are markedly nonplanar, and often display corrugations, asperities, and longer wavelength low-amplitude curvature. Such irregularities can have major effects on the mechanical behavior of the fault during slip, and a better understanding of the interplay between irregular geometry, fault kinematics, and dynamics is an important step in improving our knowledge of a variety of geological processes, including earthquake nucleation, fault reactivation, fault seal, and fluid flow in faulted reservoirs and aquifers. Despite this, most modeling of fault dynamics assumes idealized planar geometries, and no failure criteria have yet been developed for curved faults.

One of the main limitations in incorporating realistic fault geometries into dynamic fault models is the paucity of quantitative studies of three-dimensional geometry of natural fault zones. In this paper we present detailed data acquired using terrestrial laser scanning (TLS), also known as ground-based LIDAR (light detection and ranging), from well-exposed fault surfaces in an area of active tectonism, the Aegean. The main purpose of this paper is to describe methods to process, visualize, and analyze these kinds of large geospatial data set, including techniques to highlight spatial variation in different geometric parameters across a fault surface.

The study site for this work is the area of spectacularly exposed fault surfaces of the Arkitsa fault zone, reported by Jackson and McKenzie (1999), located on the southern side of the northern Gulf of Evia in central Greece (Fig. 1A). The area is currently undergoing active north-south to northeast-southwest regional extension, with rates on the order of 1–2 mm/yr (Clarke at al. 1998), in response to the complex tectonic interplay between subduction beneath the Hellenic Arc, backarc extension in the Aegean, and the westward movement of the Anatolian plate (e.g., Dewey et al., 1986; Doutsos and Kokkalas, 2001; Kokkalas et al., 2006). In central Greece, much of the extensional strain is localized into a number of west-northwest–east-southeast–trending grabens, which are mostly characterized by complex geometries in map view (Fig. 1B), and a high degree of segmentation along strike (Doutsos and Poulimenos, 1992; Ganas et al., 1998; Kokkalas et al., 2006). Most of the major normal faults in the region are thought to have been active in the Pleistocene, although the relative time of their activity is not well constrained. Some historical earthquakes have been recorded in the area, but none can be directly associated with the Arkitsa fault zone, or with any other fault with certainty, other than the 1894 rupture events of Atalanti fault (Ambraseys and Jackson, 1990).

The Arkitsa fault zone together with the Agios Konstantinos and Kamena Vourla faults form a west-northwest–east-southeast, left-stepping, north-dipping fault margin that links with the southern margin of the Almyros-Sperchios graben (Fig. 1). The Arkitsa fault zone, with a length of ~10 km, separates Late Triassic– Jurassic platform carbonates in the footwall from lower Pliocene–Quaternary sediments in the hanging wall. The amount of throw apparently varies along strike, typical of highly segmented fault zones. The minimum throw is estimated to be ~500–600 m, based on the thickness of Neogene–Holocene sediments in the hanging wall, plus the scarp height and topographic relief of the footwall block. This allows a slip rate of 0.2–0.3 mm/yr to be estimated for the Arkitsa fault zone, given that activity on faults of the Evia rift zone started within the past 2–3 m.y. (Ganas et al., 1998). Slip rates for the adjacent Atalanti fault show similar values, on the order of 0.27–0.4 mm/yr (Ganas et al., 1998). In addition to the overview given in Jackson and McKenzie (1999), further description was given by Kokkalas et al. (2007) and Papanikolaou and Royden (2007).

There is evidence of Holocene seismic activity along the fault, based on the geomorphological expression of the Arkitsa region, as well as the back-tilted terraces on the hanging-wall block of the Arkitsa fault, and the exposure of an ~1-m-high scarp of unweathered fault surface that existed prior to quarrying. In addition, the coastline in the region shows evidence of Holocene uplift and subsidence that is controlled by the motion and location of the faults (Roberts and Jackson, 1991; Stiros et al., 1992).

The Arkitsa fault is best exposed in an outcrop 450 m south of the main Athens-Lamia national road (Fig. 2), where recent quarrying activity during the past 20 yr has removed hanging-wall colluvium to reveal clean, fresh exposures of the upper 65–70 m of three large fault panels. These fault panels form the northern scarp edge of a footwall block that rises to 255 m above sea level (a.s.l.). The hanging-wall sediments form the low-lying plain at 70 m a.s.l. We have recorded the detailed geometry of the fault panels using terrestrial laser scanning, and this paper presents a description of the acquisition, visualization, and preliminary analysis of the laser-scan data.

A detailed survey of the exposed fault planes and surrounding footwall was carried out in early December 2005 using terrestrial laser scanning (TLS), also commonly termed ground-based LIDAR (light detection and ranging). TLS is a very rapid surveying method in which the precise geometry of the outcrop can be measured in detail by recording the three-dimensional (3D) position of many millions of points across the surface of an exposure. Since the method is nonpenetrative, the raw output is a surface 3D dataset (Dogan Seber, 2007, personal commun.; Jones et al., 2008), containing no subsurface data from within the outcrop. The high resolution provided by TLS allows outcrops to be analyzed quantitatively in great detail, leading to increased use in geosciences (e.g., Ahlgren et al., 2002; Jones et al., 2004; Bellian et al., 2005; McCaffrey et al., 2005; Pringle et al., 2006; and many others. See also Wawrzyniec et al., 2007, and associated papers in the Geosphere themed issue “Unlocking 3D earth systems—Harnessing new digital technologies to revolutionize multi-scale geological models”).

The equipment we used at Arkitsa was a Riegl LMS-Z420i (Fig. 3) (LMS—laser measurement system), which is a long-range scanner with an effective range of up to 1000 m, a maximum resolution (i.e., angular acuity) of 0.004° (equivalent to 7 mm at 100 m), and a specified longitudinal precision of 5 cm. The long range of the scanner was essential for enabling us to capture the morphology of the entire footwall block along the length of the hillside, and together with its high-resolution color capability and very rapid speed of data acquisition (12,000 points/s), this made the LMS-Z420i an ideal instrument for this project. However, this model of scanner may be less suitable for some other kinds of close-range project, because data acquired when the scanner is used in normal mode generally contain significantly more noise than shorter range scanners such as the Trimble GS100/GS200/GX200 or Leica HDS3000 (Lichti and Jamtsho 2006; Renard et al., 2006; Geospatial Research Ltd., unpublished tests). Therefore, when using Riegl LMS scanners (e.g., Z360, Z390, Z420i) for detailed studies of small-scale surface morphology (e.g., Renard et al., 2006; Sagy et al., 2007), it is important to improve precision by changing the mode of operation to acquire multiple scan sequences. This is a straightforward option in the RiScan Pro software that is used to drive the LMS-Z420i, and in our experience can increase longitudinal precision to 10 mm or better.

The fault zone was scanned using a standard acquisition strategy of taking multiple overlapping scans from a number of different tripod positions (scan stations) along the length of the outcrop. This aims to minimize any large gaps in the data (data shadows) caused by the irregular topography of the outcrop surface. In practice, there will always be some data shadows, but careful planning of the data acquisition can help to keep these to a minimum. In this study, 10 scan stations were used along the outcrop length, with 9 stations sited <75 m from the outcrop, and 1 station ~500 m away to provide extra coverage of the footwall hillside above the fault scarp (Fig. 2). Over a period of 2 days, 25 scans with a total of 75 × 10⁶ points were acquired from these 10 stations.

Data from each of the scan stations need to be merged together into a single model. This is achieved using two methods in combination. First, a set of highly reflective cylindrical targets is placed carefully in prominent positions spread around the survey area. Then from each scan station, all the reflectors that are visible from that station are scanned at high resolution. By matching the same reflector from different scan stations, all the different scans can be stitched together. In addition, a differential global positioning system (dGPS) was used to locate the positions of the most prominent reflectors, so that they can be georeferenced to a global coordinate system. For this we used an Ashtech Promark 2 in survey mode, capable of subcentimeter precision if the satellite configuration is suitable (McCaffrey et al., 2005).

A number of standard processing steps are needed to combine the raw scan data into a virtual copy of the real outcrop. The virtual outcrop is then suitable for further geospatial and geological analysis.

The main functionality of a terrestrial laser scanner measures only the distance traveled by the laser pulse, plus the intensity of the beam that is returned. This provides a grayscale image of the outcrop (Fig. 4A). In our experience, the capability to superimpose high-resolution color information onto the raw point cloud (Fig. 4B) can be extremely valuable during subsequent geological interpretation based on the virtual outcrop (see White and Jones, 2008, for further examples). There are a number of different solutions available to capture additional color information. Several scanner manufacturers use on-board video circuitry mounted inside the scanner. The Riegl LMS scanners use a separate off-the-shelf digital single lens reflex camera (dSLR), which is tightly calibrated so that the optical properties of the cameralens combinations are known. The camera is attached to the top of the scanner using a precision-engineered mount (Fig. 3). In this way, the pixels in the 2D images taken with the camera can be accurately mapped to the 3D point data from the scanner. Although using a separate dSLR requires more time and effort during data acquisition, the main advantage is that this solution provides very high resolution and faithful color rendition.

Scans taken from different scan stations are merged together based on common reflector targets and dGPS, as described above. We used dGPS to measure 16 tie-point positions—12 of which gave subcentimeter spatial precision following post-processing. These 12 positions were used for georeferencing (10 reflector positions plus 2 dGPS base stations).

Semiautomatic merge functionality is included in the RiScan Pro software, which iteratively minimizes the error involved in matching scan stations to each other. For this project the mean of the average errors (i.e., the mean of matching all stations to each other) was 5 mm.

Once the colored scans from different scan stations have been merged, then the data can be rendered as a virtual outcrop, either as a large model covering the entire area scanned (Fig. 5A) and/or as higher resolution models for detailed analysis of more localized parts of the outcrop (Fig. 5B,⁵⁰²; see also Animations 1 and 2 associated with Figs. 5A and 5B, respectively). Normally, at this stage of our specific work flow, the data are still in the form of a point cloud (Fig. 4B; i.e., the outcrop surface has not been meshed), as this is usually adequate to give a good impression of the nature of the outcrop prior to further processing.

The next step in the processing work flow generally depends upon the overall purpose of the project. If the main priority is to recreate a virtual copy of the outcrop that is as true to life as possible, then the point cloud data can be filtered, meshed, and draped with bitmap polygons taken from the scan photos (Figs. 4D–4F). However, for the work presented in this paper, the aim has been to use the TLS data to derive quantitative information about the geometry and kinematics of the main fault surfaces, and so the main priority was to extract the scan points from the large fault panels and isolate them from the background data. This meant that much effort was needed to clean the scan data by removing all traces of vegetation, talus, colluvium, weathered parts of the outcrop, and points reflected off airborne dust particles (Fig. 4C). This can be a very laborious process, and automated approaches to data cleaning have the potential to save considerable time and effort (Erik Monsen, Schlumberger, 2006–2008, personal communs.).

Before the geometry of the cleaned fault surfaces can be analyzed, it is necessary to filter the data to reduce the amount of points down to a more realistic size ready for computational processing. At this stage, each of the three fault panels is represented by between 10 × 10⁶ and 16 × 10⁶ points, and for computational efficiency it is ideal to reduce this to <200,000. More critically, it is very important that the data are smoothed in order to eliminate any irregularities (noise) on scales smaller than the precision of the scanner. In practice, both these aims can be met by using RiScan Pro's octree filtering functionality (e.g., Foley et al., 1997), typically using an edge threshold of 10–15 cm. The filtered result is a heavily decimated point cloud (Fig. 4D) that can then be meshed (Fig. 4E), and if necessary also image draped (Fig. 4F). This produces a processed data set that is ready for more specific visualization and geometrical and geological analysis.

Detailed analysis of fault slip surfaces can lead to better understanding of fault zone processes and fault dynamics. TLS has already proven to be an effective tool for quantifying surface roughness of faults (Renard et al., 2006; Sagy et al., 2007). Roughness is a useful parameter to characterize the overall surface irregularity of rock (or other materials), and can be measured along 1D sampling lines in specific directions (e.g., Brown and Scholz, 1985; Power et al., 1987), or as a 2D property of the surface as a whole (Renard et al., 2006). In this paper we present additional methods of analysis that allow other geometrical properties of fault planes to be studied. This approach shows that different types of analyses tend to highlight different geometrical properties of the faults, and therefore a range of methods should be used in combination when analyzing the 3D geometry of fault surfaces.

Unlike roughness, which is typically expressed as an overall bulk property of a fault surface, attributes such as orientation, curvature, and distance from the mean plane are more usefully analyzed in terms of their spatial variation across the surface of the fault. After such geometrical properties are derived for each point on the surface, the calculated values can then be mapped back onto the fault surface by recoloring each corresponding point in the virtual outcrop data. Displaying the colored fault surfaces using 3D visualization software provides a powerful tool to highlight spatial variation in fault geometry (e.g., Fig. 5C).

This is a simple method to highlight the deviation of a fault surface from a plane, in this case the mean plane through all measured points (nodes) on the fault. This method was used by Renard et al. (2006) and Sagy et al. (2007) to visualize fault irregularities. There are a number of different ways to calculate the orientation of the mean plane, including 3D regression methods and analysis of the moment of inertia (Fernández, 2005). A useful alternative is to use vector calculus to derive the orientation of the mean of all node normals. This approach is efficient, since we need to calculate node normals as an intermediate step for other types of analysis (see following for more details about the calculation of node normals and the mean plane). Once the orientation of the mean plane is known, 3D trigonometry can be used to measure the perpendicular distance of each point from the mean plane. The values for perpendicular distance are then projected back onto the fault surface and shown using a standard color ramp (Fig. 6A). To highlight more subtle geometrical variations, several cycles of the color ramp can be used (Fig. 6B).

In essence, this method of analysis is comparable to a contour plot of the fault surface. It is therefore particularly suitable for highlighting culminations and depressions on the surface, and is capable of revealing relatively subtle features of fault surface topography. In contrast, the method is less suitable for highlighting the precise orientation and wavelength of corrugations and striae, for which methods involving variation in orientation and curvature are more appropriate.

This section describes a method for visualization of the spatial variation in the orientation of a fault surface across the outcrop. Variation in orientation is calculated relative to the mean orientation of the fault panel, as described below. Once the scans from a fault panel have been merged, cleaned, filtered, and meshed (Fig. 4), the panel geometry is represented by tens of thousands of points (nodes), connected to form a mesh of small triangular polygons (Figs. 4E and 7A). For each node in the mesh, the 3D positions of adjacent nodes are used to define the node's tangent plane (this is a standard procedure in 3D computer graphics; e.g., Foley et al., 1997), and the vector normal of the tangent plane represents the pole to the fault plane at that point. Therefore, plotting all the vectors on a stereonet shows the amount of variation in the orientation of a single fault panel (Fig. 7B). In practice, not all stereonet software is capable of handling tens of thousands of points, but for panel B analyzed here, a random 10% subsample of the entire data set is sufficient to show that the orientation of the panel varies by at least 40°–50° in strike and 20°–30° in dip.

To show how this variation in orientation is distributed spatially across the fault panel, each point in the filtered point cloud is allocated a color according to its orientation relative to the mean orientation of the panel as a whole. This entire process is implemented in software and can be rapidly applied to the TLS data in a single operation, but in order to explain and illustrate the underlying concept of the method, it is described here as a series of separate steps.

The mean orientation of the fault panel is taken as the vector mean of all the normals to the tangent planes for all the nodes. Strictly speaking, each normal should be weighted according to the area of the polygons surrounding the node, but the difference is extremely small unless the meshing process has produced a large range in polygon sizes. The angle between each of the normals and the mean orientation of the fault panel is also calculated at this stage (this angle is the scalar product of the cosines of the two vectors). The strength of the color assigned to each node is dependant on the magnitude of this angle. This is easier to conceptualize if the data are rotated so that the mean orientation of the panel is at the center of the stereonet (Fig. 7C). For each point, the shade (hue) of the color applied depends on the relative azimuth of the point, with red for values to the west, cyan to the east, yellows and greens in the southwest and southeast quadrants, and magentas and blues in the northwest and northeast quadrants. The strength of the color (saturation) is proportional to the angular distance (plunge) of the point from the center of the net, with gray used for values at the center, and stronger colors with increasing angle outward. The color saturation is ramped to reach 100% at 4 standard deviations from the mean (the small circle in Figs. 7C, 7D), so that all outlying points are allocated full saturation. Using 4σ is an arbitrary amount, but using a threshold in this way (rather than the entire net or the full range of angular differences) helps to increase the sensitivity of the method. For panel B, 97.6% of the data are within 4σ of the mean.

Once a color value is calculated for every point, the color is allocated to the corresponding node in the meshed view of the fault panel (cf. Figs. 7D, 7E). Our software implementation outputs the resultant colored mesh in either VRML (virtual reality modeling language; www.web3d.org/x3d/specifications/vrml/) or VTK (Visualization Toolkit) format (www.vtk.org), suitable for display in many standard graphics packages, including Kitware's open-source program ParaView (www.paraview.org). Meshed surfaces colored in this way help to highlight surface irregularities of different sizes, particularly corrugations, fault bends, and secondary crosscutting faults (Figs. 7E and 8). Small-scale corrugations (~0.5–5 m wavelengths) appear as thin, bright bands of contrasting color, while fault bends representing larger scale changes in orientation (~10–50 m wavelengths) appear as broad belts of overall color variation. For example, Figures 7D and 7E show wide areas of predominantly northwest- southeast strike in red colors, separated by a region of west-northwest–east-southeast strike (shown in cyan). Thin traces shown in green and yellow correspond to small-scale tension fractures and thrust faults that cut the main surface of the fault panel.

The curvature of a surface is a measure of its departure from planarity. There are a number of different ways to express curvature, including normal, Gaussian, and mean curvatures (see Appendix). These were described in relation to geological surfaces by Lisle (1992, 1994), Bergbauer and Pollard (2003a, 2003b), and Pearce et al. (2006). In mathematical terms, curvature is a second-order derivative of the surface, and hence curvature values are independent of the surface's orientation (i.e., if a surface is rotated, the orientation of the surface at a given point will generally change, but curvature properties of the surface at that point will not). For folded or corrugated surfaces such as the fault panels at Arkitsa, curvature values vary across the surface, and analyzing the spatial variation of different curvature parameters is a useful step in understanding the three-dimensional geometry of the surface.

Measurement of curvature is highly scale dependant, and analysis will emphasize the shortest wavelength structures of a surface. Hence, it is important to remove any high- amplitude, short-wavelength noise from the data, as this will mask genuine features of surface geometry (e.g., Bergbauer and Pollard, 2003a). In this study we have minimized noise by resampling and smoothing the cleaned fault panels shown in Figures 4, 7,⁷⁰², and 8.

The results of the curvature analysis are shown in Figure 9 for fault panels A and B. Displaying the magnitude of normal curvature shows the location of the larger corrugations and cross faults (Figs. 9A, 9B; Animation 3), but is not particularly effective in showing more subtle geometrical features. Sensitivity of the plots can be increased (Figs. 9C, 9D) by plotting the absolute values of normal curvature, and by using an upper threshold (such that the color scale is no longer stretched to include small areas of anomalously high curvature). This method is an effective way to highlight the corrugations on the surface, and is the best plotting method to use for precise measurements of geometrical aspects such as wavelength of the corrugations. Based on peak-to-peak measurements from these plots, the average wavelength of meso-scale corrugations on panel A is 4.04 m, and on panel B is 4.43 m.

Fault-plane corrugations are believed to form from overlapping faults by two main mechanisms, including lateral propagation of curved fault tips and linkage by connecting faults (Ferrill et al., 1999). Such mechanisms operating on initially en echelon fault arrays might explain the meso-scale corrugations seen in each of the Arkitsa fault panels. Large-scale corrugations may also form by the progressive breakthrough of originally segmented fault systems (e.g., Childs et al., 1995), which is probably the most likely cause of the decametric-scale warping at Arkitsa. Other mechanisms for the formation of corrugated fault surfaces include reactivation of preexisting faults and fractures (e.g., Donath, 1962), and later folding of a preexisting fault surface (e.g., Holm et al., 1994; Krabbendam and Dewey, 1998).

Gaussian curvature can be a useful parameter to highlight domes, basins, and saddles on the fault surface. It is defined as the product of the principal curvatures, such that it is nonzero if there is curvature along both of the principal directions. Thus Gaussian curvature provides a quantitative measure of culminations on the surface. Plots of Gaussian curvature from panels A and B are shown in Figures 9E and 9F. These plots show that there is negligible Gaussian curvature across the fault panels, suggesting that there are no noticeable culminations on the panels (at least at the scale at which the surfaces are smoothed). This contrasts with the plots shown in Figure 6, which clearly show that domes and basins are present. This apparent contradiction arises because of the scale dependence of the curvature analysis: the domes and basins in Figure 6 have wavelengths of ~10¹ m, while the smoothing of the data used in Figure 9 was chosen to highlight the corrugations (with typical wavelengths ~10⁰ m). In order for the curvature analysis to detect and quantify larger scale structures, it is necessary to smooth the data more coarsely, to iron out small-scale features.

The slight mottling effect seen in parts of the plots represents very small variations in curvature values. These are caused by minor residual noise in the smoothed surface, and reflect the limits of precision of the laser scanner equipment. In these specific examples, increasing the sensitivity of the plotting method only emphasizes the small amounts of noise in the data, without highlighting any other aspects of surface geometry.

The various curvature plots shown in Figure 9 help to emphasize that the output of curvature analysis can be heavily dependant on the processing steps involved (particularly the way that the raw data are cleaned and smoothed), and the visualization parameters used for plotting the results. Quantification of curvature parameters is very useful, but because the curvature of fault surfaces is generally quite subtle (e.g., compared with the curvature seen in folded bedding), more time and effort is needed to fine-tune the variables used for analysis and visualization. For this reason it may often be ineffective to use geological software packages that have black box curvature functionality that cannot be customized on a case by case basis.

A number of previous studies have shown that the orientations of slip lineations can show significant systematic variation along the strike of a fault, particularly in parts of the fault that are close to the tip zones or in relay zones (Roberts, 1996; Maerten, 2000; Maniatis and Hampel, 2008). The use of laser scanning to make a large number of precise, spatially referenced measurements of slip lineations along the length of the fully exposed fault panels allows us to test whether systematic variations in slip direction are also present in this part of the Arkitsa fault.

At Arkitsa, fault striae that are subparallel to millimeter- to meter-scale corrugations are very prominent on the surface of the fault panels in outcrop and are also clearly visible in the TLS point clouds (Figs. 5–9,⁵⁰²). In addition, there are also other orientations of fault striae locally preserved; however, these are less well developed, and are below the centimeter resolution of our TLS data. Further work is needed to relate these secondary slip directions to fault kinematics and regional tectonics, and in the following analysis we focus only on the most prominent set of slip lineations.

Although the orientation of the fault surfaces shows significant variation within a single panel (Figs. 7,⁷⁰² and 10), the orientation of the slip vectors on each panel is remarkably consistent (Fig. 10). Along the polished fault panel surfaces, the well-defined striations concentrate along a north-northwest–south-southeast azimuth (mean azimuths 336°–344°), plunging toward the north-northwest (mean plunge values 51°–57°). This configuration is consistent with slightly oblique normal dip-slip motion with a component of left-lateral displacement. Within each fault panel, typical variation from the average plunge values is ~4°–9°, while variation from mean slip azimuth is ~7°–15°.

Our analysis is based on 902 slip vectors picked from the digital surface (124 measurements for panel A, 678 measurements for panel B, and 100 for panel C). This is a comparable though more extensive study to that presented in Kokkalas et al. (2007), for which 121 slip vector measurements were used. Statistical properties for the entire data set are given in Table 1. The eigenvector analysis we applied is a non-parametric method that is particularly suitable for analyzing directional data (Fisher et al., 1987; Mardia and Jupp, 2000). Eigenvectors are ordered such that E₁ ≥E₂ ≥E₃, with E₁ being the primary orientation of slip vectors. The normalized eigenvalues (S₁, S₂, and S₃) reflect the variance of the vector data. For example, panel A with an S₁ value of 0.9981 shows that 99% of the total data are consistent with the specified orientation of E₁, within statistical significance. Similar values of 0.9952 and 0.9971 have been derived for panels B and C, respectively. The relative size of each S value indicates the shape of the distribution of the data. Data having S₁ >> S₂ ≈ S₃ suggest a highly clustered distribution, as is the case for all three examined fault panels.

In addition, the fabric shape parameter K can be used to test for girdle distributions (K < 1.0), clustered patterns (K > 1.0), or a uniformly random distribution (K ≈1.0), while parameter C determines the strength of the shape parameter K (Woodcock, 1977; Woodcock and Naylor, 1983). If C > 3.0, a preferred orientation and shape are strongly indicated, and if C < 3.0, there are weak orientation and shape fabric of the data. Values of K and C parameters for all three panels show marked clustering with strong preferred orientation and consistent fabric shape.

The statistical analysis of populations of minor faults, usually based on measurements of fault orientation, slip direction, and sense of slip, is commonly used for areas of brittle deformation. Results from this type of fault slip analysis are typically used to make inferences about the stress field at the time of deformation (assuming that the applied stress is homogeneous). By definition, fault-slip data provide more direct information about kinematics than the mechanics of faulting, with slickenlines used as direct indicators of displacement directions on fault planes. In general, this analysis is applied using individual fault-slip measurements from a large population of striated fault surfaces, rather than many measurements from a single nonplanar fault surface. Using TLS data to make precise measurements of large numbers of slip lineations and fault orientations from a single fault surface also gives a unique opportunity to determine the local stress (strain) tensor for each fault panel, and hence to test the consistency of slip direction between the different panels. Stress analysis using 100–150 fault slip data from each panel suggests that the Arkitsa fault zone is slipping in response to an extensional stress field with a north-south orientation of minimum principal stress (σ₃), with a small deviation toward the north-northeast orientation calculated only for panel B (for further details, see Kokkalas et al., 2007).

Although the stereonet plots (Fig. 10) and statistical analysis (Table 1) show that the slip data are very tightly clustered, we have also analyzed the data to try to detect small, statistically significant spatial variations in the orientation of slip vectors across the fault panels. This is because, as noted by Jackson and McKenzie (1999, p. 5), some striae appear to curve toward the top of the exposed panels. Although we concur with their field observation (some parts of panel B do show gradual variation in orientation of ~1°–2° between the top and bottom of the outcrop), this variation appears to be of localized extent, and apparently does not reflect the whole of panel B, or the other exposed panels. Note, however, that we have not yet been able to make quantitative measurements from the TLS data that are precise enough (±0.5°) to refine the estimated variations in slip vector orientation made by Jackson and McKenzie (1999). Further work is needed to optimize the processing of the TLS data prior to repicking of the fault striae.

The present-day regional stress field, derived from fault-slip data, focal mechanisms of recent earthquakes, and several GPS campaigns, represents north-south to north-northeast–south-southwest extension of central Greece (Roberts and Jackson, 1991; Armijo et al., 1996; Roberts and Ganas, 2000; Hollenstein et al., 2008). Recent GPS velocities and strain rate tensor calculations, derived from the velocity field, indicate that velocity differences observed across the North Evia Gulf are small (2–4 mm/yr; Hollenstein et al., 2008). This is consistent with the low seismic activity of this rift structure. Regional strain rates can be resolved into normal and shear components relative to the west-northwest–trending Arkitsa–Kamena Vourla fault zone (Figs. 12 and 13 of Hollenstein et al., 2008). Normal strain rates represent low extensional strain values in a north-northeast orientation, while shear strains indicate sinistral motion along the Arkitsa–Kamena Vourla fault zone. These geodetic calculations fit closely with the north-south orientation of the stress tensor we derive from the Arkitsa fault slip data, indicating that the stress/strain tensor has not changed significantly in orientation since the last major slip event recorded on the Arkitsa fault surface. Despite the close match between our fault data and regional strains inferred from GPS, we still cannot be certain whether the precise orientation of oblique slip we have measured at Arkitsa is entirely indicative of regional strain, or whether there are also additional, more localized effects close to the eastern tip of the Arkitsa–Kamena Vourla fault zone (cf. Roberts and Ganas, 2000).

Terrestrial laser scanning (i.e., LIDAR) is an effective tool for acquiring detailed geometrical data from exposed outcrop surfaces, and is therefore ideal for studying surface irregularities and curvature of nonplanar faults. Field acquisition of LIDAR data is rapid, though the millions of data points gathered will subsequently need to be merged, cleaned, and filtered prior to further geological analysis. Various methods of analysis can be used to highlight different geometrical properties of the scanned surfaces. Roughness provides a bulk measure of the overall surface irregularity of a fault, while other methods allow the spatial variability of fault geometry to be studied. These methods include measuring the perpendicular distance of points from the mean plane, plotting variation in orientation relative to the mean plane, and curvature analysis (including normal, mean, and Gaussian curvature). Studying curvature allows quantitative measurements of fault bends and corrugations, although the measurements are highly scale sensitive and require care during data preparation and analysis. Spatial variation in the various geometrical properties can be highlighted by projecting the resultant analytical values back onto the virtual fault planes (Animation 3). 3D visualization software is particularly effective in making it easy to compare the derived parameters with the underlying original fault geometry.

Analysis of three well-exposed fault panels from Arkitsa (central Greece) shows that each of the faults contains surface irregularities, corrugations, fault bends, and crosscutting minor faults, and deviates significantly from planarity. Each panel shows a spread of at least ±20° in the orientation of the surface perpendicular to the main slip vector. Normal curvatures delineate clearly the location of corrugations on the surfaces of the panels. Gaussian curvature is very low at this scale of analysis. Corrugations have axes orientated parallel to the main slip direction, and have average wavelengths of ~4 m. The orientations of the slip vectors are tightly clustered, plunging ~55° north-northwest. This is consistent with the present-day regional velocity field derived from GPS measurements. These results can be compared with measurements given by Resor and Meer (2006).

The high resolution of the LIDAR data helps to illustrate that the detailed geometry and kinematics of the fault zone are highly complex, and that different complexities exist at different scales. It is not yet clear the degree to which this part of the Arkitsa fault zone is typical of faults in general, or whether the high structural complexity is related to the regionally complicated rotational strains seen throughout the Aegean. In either case, much further work is needed to unravel these complexities on the outcrop scale, and to relate the spatial variation in fault geometry to larger scale fault kinematics and dynamics.

We thank Phil Clegg and Tim Wright for assistance during field work, Jonny Imber for discussion on faulting, Mark Pearce for collaboration on curvature, and Steve Waggott and Ruth Wightman for collaboration related to laser scanning. Thomas Dewez kindly commented on the manuscript, and Emily Brodsky, Ioannis Papanikolaou, and Ramon Arrowsmith provided useful reviews and editorial comment.

The normal curvature is the magnitude of curvature corresponding to the maximum and minimum principal curvature directions, k₁ and k₂. These values help to define the 3D traces of hinge lines along folds or corrugations. It can sometimes be useful to visualize absolute (i.e. modulus) values of normal curvatures (|k_n|), as shown in Figures 9C and 9D.

Mean curvature, M, is the average of the two principal curvature values (Roberts, 2001; Bergbauer and Pollard, 2003a, 2003b):

Since curvature is a function of the second derivative (equations 1, 3, and 7), it is an invariant property of the curve. This means that curvature values are independent of the dip of the fault surface. In contrast, the positions of crests and troughs are not invariant features, but are closely dependent on the orientation of the fault. For this reason, the method of analysis of spatial variation of fault surface orientation shown in Figures 7,⁷⁰² and 8 is done relative to the mean fault plane.