Extensive, very high quality outcrops on the southern flank of Jabal Shams, Oman Mountains, expose the complex fracture and fault network of an exhumed high-pressure cell. Veins are filled with white calcite in grey host rock, allowing 0.7 meter-resolution mapping based on satellite image interpretation, followed by ground-truthing with detailed field observation in selected areas. From Quickbird and Landsat data, a model was developed using interpreted 562 faults and 145,000 fractures in an area of 31 km2; the area of a typical 3-D seismic survey but with one hundred times higher resolution. Four generations of veins, with strikes of 130°, 000°, 090° and 045°, are overprinted by normal faults with wide cemented damage zones and offsets up to 500 m. Veins are generally a few 10s of meters long, and lack signs of strong mechanical interactions such as curving or abutting. This suggests a restoration of the bulk strength during the formation of the joints, and an anticlockwise rotation of permeability anisotropy during progressive fracturing. The spatial density of all cemented fractures is rather homogenous, but the spatial density within individual joint sets is heterogeneous, suggesting patchy fracturing during the different evolutionary stages. A weak mechanical anisotropy, caused by the veins, affects the nucleation and evolution of the anastomosing normal fault system that overprints the cemented joints.


High-pressure cells are common in the Earth’s upper crust (e.g. Bradley, 1975; Powley, 1990; Bradley, 1994; Ortoleva, 1994; Ortoleva et al., 1995; Finkbeiner et al., 2001; Nollet et al., 2005; Hilgers et al., 2006a; Hilgers et al., 2006b). Understanding the cause, sequence and magnitude of the fracturing events is important for understanding the fluid flow in such a system and is therefore of interest for many basic and applied fields of geoscience. Models of high-pressure cells are derived from volume data (3-D seismic) with a resolution of around 10 m, often combined with outcrop studies and information from wells. Fracture networks are a prominent feature in high-pressure cells due to low effective stresses in the system. Knowledge of their spatial distribution and temporal evolution is, however, limited because they are not adequately resolved by the 3-D data, and cannot in general be mapped in outcrops of limited extent (Bradley, 1975; Noorishad et al., 1984; Patton and O’Conner, 1988; Powley, 1990; Bradley, 1994; Dewers and Ortoleva, 1994; Ortoleva, 1994; Ortoleva et al., 1995; Nguyen and Selvadurai, 1998; Finkbeiner, 1999; Ingram and Urai, 1999; Petit et al., 1999; Townend and Zoback, 2000; Finkbeiner et al., 2001; Abrams, 2003; Hilgers et al., 2006a; Nollet et al., 2006).

This study integrated high-resolution remote sensing and field data (Holland et al., 2009) to investigate how the meso- to microscale observations in outcrops can be laterally extrapolated. It discusses the required scale to describe the governing structures (Gillespie et al., 1993; Bonnet et al., 2001; Shipton and Cowie, 2001; Ben-Zion and Sammis, 2003; de Keijzer et al., 2007), and characterizes the evolving anatomy of a high-pressure cell in carbonates at high spatial resolution (Willemse et al., 1997; Rawnsley et al., 1998; Gillespie et al., 2001; de Keijzer et al., 2007; Holland et al., 2009). Multi-scale data from field observations and image interpretations were integrated in a geographic information system (GIS). The database contains georeferenced field photographs, measurements, sketches and notes, as well as a collection of remote sensing data sets. The studied system comprises Mesozoic carbonates located in the Oman Mountains (Figure 1) (Hilgers et al., 2006a).

In the study area, multiphase deformation under close-to lithostatic fluid pressure conditions led to the formation of four joint sets overprinted by bedding-parallel shear and normal faulting (Loosveld et al., 1996; Breton et al., 2004; Glennie, 2005; Al-Wardi, 2006; Hilgers et al., 2006a; Searle, 2007; Holland et al., 2009). The ubiquitous joints are predominantly cemented with bright white calcite in gray carbonate host rock. This high optical contrast allows for the interpretation and mapping of fractures on the basis of remote-sensing data. It was carried out in an area of 31 km2 with a spatial resolution of up to 0.7 m. The study produced a high-resolution regional fault and fracture map, and analyzed the first-order characteristics of the spatial distribution of the structures. Field observations of overprinting relationships, and of structures that are smaller than the resolution of the satellite image, guided the image interpretation. This approach provided reliable constraints on the applicability of the results (Abrams, 2003; de Keijzer et al., 2007; Holland et al., 2009).


