2.a. The Niagara escarpment

The Niagara Escarpment is an erosional cuesta, formed in part by differential erosion of alternating layers of shale, sandstone and carbonate-rich lithologies (Fig. 4; Hewitt, Reference Hewitt1971). Locally, the easily weathered Queenston and Rochester Formations, composed of shale, undercut the sub-vertical slopes of the more resistant Whirlpool sandstone and dolostone-rich Irondequoit Formation and Lockport Group, which form overhangs (Hewitt, Reference Hewitt1971; Brunton & Brintnell, Reference Brunton and Brintnell2020). This differential erosion yields an abrupt termination of strata at the escarpment face (Fig. 5; Chorley, et al. Reference Chorley, Dunn and Beckinsale1964: 336–7; Armstrong & Dodge, Reference Armstrong and Dodge2007; Hayakawa & Matsukura, Reference Hayakawa and Matsukura2010). Areas of the escarpment protected from basal undercutting have lower slope angles, whereas areas subject to active basal undercutting are characterized by steeper slopes; this highlights the importance of differential erosion processes in the development of extremely steep escarpment slopes in the Hamilton area (Fig. 5; Chorley et al. Reference Chorley, Dunn and Beckinsale1964: 336–7).

Fig. 5. A local north-facing outcrop of the Niagara Escarpment observed from the Chedoke Radial Trail (Fig. 2) showing geological formations (Fig. 4) including the Reynales Formation, the Irondequois Formation and the Rochester Formation. Formation boundaries shown by solid red lines.

2.b. Studied lithological units

The Niagara Escarpment exposes a series of distinctive lithological units grouped into a number of regionally extensive formations (Armstrong & Dodge, Reference Armstrong and Dodge2007; Brunton & Brintnell, Reference Brunton and Brintnell2020; Fig. 4). This study focuses on sub-vertical fractures within formations exposed in the upper part of the Niagara Escarpment, above and including the Reynales Formation (Figs 3c, 5a); these upper units were selected for study as they are easily accessible in the city of Hamilton and have made the greatest contributions to significant rockfall events (Van Dongen, Reference Van Dongen2016). The uppermost formations, including the dolostone-rich Gasport Formation and Ancaster Member of the Goat Island Formation, are included in the Lockport Group (Brunton & Brintnell, Reference Brunton and Brintnell2020), and determine the overall stability of the uppermost part of the cuesta in Hamilton. The Lockport Group has experienced significant rockfalls due to its position overlying the Rochester Formation, which is composed of shale and tends to undercut overlying units (Hewitt, Reference Hewitt1971; Brunton & Brintnell, Reference Brunton and Brintnell2020). The following are brief descriptions of the lithological and structural characteristics of the geological units examined in this study (Fig. 4).

The Reynales Formation is composed of medium beds of laminated fine-grained dolostone with shale interbeds (Hewitt, Reference Hewitt1971). The proportion of shale increases upwards in the formation through increasing thicknesses of shale interbeds and an increase in the amount of shaly dolostone (Hewitt, Reference Hewitt1971; Brunton & Brintnell, Reference Brunton and Brintnell2020). The overlying Irondequoit Formation consists of a single massive unit of fossiliferous and highly bioturbated vuggy dolostone that has a relatively consistent thickness of between 1 and 2 m in the Hamilton area (Hewitt, Reference Hewitt1971; Figs 4, 5). This unit is overlain by the Rochester Formation, a relatively thick unit of grey, laminated shale with minor interbeds of dolostone and limestone which become more prominent upwards in the stratigraphy (Hewitt, Reference Hewitt1971; Brunton & Brintnell, Reference Brunton and Brintnell2020). The Lockport Group forms the resistant cap rock of the Niagara Escarpment, which consists of two distinct units, the Gasport Formation and the Ancaster Member of the Goat Island Formation, which are prominent in the Hamilton area (Armstrong & Dodge, Reference Armstrong and Dodge2007). The Gasport Formation lies at the base of the Lockport Group, and is composed of medium to thickly bedded fossiliferous dolostone with minor chert nodules and discontinuous thin shale interbeds and lenses (Brunton & Brintnell, Reference Brunton and Brintnell2020). The Ancaster Member of the Goat Island Formation has a gradational contact with the Gasport, and is composed of highly irregular and fractured thin to medium beds of dolostone with abundant chert nodules, typically elongated parallel to the bedding (Hewitt, Reference Hewitt1971; Brunton & Brintnell, Reference Brunton and Brintnell2020).

2.c. Fractures in the Niagara Escarpment

