Tectonic processes, flexure due to crustal loading, and dynamic mantle flow each impart a unique imprint on topography and geomorphic responses over time scales of 104 to 106 yr. This paper explores the mobility of regional drainage divides as a key geomorphic metric that can distinguish between the various processes driving crustal deformation in the greater Yellowstone region of the northwestern United States. We propose a new analysis that quantifies the differences between the location of the present-day drainage divide from divides synthetically generated from filtered topography to determine the relative impact of tectonic and dynamic mantle influences on landscape development. The greater Yellowstone region is an opportune location for this investigation because contrasting models have been proposed to explain the parabolic shape of elevated topography and active seismicity that outline the imprint of hypothesized hotspot activity. Drainage divides synthesized from topography filtered at 50, 100, and 150 km wavelengths within the greater Yellowstone region show that the locations of the actual and synthetic Snake River drainage divides are controlled by both dynamic and flexural mechanisms in the eastern greater Yellowstone region, but by flexural mechanisms only in the western greater Yellowstone region. The location of the actual divide deviates from its predicted position in the filtered topography where tectonic controls, such as active faults (e.g., Centennial and Teton faults), have uplifted large footwall blocks. Our results are consistent with the notion of a northeastward-propagating greater Yellowstone region topographic and seismic parabola, and suggest that Basin and Range extension follows from, rather than precedes, greater Yellowstone region dynamic topography. Furthermore, our analysis suggests that eastward migration of the Snake River drainage divide lags behind the continued northeastward propagation of high-standing topography associated with the Yellowstone geophysical anomaly by 1–2 m.y.
Watershed divides are topographic features that provide a potential metric by which to assess crustal deformation at a variety of length scales. In tectonically active areas topography and drainage divides are mobile elements of the landscape that can migrate kilometers to tens of kilometers on 105 to 106 yr time scales (e.g., Anderson, 1947; Taylor et al., 1987; Repenning et al., 1995; Harbor, 1997; Willett, 1999; Stüwe and Hintermueller, 2000; Marshall et al., 2003; Beranek et al., 2006; Lock et al., 2006). Many studies have demonstrated the relative influence of tectonic processes on landscape evolution by examining parameters derived from fluvial networks, such as drainage patterns, basin asymmetry, river long profiles, drainage evolution via reconstructed terrace longitudinal profiles, indices of channel concavity, stream power distribution, and channel slope–basin area relationships (e.g., Wells et al., 1988; Pazzaglia et al., 1998; Wegmann and Pazzaglia, 2002; Tomkin et al., 2003; Molin et al., 2004; Lock et al., 2006; Wobus et al., 2006). The spatial position of present-day and synthesized regional-scale drainage basin divides is an underutilized geomorphic metric that can be used to identify various overlapping geodynamic processes that drive landscape evolution. The approach that we present here, of generating synthetic divides on filtered digital topography, is a relatively rapid way to analyze the topographic imprint of different tectonic and dynamic processes, and has the advantage of not being limited, as are many of the traditional metrics, to any particular length scale of topography. Because topography is the integrated result of processes operating at a variety of spatial scales, the potential exists to interpret the topographic fingerprint of underlying geodynamic processes by analysis of synthetically generated drainage divides from filtered topography at equivalent length scales. The first goal of this study is to test the applicability of using synthetically generated drainage basin divides as a means of identifying the crustal deformation processes proposed for the greater Yellowstone region in the western United States. The second goal is to assess how different deformation processes collectively contribute to the building of the Yellowstone topography. We demonstrate that synthetically generated drainage divides are a portable metric that can be used to assess the geodynamics of regions where the underlying tectonic and dynamic processes are poorly understood.
Studies have shown that over geologic time the spatial positions of regional-scale drainage divides are sensitive to crustal deformation and the uplift of rocks. In considering the stability or mobility of regional-scale drainage basin divides, two end-member scenarios may be envisioned. The first is where the divide remains spatially fixed over long periods of time during which an overthickened orogenic crustal root is passively consumed via erosional unloading and its isostatic response. Examples of these fixed divides include those of great escarpments such as the eastern continental divides of Australia (Young, 1989) and North America (Hack, 1979). The second is where divides actively migrate in response to dynamic support of topography, which is the component of topography deriving from buoyancy and flow variations beneath the lithosphere, such as the axial divide of the Italian Apennines (D'Agostino et al., 2001). Knowing that active drainage divide migration can be the key feature in distinguishing between passive and active geodynamic settings, we investigate drainage divide migration potential in the greater Yellowstone region, a geodynamically active area where the processes influencing the present-day topography are fairly well defined. We utilize the spatial position of synthetically generated drainage basin divides from topography filtered at specific wavelengths as a geomorphic fingerprint to elucidate the contributions of dynamic, tectonic, and flexural processes to the integrated topography of the greater Yellowstone region.
We model drainage divides of the eastern Snake River, which drains a large portion of the greater Yellowstone region and forms a portion of the continental divide (Fig. 1), by filtering digital topography at three separate wavelengths (50, 100, and 150 km) in order to reveal the relevant length scales of the regional topography and the underlying geodynamic processes that control drainage divide location. This multiscale approach may prove useful in identifying the transition between long-wavelength dynamic (≥150 km) and shorter wavelength tectonic (faulting) and flexural (≤50 km) crustal deformational processes. We are also able to predict the geomorphic lag or landscape-response time associated with drainage divide migration in areas of active dynamic deformation, such as the greater Yellowstone region.
The greater Yellowstone area is an ideal location to test this approach because it is a geologically complex region in which dynamic, tectonic, and flexural features associated with the eastward propagation of the loci of Yellowstone volcanism and crustal deformation are superimposed on Sevier and Laramide compressional and Basin and Range extensional structures (Suppe et al., 1975; Pierce and Morgan, 1992; Smith and Braile, 1994; McQuarrie and Rodgers, 1998). We hypothesize that the Snake River drainage divide is a quickly evolving and sensitive indicator of active deformation associated with the Yellowstone hotspot (e.g., Link et al., 2005; cf. Young, 1989) and that it responds to deformation at a variety of length scales (Beranek et al., 2006).
Over the past 10 m.y., a steady northeast-trending progression of magmatic activity, currently centered at Yellowstone, has significantly changed the regional topography of the Eastern Snake River Plain, which today is a broad lowrelief depression circumscribed by an arcuate belt of high-standing topography (3000–3500 m) (Fig. 1) (Leeman, 1982; Anders et al., 1989; Malde, 1991; Pierce and Morgan, 1992). Although discussions continue on the precise geodynamic nature of the Yellowstone region, it is generally accepted that high heat flow, magmatism, and a parabolically shaped region of high-standing topography are manifestations of dynamic support resulting from anomalous buoyancy contrasts and convection in the upper asthenospheric mantle (Humphreys et al., 2000; Lowry et al., 2000; Burov and Guillou-Frottier, 2005).
Analysis of topographic, crustal thickness, and Bouguer gravity data (Fig. 2) suggests that topography in the greater Yellowstone region is not uniformly supported by low-density crustal roots (Airy isostasy), which in turn is an indication that at least a component of the long-wavelength topography is dynamically supported by active lithospheric and/or sub-lithospheric processes (e.g., Lowry et al. 2000). Perhaps the strongest observation to suggest a dynamic component to the topography in the greater Yellowstone region is that some of the thinnest crust (∼35 km) is observed in regions of high-standing topography (e.g., north of the Eastern Snake River Plain), and thus no correlation exists between regional topography and crustal thickness (Figs. 2A, 2C), implying the lack of a crustal root. Regionally slow seismic P-wave velocities in the upper mantle (Fig. 2D) (Yuan and Dueker, 2005), suggestive of anomalously warm and buoyant mantle at the base of the crust (Burov and Guillou-Frottier, 2005), and a long-wavelength positive geoid anomaly centered on Yellowstone (Roman and Wang, 2003) provide further evidence for dynamic support of the greater Yellowstone region.
Dynamic models, where the swell of high topography surrounding Yellowstone is supported by sublithospheric processes (e.g., Anders et al., 1989; Christiansen et al., 2002; Pierce and Morgan, 1992; Humphreys et al., 2000; Pysklywec and Mitrovica, 1997), and flexural models (McQuarrie and Rodgers, 1998), where flexural subsidence of the Eastern Snake River Plain is driven by the emplacement of a lower-crustal load resulting in uplift of the shoulders of orthogonally flanking mountain ranges, have been proposed for the Yellowstone–Eastern Snake River Plain region to account for the observed crustal deformation and belt of arcuate high-standing topography. In a related study of a landscape shaped by dynamic support, D'Agostino et al. (2001) documented short- and long-wavelength topographic features in the central Apennines of the Italian Peninsula and concluded that the short-wavelength relief is supported by flexural and tectonic stresses, while the longer wavelength (≥150 km) topography of the Apennines is an expression of dynamic support by mantle convection and is uncompensated by crustal roots. Similarly, in the Yellowstone region both short- and long-wavelength features are observed in the topography (Fig. 2A) and gravity anomalies (Fig. 2B) that are attributable to tectonic (individual fault-bounded uplifts), flexural response (crustal loading beneath the Eastern Snake River Plain), and dynamic support mechanisms. The use of drainage divides generated from filtered (50, 100, and 150 km) and non filtered digital elevation models (DEMs) allows for the discrimination of contributing geodynamic mechanisms to the present-day position of the eastern portion of the Snake River watershed, as well as to an estimation of the landscape-scale lag times linked to the continuing north-eastward migration of the Yellowstone hotspot.
The geomorphic and topographic signatures of dynamic support in the greater Yellowstone region were sought through modeling of drainage basin divides to topography filtered at short (50 km), medium (100 km), and long (150 km) wavelengths. A DEM was produced for the greater Yellowstone region by combining a 90 m gridded Shuttle Radar Topography Mission (SRTM) elevation data set with resampled (90 m) GTOPO30 elevation data, utilized to remove data gaps in the original SRTM product. The resulting continuous DEM of the greater Yellow stone region was resampled to a cell size of 0.004° (∼400 m). The resampled DEM was filtered at wavelengths of 50, 100, and 150 km in the spectral domain using the two-dimensional forward Fast Fourier Transform (FFT) function of the Generic Mapping Tools (Wessel and Smith, 1991) by applying a band-pass filter with a 10-km-wide cosine taper (Fig. 3), a standard consideration for maintaining a smooth, continuous function in the frequency domain of the FFT filter. For example, smoothing of the SRTM data set with a 50 km filter resulted in the removal of all topographic features with a spectral dimension of <50 km, preserving all features with a spectral dimension >60 km, and tapering all features between 50 and 60 km. The spatial area of the DEM was increased well beyond the greater Yellowstone region to the point at which edge effects produced during the filtering process were nonexistent within the region of interest.
Drainage Divide Generation
Watersheds and drainage basin divides were generated for both the unfiltered and filtered DEM data sets by using the FLOWDIRECTION and WATERSHED routines within ArcInfo (Fig. 3). Positional errors in the location of the generated drainage divides, both for the unfiltered and filtered topography, are a function of the horizontal and vertical measurement errors associated with the original SRTM data, considered to be ≤15 m for each 90 m grid cell (Rodriguez et al., 2006). For this study horizontal positional errors are considered to be <400 m for the synthetic divides generated from the resampled DEM with a grid node spacing of ∼400 m. Potential erroneous horizontal shifts of as much as 400 m for the synthetic drainage divides with respect to the actual divide are insignificant at the spatial scale of analysis, and are thus ignored.
Sensitivity analyses were performed to test the degree to which local topographic features are capable of producing farfield effects on the location of the synthetic 50 km divide on the resampled ∼400 m DEM. As a test of the sensitivity of the position of the 50 km filtered divide location to individual mountain ranges, two additional DEM data sets were generated where the topography underlying the Centennial and Teton Ranges (Fig. 1) was removed from the DEM and replaced with the average elevation of adjacent basin floors (∼2000 m). These two ranges were selected because of their proximity to the actual Snake River drainage divide and because they represent youthful footwall uplifts of active Basin and Range faults that may locally control divide location.
Drainage Basin Comparison
Swath Topographic Profile
A swath topographic profile, 500 km long, 10 km wide, and oriented northwest-southeast across the Eastern Snake River Plain (Fig. 4), was collected from the 90 m DEM to better illustrate the similarities and differences in the location of actual and synthetic drainage divide locations (Fig. 5). The swath topographic profile is the combination of maximum, mean, and minimum elevation data collected from six parallel individual topographic profiles evenly spaced 2 km apart. The sampling frequency for individual profiles was 2 km.
Spatial Relationship between Actual and Synthetic Drainage Divides
The actual Snake River drainage divide roughly follows the highest topography within the general confines of the seismic parabola of Anders et al. (1989)(Fig. 4). Exceptions to this are in the apex of the parabola, centered on the active Yellowstone hotspot, and adjacent to the Centennial Range (northeast quadrant), where the highest topography is ∼50 km north of the divide. Conversely, all three of the synthetic drainage divides (50, 100, and 150 km) roughly follow the highest topography around the apex of the seismic parabola. At the latitude of Yellowstone Lake and to the north, all three synthetic divides deviate significantly (55 ± 5 km) from the actual divide, following the trend of highest topography (Fig. 4). This deviation continues to the north and west of Yellowstone for ∼100 km. West of the Centennial Range, the synthetic divides begin to converge toward the actual divide, such that by long 113°W, all divides on the north side of the Eastern Snake River Plain are within 10 km of each other (Fig. 4). Within the southwest quadrant along the southern flank of the Eastern Snake River Plain, the actual and synthetic divides are similarly located, within 15 km of each other, with the notable exception of two basinward deviations of the 50 km synthetic divide centered on long 114°W and 111.5° W, respectively (Fig. 4). Within the southeast and northeast quadrants, the location of the 50 km synthetic divide follows the crest of the Teton Range, 65–95 km farther west than the actual Snake River divide and the 100 km and 150 km synthetic divides (Fig. 4). The medium- and long-wavelength 100 km and 150 km synthetic divides follow similar paths around the seismic parabola, alternating between which divide is located farther from the axis of the Eastern Snake River Plain (Fig. 4). A striking observation is that all of the divides, synthetic and actual, converge along the high topography at the western end of the Eastern Snake River Plain (Fig. 4; 01Table 1, northwest and southwest quadrants).
Deviation in Drainage Basin Area between Actual and Synthetic Basins
The percent deviation (D) in area (km2) between the actual Snake River watershed and modeled watersheds numerically illustrates the relationships described in the previous paragraph and visible in Figure 401(Table 1). For the entire eastern Snake River basin, the 50 km synthetic divide deviates the most (36%) from the aerial extent of the actual basin 01(Table 1). The greatest individual and mean deviations occur in the northeast quadrant, the focus of active volcanism at Yellowstone today. In the southeast quadrant, the 50 km synthetic drainage basin has a deviation of 80% compared to the actual basin area, in large part due to control of the 50 km divide location by the Teton Range. The smallest individual and mean deviations occur in the northwest and southwest quadrants, at the western end of the Eastern Snake River Plain, where the actual and synthetic divides converge at all wavelengths of topography (Fig. 4; 01Table 1).
Deviations exist between the actual and 50 km synthetic divides near the Centennial and Teton Ranges (Fig. 4). The extent to which these two ranges are capable of controlling the position of both the actual and synthetic 50 km drainage divide was investigated via sensitivity analyses. The Centennial Range strongly localizes a portion of both the actual and 50 km synthetic drainage divides along the Montana-Idaho border, but this control disappears when the range is removed (Fig. 4, comparison between the red, orange, and purple lines). In the absence of the Centennial Range, the 50 km synthetic divide is aligned with the 100 and 150 km synthetic divides 30–35 km north of the actual divide. In contrast, the actual drainage divide is ∼75 km east of the north-south–trending Teton Range, a position well predicted by the 50 km filtered topography when the range is removed (Fig. 4, comparison between the red, orange, and yellow lines). This eastward shift of the 50 km synthetic divide in the absence of the Teton Range brings it into line with the actual drainage divide, at the latitude of the range; however, it is important to note that the removal of the Teton topography from the DEM does not affect the location of the distended 50 km divide southwest of the Teton Range, centered on long 111.5°W (Fig. 4).
Swath Topographic Profile
The swath topographic profile resolves two swells in topography on the southeastern side of the Eastern Snake River Plain (Figs. 4 and 5). The 50 km divide approximately follows the lower of the two swells with all other divides located along the crest of the higher swell. Two swells are also observed along the northwestern side of the Eastern Snake River Plain, but are not as distinct as across the Eastern Snake River Plain. All drainage divides follow the axis of the highest topographic swell north of the Eastern Snake River Plain.
The position of the Snake River watershed divide as it wraps around the greater Yellowstone region is controlled by short, intermediate, and long topographic length scales, as revealed in the mismatch between the actual and synthetic divides generated on the filtered topography. The most obvious mismatch in the actual and synthetic divides occurs at Yellowstone National Park (Fig. 4). Here the actual divide is controlled by small scale-length (<50 km) features and is positioned west of its predicted location (Figs. 3 and 4). These observations are consistent with previous studies (e.g., Lowry et al., 2000) that call upon dynamic mantle support for Yellowstone, based upon slow mantle P-wave velocities, observed crustal thickness, and Bouguer gravity anomalies (Fig. 2). P-wave tomography from the Yellowstone area suggests that the long-wavelength topographic highs controlling the location of the synthetic drainage divides are the result of upwelling mantle beneath the northeastern part of Yellowstone National Park (Yuan and Dueker, 2005).
The mismatch between the actual and predicted divide in Yellowstone National Park may reflect the local influence of short-wavelength topographic features such as drainage rearrange ment by effusive Quaternary volcanism at Yellowstone, active faults, highlands supported by resistant rock types, or valleys sculpted by fluvial and glacial processes. If we take the synthetic divides as accurately locating the apex of dynamic support, the mismatch represents the spatial length scale across which surface processes operate in keeping the divide astride the highest standing part of the landscape. These surface processes ultimately drive drainage capture, a process that is both unsteady and beholden to complex interactions of fluvial and glacial erosion, hillslope processes, and groundwater sapping. Given that the greater Yellowstone region geophysical anomaly has propagated northeastward at ∼50 km/m.y. based upon the time-transgressive locations of earlier silicic volcanic centers (Fig. 1), and that the predicted synthetic divides are 50–100 km northeast of the actual divides, we argue for a 1–2 m.y. time lag between the dynamic mantle processes that uplift the landscape and the surface processes that respond to that uplift by working to center the drainage divide astride the apex of the uplift.
Another marked deviation in the locations of synthetic and actual divides exists southwest of the Teton Range, where the 50 km synthetic divide distends significantly toward the Eastern Snake River Plain (Fig. 4). Embedded in the swath profile in this location is a bimodal distribution of elevations that are most prominent south of the Snake River Plain (Fig. 5). At the 50 km filtered wavelength the synthetic divide is coincident with the approximate crest of the lower amplitude topographic swell. This shorter wavelength topographic feature is interpreted to be the result of lithospheric flexure resulting from the high-density load of late Miocene lower-crustal basaltic intrusives beneath the Eastern Snake River Plain onto relatively weak lithosphere (e.g., Sparlin et al., 1982; McQuarrie and Rodgers, 1998; D'Agostino et al., 2001).
Whereas the modeled and actual divides diverge in the area surrounding Yellowstone, they converge farther west along the margins of the hotspot track (Fig. 4; 01Table 1), a visual trend supported by the drainage basin area percent deviation analysis 01(Table 1). Westward cooling of the lithosphere, as evidenced by the progressive increase in age of volcanism away from Yellowstone (Fig. 1), a westward deepening of the Eastern Snake River Plain (Fig. 2C), and a more pronounced Bouguer anomaly underlying the western end of the Eastern Snake River Plain (Fig. 2B) argue for an increase in effective lithospheric elastic thickness and a topographic expression consistent with lithospheric flexure. If this is true, then the convergence of the actual and synthetic drainage divides at all wavelengths of filtering on both sides of the western end of the Eastern Snake River Plain is likely the result of flexural focusing of topographic axes parallel to the downwarped Snake River Plain. In light of the time-transgressive eastward migration of silicic volcanism (Morgan, 1972) along the Eastern Snake River Plain (Fig. 1), this interpretation enforces the diachronous nature of topographic growth and decay that is likely to result from dynamic mantle processes.
Local deviations in the position of the actual and synthetic divides result from short wavelength features such as individual tectonically active fault-bounded ranges. For example, a local deviation exists along the north flank of the Eastern Snake River Plain where the Centennial Range (Centennial fault) controls the location of the actual divide (Fig. 4). The tectonically active east-west–oriented Centennial Range (Pierce and Morgan, 1992) has retarded north to northwest expansion of the modern Snake River drainage divide, even though the region of highest standing topography, and thus the most likely location for divide development, is located 30–35 km north of the range, as attested to by the collocation of the 50, 100, and 150 km synthetic divides due north of the range, and by the 65 km northwest shift in the 50 km synthetic divide location upon digital removal of the Centennial Range within the DEM (Fig. 4). Likewise, the digital removal of the Teton Range topography indicates that individual ranges are capable of controlling the location of the synthetically generated divides, in this case the 50 km synthetic divide, as it jumps 60–65 km to the east upon removal of the range within the DEM (Fig. 4). However, the spatial influence of the Teton Range is limited, as the location of the 50 km synthetic divide north and south of the latitude of the Tetons is unaltered upon digital removal of the range from the DEM (Fig. 4).
A consideration in how to interpret the synthetic divides from the actual divide arises in terms of the possible impact that preexisting topography may have on the location of the synthetic divides. The Yellowstone plume is currently pushing northeastward into major Laramide uplifts, including the Beartooth, Absaroka, Gros Ventre, and Wind River Ranges (Figs. 1 and 4). We conclude that the influence of preexisting topography associated with these ranges must be small given that the synthetic Eastern Snake River Plain divides follow paths that take them across the strike of these ranges at oblique angles, a characteristic not shared by the actual Eastern Snake River Plain divide, which comprises a component of the continental divide that is located significantly west of the axis of high-standing topography in the vicinity of Yellowstone National Park. Furthermore, topographic expression of these Laramide Ranges has been enhanced by dynamic support processes at the greater Yellowstone region and elsewhere (Dethier, 2001), so we postulate that the effect of preexisting topography of the scale of an individual Laramide range is a second-order consideration, subordinate to dynamic mantle processes responsible for whole-scale uplift of large regions.
Synthetic drainage divides, produced by band-pass filtering topography at short (50 km), medium (100 km), and long (150 km) wavelengths, illustrate the competing effects of tectonic, flexural, and dynamic mantle flow processes, respectively, on the collective topographic expression of the greater Yellowstone region. At Yellowstone National Park and along the northern arm of the seismic parabola proximal to the park, the actual divide appears to be controlled by topographic wavelengths shorter than 50 km demonstrated to be locally linked to actively uplifting footwall blocks along active Basin and Range normal faults. These foot-wall uplifts are all in the wake of the greater Yellowstone region dynamic topography, suggesting that Basin and Range extension in the greater Yellowstone region results from, rather than driving, the greater Yellowstone region geophysical anomaly. Other short-wavelength topographic effects resulting from Quaternary volcanism, rock type, or glacial and fluvial sculpting of valleys likely also play a role in localizing the actual divide within Yellowstone west of its predicted dynamically supported location. There is also almost certainly a lag time in the surficial processes working to localize the actual divide atop its eastward predicted position from its current westward position within the park. Such surface processes relocalize the divide through drainage captures and are likely to be unsteady and influenced by alternating glacial versus fluvial erosion, and perhaps extrusive volcanism, at least during the Quaternary.
West of Yellowstone National Park the three filtered synthetic divides are all more or less located along the crest of the arcuate belt of high-standing topography and coincident with the seismic parabola that we interpret to be dynamically supported. Positioning of the 50 km synthetic divide southwest of the Teton Range is controlled by the strong local effects of lithospheric flexure of the Eastern Snake River Plain flanks in response to flood basalts loading a relatively weak lithosphere (Fig. 5). Proximal to Yellowstone National Park, these short-wavelength flexural effects are superimposed on the broader effects of dynamic support. However, in the western portion of the greater Yellowstone region, all divides converge along the highest topography paralleling the Eastern Snake River Plain; we interpret this to indicate a transition to topography shaped only by lithospheric flexure of a broader wavelength.
This analysis has demonstrated that geodynamic processes such as dynamic support of topography by upwelling mantle and lithospheric flexure impart decipherable signatures upon the gross topography of a region. The generation of synthetic drainage basins from digital topography filtered at specific wavelength allows for the interpretation of the geodynamic processes important to the shaping of landscapes at varied spatial scales. This geomorphic analysis tool should be applicable to other regions characterized by the varying importance of different geodynamic processes.
We thank Skye Cooley for a review of an early draft of the manuscript, and Cameron Wobus and Paul Link for constructive and beneficial formal reviews. This research was the outgrowth of a graduate course project supported in part by the Department of Earth and Environmental Sciences at Lehigh University. The research received no formal external funding.