A number of digital and printed data sets were used. Printouts were used in the field together with data sets stored on a hand-held computer linked to a global positioning system (GPS) receiver. This allowed on-site comparison of the field observations to verify the relationship between remote sensing data and field observations. The map data consist of a geological map (Beurrier et al., 1986), radar-derived elevation data, air photographs together with two multi-spectral satellite images.

Geological Map

The study area was mapped at 1:100,000 scale by Beurrier et al. (1986). The geological map was used to distinguish between the Albian Nahr Umr Formation (lower Wasia Group), Albian – Turonian Natih Formation (upper Wasia Group) and Aptian Shu’aiba Formation (upper Kahmah Group).

First-order faults are also shown on this map. Figure 2 shows data from this map as geo-referenced against a Landsat image, and draped on a coarse-resolution digital elevation model (DEM). The good correspondence of the main geological units with morphology is apparent.

Digital Elevation Model

Terrain morphology data was derived from two DEM data sets. The first is a SRTM file (Shuttle Radar Topographic Mission) obtained from the Global Land Cover Facility (USGS/GLCF 2000). The one-degree tile offers a spatial resolution of 90 m (Rabus et al., 2003; Farr et al., 2007). The large off-nadir angle (>30°, near-vertical angle), during the acquisition of SRTM data (Farr et al., 2007), limits its collection near to steep morphological elements because the subvertical wadi walls and steep cliff sections cast shadows (Figure 3).

To fill these gaps we used Aster files (Advanced Spaceborn Thermal Emission and Reflection Radiometer) as a second data set obtained from the Earth Observing System (EOS) Data Gateway. Aster DEMs are generated optically by stereoscopic images made from a nadir viewing sensor and a backward viewing sensor. The two data sets were processed to a DEM with a spatial resolution of 30 m (Abrams, 2003). The combined ASTER and SRTM data sets offer a good representation of the regional morphology (Figure 2). The spatial resolution of about 30 m does not, however, fully display the structure of the deep wadis.

Landsat Enhanced Thematic Mapper

A Landsat scene (NASA Landsat Program 2000) from the Global Landcover Facility was used to derive data on lithology. Landsat images (ETM+, Enhanced Thematic Mapper) provide, in total, seven bands of spectral information and an additional panchromatic band with high-resolution data, acquired at nadir. Most bands (1 to 5 and 7) have a spatial resolution of 30 m. The thermal bands 6.1 and 6.2 (high- and low-gain bands) have a resolution of c. 60 m. All bands can be pan-sharpened against the 15 m resolution panchromatic band 8 to obtain a higher-resolution image. Combined bands 1 to 3 from the visual spectrum produce an image similar to human vision (Figures 1 and 4). The relatively low optical contrast among the carbonates on the image is a result of the high correlation of DN (data number, measure of the signal) within the first three bands. Using the near-IR (Infra-Red) and thermal-IR bands produces higher visual contrasts due to differences in the thermal properties of the rock types.

The combination of the bands 1 (0.45–0.515 µm), 4 (0.75–0.9 µm) and 7 (2.09–2.35 µm), displayed as an RGB image, can be compared with the geological map (Figures 4b and 2). In this image, the optical difference between the upper Kahmah Group (predominantly Shu’aiba Formation) and the upper Wasia Group (Natih Formation) become visible (Figure 3b).

To intensify the spectral differences of the exposed rocks, the information of all available bands can be incorporated. A subset of the Landsat scenes was created excluding areas of other rock exposures to optimize the analysis to the region of interest. The seven bands were processed in a principle component analysis, which reduces the amount of bands (or dimensions) within highly correlated sets (Pearson 1901; Fukunaga 1990). This procedure intensifies the contrast between the lithologies. Bands 2, 3 and 4 are displayed as an RGB image in Figure 4c; the significant increase in detail allows tracing the boundaries of several formations and also members (e.g. SW-flank).

However, because the small-scaled vertical changes in lithology create small and complex patterns (Holland et al., 2009), the resolution of the data set proved insufficient to consistently map lithology at this scale.

Aerial Photograph and High-Resolution Digital Elevation Model

Petroleum Development Oman (PDO) provided an high-resolution aerial photograph with a corresponding high-resolution DEM, which covers a section of the central part of the field area and its southern extension down to Wadi Ghool (Figures 2 and 5). It was acquired approximately at nadir and geo-referenced, and provides spectral information with three visual bands at a spatial resolution of 1 m. The corresponding DEM is seamless and provides elevation data with 10 m spatial resolution. The two combined data sets provide an accurate topographic model of the geometry of the strata on this southern flank. While the Landsat offers coarse information and the ASTER and SRTM data only a first-order morphology, this data pair offers insight into a number of different features related to the interplay of erosion and lithology: the differences of the weathering characteristics are reflected by small-scaled morphology as terrain steps.