Within the studied outcrops in the Hamilton area (Fig. 2), discontinuities are present in two primary forms: (1) horizontal bedding discontinuities, separating strata and bedding habits of lithological units, and (2) sub-vertical, primarily stratabound joints (Fig. 6; Gross & Engelder, Reference Gross and Engelder1991; Priebe et al. Reference Priebe, Brunton, Rudolph and Neville2019). Slope stability is intrinsically linked to the characteristics and extent of these discontinuities, as fracturing controls the formation of blocks of rock that are subject to rockfalls from exposed rock faces (Priebe et al. Reference Priebe, Brunton, Rudolph and Neville2019, Reference Priebe, Frape, Jackson, Rudolph and Brunton2021).

Fig. 6. Geological fractures in (a) the lower Gasport Formation; (b) the middle Gasport Formation; and (c) the Ancaster Member of the Goat Island Formation of the Lockport Group. Horizontal bedding planes (grey) and sub-vertical joints (red) are indicated on the photographs.

In addition, Priebe et al. (Reference Priebe, Brunton, Rudolph and Neville2019, Reference Priebe, Frape, Jackson, Rudolph and Brunton2021) examined hydraulic conductivity values associated with the Palaeozoic carbonate groundwater systems in southern Ontario, concluding that high hydraulic conductivity correlates with dense fracture networks and karstic features (Priebe et al. Reference Priebe, Brunton, Rudolph and Neville2019). Flow is especially notable at stratigraphic breaks in the rocks compared to carbonate rock matrices, highlighting the importance of bedding planes, which connect sub-vertical joints, in controlling fluid-flow properties (Priebe et al. Reference Priebe, Brunton, Rudolph and Neville2019).

3.a. Field data acquisition

Fractures have a quantifiable orientation and aperture (Schultz, Reference Schultz2019) which can be documented for a range of purposes including the understanding of hydrocarbon reservoirs (Massaro et al. Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018, Reference Massaro, Corradetti, Tramparulo, d’Assisi, Prinzi, Iannace, Parente, Invernizzi, Morsalnejad and Mazzoli2019; Giuffrida et al. Reference Giuffrida, La Bruna, Castelluccio, Panza, Rustichelli, Tondi, Giorgioni and Agosta2019, Reference Giuffrida, Agosta, Rustichelli, Panza, La Bruna, Eriksson, Torrieri and Giorgioni2020), structural inheritance and fluid-flow properties (Gross & Engelder, Reference Gross and Engelder1991; Eyles & Scheidegger, Reference Eyles and Scheidegger1995; Peace et al. Reference Peace, Dempsey, Schiffer, McCaffrey, Imber and Phethean2018; Samsu et al. Reference Samsu, Cruden, Micklethwaite, Grose and Vollgger2020). Semi-automated methods of fracture quantification have also been developed which have been shown to be useful in the acquisition of large data sets, in addition to limiting observer biases (Thiele et al. Reference Thiele, Grose, Samsu, Micklethwaite, Vollgger and Cruden2017; Lee et al. Reference Lee, Maclachlan and Eyles2020; Palamakumbura et al. Reference Palamakumbura, Krabbendam, Whitbread and Arnhardt2020; Samsu et al. Reference Samsu, Cruden, Micklethwaite, Grose and Vollgger2020). While digital quantification methods minimize observer biases, manual observations allow measurements to be taken on a case-by-case basis and are most appropriate for documenting complex fracture patterns (Thiele et al. Reference Thiele, Grose, Samsu, Micklethwaite, Vollgger and Cruden2017; Andrews et al. Reference Andrews, Roberts, Shipton, Bigi, Tartarello and Johnson2019). In this study, we use a manual method for fracture quantification (Mauldon et al. Reference Mauldon, Dunne and Rohrbaugh2001; Watkins et al. Reference Watkins, Bond, Healy and Butler2015; Andrews et al. Reference Andrews, Roberts, Shipton, Bigi, Tartarello and Johnson2019), rather than a digital method (Thiele et al. Reference Thiele, Grose, Samsu, Micklethwaite, Vollgger and Cruden2017; Lee et al. Reference Lee, Maclachlan and Eyles2020; Palamakumbura et al. Reference Palamakumbura, Krabbendam, Whitbread and Arnhardt2020; Samsu et al. Reference Samsu, Cruden, Micklethwaite, Grose and Vollgger2020), as it allows more detailed examination of the relationship between fracture strike direction and the orientation of the rock face (Thiele et al. Reference Thiele, Grose, Samsu, Micklethwaite, Vollgger and Cruden2017; Andrews et al. Reference Andrews, Roberts, Shipton, Bigi, Tartarello and Johnson2019), a relationship that is highly variable in the study location. To minimize observational biases in measuring fracture orientation, ≥30 measurements by a single observer were made in each survey (Mauldon & Dershowitz, Reference Mauldon and Dershowitz2000; Andrews et al. Reference Andrews, Roberts, Shipton, Bigi, Tartarello and Johnson2019).

