Abutting faults: a case study of the evolution of strain at Courthouse branch point, Moab Fault, Utah

. Slip planes and slip directions of subsequent generations of faults were measured and analysed in the interaction damage zone of two abutting faults in porous sandstones, to understand the paleostress/paleostrain 10 evolution. The Courthouse branch point of the Moab Fault in SE Utah (USA) is a much-studied, spectacular outcrop of two abutting faults, located in the footwall block of the main fault and in the hanging wall block of the abutting fault. The abutting fault is synthetic to the main fault. The outcrop shows a wide range of deformation structures and fault related diagenesis: striated slip planes, deformation bands, veins, Liesegang bands and copper-rich 15 mineralization. By combining our own measurements with published data on the relative age of these structures, we classified the data in four sets. Using the Numeric Dynamic Analysis (NDA) to calculate the orientation of the kinematic axes we found three different paleo-extension directions in the four sets, recording the evolution of stress/strain axes during the abutting process. 20 The first phase of deformation is regional extension in the NE-SW direction. As the second fault approached the main fault from its footwall side and the two faults start to become kinematically linked, the extension direction changed so that the overall extension became perpendicular to the approaching fault (NW-SE). Finally, the extension direction changed back to perpendicular to the first segment (NE-SW), when the two faults become geometrically-linked and regional extension became