Structural elements can be distinguished from morphology. One example is the stair-case geometry of the Nahr Umr Formation (Figure 4). The Formation is exposed on the footwall block of a major fault and offset against the massive Natih Formation (Holland et al. 2009). The latter protects the exposure of the argillaceous Nahr Umr Formation from erosion. The latter section shows morphology distinct from the overall dip-slope style in the area. The DEM provides vital information to understand this geomorphology (Figure 4). The same applies for a ramp structure exposed on the southern flank (Holland et al., 2009). The high-resolution data pair shows the lateral extent of the ramp as well as its dip angle.

Quickbird Satellite Image

For synoptic high-resolution mapping, a data set with constant image properties is required. For this purpose the study used a 75 km2 Quickbird scene covering the central part of the study area. The data set is c. 5.5 km × 17.5 km, geo-referenced and geo-corrected by the vendor. The scene covers the main parts of the southern slope of Jabal Shams (Figures 1, 2 and 5) and includes Wadi Nakhr – Oman’s “Grand Canyon”. The satellite acquired data with four spectral bands and one panchromatic band. The spectral bands encompass the visual spectrum and a near-IR band (0.76–0.90 µm) with a spatial resolution of approximately 2.8 m. The panchromatic band offers a significantly higher spatial resolution of 0.7 m. The used scene was acquired with 23° off-nadir angle. This leads to some distortion and data loss since the oblique view of the satellites’ sensor misses steep sections not facing the satellite (Figures 3 and 6). Although geo-corrected, the satellite image is not a full map view of the area. Consequently sections of the steep wadis were excluded from the data set; in total c. 15 km2.

For the purpose of mapping, bands 1, 2 and 4 were combined to an RGB image. The use of the near-IR band 4 (0.76–0.90 µm) allows spotting sparse vegetation, which appears as red (Figures 6 and 7). The RGB image was multiplied with the pan-chromatic image to synthetically raise the resolution of the color image to 0.7 m. The high-resolution data set resolves a great amount of detail on the shallow-dipping flank of Jabal Shams (Figures 6 and 7). At 1:10,000 and 1:5,000 scales, it shows houses, roads, trees and bushes, which can be used to accurately determine the position in the field even without the use of a GPS receiver.


The data sets were uploaded and organized in a GIS system (ArcGIS 9.0) before processing. Subsets of the data, prepared from the Landsat and the Quickbird scenes, were taken into the field as map printouts as well as in digital form uploaded to a hand-held computer. Attached to a GPS, the latter device allowed instant navigation on the different data sets and allowed convenient collection of geo-referenced information. The GPS receiver has an RMS error of less than 15 m. Correlating recognizable ground points in the field, with the GPS-location on the data sets, validated an accuracy of a few meters.

Ground Truthing

Several prominent outcrops were documented in more detail to verify the quality and the limitation of the remote sensing data. The first area is located a few 100 m east of Wadi Nakhr with a well-developed set of veins (Figure 8). The GPS system was used to navigate to single veins on the map where a strike reading was taken with a structural compass along with a GPS waypoint and photograph. This was done for thirteen veins. The veins that were photographed within this small survey all showed bright white calcite cement. The aperture of these cemented fractures varied from a few to 10s of centimeters. Thus, although the aperture is significantly smaller than the spatial resolution of 0.7 m, the brightness of the pixels is sufficiently enlarged by the vein to make it visible on the image. The visibility of cemented fractures on an image is related to a number of parameters:

  • (1) The ratio of area of cement and bedrock within a single pixel.

  • (2) The contrast of the DN (data number, quantitative measure of the signal) of the vein and the DN of the wall rock within a pixel.

  • (3) The contrast of the vein pixel to a neighboring pixel surrounding the vein.

Figure 9 illustrates these effects with a synthetic example. A dark area and a bright area have an isolated vein and a vein array. The area of the cement for both vein systems is the same. Applying a mosaic filter that averages the DN of the image within a rectangular area simulates a high- and a low-resolution image. The example illustrates that bright veins in darker areas are much better to resolve due to the higher contrast. Changes in wall rock color or brightness can therefore affect the density of interpreted joint patterns. The nature of the vein is not recognizable anymore when the spatial resolution is of the same order as the aperture. In this case, the single-strand vein and the vein array produce an identical signal (White, 1979; Adelson, 2000; Gilchrist and Annan, 2002).