Manual fracture quantification was completed using both scanline and areal (also called ‘window’ (Watkins et al. Reference Watkins, Bond, Healy and Butler2015)) surveys, both of which are standard methods used in similar studies (Gross & Engelder, Reference Gross and Engelder1991; Watkins et al. Reference Watkins, Bond, Healy and Butler2015; Massaro et al. Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018, Reference Massaro, Corradetti, Tramparulo, d’Assisi, Prinzi, Iannace, Parente, Invernizzi, Morsalnejad and Mazzoli2019; Andrews et al. Reference Andrews, Roberts, Shipton, Bigi, Tartarello and Johnson2019; Giuffrida et al. Reference Giuffrida, La Bruna, Castelluccio, Panza, Rustichelli, Tondi, Giorgioni and Agosta2019, Reference Giuffrida, Agosta, Rustichelli, Panza, La Bruna, Eriksson, Torrieri and Giorgioni2020). Surveys of both types were completed initially to capture vertical and horizontal variations in fracture distribution. Strike and dip of fractures were measured using a compass–clinometer, while aperture was recorded using a digital caliper in the centre of each fracture within the survey areas, or where the scanline intersected the fracture. Bedding thickness was also recorded.

Scanline fracture quantification involved establishing a level datum on an outcrop and recording fracture attributes (strike, dip and aperture) for all fractures intersecting the line (Gross & Engelder, Reference Gross and Engelder1991; Watkins et al. Reference Watkins, Bond, Healy and Butler2015). Scanlines were selected to ideally be 10 m in length, although shorter scanlines were necessary for the Ancaster Member due to the intensity of fractures in these strata for feasibility. When a greater length of outcrop was available, longer scanlines were surveyed. The scanline methodology allows fracture distribution to be measured laterally across an outcrop, but is subject to orientation bias, as the scanline does not intersect sub-horizontal bedding planes, which were also quantified as discontinuities in this study (Watkins et al. Reference Watkins, Bond, Healy and Butler2015).

Area surveys were completed by marking an area (window) on the outcrop of interest and then quantifying each fracture within the area, together with those intersecting its borders, to capture multiple stratigraphic layers (Watkins et al. Reference Watkins, Bond, Healy and Butler2015). This method should involve consistent area size between study sites, and is associated with minimum orientation bias (Watkins et al. Reference Watkins, Bond, Healy and Butler2015). Areas were selected to be ideally 1 m × 1 m in dimensions for most surveys, although some surveys of the Ancaster Member of the Goat Island Formation were slightly smaller (0.5 to 0.75 m by 0.5 to 0.75 m).

Outcrops of interest were located along the exposed, high-angle cliff face of the Niagara Escarpment in Hamilton (Fig. 5) with the help of satellite imagery. The following criteria were used in selecting study sites: (1) a minimum of 10 m in length and 1 m in height, (2) inclusion of at least one of the lithological units identified for the study, and (3) safe accessibility for the observer.

Fracture characteristics (strike, dip and aperture) were documented at 11 sites across Hamilton (Fig. 2) using both scanline and area surveys. Data from seven sites (Rock Chapel, Sydenham Cut, the Sherman Falls Trail, the Chedoke Radial Trail, the Bruce Trail, Jolly Cut and Mountain Brow Trail; Fig. 2) were selected and used for further analyses due to their visual lithological homogeneity. These sites were used to produce a representative model of fracture characteristics for the uppermost units of the Lockport Group. Data from other study locations (Fig. 2) were excluded from statistical analysis and DFN modelling because they did not pertain to the Lockport Group (e.g. Smokey Hollow), did not reflect primary local lithologies in Hamilton (e.g. Rattlesnake Point and Mt Nemo) and/or lacked sufficient measurements (e.g. Devil’s Punchbowl). However, the data from these excluded sites were used to document overall fracture and rockfall characteristics, and regional context.

Subsequent analyses to define and model fracture sets were conducted separately for the Gasport Formation and overlying Ancaster Member of the Goat Island Formation as these units continuously crop out at each study site and are associated with most rockfall events. The collected data were divided into four groups based on the method of acquisition and sampled stratigraphic unit as follows:

1. (1) Area survey data from the Gasport Formation, comprising 11 surveys (n = 313);

2. (2) Scanline survey data from the Gasport Formation, comprising 11 surveys (n = 404);

3. (3) Area survey data from the Ancaster Member, comprising 9 surveys (n = 310);

4. (4) Scanline survey data from the Ancaster Member, comprising 7 surveys (n = 249).

These data groups are analysed and discussed in detail in the following subsections. Statistical analyses were conducted on each individual data set listed above, while DFN modelling combined the scanline survey data sets and the area survey data sets to model each unit (Fig. 4).

3.b. Statistical analysis and data preparation