Although at the regional scale the stress field and kinematics of extensional tectonics are well established (Anderson, 1951), less is known of the interaction of faults at the local scale (Walsh and Watterson, 1991;Peacock, 2002;Nicol et al., 2010;Duffy et al., 2015;Delogkos et al., 2017;Peacock et al., 2017a), though this is where damage zones and the associated flow of geofluids are most important Rotevatn et al., 2009a;Yale, 2003;Eichhubl et al., 2009;Duffy, 2015). Peacock et al. (2017a, b) suggested that flow was already affected in the "interaction damage zone" as the faults approached each other before becoming kinematically linked. The partitioning of strain at this scale can compartmentalize reservoirs or provide conduits for fluid flow, depending on the mode of failure (Aydin and Johnson, 1978;Ingram and Urai, 1999;Johansen et al., 2005;Rotevatn et al., 2007;Bozkurt Çiftçi and Bozkurt, 2007;Urai et al., 2008;Eichhubl et al., 2009;Ferrill et al., 2009;Rotevatn et al., 2009b;van Gent et al., 2010a, Vrolijk et al., 2016. In this study we focus on the evolution in the stress and strain field in a single-tip fault interaction , in which one fault grows towards another. Peacock et al. (2017a) called this an "abutting fault interaction", where the synchronously growing abutting fault is synthetic to the main fault until the two faults become hard linked and slip together, without crosscutting.
We use the orientation and slip direction of subsequent generations of minor faults as input for palaeostress and palaeostrain analysis. We document changes in the extension direction and related local stress field over time, which we correlate with the transition from isolated faults, (stage i), through interacting, kinematically linked faults, (stage ii), to geometrically linked abutting faults, (stage iii). Understanding the stress and strain evolution of such a structure is a key step towards understanding the fluid pathways in the subsurface. Peacock et al. (2017a) described the interaction of abutting faults as a continuous process in a regional stress field with three evolutionary stages. First the faults are far from each other and do not interact. The faults start to interact as they approach each other and become kinematically linked. Interaction can be indicated by deformation structures in the acute bisector. In the final stage the approaching fault abuts to the main fault, and they become geometrically linked along a branch line. From this point on, both faults slip parallel to this line; the faults become hard linked.

Regional geology
The study area is located on the Colorado Plateau near Moab, Utah (USA), where evaporites of the Paradox Formation were deposited on basement in the Carboniferous (Doelling, 2002;Nuccio and Condon, 1996;Foxford et al., 1996). The Lower Jurassic aeolian Navajo sandstone is overlain by the Middle Jurassic Dewey Bridge Member, a calcareous, medium-to fine-grained sandstone. The Slick Rock Member is a part of the Middle Jurassic Entrada Fm. This is a reddish to orange, fine-grained and cross-bedded aeolian sandstone, which is overlain by the greyish-yellow sandstones of the Moab Member (Curtis Fm.). Directly overlying the Moab Member is the Late Jurassic terrestrial Morrison formation, consisting of siltstones, fine-grained sandstones, nodular limestones and conglomerates, with occasional dinosaur bones.
In the Jurassic, major salt-cored anticlines (Fig. 1a) developed in the Moab area (Doelling, 2002;Foxford et al., 1996), with numerous salt-related faults and salt welds in the overburden. The salt-related Moab Fault is described in detail in Foxford et al. (1998) and in many other publications. With a maximum throw of ∼ 950 m near the entrance of the Arches National Park and a length of around 45 km (Foxford et al., 1996(Foxford et al., , 1998, the Moab Fault likely initiated as salt weld but shows clear characteristics of brittle normal faulting in the study area. The fault offsets rocks from Late Carboniferous to Cretaceous. It was active from the Triassic to mid-Jurassic during salt tectonics and again from mid-Cretaceous to Early Tertiary when the salt dissolved and the salt structures collapsed (Foxford et al., 1996(Foxford et al., , 1998. Maximum burial depths of 2-3 km and temperatures of around 80 • C were reached in the Entrada sandstone (Nuccio and Condon, 1996;Fossen, 2010).
In aerial and satellite photos the exposed sandstones of the Colorado Plateau display a dense network of subvertical joints in many locations, which formed due to the uplift of the Colorado Plateau in the Tertiary (Cruikshank and Aydin, 1995;Foxford et al., 1996;. Similar structures are visible in the footwall of the Moab Fault in the Moab Member. They postdate the activity of the Moab Fault (Cruikshank and Aydin, 1995;Fossen, 2010).

Outcrop description
We focus here on an outcrop near the northern termination of the Moab Fault ( Fig. 1) with a number of hard-linked strands ( Fig. 1), which likely developed as breached relay ramps (Foxford et al., 1998;Fossen and Rotevatn, 2016). This is opposed to a secondary fault splaying from the main fault. This implies that the growing abutting fault (Segment B, Figs. 1b and 2) arrived at the junction after the main Moab Fault (Segment A). A third segment (Segment C) is exposed further to the west. At the intersection of Segments A and B, a spectacular outcrop of abutting faults is found, Courthouse Junction branch point (CHJ) Figs. 1, 2 and 3;Foxford et al., 1996;Davatzes and Aydin, 2003;Davatzes et al., 2005;Eichhubl et al., 2009;Johansen et al., 2005). The two faults abut at an angle of 64 • . The intersecting damage zones are exposed in the footwall of Segment A and the hanging wall of Segment B (see box in Figs. 1b and 2a).
The triangular CHJ outcrop exposes grey-yellow sandstones of the Moab Member (with patches of Morrison Fm on top), while the Courthouse Canyon to the south exposes older rocks, including the Navajo sandstone. The Moab Member typically has a porosity of around 20 %-25 % (Antonellini and Aydin, 1994), which is strongly reduced, to between 17 % and 1 % at the CHJ (Eichhubl et al., 2009;. Foxford et al. (1998) constrained the throw along Segment A to about of 370 m just south of Courthouse Rock (see Fig. 2 for location) and between 100 and 20 m in Mill Canyon, just north of Courthouse Rock (Fig. 1b). Eichhubl et al. (2009) estimated a throw on Segment A of about 200 m at the CHJ and on Segment B on the order of 100 m. Johansen et al. (2005) published a detailed map of the outcrop, showing the distribution, location and interaction (crosscutting relationships) between the different structure types mapped. These include joints and three types of deformation bands, including thick deformation bands (ThickDB) and thin deformation bands (ThinDB), occasionally with a slip plane (ThinDB-S). All of these types occur in bundles or clusters. While the ThickDB (∼ 1-1.5 mm wide) are domi-  Nuccio and Condon, 1996;Doelling, 2002;Gutiérrez, 2004;Fossen et al., 2005;Johansen et al., 2005;and Eichhubl et al., 2009, with  nated by large grain size (∼ 50-150 µm) and occasional undeformed grains, the ThinDB are truly cataclastic with grain size < 5 µm. Johansen et al. (2005) used crosscutting relationships to demonstrate that the ThickDB are older than the ThinDB. In addition, ThinDB were subdivided into seven separate evolutionary sets ( Fig. 12 in Johansen et al., 2005). Our own observations are consistent with this, and we will rely on the analysis of Johansen et al. (2005) in the remainder of the paper.
Occasionally, ThinDB in a bundle display smooth, corrugated surfaces (Fig. 2b-e), which are interpreted as slip planes (ThinDB-S) by Johansen et al. (2005). These are the slip surfaces that are measured in this study. Because they are corrugated slip surfaces, rather than cataclastic bands, we assume that they follow the Mohr-Coulomb failure criterion. These slip surfaces are relatively small compared to the entire deformation band bundle. Eichhubl et al. (2009) also reported joints in CHJ. Some occur along deformation band bundles, indicating that some of the jointing postdates cataclasis, in agreement with the regional trend. Davatzes et al. (2005) used a static-elastic calculation of stress at CHJ to explain the formation of these joints, controlled by changes in remote tectonic stress, burial depth, fluid pressure and rock properties.
Palaeostress and palaeostrain reconstruction algorithms share many assumptions, although they differ strongly in their approach. Only the orientations and the relative sizes (expressed in the ratio R) of the principal stresses are reconstructed, not the full tensor. This produces (in the case of palaeostress using the "dynamic fault slip analysis") a "reduced stress tensor" and for the "kinematic palaeostress analysis" a reduced strain tensor. In the case of coaxial deformation these can be considered to coincide with the stress axes and the reduced stress ratio (Sperner et al., 1993;Sperner, 1996;Ilic and Neubauer, 2005), but the differences in infinitesimal strain (earthquake analysis) and finite strain (fault slip data) have not yet been fully addressed. For a more complete discussion of these issues see Angelier (1994), Ramsay and Lisle (2000), Celerier et al. (2012) and Simon (2019).
The kinematic approach is based on the Mohr-Coulomb criterion (Coulomb, 1776;Mohr, 1900) and assumes a high shear-stress-to-normal-stress ratio. An example is the Numeric Dynamic Analysis (NDA), which calculates the orientation of the kinematic axes and the kinematic reduced ratio. These methods treat all faults as being newly formed. A further requirement is that the contraction and extension axes lie in the plane defined by the slip direction and the fault plane normal. This makes this method unsuitable for the use with faults with oblique striae and reactivated faults (Sperner, 1996). These methods also require the assumption of an angle of internal friction (generally assumed to be 30 • ), neglecting the natural variability of this parameter (Sperner, 1996;Sperner et al., 1993).
The dynamic method assumes that a slip on a fault plane will be parallel to the direction of maximum-resolved shear stress on the fault plane (Wallace, 1951;Bott, 1959). The Direct Stress Inversion method (DSI) (Angelier, 1979(Angelier, , 1984(Angelier, , 1990 is an example of a Wallace-and-Bott-type palaeostress analysis, in which this argument is inverted. Important assumptions of stress inversion are that the stress field is homogeneous, each fault is independent of the others, they are not coeval and they do not interact (Simon, 2019). Celerier et al. (2012) and Simon (2019) discuss the similarities between palaeostress and palaeostrain analysis and the simple geometric relationships that are sometimes observed between strain and stress axes when the fault and slip orientations are clustered. However, in cases when these orientations are more scattered, the two approaches yield different results. Celerier et al. (2012) argue that "by summing each fault contribution to obtain global strain tensor, strain tensor may be well suited for averaging strain within a volume". The inverse methods of stress analysis on the other hand may be more suited to separate heterogeneous data.
Separating heterogeneous fault orientation-slip direction pairs into homogenous subsets, which can reasonably be considered to originate from a single phase, is one of the major challenges in palaeostress calculations (see for example Sippel et al., 2009, and references therein), and individual faults that are wrongly assigned or are overprinted by a later stress state can cause the inversion algorithm to fail. We assume that in CHJ the structures evolved under the same regional stress field but are interacting with each other, causing overprinting and local stress rotations. For this reason, a palaeostrain analysis is in our opinion more suited to analyse this outcrop.
We used the implementation of NDA in the program Tectonics FP (Ortner et al., 2002, http://www.tectonicsfp.com/, last access: January 2011). For completeness and discussion, the palaeostress results of the DSI method are also included in the figures. We indicate cases where palaeostress and palaeostrain axes are very similar.

Field observations
We measured the orientation of both the slip planes with polished surfaces as well as the undulations, which are interpreted to represent the slip direction (see Angelier, 1994;Doblas, 1998, Fig. 2b and c) in the different generations of deformation band bundles, following the set numbering and evolutionary sequence of Johansen et al. (2005). Macroscopic observations on samples cut perpendicular to the slip plane show an increase of deformation band frequency towards the slip plane, in addition to a bleaching of the matrix (Fig. 2d). Observations in thin sections show an anastomosing network of cataclastic deformation bands with a significant reduction in grain size (Fig. 2e), in agreement with Aydin and Johnson (1978) and Fossen et al. (2007). This demonstrates that there is a genetic relation between bundles and slip planes and suggests that these slip planes form under the same stress conditions as the bundle. The slickensides do not react with HCl. We note that if the slip planes had formed late and in a single event, a palaeostress/strain analysis on the different sets would yield the same result. More data from detailed optical and broad-ion-beam-polished scanning electron microscope (SEM) micrographs of samples of deformation bands from this outcrop can be found in Komoroczi (2014).
In Fig. 3, we augment the different sets mapped by Jonhansen et al. (2005) with our own observations. We have only included deformation band bundles longer than 5 m. GPS stations not associated with a drawn line are of shorter bundles or isolated slip planes. The "forking geometry" (Simón, 2019), formed by Segment B and the long deformation band bundle of Set 3 and 4 (see Fig. 3a) against Segment A, supports the interpretation that Segment B abuts against Segment A. We collected 207 separate fault orientation-slip direction measurements from the larger CHJ area (see Fig. 3). The bulk of these measurements (n = 123) are collected at the CHJ outcrop (Fig. 3a). The remaining data come from both inside the fault zone and from the footwall of Segments B and C (Fig. 3b) where very similar-looking slip planes were found, associated with deformation band bundles. The observed slip planes often form conjugate sets. Measurements at the CHJ and those on top of the unnamed tower directly west of the Courthouse Rock are in the Moab Member, while those at the base of the Courthouse Rock (Fig. 3b) are in the Slick Rock and Dewey Bridge members. Those inside the fault zone of Segments B and C are in either Slick Rock or Moab members.
Where ThinDB-S have undulations on the metre and smaller scale (see Fig. 4), the different slip directions are generally parallel (Fig. 4d). The observation that slip planes typically occur in the core of a deformation band bundle is illustrated by this example (Fig. 4b and c). To better capture the variability of slip direction and reduce the effect of individual outliers due to local spurious orientations or measurement errors, we measured the fault plane and striation orientation up to three times on a single structure, about a metre apart. The authors strived to include all visible slip planes, but because a reasonable outcrop quality is required to measure both the orientation of the surface and the corrugations, not all of the visible slip planes were included in the final analysis.
Like Johansen et al. (2005), we did not observe slip planes in the set consisting of ThickDB, which are associated with Segment A, but the seven sets associated with ThinDB can all be recognized in our data set of fault orientation-slip direction pairs.
Fault orientation-slip direction pairs measured outside of the Branch Point were grouped into four locations (I-IV) from east to west (see Fig. 3b), based on their geographical spread. Data from location I contain many observations along the base of the Courthouse Rock in the footwall of Segment B, directly south of the branch point (see Figs. 1 and 5).
No polished slip planes and slickensides are observed on Segments A and B, but observations of the metre-scale undu-lations of the main fault structures show slip on these structures was dip slip (at least in the final stages of movement; Fig. 3a).

Separating the data from Courthouse Junction
We separated our 123 fault orientation-slip direction pairs from CHJ into the seven evolutionary sets, following the subdivision of Johansen et al. (2005), using the average orientations and published locations to differentiate them. Our own observations were classified into their scheme, using the orientation only. Since Set 1 and Set 2, Set 3 and Set 4 are conjugates (the average poles are at an angle of around 60 • ), we combined these into "Set 1 and 2", "Set 3 and 4". The remaining data were grouped into "Set 5, 6 and 7" when we demonstrated this combined set produces kinematically consistent results (see Fig. 6).

Palaeostrain results at the Courthouse Junction
The calculated palaeostress/palaeostrain results at CHJ are shown in Fig. 6. As discussed above, we focus on strain axes calculated using NDA. Results show extension in the NW-SE direction for Set 1 and 2 and Set 3 and 4. For Set 1 and 2 the extension axis has the direction 334 • /03 • and for Set 3 and 4, 336 • /10 • , with R values of 0.48 and 0.46 respectively. The maximum compressive strain axis is near vertical for Set 1 and 2 but is rotated slightly towards the west in the lower hemisphere stereoplot for Set 3 and 4. The results for Set 5, 6 and 7 show extension in the NE-SW direction (orientation 219 • /03 • ), with an R value of 0.47. We note that the DSI palaeostress axes for this particular set are broadly similar to the palaeostrain axes, although the R values for both methods are very different. Figure 7 shows the palaeostrain/palaeostress results outside CHJ. At location I (Figs. 1, 5), 37 measurements were collected in the footwall of Segment B (Fig. 5). Based on the distribution of slip directions, we separated the data in two subsets (Set Ia and Set Ib). Since we did not observe crosscutting relations at this locality, we cannot determine relative age relationships. By separating the data from location I into Sets Ia and Ib, the variance reduces for both methods. The histogram of Fig. 7 shows the angular mismatch between the measured orientation of slip on a fault plane and the calculated slip direction on that plane, providing a quality control. It can be seen that the mismatch is quite small. Both Set Ia and Set Ib as well as the full data set in location I show a vertical maximum compressive strain axis and a north-south In location II and location III, the palaeostress and palaeostrain axes are very similar. As these locations are in the footwall of Segment B, away from major fault interactions, it seems reasonable that here the main palaeostress requirements of a homogeneous stress field and independency between the faults are met. As a result, NDA and DSI produce results that are similar, with a vertical maximum compressive stress/strain axis and N-S extension direction.

Palaeostrain results for the study area outside the CHJ
The data at location IV in Fig. 3 are insufficient to calculate a stress or strain tensor; at least four but preferably 15 independent measurements are required (Sperner, 1996;Sperner et al., 1993;Simon, 2019). The sample location east of the stream in location IV (red great circle in Fig. 3) has an anomalous orientation compared to the north-dipping fault at the western location (blue great circles), but slip directions are near parallel, suggesting reactivation of an earlier structure. The slip directions do not fit the expected north-south extension. This might be due to the perturbation of the stress  field by the nearby branch point on Segment C with the minor fault structure (highlighted with an * in Fig. 1b). . Fault plane orientation-slip direction and palaeostress (DSI) and palaeostrain (NDA) results (for the three sets) from CHJ (Fig. 3a). The separation is based on Johansen et al. (2005), with kinematically related sets combined into a single data set.

Interpretation and discussion
Our study shows that Set 1 and 2 and Set 3 and 4 formed under similar extension directions, supporting the interpretation of Johansen et al. (2005) and suggesting a more complex history with gradually rotating deformation axes. Our Sets 5, 6 and 7 are kinematically consistent (but not conjugate); thus we interpret them to have formed roughly at the same time. In this set, the palaeostress axes are consistent with the palaeostrain axes, while this is not the case for combined Set 1 and 2 and Set 3 and 4. Similar to the situation described above for locations II and III, Set 5, 6 and 7 is interpreted to be the only set in CHJ where the requirements of a homogeneous stress field and independency between the faults are met. Thus, both methods produce similar results here. The same can be observed for the palaeostress results of location Ia, outside of CHJ.

Deformation during the transition from kinematically to geometrically linked abutting faults
In this study we used both kinematic and dynamic analysis of fault slip data. Using the analysis of Johansen et al. (2005) to separate data from different generations, the strain-based NDA method, which averages strain in a volume, can analyse this evolution. The NDA method yields consistent results over the study area, which is a reasonable result for this system (Celerier et al., 2012). In those locations where we can reasonably assume a homogenous stress field and no fault interaction, we found similar stress and strain axes. In the other locations where the relation between calculated stress and strain axes diverges, we infer that this is because DSI cannot separate the effects of the gradually rotating stress field as Segment B approaches Segment A. Figure 8 summarizes the evolution of the stress state in three stages, assuming coaxial deformation, as (i) Segments A and B are isolated, (ii) kinematically linked and (iii) geometrically linked.
For stage (i), our observations are consistent with the interpretation of Johansen et al. (2005) that Segment A formed first. Since no striated slip planes are observed here, we assume a vertical σ 1 and estimate σ 3 at 40-45 • , perpendicular to Segment A. This is interpreted as the regional stress field. ThickDB also formed during this phase .
In stage (ii) (Fig. 8b), the two faults start to interact. Set 1 and 2 start to form in the damage zone around the fault tip of Segment B and develop further in the interaction damage zone, in a heterogeneous rotating stress field. This stress field eventually causes a second group of deformation structures (set 3 and 4) to form. The difference in orientation of the extension axis between the results of Set 1 and 2 and Set 3 and 4 is small, but the maximum principle strain (λ 1 in Fig. 6) rotates slightly away from vertical.
In stage (iii) the Segments A and B connected and became geometrically hard linked. The stress state in the Courthouse branch point rotates back to regional extension direction of NE-SW tension as seen in the NDA analysis of Set 5, 6 and 7 (Fig. 8c).
The relative timing of the observations outside CHJ is poorly constrained, but extension at locations I, II and III was roughly N-S (Fig. 8b). If we assume that Set Ia formed before Set Ib, the rotation caused by the linkage of Segments A and B corresponds to the minor rotation of axes at this location. In locations II and III we infer that the stress states were not influenced by the connection of Segments A and B, and deformation axes are assumed to be relatively constant here. At location IV our data suggest that extension was with a slight anticlockwise rotation.
The changes in orientation of the principal stresses over time in CHJ are interpreted to show the local, third-order stress evolution. The northern termination of the Moab Fault (Segments A, B and C but before the relay ramps between them were breached) initially formed in a regional stress field, an overall NE-SW regional extension. The approach of Segment B as the relay ramp was breaching caused local stress rotations around the fault tip, causing new structures to form. Thomas and Pollard (1993) showed that neartip stresses in models of en échelon mode I fractures in a similar geometry dominate over the regional stresses and these stresses change before oblique intersection. Davatzes et al. (2005) modelled the stress state at a branch point by combining observations from the Courthouse branch point and the branch point close to location IV (Fig. 3). Their models are iterative, and the calculated stress orientations are similar to our results. Peacock at al. (2017a) discussed how the interaction damage zone can influence fluid flow before the faults become geometrically linked. Eichubl et al. (2009) inferred that most diagenetic cement deposits in the CHJ are associated with Set 5, 6 and 7, which formed after the faults became hard linked. They identify the structures associated with this set as joints and not as deformation bands. Johansen et al. (2005) however shows that joints in the CHJ formed inside a deformation band (an observation shared by us), and thus jointing  Peacock et al., 2017b). Blue arrows point in the direction of extension; light blue symbols are inferred. (a) The faults are isolated. The Moab Fault (Segment A) forms under the influence of the regional NE-SW extension direction (based on fault undulations on Segment A and the thick deformation bands set, described by Johansen et al., 2015). (b) When Segment B, initially formed as a parallel fault, started to rotate and approached Segment A, the segments became kinematically linked, and the minimum principal stress rotated at CHJ to the NW-SE (deformation band Sets 1 and 2 and 3 and 4). (c) When Segments A and B connected and became hard linked or geometrically linked, the minimum principal stress at CHJ rotated back to the regional stress state (deformation band Set 5, 6 and 7). The outline of these sketches corresponds roughly to the area shown in Fig. 3b. postdates the formation of deformation bands. This suggests that the diagenetic overprint in the CHJ is later than the abutting of the fault.

Conclusions
We analyse subsequent generations of striated minor fault planes in the overprinting damage zones of interacting normal faults in porous sandstone, at Courthouse Junction, near Moab, UT.
By combining our own measurements with published data on the relative age of these structures, we classified the observations in four sets. Using numeric dynamic analysis (NDA) the orientation of the kinematic axes was calculated for three of these sets, and inferences were made on the first set. We found three different palaeo-extension directions in the sets, recording the evolution of stress/strain axes during the abutting process.
Stage 1 of deformation direction is a regional extension in the NE-SW direction. In this case study of a single-tip interaction of abutting faults, we find a near-90 • rotation of the extension direction (from NE-SW to NW-SE) as the stress fields around the main fault and the abutting faults start to interact when they are kinematically linked in stage (ii). Once the faults abut and become geometrically linked, stage (iii), there is a rotation back to the regional extension direction. This evolution represents a third-order local change in the regional stress field, caused by relay ramp breaching and abutting faults.
Data availability. Data of this study are available on reasonable request from the author.
Author contributions. HvG performed the field investigation, data analysis and conceptualized, and wrote the first draft. JLU provided supervision, in-depth discussions and validation throughout, including extensive revisions of the original draft.