For the discussed reasons the satellite’s sensor cannot resolve whether the signal is caused by a single vein or a sub-parallel vein array with the same cumulative area. Nor can it reveal kinematic indicators like splays or en-échelon segments for which field observations are necessary (Holland et al., 2009). Dark host rock surfaces without major accumulation of soil and debris consequently deliver high-quality outcrops, whereas brighter surfaces with homogenously scattered veins may not provide the sufficient optical contrast for proper image interpretation.

These optical effects are clearly shown in the second outcrop example (Figure 10), which compares outcrop and image data. The outcrop surface is slightly brighter as in the previous example and shows a dense network of veins striking in different directions. Since the DN is integrated over the entire area of a pixel (0.7 m × 0.7 m), the signal lacks a sharp contrast (compare with Figure 9). In this outcrop only the thick veins and vein arrays are resolved. The macroscopic outcrop observations show information invisible to the satellite image: the morphology of the veins incorporates en-échelon segments, braided systems and splays at vein terminations (Holland et al., 2009). Overprinting relationships and observations on the fractures cement that document repeated crack/seal events are vital for the understanding of the fracture network. This information can be obtained in outcrop only. Keeping these well-constrained limits of resolution in mind, the remote sensing data are well-suited for regional mapping of a very large number of veins.


The high-resolution data set provides a remarkable view of the spatial distribution of faults and fractures. The resolution is one order of magnitude greater than the scale of seismic data, which is of the order of 25 m. This study focused on the high quality and extremely large outcrop surfaces on the southern flank of Jabal Shams, with surface exposures of the Natih Formation in the western part and upper Kahmah Group in the eastern part (Figure 4). Bound by two major fault zones, this central region of interest exposes abundant medium-throw faults as well as a large number of veins and uncemented fractures (Holland et al., 2009).

The interpretation and mapping was done manually, guided by the data of the field campaigns, using field observations and the database of georeferenced photographs as a guide the interpreter. We mapped the following elements:

  • (1) A polygon shape file masking the steep wadis, which were excluded from the interpretation. Of the 75 km2 satellite image, approximately 15 km2 are masked.

  • (2) A polyline shape file with interpreted fault lineaments.

  • (3) A polyline shape file with interpreted fracture lineaments.

The decision whether a lineament on the image is interpreted as a fault or fracture is based on a number of criteria. Criteria for faults are:

  • (1) Ground truth information confirms the structure to be a fault with offset (Figure 11a).

  • (2) Offset is visible on the satellite image in terms of juxtaposition of lithology or morphological elements (Figure 11b).

  • (3) The width of bright zones exceeds several pixels, corresponding to a meter to a several meter-wide zone (Figures 11a, c).

  • (4) Bright zone with curvature, branching or anastomosing patterns on a scale of 10s of meters or more (Figure 11b, d).

Using these criteria, a clear distinction is not always possible between a fracture and a fault. In a limited number of cases field observations showed features to be massive veins with apertures of up to 20 cm. These so-called “mega-veins” showed no offset although being several 100 m long (Figure 11e). A misinterpretation of similar fractures as faults is possible in a number of other cases where no ground truth information is available.

Thinner lineaments in regular sub-parallel patterns are in general interpreted as fractures. Although most fractures and faults are cemented, some fractures are uncemented (joints). Although common in the field (Figure 10), only a few uncemented joints are well imaged in the satellite data. The interpretation of the faults was done on a scale of 1:3,000. The interpretation of the fractures (joints, veins) was done at a higher magnification with a scale up to 1:1,500. The interpreter refrained from interconnecting aligned lineaments to maintain the geometry of segmented structures. The fracture vertices were picked depending on the shape of the fracture. Equal spaced vertices were not picked but are recommended for robust statistics based on unit length.

In branching networks of faults the most prominent strand was picked as the continuous segment. The interpretation process was carried out over a period of three months, after a training phase. To ensure consistency, the interpreter occasionally returned to an area already interpreted and reinterpreted it to check that the results are the same and there was no drift in the picking criteria.


Based on the raw interpretation, further processing to secondary data sets was done within the GIS system using the standard toolboxes. Additional statistics were calculated within the Matlab environment. The data conversion of the shape file material was done using the ARC_MAT toolbox (LeSage and Pace 2004). In addition, functions to analyze and visualize the data were programmed.

Fault Interpretation Map and Directional Statistics

A total of 562 faults were interpreted on the Quickbird scene. Combined with throw estimates made during fieldwork (Holland et al., 2009), the results are displayed in Figure 12. Bound by two large fault zones, with throws exceeding several hundreds of meters, the majority of the interpreted faults show offsets up to a few 10s of meters. Isolated fault segments that do not abut-against or merge-with other faults are rare. The faults typically branch to form interconnected anastomosing networks.