Data collected from field surveys of different study sites along the escarpment (Fig. 2) were analysed to identify statistical patterns (Fig. 7). The data were prepared for analysis by defining sets of sub-parallel fractures from the collected field data based on observed fracture strike. In order to consider such linear features in a three-dimensional (3-D) vector space, trend and plunge values of poles for each measurement were entered into Orient^TM spherical data analysis software (version 3.12.0; Vollmer, Reference Vollmer2015) in which analysis was conducted. Statistical analysis of the fractures involved calculating the mean fracture aperture, identifying strike orientation clusters and calculating the respective maximum and minimum eigenvectors and associated confidence intervals for each fracture set in each data group (Vollmer, Reference Vollmer1990, Reference Vollmer1995). Apertures measured as <1 mm were assigned representative values of 0.5 mm in the statistical analyses. In addition, some fractures had apertures that were not measurable, and thus were excluded from average calculations. As most beds identified in this study had sub-horizontal orientation, fractures parallel to beds were treated as horizontal discontinuities and excluded from the generated spherical projections (Fig. 8). In Orient, cluster analysis and cones of confidence were used to define sub-vertical joint sets (Vollmer, Reference Vollmer1990, Reference Vollmer1995).

Fig. 7. PGR (point-girdle-random) plots of the N–S and E–W joint sets for (a) the Gasport Formation area surveys; (b) the Gasport Formation scanline surveys; (c) the Ancaster Member area surveys; and (d) the Ancaster Member scanline surveys. All joint sets plot closest to the P (point) index, indicating that they represent a point distribution.

Fig. 8. Fracture orientations plotted as poles from field data shown on equal-area spherical projections produced in Orient (Vollmer, Reference Vollmer2015) for (a) the Gasport Formation area surveys; (b) the Gasport Formation scanline surveys; (c) the Ancaster Member area surveys; and (d) the Ancaster scanline surveys. The N–S and E–W sets were defined using the cluster analysis function in Orient with the respective minimum and maximum eigenvectors representing the locations of minimum and maximum pole density, respectively (Vollmer, Reference Vollmer1990).

3.c. Discrete fracture network (DFN) modelling in MOVE^TM

The escarpment strata and fractures documented in the field were modelled in MOVE^TM, a discrete fracture network (DFN) modelling program successfully applied in many previous studies (Panza et al. Reference Panza, Agosta, Zambrano, Tondi, Prosser, Giorgioni and Janiseck2015, Reference Panza, Agosta, Rustichelli, Zambrano, Tondi, Prosser, Giorgioni and Janiseck2016; Massaro et al. Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018, Reference Massaro, Corradetti, Tramparulo, d’Assisi, Prinzi, Iannace, Parente, Invernizzi, Morsalnejad and Mazzoli2019; Giuffrida et al. Reference Giuffrida, La Bruna, Castelluccio, Panza, Rustichelli, Tondi, Giorgioni and Agosta2019, Reference Giuffrida, Agosta, Rustichelli, Panza, La Bruna, Eriksson, Torrieri and Giorgioni2020; Rimando & Peace, Reference Rimando and Peace2021). DFNs are computer-based, stochastic models that represent the geometrical characteristics of fracture networks, from which porosity and hydraulic conductivity can be estimated (Panza et al. Reference Panza, Agosta, Zambrano, Tondi, Prosser, Giorgioni and Janiseck2015, Reference Panza, Agosta, Rustichelli, Zambrano, Tondi, Prosser, Giorgioni and Janiseck2016; Lei, Latham & Tsang, Reference Lei, Latham and Tsang2017; Massaro et al. Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018, Reference Massaro, Corradetti, Tramparulo, d’Assisi, Prinzi, Iannace, Parente, Invernizzi, Morsalnejad and Mazzoli2019; Giuffrida et al. Reference Giuffrida, La Bruna, Castelluccio, Panza, Rustichelli, Tondi, Giorgioni and Agosta2019, Reference Giuffrida, Agosta, Rustichelli, Panza, La Bruna, Eriksson, Torrieri and Giorgioni2020).

Massaro et al. (Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018) used the MOVE^TM software to model a series of DFNs in dolostones in Italy, linking the fracture properties of individual beds to determine resultant characteristics in the rock mass. A similar approach was used herein to model strata of the Niagara Escarpment based on data collected in the field. It should be noted that the model size and number of modelled beds were limited by the available computational power.

In this study, beds from the Ancaster Member of the Goat Island Formation were modelled, forming a total model domain volume of c. 5 m × 5 m × 5 m (Fig. 9). We modelled horizontal bedding planes as DFNs since bedding discontinuities contribute to fluid flow. As the regional bedding in the study area is very close to horizontal according to our observations and previous work (dipping ∼0.3° towards the south; Formenti et al. Reference Formenti, Peace, Waldron, Eyles and Lee2021), bedding discontinuities were considered horizontal in both models. Two separate models were produced, each comprised of beds from both the Gasport Formation and the Ancaster Member, one using data from area surveys and a second using data from scanline surveys. This was done to examine the differences, if any, between area and scanline data and to test the reproducibility of the modelling process. The model-building process is summarized in Fig. 9.

Fig. 9. Workflow from field data acquisition through modelling in MOVE^TM created to model the sub-vertical joints and bedding planes of the Niagara Escarpment in the Hamilton area.

To produce DFNs in MOVE^TM, input parameters were derived from field data (Tables 1–3). The models were built by defining sedimentary beds, separated by horizontal, manually placed fractures representing bedding discontinuities, and assigning each bed two DFNs representative of the two sub-vertical joint sets observed in the field (N–S and E–W striking sets; Fig. 10). Parameters were defined using mean values from field measurements and their associated confidence intervals (Table 1). Bed thicknesses, fracture apertures and the strike of joint sets were derived from fracture set characteristics selected for each bed by generating semi-random parameters from a Gaussian distribution of parameters measured in the field. The values were used to define bed thicknesses, apertures and strikes of the joint sets.

Table 1. Average orientations and confidence intervals of fracture sets, apertures and bedding thicknesses used to derive model parameters. These values were calculated from field data

Table 2. Final parameters for the area survey model derived from field data and the literature. These parameters were used in the fracture modelling module in MOVE^TM to generate unique pairs of DFNs (composed of a N–S set and an E–W set) which were repeated three times for a total of nine beds in each formation. An additional DFN was created to represent the beds

Table 3. Final parameters for the scanline survey model derived from field data and the literature. These parameters were used in the fracture modelling module in MOVE^TM to generate unique pairs of DFNs (composed of a N–S set and an E–W set) which were repeated three times for a total of nine beds in each formation. An additional DFN was created to represent the beds

Fig. 10. Two different views of 3D DFN models produced in MOVE^TM for the Lockport Group from area survey data (a, b) and scanline survey data (c, d). The DFNs for the Lockport Group comprise two units in Hamilton: the Gasport Formation joints (blue and teal) and the Ancaster Member joints (red and orange). Each unit comprises nine beds for each Member containing two stratabound joint sets generated with the fracture modelling module in MOVE^TM. Both models have a width and length of 5 m, while the height varies based on the bedding thicknesses, totalling to 5 m for the area survey model and 4.8 m for the scanline survey model. The bedding planes (grey) are modelled as a DFN representative of their role as a real-world geological discontinuity.

The fracture aspect ratio applied in the DFNs was selected based on values from within the literature (Odling, Reference Odling1997; Olson, Reference Olson2003; Massaro et al. Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018, Reference Massaro, Corradetti, Tramparulo, d’Assisi, Prinzi, Iannace, Parente, Invernizzi, Morsalnejad and Mazzoli2019). To facilitate the construction of the DFNs, the dip of the bedding-confined fractures was set constant at 90° as the majority (87 %) of joints measured from both the Gasport and Ancaster were sub-vertical (dip above 75°).

Once the DFNs based on the field data obtained by area and scanline acquisition methods were constructed, calculation of the porosity and permeability resulting from the synthetic fractures across the entire domain volume was undertaken using a resolution defined by 5 cm × 5 cm × 5 cm cubic cells. Giuffrida et al. (Reference Giuffrida, La Bruna, Castelluccio, Panza, Rustichelli, Tondi, Giorgioni and Agosta2019) successfully modelled stratabound fractures with DFNs in MOVE^TM using a 20 cm × 20 cm × 20 cm cell size for a 7 m × 7 m × 7 m model volume. This suggests that the 5 cm × 5 cm × 5 cm cell size utilized herein is more than sufficient to quantify and predict the hydraulic properties of interest. Two DFN models (Fig. 10) were produced, in which all fluid flow was assumed to occur via geological fractures; no further lithological parameters are included in these models as fractures are considered to have the dominant influence on fluid-flow properties (Allan & Sun, Reference Allan and Sun2003; Massaro et al. Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018, Reference Massaro, Corradetti, Tramparulo, d’Assisi, Prinzi, Iannace, Parente, Invernizzi, Morsalnejad and Mazzoli2019; Giuffrida et al. Reference Giuffrida, La Bruna, Castelluccio, Panza, Rustichelli, Tondi, Giorgioni and Agosta2019, Reference Giuffrida, Agosta, Rustichelli, Panza, La Bruna, Eriksson, Torrieri and Giorgioni2020). Modelled permeability is considered to be representative of migration pathways for meteoritic water in the fracture networks.

4.a. Field observations