The fault density subjectively increases towards the eastern part of the study area where the throws of the bounding northern and southern fault zones also increase. This indicates that the total extension increases to the east of the study area.

The length of the interpreted faults ranges from a few meters to 1,600 m (Figure 13). The length of a fault is in most cases difficult to determine exactly (e.g. Rawnsley et al., 1992; Yielding et al., 1992; Gillespie et al., 1993; Mansfield and Cartwright, 1996; Needham and Yielding, 1996; Ouillon et al., 1996; Watterson et al., 1996; Cello, 1997; Odling 1997; Rawnsley et al., 1998; Aydin, 2000; Cello, 2000; Bonnet et al., 2001; Gillespie et al., 2001; Ben-Zion and Sammis, 2003). Many faults extend beyond the study area. The mean fault length is 422 m within a skewed population of 562 faults. The most common interpreted fault length is 295 m.

When plotted in a log/log histogram with logarithmic binning, a power law distribution is apparent with an exponent of m = -2.68. Reported exponents for cumulative distributions of the length of faults and fractures are reported in literature in the range of -1.5 to -3.0 (Reches, 1986; Scholz and Cowie, 1990; Yielding et al., 1992; Castaing et al., 1996; Ouillon et al., 1996; Odling, 1997). The population shows a lower limit at c. 316 m, suggesting that smaller faults are under-represented with respect to a power-law distribution. This is typical in many studies because smaller structures are commonly overlooked in the interpretation.

To characterize the directional distribution, two measures were considered: (1) the direction of the entire fault or fault strand as defined by the fault tips (start and end nodes), and (2) the direction of each fault element (element defined by two successive fault nodes). Not weighted or corrected for the mean direction, these two measures show the bulk orientation (regional) and potential orientation of anisotropy or local propagation. The left rose plot in Figure 13c shows strike values of every interpreted fault element. The right plot (Figure 13c) takes the strike values of the entire fault strands.

The rose diagram and the Cartesian projection of the same data (Figures 13c, d) show peaks in the direction of 090° for the fault elements. The mean direction of the strike of the whole faults is oriented 100°.

Fracture Interpretation Map and Directional Statistics

The fracture interpretation map consists of more than 145,000 interpreted lineaments (Figures 14 and 15c). These incorporate both joints and cemented joints (veins). This massive data set allows statistical analysis related to location, strike direction and the length distribution of the elements.

Most interpreted fracture elements are straight (Figure 14) and defined by two nodes only (145,000 fractures versus 217,000 elements). Plotted in a rose diagram (Figure 15c) the main strike directions of the system are NW/SE (c. 130°), NE/SW (c. 045°), N/S (c. 000°) and E/W (c. 090°). The 130° direction is the predominant strike direction of the system. The maximum length of the interpreted fractures is c. 400 m with a mean length close to 25 m. The most common length (mode) in the interpreted set is 12 m (Figure 15a).

Figure 15b shows the cumulative frequency plotted against the length of fractures using logarithmic binning. The distribution produces a graph with a curved and a more linear segment. The linear segment has a slope of m = c. 4.2, being considerably steeper than the corresponding plot of the fault population. The length that separates the curved from the linear segment is c. 56 m as opposed to the 316 m of the fault population (Figures 13b and 15b).

The magnitude of the mean length and the maximum length was examined to test for preferential orientations (Figure 15). The mean length of the fractures is highest in a quadrant from 090° to 140°, exceeding the minimum values by more than 50%. The variance of the mean fracture length in a directional sense is rather small.

The maximum lengths are less robust as this measure is based on only one value. Trends are however visible and for example, show that the longest interpreted features strike c. 140°, 10-15° off the most common strike direction. The large scatter in this plot inhibits defining groups or correlating the values with the frequency or mean length plots. To visualize dependencies we used three scatter plots of the values of frequency, mean length and maximum length against each other (Figure 16). For this purpose we divided the main population into 90 groups based on their strike direction (2° bin size). From every group the frequency, mean length and maximum length was measured and plotted as a data point in the scatter plots. Labeled and marked with a gray value the plots show distinct groups and trends.

The mean length and the maximum length are correlated in a linear fashion (Figure 16a, b and c). The directional distribution of the maximum length is subject to some scatter. The largest lengths are interpreted in the 140° direction. These data points deviate from the neighboring data points. This variation can be based on a misinterpretation or by reactivation. The latter process could cause veins to grow in length. Although technically being faults, these reactivated features may be interpreted as fractures.