Observations from the field sites show two distinct sub-vertical fracture sets with different orientations, along with bed-parallel sub-horizontal discontinuities, making up three fracture sets in total (Figs 2, 4, 5). Bed thickness varies in each lithological unit. All sites display sub-vertical, primarily stratabound joints, formed without shear offset, representing mode I fractures. While most of these fractures abut horizontal beds, some throughgoing fractures were observed. Visually, chert nodules, which typically form horizontal lenses of c. 10 cm in diameter, appear to facilitate the propagation and deflection of irregular joint patterns, often changing the orientation of the joints, which run along the edges of the chert nodules. Numerous micro-fractures were observed within the chert nodules but were not included in this analysis. Chert nodules were observed to increase in frequency upwards through the studied sections, with relatively few in the Gasport Formation (Fig. 6a, b) and many in the Ancaster Member of the Goat Island Formation (Fig. 6c). The resulting joint deflections made the measurement of fracture orientation in the Ancaster Member particularly challenging.

Shale interbeds, such as those separating dolostone beds in the Reynales Formation (Figs 2c, 3) at the Chedoke Radial Trail (Fig. 2), terminate sub-vertical joints. Joints cross-cutting sub-horizontal bedding planes, observed in other formations where bedding discontinuities occur between two carbonate lithologies, are not observed in shales, suggesting that deformation in shale was mesoscopically ductile, and therefore did not develop visible fractures (Sone & Zoback, Reference Sone and Zoback2014; Massaro et al. Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018).

4.b. Statistical analysis of joint orientations

Statistical analysis of joint orientation data collected in the field demonstrated that the sub-vertical joints make up two distinct sets striking roughly N–S and E–W respectively (Figs 7, 8). Both sets were not visible simultaneously at individual field sites as each site consisted of a single two-dimensional rock face (Fig. 6). One exception was at a curved section of the rock face on the Chedoke Radial Trail (Fig. 2) which exposed the two sub-vertical joint sets later identified through analysis of field measurements. Plotting data collected at sites with different aspects on equal-area spherical projections clearly shows the two fracture sets with two dense clusters of poles (Fig. 8).

The equal-area spherical projections (Fig. 8) paired with point-girdle-random (PGR) distribution plots (Fig. 7) show that both joint sets form point distributions (Vollmer Reference Vollmer2015), corresponding to clusters of similarly striking joints, which guided cluster analysis. Two clusters were identified and summarized by calculating the maximum and minimum eigenvectors for each cluster (ϵ ₁ and ϵ ₃). The average strikes from ϵ ₁ for each joint set in each data group are presented in Table 1 with the respective confidence limits derived from the elliptical cones of confidence.

By definition, ϵ ₃ plots at the point of lowest fracture density and is oriented 90° from the point of highest density (Vollmer, Reference Vollmer1990). Orientations of ϵ ₃ for each of the two sub-vertical joint sets (Fig. 8) lie near the origin of the projection, corresponding to planes with dip ≈ 0°. If all sub-vertical joints making up the fracture sets had a dip of 90°, as assumed in the models, the ϵ ₃ values would be representative of the bedding planes, being oriented 90° from the joints. The data do not consist of exclusively vertical joints; however, since the ϵ ₃ values plot as poles representative of planes with a dip of 0°, our inference that the high-angle joints are bounded by horizontal or sub-horizontal discontinuities is supported mathematically.

4.c. Discrete fracture network modelling results

The two DFNs produced in this study, generated from two different sets of data (i.e. scanline- and area-survey-derived data), are shown in Figure 10, and the resultant corresponding porosity and meteoritic water migration pathway models are shown in Figure 11. The models show that similar patterns are displayed by either set of input data (from area- or scanline-derived measurements), and that a clear distinction between the hydraulic and fracture parameters of the Gasport Formation and Ancaster Member can be identified. The porosity and fluid migration pathway models in both the area and scanline DFNs show that the longer and deeper fractures present in the Gasport Formation compared to the Ancaster Member lead to a higher degree of horizontal connectivity in the Gasport.

Fig. 11. The permeability (red scale), calculated in darcys (Oda, Reference Oda1985), and porosity (blue scale), calculated as the ratio of fracture volume to cell volume, of each model are displayed with a 125 cm³ cell size. Three faces of the 3D model are shown: the north face, the bottom face (showing joints in the Gasport Formation) and the top face (showing joints in the Ancaster Member). Note that the permeability and porosity scales are different for the area survey model (a, b) and the scanline survey model (c, d) due to variations in aperture. The larger fracture apertures in the scanline model drastically increased the permeability and porosity values, rendering the need for a different scale. Cells with calculated values of zero (inactive cells) are excluded from the modelling for viewing ease. These areas are representative of solid blocks of rocks containing no joints.

In addition, all three modelled fracture sets show a high degree of connectivity in three dimensions, in both DFNs created from the scanline and area data (Fig. 11). The three fracture sets are all oriented approximately perpendicular to one another, forming discrete rectangular blocks. This is consistent with the field observations of displaced blocks from rockfalls (Fig. 3). The DFNs illustrate how smaller blocks are formed in the Ancaster Member because of the abundance of shorter, more densely clustered joints, and thinner beds, in comparison to blocks released from the Gasport Formation (Fig. 6).

Fluid-flow magnitude (Fig. 11) is proportional to the aperture of the joint sets, and thus the scanline model yields substantially greater porosity and permeability values due to larger apertures provided as input parameters (Table 3).

5.a. Escarpment rockfalls and implications for geohazards

This study identifies three sets of discontinuities in the studied sites along the Niagara Escarpment, striking at c. 90° to one another, which facilitate the release of rectangular blocks from the exposed rock face (Fig. 3). From the results of this study, we hypothesize that geological fractures present a primary control on the size of blocks included in rockfalls generated along the face of the escarpment. The size of blocks corresponds to fracture density and bedding thickness, block size being proportional to stratigraphic bed thickness and inversely proportional to fracture density. Quantifying block size in the field was beyond the scope of this study. However, the three-dimensional models produced (Fig. 9) illustrate rock compartmentalization by fractures similar to that described by Massaro et al. (Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018, Reference Massaro, Corradetti, Tramparulo, d’Assisi, Prinzi, Iannace, Parente, Invernizzi, Morsalnejad and Mazzoli2019). The modelled block size can be compared with qualitative observations (Fig. 3) from rockfall events, yielding a practical application of the workflow presented herein.

The DFN results demonstrate how the Gasport Formation generally produces fewer but larger blocks compared to the overlying Ancaster Member of the Goat Island Formation, due to differences in bedding thickness (Formenti et al. Reference Formenti, Peace, Waldron, Eyles and Lee2021) and fracture density. The larger blocks released from the Gasport Formation present the potential for greater damage and risk to urban infrastructure around Hamilton (Fig. 2). However, the Ancaster Member may have a higher frequency of dispersal of lower-mass blocks. Both of these units of the Lockport Group are commonly found in close proximity to urban infrastructure around Hamilton, and their fracture characteristics thus require consideration from a geohazard risk and mitigation perspective. The workflow used in this study could practically be applied in the characterization (in terms of block size), and management (in terms of risk analysis), of rockfall events. It may be used to enhance understanding of the 3-D nature of fracture networks in local rock outcrops and the creation of different block volumes via the intersection of the three perpendicular fracture sets. However, it should be noted that variations in lithology will affect fracture and rock block geometry such that the model presented here is only applicable in the present study area.

5.b. Fluid flow through rocks of the Niagara Escarpment

Field observations show that sub-vertical joints are consistently connected via bedding planes, thereby yielding high porosity and permeability. Permeability was calculated in MOVE^TM, but since the models are based on DFN characteristics and the nature of the rock matrix was not considered, the results more accurately describe fluid migration pathways through the rock body. Erosion is commonly facilitated by fluid flow (Collins & Stock, Reference Collins and Stock2016; Lebedeva & Brantley, Reference Lebedeva and Brantley2017; Gage et al. Reference Gage, Eyles and Lee2021) and this is especially true for the carbonate strata of the Niagara Escarpment where interbedded shales with low permeability force fluids into the carbonates (Priebe et al. Reference Priebe, Brunton, Rudolph and Neville2019, Reference Priebe, Frape, Jackson, Rudolph and Brunton2021). The fracture apertures used in the Gasport scanline survey model (Table 3) were larger than those used in the area survey model (Table 2). This was unexpected, but possible given the input parameters and semi-random parameter selection using a Gaussian distribution. This yielded far higher porosity and permeability values for the scanline survey model of the Gasport Formation compared to the area survey model, so much so that the resolution of the model domain itself appears coarser due to the greater range of values between cells with and without fractures. Although the large differences in calculated fluid-flow values cause visual presentation challenges, this demonstrates a worthwhile point. The aperture of joints has a significant effect on fracture-induced fluid-flow properties. Since a wide range of apertures were observed at the studied sites on the Niagara Escarpment (<1 mm to >40 mm wide), fluid-flow properties should also be considered to be highly variable. This study also supports previous works (Allan & Sun, Reference Allan and Sun2003; Massaro et al. Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018, Reference Massaro, Corradetti, Tramparulo, d’Assisi, Prinzi, Iannace, Parente, Invernizzi, Morsalnejad and Mazzoli2019; Giuffrida et al. Reference Giuffrida, La Bruna, Castelluccio, Panza, Rustichelli, Tondi, Giorgioni and Agosta2019, Reference Giuffrida, Agosta, Rustichelli, Panza, La Bruna, Eriksson, Torrieri and Giorgioni2020) that suggest fracture aperture (Fig. 11) exerts a first-order control on fluid-flow properties.