In Figure 16b, the mean length is plotted against frequency showing clustered data. This is a result of the distinct directional classes. Two data points deviate from the groups. The first is the 02 point (000–002°) and the second is the 92 data point (090-092°). Both data points show shorter lengths and higher frequencies as opposed to the data points of adjacent directions. This could be an effect of alias artifacts that biased the interpretation. Such an effect is created when structures that have a width comparable to the image resolution are displayed with different orientations in respect to the image axes.

Figure 17 is a two dimensional histogram showing strike, length and frequency in a single plot. On the horizontal axis the strike direction of the fractures is plotted using 2° bins. The vertical axis is joint length with a 1 m bin increment. The number of elements in a specific strike bin and a specific length bin is represented by different colors making this plot a directional length histogram. The plot is clipped at a fracture length of 100 m, since less than 1% (about 1,000 elements) exceeds this length. The frequency of the elements in the individual bins is plotted with a linear and a logarithmic color bar. The two dimensional histogram shows the directional distribution as well as the mode and skewness of the joint population at the same time. It becomes apparent that the lower length limit of the joint interpretation is constant among the entire population with a little less than 10 m (Figure 17 b). The Figure shows – similar to the rose diagrams – four major strike groups at c. 000°, 45°, 90° and 130°. The populations themselves differ in a number of aspects: the 130° group is the one with most elements. A large variance with the majority of the joints occurs at approximately 130°. In comparison to the other groups, the interpreted length of these fractures is much longer. The asymmetric color distribution of this group indicates that the length distributions within the strike bins are skewed.

The second major group strikes north/south (000°). The directional distribution is much narrower as compared to the 130° group. The joints within this group appear to be longer as opposed to the next larger group that strikes at c. 045°. The 045° group shows a diffuse pattern. The skewness is not as large as the previous two groups resulting in a more circular color pattern.

The fourth group strikes c. east/west (090°) and essentially falls into two 10 degree bins. Field work shows the existence of this group, however an image analysis bias parallel to the pixels boundaries (0°, 90°) and diagonal to the pixel boundary (45°, 135°) is possible. A close look at Figure 17b confirms minor distortions in the population at exactly these directions.

Fracture Density Maps

In addition to orientation distribution, our dataset offers spatial information. The interpretation covers primarily the Natih Formation on the flank west of Wadi Nakhr and formations of the upper Kahmah Group (primarily the Shu’aiba Formation) on the east flank of the canyon (Figures 2 and 4), both with many clearly distinguishable beds.

The density of the interpreted fractures was analyzed by box counting with a kernel of 25 m × 25 m, calculating fracture length per unit area (m/m2) (Figure 18). The fracture density ranges from 0.0 m/m2 to 1 m/m2. The lowest densities (Figure 18) are found at the boundaries of the interpretation area and close to some of the wadis where no data is available. Variations in fracture density occur on the western flank within the Natih Formation (Figure 18, location 1).

Virtually the entire exposure of the Nahr Umr Formation shows low line densities (Figure 18, location 2) as compared to the adjacent Natih or Shu’aiba formations. Lower joint densities are also found close to the major fault zones in the area (Figure 18, location 3).

Fracture Density Maps for Strike Populations

The line density of the fractures was also calculated individually for the four different strike groups to check for the spatial distribution of the different sets within the field area. Figure 19 shows the distribution of the fracture sets with a relative color bar. From top to bottom the densities cover the groups 160°–020°, 020°–060°, 080°–100° and 115–145°, displayed by the gray sectors in the rose diagrams. The plots show that the joint distribution in the sub-populations have much larger variations in the distribution as opposed to the entire population in Figure 18. The densities within a strike population seem to form patches in the study area. Although the bulk joint density is a function of the lithology, the density of the individual sets is a more a function of location. The spatial variation forms domains in which individual sets predominate. These domains are in some cases bound by normal faults.

DISCUSSION Limitations and Possibilities of the Data Sets

Landsat and Quickbird data were used for mapping the fracture and fault system at 0.7 m ground resolution (Figure 11). The resulting dataset contains 145,000 fractures and 562 faults in an area of approximately 31 km2 of outcrops of predominantly the Natih and Shu’aiba formations. Ground-truthing showed that large aperture vein systems (a few centimeters aperture and more) and cemented faults in the predominantly dark grey carbonates are consistently imaged on satellite data (Figures 8 and 10), but the visibility of the data set is a function of color of the wall rock (Figure 9). This causes some bias in the data set and its full quantification requires more fieldwork. In the studied areas, nearly all joints visible on the satellite image were confirmed, thus validating the interpretations. It is estimated that more than 95% of the interpreted fractures and faults correspond to the correct features in outcrop.

Since the aperture of the fractures is smaller than the image resolution, internal structures (segmentation, en-échelon structures) were not resolved. A line drawn on the interpretation can be a bundle of thin veins in the field, or a single thick vein. Moreover, a large number of veins, which are too thin or insufficiently clustered to produce significant color contrast, may not be mapped (Figure 10). Exceptionally long veins with large apertures are difficult to distinguish from low-offset faults.

Keeping these limitations in mind, the results of the satellite image mapping are in excellent agreement with field observations. Integrating overprinting relationships and mm-scale structures from outcrop observations produced a well-defined, high-resolution model of faults and fractures over a very large area.

The present study offers insights into the detailed high-resolution anatomy of fault and fracture sets in massive carbonate sequences buried to several kilometers in depth. The stress history and fracture sets of these rocks differs from those in adjacent Interior Oman studied in the subsurface by Filbrandt et al. (2006). The latter study is primarily based on the interpretation of faults from 3-D seismic, whereas we incorporated a large number of fractures that formed prior to the faults.

Mechanical Layering

Fracture spacing is commonly related to the thickness of mechanical layers in systems where the fractures are not re-sealed (Bai and Pollard, 2000; Gillespie et al., 2001; Cooke et al., 2006; de Keijzer et al., 2007). In the present study, however, a simple correlation between vein spacing and layer thickness was not apparent. The absence of such a relationship is surprising because the thickness of the mechanical layers was sometimes smaller than one meter. This observation may be explained by a crack-seal model (as proposed in Holland and Urai, in review), where fractures are quite effectively sealed by vein cement before new fractures form. Vein density is therefore determined by the mechanical properties of the (re-sealed) host rock, and by the properties of the plumbing system in which the fluid pressures are built-up and locally released (Petit et al., 1999; Townend and Zoback, 2000). Although the spacing of fractures among small-scaled layers can show clear changes, the larger scale density (scale of the interpretation) is rather constant over large sections (Figure 18).

Fault Distribution

An anastomosing network of in total 562 faults is interpreted on the satellite image. The main strike direction is 100°. The lengths of interpreted faults show a power-law distribution down to 300 m. Faults shorter than this are under-represented with respect to the theoretical distribution. This may relate to a real property of the fault system, but it is more likely that faults with a length of less than 300 m are not recognized due to several reasons:

  • (1) The fault does not show significant throw to produce a visible juxtaposition.

  • (2) The fault does not show significant throw to produce a wide enough damage zone to be classified as such, and is interpreted as a fracture. For example, the low-offset systems can consist of a flexure, with a minor localized fault zone. The pattern of small cemented fractures distributed over such a wide zone may not produce a distinct enough signal for remote sensing interpretation.

  • (3) Some normal faults evolve by reactivation of pre-existing veins (Holland et al., 2009), and at small offset these again will be interpreted as fractures. This could explain the long veins striking 130° (Figure 16a).

Fracture Distribution

As discussed above, not all of the fracture population can be interpreted in the high-resolution Quickbird satellite image because the detection limit is a function of fracture aperture and optical material contrast (Figure 9). This effect is enhanced close to the major fault systems in the mapping area. The wide damage zones of some fault systems show dense networks of cemented fractures below the optical resolution (e.g. Figure 10). Although not resolved by the satellite, the fractures increase the overall brightness and inhibit the interpretation of veins that formed prior to the faulting (Figure 18a, b). A similar bias is caused by uncemented fracturing that locally increases erosion and weathering, or promotes the accumulation of rubble and soil. These processes influence the optical properties of the mapped rock surface and lead to a virtual lower fracture density.

From the population of more than 145,000 interpreted veins and joint lineaments a number of plots can be made to characterize the fractures according to their directional and spatial distributions (Figure 17). As in the field, abutting or curvature of interacting fractures are not present. The most prominent group of fractures strikes at 130° and shows a large directional variance (Figure 17). Confirmed by ground truthing the other groups strike 000°, 090° and 045°. Minor distortion due to aliasing effects are visible close to these peak directions (Figure 17).