The variability in fracture aperture recorded in this study, paired with the associated variability in fluid-flow properties, suggests that erosion potential is highly variable along the escarpment. Examination of fractures and block formation from individual sites would be a useful means of predicting the highly variable geohazards associated with rockfalls around the City of Hamilton.

5.c. Fracture formation in the Niagara Escarpment

The results of this study, in conjunction with previous works, have allowed us to consider the characteristics of, and mechanisms governing the origin of, the sub-vertical joint sets. More research is needed to confirm the origin of these joints, although the following discussion introduces some of the possible factors contributing to their formation.

Along the Niagara Escarpment, fractures may be related to: (1) horizontal regional tectonic stress (Heidbach et al. Reference Heidbach, Rajabi, Reiter and Ziegler2016 and references therein; Gross & Engelder, Reference Gross and Engelder1991), (2) structural inheritance from underlying bedrock joints (Eyles & Scheidegger, Reference Eyles and Scheidegger1995; Eyles et al. Reference Eyles, Arnaud, Scheidegger and Eyles1997), (3) non-tectonic stresses due to post-glacial rebound from the retreat of the Laurentide Ice Sheet (Quinlan, Reference Quinlan1984; Adams, Reference Adams1989), and (4) post-glacial freeze–thaw action (Moss & Nickling, Reference Moss and Nickling1980; Fahey & Lefebure, Reference Fahey and Lefebure1988). These hypotheses all relate to the sub-vertical joint sets, which abut at bedding planes formed during deposition and are discussed in the following subsections. Throughgoing fractures were found to be less common and are not the focus of the present study.

5.d. Field data acquisition: area vs scanline methods

Recording data from both area and scanline surveys was initially done to capture the horizontal and vertical variation in fracture attributes at the selected study sites (Watkins et al. Reference Watkins, Bond, Healy and Butler2015). The data and resulting models (Table 1; Fig. 9) were found to be nearly identical for both survey types, showing only minor differences in average joint orientation and aperture values (Fig. 8). The similarity of data collected in each survey type, the redundancy of measuring sub-horizontal bedding planes, and the overall ease of scanline surveys compared to area surveys, suggest that scanline surveys provide a more efficient methodology for the documentation of fracture characteristics in this area, specifically the sub-vertical joint sets.

The selected study sites spanned almost 50 km along the escarpment in Hamilton (Fig. 2), which proved to be valuable, as the variability in site location and aspect allowed identification of two distinct sub-vertical joint sets. Most individual sites (except a curved portion of the Chedoke Radial Trail; Fig. 2) consisted of two-dimensional exposures, revealing only one sub-vertical joint set, striking either N–S or E–W, and only through examination of multiple exposures with different orientations (Fig. 2) could a three-dimensional understanding of joint orientation be gained. Moreover, the consistency of the documented joint orientations, shown by dense clusters of poles on equal-area spherical projections (Fig. 8), indicates that these trends are observed across the study area and that the models produced (Fig. 9) are representative of the study area.

5.e. Limitations of geological fracture modelling workflow

Several limitations of using DFNs to model fracture sets in the Lockport Group became apparent during this study. Most notably, modelling stratabound joints and beds that acted equally as discontinuities was not an inbuilt function within the MOVE^TM workflow. Facilitating stratabound joints in the DFNs required editing the stochastic nature of the DFNs by manually placing individual sub-vertical joints such that they were situated at the centre of each bed, the boundaries of which were defined by manually placed horizontal bedding planes. In addition, variations in dip angle would not allow all joints to be absolutely bounded by the beds, meaning that each joint set in the DFN was given a dip of 90°. In other words, within the MOVE^TM workflow it was not possible to both vary the dip and terminate joints at bedding planes. However, the majority of recorded joints are sub-vertical (87 % have dips >75°), and thus, in the context of the field observations, these simplifications are reasonable approximations of a complex real-world system.

Massaro et al. (Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018, Reference Massaro, Corradetti, Tramparulo, d’Assisi, Prinzi, Iannace, Parente, Invernizzi, Morsalnejad and Mazzoli2019) present similar DFN models for stratabound fractures in carbonate reservoirs for the purposes of estimating fracture connectivity, fracture-induced porosity and fluid-migration pathways. The results of the current study support their conclusion that fractures have a significant effect on fluid-flow properties such that those of matrix porosity can effectively be ignored due to its low contribution to fluid flow when fractures are present (Massaro et al. Reference Massaro, Corradetti, Vinci, Tavani, Iannace, Parente and Mazzoli2018, Reference Massaro, Corradetti, Tramparulo, d’Assisi, Prinzi, Iannace, Parente, Invernizzi, Morsalnejad and Mazzoli2019).