The orientation distribution of four major groups was produced in an anticlockwise rotating stress field prior to the formation of the normal faults (based on field observations presented here the model (Holland et al., 2009). The fracture formation is characterized by the complete re-healing of fractures before the formation of a new generation implying a major rotation of permeability anisotropy in this system over geologic time. The anisotropy must have been very strong, as only one fracture set was open at any given time, so that connectivity by cross-cutting open fractures was not formed.

The bulk fracture density appears to be rather constant among different lithologies. Densities of the individual fracture sets (Figure 19) have in contrast a heterogeneous distribution. These patches with higher densities have a characteristic size of about one square kilometer. Some patches are bounded by normal faults, which juxtapose different layers with apparently slightly different fracture densities. In other areas, the boundaries of the fracture patches are not clearly related to layer transitions, and are interpreted to represent true lateral variations, similar to the formation of local patches in sandstone (Cruikshank and Aydin, 1995). Ongoing work is aimed at quantifying this pattern more accurately.

No direct indications were observed for the mechanical interaction of the joints, such as abutting or curving (e.g. Rawnsley et al., 1992, Rawnsley et al., 1998). Nevertheless, a weak mechanical interaction is inferred from the detailed study of some of these vein systems (Holland and Urai, in review), which affects the evolution of the normal faults and reactivation of the veins (Holland et al., 2009).


High-resolution mapping of remote sensing data on the southern flank of Jabal Shams in an area of 31 km2 identified more than 145,000 fractures and 560 faults (Figures 12 and 18). Ground truthing the satellite image interpretations with field observations confirms that the fractures are isolated veins or vein bundles (Figures 8 and 10) (Holland et al. 2009). The high data resolution proves however to be insufficient to resolve small-scale structures that allow observations on kinematic and overprinting features like splays, en-échelon systems, crack-seal evidence and cross-cutting relationships for which field observations are required.

Four generations of veins were observed (Figure 15) that lack a strong mechanical interaction as abutting and curving is not observed in the study area (Holland et al. 2009). This suggests rapid cementation and the restoration of the rock strength. These fracture generations show high and reasonably constant densities over the whole study area, in the Natih and Shu’aiba formations.

Fracture density of the individual sets within a major carbonate unit is laterally heterogeneous (Figure 19). This suggests phases of jointing in patches, followed by rapid healing of the fractures with calcite cement.

An anastomosing network of normal faults overprints the joints. Incipient faults locally develop by reactivation of pre-existing veins (Holland et al. 2009) reflected by a small discrepancy of the orientation of fault elements versus fault strands as well as extraordinary long veins. The presented study indicates that the spatial distribution of faults and fractures in such a high-pressure resealing system is difficult to characterize and predict as an effect of the variety of scale and mechanical interaction.


We would like to thank Petroleum Development Oman for providing us with the high-resolution DEM and aerial photograph of the field area. The support of Erik Wanningen on ArcGIS is appreciated. We also acknowledge the constructive comments of Joerg Mattner and two anonymous reviewers that helped to improve the paper. This project was funded by Shell International Exploration and Production B.V. The final design and drafting by GeoArabia Nestor Niño Buhay is appreciated.


Marc Holland studied Geology at the University of Hawaii at Manoa (UHM, USA) and the RWTH Aachen University in Germany, receiving his diploma with honours in 2004. In his subsequent PhD project at the department “Structural Geology, Tectonics and Geomechanics” at the RWTH Aachen University, he worked on several industry funded projects on multiscale fault and fracture properties, integrating field work, analogue modelling and satellite image interpretation. In April 2008 Marc started a position in the technology group of GeoMechanics International (GMI).


Nishank Saxena is currently working as a Geophysicist for Schlumberger. He graduated from the Indian Institute of Technology, Kharagpur, with a five-year integrated MSc in Exploration Geophysics. An Institute Silver Medalist and P.K Bhattachrya Award Winner in 2008, he also received the Schlumberger Indian University Handshake Scholarship and was selected to represent India at the International Petroleum Technology Conference held in 2007in Dubai. In 2006, he was funded by DAAD Germany to work on a Shell-sponsored industry project at the department of Structural Geology, Tectonics and Geomechanics of the RWTH Aachen University.


Janos L. Urai is Dean of the Faculty of Science and Head of the Department of Applied Geoscience at the German University of Technology GUtech in Muscat, Oman, and Professor of Structural Geology, Tectonics and Geomechanics at RWTH Aachen University in Germany. He was C&C Huygens Fellow of the Netherlands Organisation for Scientific Research and Senior Research Geologist at Shell Research. Janos’ current projects aim at understanding salt deformation and salt tectonics, evolution and resealing of fault and fracture systems in carbonates, the morphology of pore networks in fine-grained rocks, and the development of state-of-the-art geoscience teaching in the Middle East.