The Significance of Ancient Buried Landscapes as Natural Geomorphic Experiments

There is considerable interest in developing quantitative methods for analyzing present‐day fluvial landscapes with a view to extracting information about tectonic forcing and drainage evolution, together with the influence of lithologic substrates and of paleoclimatic variations. In view of the multifactorial nature of this complex problem, it has previously been proposed that natural geomorphic experiments could play a significant role in developing a quantitative understanding of landscape growth and decay. Here, we describe and analyze a stacked sequence of five buried transient landscapes that punctuate marine strata along the fringes of the North Atlantic Ocean. We propose that these landscapes constitute a suite of natural experiments, which illuminate significant aspects of quantitative fluvial geomorphology. Our preliminary analysis of four of these buried landscapes suggests that the amplitude of external tectonic forcing plays a significant role in fluvial landscape evolution. In future, we hope that this suite of natural experiments will be further exploited by the fluvial community with a view to identifying the most appropriate analytical techniques.


Introduction
There is considerable interest in developing quantitative methodologies for analyzing fluvial landscapes on Earth.The primary reasons for this interest are threefold.First, fluvial landscapes reflect, and therefore contain indirect information about, regional tectonic forcing as functions of space and time.There is increasing evidence that this forcing is often triggered by mantle convection.An important corollary is that fluvial landscapes can potentially inform otherwise inaccessible sub-plate processes.Secondly, fluvial landscapes can be exploited to gain a better understanding of erosive processes, which depend upon a combination of paleoclimate and lithology.Finally, there is a hydraulic engineering imperative in understanding the relationship between precipitation, channel flow and deposition.
Fluvial analysis is underpinned by a generalized version of the empirical stream power formulation, which relates the rate at which height along a river channel, z, varies with time, t, to some combination of tectonic uplift, U, and erosion, E, such that ∂z ∂t = U(x,t) E(x,t). (1) Erosion can be parametrized using where A is measured upstream area and P is precipitation.v, m, n, and κ are erosional parameters, which are expected to vary as a function of space and time (Howard & Kerby, 1983;Whipple & Tucker, 1999).It is generally accepted that channel erosion is predominantly an advective process, which means that κ can be disregarded within fluvial channels themselves.However, there is ongoing debate about the values of the other three erosional parameters which are usually determined either by slope-area analysis or by calibrated inverse modeling (O'Malley et al., 2021;Whipple & Tucker, 1999).The value of v sets the pace of fluvial erosion.A major difficulty is reliably determining how v varies as a function of time and space.The value of m is accepted to be fractional if n = 1 (i.e., m/n is 0.35-0.6;Whipple & Tucker, 1999).However, the value of the empirical constant n is subject to debate.A parsimonious approach suggests that since there is no reliable evidence for kinematic shock waves within modern fluvial landscapes that n = 1 is a satisfactory assumption (Pritchard et al., 2009).Furthermore, continent-wide inventories of river profiles can be most accurately matched provided that n = 1 (Rudge et al., 2015).We note here that alternative schemes have been used to argue that n > 1 or even that n ≫ 1 (Gallen & Fernández-Blanco, 2021;Lague, 2014).
The reason for these ongoing and unresolved disagreements undoubtedly stems from the fact that fluvial landscapes are a snapshot of complex and multifactorial processes.There are at least three sources of complexity.First, landscapes are carved into a heterogeneous lithologic framework that itself evolves as a function of space and time.This intrinsic and changing heterogeneity undoubtedly influences the values of v, m and n, which in turn affect channel erosion.Secondly, the value of P, which has a controlling influence upon channel incision, is indirectly moderated by climate.Spatial variations of climate and precipitation are largely determined by latitude but temporal variation is problematic because it is so difficult to measure.The relationship between climate per se and fluvial discharge is also not straightforward.Thirdly, and perhaps most importantly, landscapes are externally forced by regional tectonic uplift whose functional form is poorly known and usually ignored by fluvial geomorphologists.
A range of strategies have been employed to circumvent these difficulties.Carefully designed and circumscribed field studies of individual monolithologic catchments where tectonic forcing is reasonably well known are often employed (e.g., Castillo et al., 2013;Hayakawa & Matsukura, 2003;Whittaker et al., 2007).Unfortunately, these studies are limited in terms of spatial and temporal scalings.They also suffer from inadequate information about tectonic forcing.An alternative approach exploits numerical experiments which simulate the development of fluvial landscapes (Nones, 2020).The underpinning algorithms are usually predicated upon the stream power formulation but they are also capable of emulating critical aspects of drainage development, including catchment growth, drainage divide migration, and stream piracy (e.g., Badlands: Salles and Hardiman (2016); Landlab: Hobley et al. (2017)).The principal drawback of these algorithms is that it is not yet possible to adapt them as the basis of inverting natural landscapes, although closed-loop modeling has been exploited in an attempt to bridge this gap (O'Malley et al., 2021).Scaled laboratory experiments of landscape evolution have also been carried out (Bonnet, 2009;Hasbargen & Paola, 2000).These experiments utilize, for example, a paste consisting of pure silica grains to ensure that vertical angle of rest and water infiltration are negligible.Controlled tectonic forcing is imposed and precipitation is generated with water droplets which have diameters <10 μm to minimize splash dispersion.Carefully controlled large-scale landscape experiments have been designed and analyzed with a view to addressing these shortcomings (Paola et al., 2009).Despite their qualitative usefulness, scaling distortions mean that both forms of experiment are not accurate analogs of nature (e.g., Densmore et al., 1997;Holland & Pickup, 1976;Lague et al., 2003;Paola et al., 2009;Schumm et al., 1987).
For these reasons, Tucker (2009) proposes that "natural experiments provide a means of testing landscape evolution theory on the large space and time scales to which that theory applies."His motivation is the development of a strategy which can address the difficulties described above "by taking advantage of experiments that nature has run for us."Tucker (2009) proposes four criteria for a natural experiment: (a) minimal variation in all but one factor (e.g., regional uplift, paleoclimate); (b) good temporal constraints; (c) well-characterized processes; and (d) appropriate resolution of topographic data.Here, we describe and analyze five buried transient landscapes which we believe fulfill these Tucker criteria.Our purpose is to demonstrate that these uniquely preserved surfaces share many of the characteristics of fluvial landscapes and that they can be exploited, as described by Tucker (2009), to yield quantitative insights into fluvial geomorphic processes and their drivers.We emphasize application of an inverse modeling strategy based upon solving a generalized stream power formulation to explore these processes and drivers but we acknowledge that many other forms of fluvial geomorphic analysis can, and should, be applied to these natural experiments.

Regional and Stratigraphic Context
During the Cenozoic Era, the Icelandic mantle plume has played a dominant role in controlling vertical motions throughout the North Atlantic region (Figure 1; Parnell-Turner et al., 2017).At the present day, residual depth measurements of 1-2 km in conjunction with surface wave tomographic models and long wavelength free-air gravity anomalies demonstrate that this plume is centered upon Iceland and has a radius of ∼1,000 km (Hoggard et al., 2017).
There is burgeoning evidence that plume activity fluctuates during the Cenozoic Era.The most convincing observations concern the diachronous geometry of Neogene V-shaped ridges and troughs that straddle the Reykjanes Ridge (Vogt, 1971; N. J. White & Lovell, 1997).Combined geophysical and geochemical observations at, and adjacent to, the zero-age ridge axis suggest that these V-shaped ridges and troughs are generated by temperature fluctuations of ±25-30°C within the plume head that have a periodicity of 3-5 Myrs (Parnell-Turner et al., 2017).Magmatic and sedimentary records imply that the Icelandic plume initiated at ∼64 Myrs.During the Paleogene period, it is inferred that temperature fluctuations also occurred since post-rift marine strata along the continental margins that fringe the North Atlantic Ocean are repeatedly punctuated by erosional surfaces  et al., 2003;Dam, 2002;Plink-Björklund & Steel, 2006;Stoker et al., 2012;Stucky de Quay et al., 2017); white circle = putative center of Icelandic plume at 55 Myrs (R. S. White & McKenzie, 1989).indicative of regional uplift and subsidence (Figure 1; Mudge & Jones, 2004).These fragmentary surfaces occur synchronously over a wide region and are indicative of periodic uplift and subsidence that is attributable to Icelandic plume activity.During the Neogene period, there is also evidence for regional uplift and subsidence events which may also be triggered by the Icelandic plume (Schoonman et al., 2017).
Within the Faroe-Shetland basin, post-rift marine strata are repeatedly punctuated by irregular reflective surfaces that have been mapped throughout a three-dimensional (3D) seismic reflection survey that extends over 2 × 10 4 km 2 (Figure 2).This merged survey is calibrated using 36 boreholes such that lithologic, biostratigraphic and depositional constraints can be carefully integrated with reflective sequence mapping.The geodynamic and paleoclimatic significance of these spectacular ephemeral landscapes is discussed in detail by Conway-Jones and White (2022).Given their unique preservation, such landscapes represent useful natural experiments that fulfill the four criteria laid down by Tucker (2009), enabling key aspects of fluvial geomorphology to be examined.In the Faroe-Shetland basin, our specific approach is threefold.First, we present and describe five transient landscapes.Secondly, we carry out a preliminary geomorphic analysis of four of these landscapes with a view to demonstrating that realistically complex channel and catchment morphologies can exist on short timescales, despite the uniformity of substrate lithology.Finally, we show that spectral analysis, slope-area measurements and inverse modeling of extracted river profiles yield helpful geodynamic insights.In this way, geomorphic processes that are otherwise difficult to study by analyzing present-day landscapes alone can be investigated.

Ancient Buried Landscapes
The chronostratigraphic context of these five buried landscapes is shown in Figure 3 (Conway-Jones & White, 2022;Shaw-Champion et al., 2008).Each landscape is carefully traced across the Faroe-Shetland 3D MegaSurvey, which was generously provided by Petroleum Geoservices (PGS) Ltd (Figure 2).Details of seismic acquisition and signal processing are summarized by Conway-Jones and White (2022).The dominant acoustic frequency is ∼50 Hz which means that the vertical and horizontal resolution of the migrated seismic reflection volume is ∼10 m.In comparison, the spatial resolution of digital topography recovered from SRTM observations is about 30 m.
On any given two-dimensional vertical slice through this seismic reflection volume, different rugose reflective surfaces are visible that cut across and truncate underlying strata (Figure 4).These subtle and irregular surfaces have variable reflectivity, which means that they are difficult to automatically track by algorithm.Instead, accurate mapping is best achieved if these surfaces are hand-picked on every in-line and cross-line of the survey (Conway-Jones & White, 2022).The five mapped landscapes are labeled Q to U (Figures 3, 5, and 6).

Corrections
Once a given surface has been mapped, four corrections of decreasing importance are carried out to ensure that detailed landscape relief is accurately recovered.The technical details of each correction are described and implemented by Shaw-Champion et al. (2008) and by Conway-Jones and White (2022).
First, post-depositional folding of Cenozoic strata is geometrically removed.For the two oldest landscapes Q and R, the simplest strategy is to flatten the seismic reflection volume on the top of the Balder Formation, which largely consists of tuffaceous material.The upper half of this formation is inferred to have been deposited uniformly on top of post-rift strata with minimal bathymetric relief (Figure 4b).For landscapes S-U, the volume was flattened on local deltaic topset surfaces.Secondly, surface relief is converted from two-way travel time, t, into depth, z.This time-to-depth correction is the most important one.Checkshot surveys, which yield t as a function of z, are available for 45 boreholes that penetrate the seismic reflection volume.A compilation of measurements is used to constrain an optimal relationship between t and z (Shaw-Champion et al., 2008).In this way, each landscape was accurately converted from time into depth.Thirdly, the flattening effect caused by compactional loading of the stratigraphic pile is removed.Compaction affects the amplitude of topographic features but it does not materially affect the shape of a recovered landscape.This correction is carried out using a standard porositydepth relationship that is calibrated using information from regional borehole logs (Shaw-Champion et al., 2008).In this way, original thicknesses are calculated from the seismically observed compacted thicknesses.Finally, the sediment fill of each landscape is isostatically converted into an equivalent air-loaded layer.It is difficult to precisely quantify uncertainties.However, given the scale of these landscapes and the nature of each correction, uncertainty in depth is dominantly at long wavelengths which means that the geomorphologic details on length scales of tens to hundreds of meters are well resolved.

Landscape Geomorphology
The stratigraphic position and regional fragmentary evidence for these different landscapes are presented by Shaw-Champion et al. (2008) and by Conway-Jones and White (2022).Biostratigraphic constraints confirm that each landscape is sandwiched between marine strata with evidence of terrestrial exposure (e.g., lignite deposits, pollen spores, carbonaceous material).Note that an older Landscape P exists which is dramatically exposed in two dimensions along the coast of West Greenland (Figures 1 and 3).A suite of animated aerial flights through each of the five landscapes is provided in Appendix A.

Landscape Q
Landscape Q is the oldest recovered surface within the Faroe-Shetland basin (Figure 5a).This landscape is the most difficult to map with accuracy, partly due to its depth and partly due to the fact that it is overlain by brightly reflective and irregular surfaces that have been interpreted as lava flows.Nevertheless, Landscape Q, which has a relief of up to 310 m, is characterized by significant and easily identifiable geomorphic features.A prominent east-west channel cuts deeply into the underlying Vaila Formation, which consists of marine mudstones and sandstones.The flanks of this sinuous channel are incised by a series of steep north-south tributaries which have interlocking spurs.Both the main channel and these tributaries exhibit both short (i.e., tens of meters) wavelength steps in slope (i.e., knickpoints) and longer (i.e., hundreds of meters) wavelength excursions (i.e., knickzones).This fluvially eroded landscape is infilled by the Kettla Tuff Member of the Lamba Formation, which contains lignite fragments (Knox et al., 1997).Biostratigraphic analysis demonstrates that sub-aerial exposure initiated at 59.15 ± 0.45 Ma and lasted ≤2 Ma (Conway-Jones & White, 2022;Speijer et al., 2020).Since it is more deeply buried and thus less well resolved, Landscape Q is not considered any further in this contribution.

Landscape R
The bulk of the Lamba Formation was deposited during Thanetian times.This formation consists of a prograding deltaic sequence that contains benthic foraminifera, corroborating the transient nature of Landscape Q. Divergent sets of clinoformal surfaces that are ∼400 m high indicate the existence of two deltaic lobes (Shaw-Champion et al., 2008).There is no evidence for topset surfaces.Instead, clinoformal surfaces are truncated by an irregular unconformity, referred to as Landscape R which lies above Landscape Q (Figure 5b; Shaw-Champion et al., 2008).This striking landscape is the most extensive and best preserved that we have encountered in this basin.Landscape R was sub-aerially exposed at 55.80 ± 0.8 Ma for a duration of <2 Ma, based upon the first and last appearance datums of Apectodinium augustum (Knox et al., 1997;Mudge & Bujak, 2001).It was previously mapped by Smallwood and Gill (2002), Shaw-Champion et al. (2008), andHartley et al. (2011).Conway-Jones and White (2022) revised and augmented mapping of this landscape using additional seismic reflection surveys (Figure 5b).Walker et al. (2022) show that Landscape R continues further northeast within the Faroe-Shetland basin.A coeval landscape has also been mapped along the northwestern flank of the North Sea basin (Stucky de Quay et al., 2017).
Detailed mapping of Landscape R reveals a dendritic drainage pattern, which is dominated by a large-scale meandering channel up to ∼60 km long.
An extensive network of at least eight steep-sided tributaries drains into the main channel.Along the northwestern edge, a broad and flat surface, interpreted as a flood plain, is visible.It culminates in a steep slope, interpreted as a sea cliff, with a height of ∼250 m.Spatial distribution of knickzones along the main channel and its tributaries reveals the existence of three tiers of interpreted peneplanation at elevations of 200, 450, and 750 m measured with respect to the base of the sea cliff.Landscape R has an average relief of 500-750 m.It is mostly infilled by the Flett Formation, which contains pollen assembles and lignite horizons (Mudge & Bujak, 2001;Shaw-Champion et al., 2008).This terrestrial formation is overlain by tuffs and mudstones of the Balder Formation, indicating a rapid return to marine conditions (Figure 3).

Landscapes S and T
The Balder Formation is overlain by the Stronsay Group, which was deposited in Ypresian times.The lower half of this chronostratigraphically undivided group consists of northward prograding clinoformal surfaces that are up to 100 m high and dominately comprise marine sandstones and mudstones (Figure 4b; Knox et al., 1997).The upper surfaces of two distinct deltaic sequences are erosional unconformities, designated as landscapes S and T (Conway-Jones & White, 2022).Both landscapes are infilled with carbonaceous-rich and tuffaceous deposits that contain lignite fragments and pollen, which are indicative of sub-aerial exposure with coeval volcanism (Figure 3).Landscape S was incised between 54.3 and 53.0 Myrs.Landscape T was incised between 53.0 and 51.7 Myrs.
Landscape S is dominated by a single large canyon with an asymmetric planform that curves southwestwards away from its mouth.This canyon is 25 km long and is characterized by steep slopes (Figure 5c).It has numerous short steep tributaries that emanate from a regionally extensive surface, which is interpreted as a peneplain.The average topographic relief is ∼200 m.Directly behind a minor sea cliff, this canyon is deeply incised into underlying deltaic deposits, generating prominent eroded spurs.The northfacing slope of this landscape has a less eroded appearance and is characterized by short steep catchments.Two or three tiers of flat surfaces are discernible, although they are much less prominently manifest than for Landscape R.
The frontal aspect of Landscape T is characterized by numerous steep catchments that are uniformly distributed along a 25 km-long northeast-facing slope (Figure 5d).Its northeast-facing slope has a more equilibrated geomorphic expression with larger, more eroded out, catchments.

Landscape U
During Early to Middle Eocene times, numerous submarine unconformities occurred within the upper half of the Stronsay Group (e.g., Stoker et al., 2018).At least some of these unconformities comprise buried landscapes.
Here, we focus on the most significant and continuous sub-aerially exposed surface, Landscape U, first identified by Robinson et al. (2004) who refer to this landscape as the Meoc2 surface (see their Figure 6c).It is carved into a clastic deltaic deposit and its infill contains lignite fragments.The precise age of this surface is debated but we infer that sub-aerial exposure occurred at 38.8 ± 1 Myr, based upon the appearance of Heteraulacysta?leptalea, Cerebrocysta bartonensis and Phthanoperidinium regalis as well as other dinoflagellate cysts and foraminifera (Speijer et al., 2020;Stoker et al., 2013).
Landscape U is characterized by a steep-sided canyon that curves toward the west (Figure 6).Its well-preserved drainage network is divided into three prominent tributaries that deeply incise into an extensive surface of low gradient.Landscape U has a relief of ∼250 m and it is overlain by definitively marine strata of the Upper Stronsay Group.

Natural Geomorphic Experiments
These transient landscapes record discrete episodes of regional uplift, sub-aerial exposure, fluvial erosion and rapid subsidence that are plausibly triggered by sub-plate temperature fluctuations which horizontally advect away from the central conduit of the Icelandic plume (Conway-Jones & White, 2022).We propose that the degree of preservation recorded by these landscapes is both remarkable and geomorphologically significant.The causes of this quality of preservation are threefold.First, these landscapes developed within a paleoenvironmental "Goldilocks Zone" which enhanced their preservation potential.If pre-existing water depths had been greater than several hundred meters, it is unlikely that this continental margin would have been sub-aerially exposed at all.If, on the other hand, pre-existing water depths had been negligible, prolonged sub-aerial exposure would have given rise to wholesale excision.Secondly, regional uplift was short-lived, which means that these transient landscapes were rapidly submerged.Thirdly, subsequent sedimentation rates were rapid so that submarine erosion did not occur.This particular combination of water depth, transient uplift and sedimentation rate has ensured preservation of these delicately sculpted surfaces.Figure 1 identifies other locations throughout the North Atlantic region where similar Cenozoic transient buried landscapes have been recorded.
It is recognized that there is a pressing need for geomorphologic studies that combine the control and functional simplicity of analog experiments with more physically realistic spatial and temporal scaling.Tucker (2009) proposes four criteria for natural experiments: sufficient control of different elements (e.g., lithology, precipitation); temporal limits; process characterization; and well-resolved topography.A number of published case studies fulfill one or more of these conditions but no single study satisfies all four criteria.Here, we propose that the five buried landscapes from the Faroe-Shetland basin be regarded as suitable, and therefore useful, natural experiments.
The first criterion can be more strictly defined as an ability to control all but a single element.A significant advantage of 3D seismic imagery calibrated by boreholes is that the stratigraphic substrates into which each landscape is incised are known.For example, Landscape R is cut into the Lamba Formation which consists of a deltaic sequence dominated by mudstones with subordinate siltstone and occasionally sandstone layers.Landscapes S, T, and U are all cut into the Stronsay Group, which is also dominated by mudstones with occasional sandy intervals.These borehole-based observations, combined with the general absence of bright seismic reflections, demonstrates that both the Lamba and Stronsay substrates can be regarded for the purposes of this study as being lithologically homogeneous.
Biostratigraphic constraints show that all five landscapes are exposed for short durations of less than O(1) Myrs, during which the variation of precipitation as a function of time appears to be dominated by Milankovitch orbital cyclicity (Conway-Jones & White, 2022; Westerhold et al., 2020).Barnet et al. (2019) analyzed a 14.75 Myrs long orbitally tuned isotopic record between 67 and 52 Myrs, a period of time that encompasses the buried landscapes considered here.They conclude that eccentricity, through modulation of precession, is the dominant pacemaker of climatic variation during Paleogene times.The coolest conditions occur between 61 and 58 Myrs, which is followed by a long-term warming trend, triggered by rising pCO 2 levels.All five buried landscapes occur during this warming trend when it is reasonable to assume a weighted average of climatic conditions (Tucker, 2009).Significantly, it has been quantitatively demonstrated that fluvial drainage patterns and the geometries of longitudinal river profiles are insensitive to temporal precipitation variations that have periodicities of <1 Myr (O' Malley et al., 2021).This insensitivity is unsurprising since the solution to the stream power equation is an integral, which means that average (i.e., not instantaneous) discharge along a river channel is more important on long timescales (Goren, 2016;Paul et al., 2014).The landscapes presented here also have a sufficiently small aerial extent (i.e., ≤100 km) such that biota and climate would not be expected to spatially vary.
Secondly, temporal control is an essential criterion.The landscapes presented here satisfy this requirement in two important ways.The initial conditions at t = 0 are known.Each landscape has developed on gently dipping topset deposits at the top of a deltaic system (Figure 4).In each case, biostratigraphic observations provide well-resolved constraints on the duration and absolute timing of sub-aerial exposure (Conway-Jones & White, 2022).As a result, it is possible to estimate the rate of uplift and erosion.This high degree of temporal control is akin to analog experiments.
Thirdly, active processes (i.e., external forcing) must be characterized.Previous studies demonstrate that transient tectonic uplift is the primary control (Hartley et al., 2011;Shaw-Champion et al., 2008).For each landscape, Conway-Jones and White (2022) used a combination of stratigraphic constraints and geomorphic relief to calculate the amplitude of regional uplift.The largest amplitude is recorded by Landscape R which experienced 400-600 m of regional uplift within ~1 Myr, which yields an average uplift rate of 0.5 mm yr 1 (Hartley et al., 2011).Landscapes Q and U each underwent about 200 m of regional uplift.Landscapes S and T were both regionally uplifted by about 100 m.The detailed biostratigraphic constraints used to determine the amplitude and duration of uplift for each buried landscape are summarized by Conway-Jones and White (2022).These inferred uplift rates are significant, suggesting that regional tectonic forcing is the dominant forcing process in terms of the modified stream power formulation.
Finally, Tucker (2009) highlights the requirement for well-resolved digital elevation models.The spatially migrated 3D seismic reflection volume exploited in our study has a vertical and horizontal resolution of about 10 m, which is principally governed by frequency content of the acoustic source (Conway-Jones & White, 2022).
Given that these ancient landscapes are buried beneath ∼1 km of sedimentary rocks in water depths of ~2 km, it is remarkable that these landscapes can be mapped at resolutions that are comparable to those of digital elevation models exploited in geomorphic studies of the Earth's surface.We conclude that each of these buried landscapes represents a natural experiment which should enable different geomorphic hypotheses and models to be tested on appropriate timescale and length scales.Each landscape is carved into very a similar substrate but tectonically forced by a different amplitude of uplift, which helps to shed light on the quantitative relationship between tectonic forcing and geomorphic expression.

Geomorphic Analysis of Landscape R
The exceptional degree of control afforded by these natural experiments suggests that they can be used to test fundamental, well-established geomorphic principles.Here, we show that buried landscapes adhere to wellknown rules of channel geometry such as width, degree of incision, sinuosity and curvature (Fisher et al., 2013).The initial focus of our attention is Landscape R, which is the most spatially extensive and best resolved landscape with a clearly expressed drainage network.The other landscapes are exposed on a smaller scale and are less dramatic but have been similarly analyzed with equal success.
The fluvial drainage network of Landscape R can be visualized in different ways (Figure 7).The planform of this landscape demonstrates that drainage is dominated by a large-scale meandering channel with eight significant tributaries (Figure 7c).The detailed geometry of this channel is easily observed by calculating the spatial gradient, which reveals varying patterns of incision along major and minor channels (Figure 7b).Along the meandering channel, the degree of incision systematically varies along each meander loop with evident slope asymmetry at meander apices.This channel narrows dramatically upstream.In contrast, minor tributaries are generally flatbottomed and less dramatically incised.These tributaries have characteristic branching and lobate minor catchments.Five major drainage basins together with about 13 minor catchments are identified on Landscape R (Figure 7c).
It has been argued that rivers are effectively "tape recorders" of external forcing since the geometry of longitudinal river profiles contains information of regional uplift rate, preserved as knickzones which retreat upstream (Whipple & Tucker, 1999).The time taken for an uplift signal that is inserted at the mouth of a river to advect upstream to a given position is given by the Gilbert Time, τ G , where A is upstream drainage area, v is the velocity of knickzone retreat, and m = 0.3 is an empirical erosional parameter (Roberts et al., 2012).Thus the speed at which an uplift signal travels depends upon a combination of v and A m where the latter is a proxy for stream power which gradually diminishes upstream.Figure 7d shows how τ G varies for Landscape R. In this example, a value of v ∼ 10 is chosen to ensure that the maximum value of τ G approximately equals the maximum exposure duration of the landscape.This value of v can be regarded as a minimum.Different values of m do not materially affect our results.Provided that the regional uplift signal is inserted at the putative paleocoastline, short rivers in proximal catchments can only constrain youthful uplift events whereas longer rivers can constrain uplift events that are as old as landscape exposure.At the major drainage divides, τ G → ∞ as A m → 0 (Paul et al., 2014).The spatial pattern of τ G is geomorphically self-consistent and realistic, which implies that digital recovery of these ancient river channels is robust.Note that comparative values of τ G are less reliable toward the edges of the landscapes as would be the case for metric analysis of present-day landscapes.The principal drawback of τ G analysis is that U is assumed to be spatially invariant.We note in passing that other popular metric measures, such as χ, suffer from the same serious limitation.It is considerably better to pose and solve the general inverse problem which does not make limiting assumption about the spatial and temporal variance of regional uplift.
A sequence of four vertical cross-sections that traverse the main channel visible on Landscape R is presented in Figure 8.As expected, channel width increases downstream (Fryirs & Brierley, 2012).Channel depth also increases downstream but to a proportionally lesser extent than channel width in agreement with global observations (Leopold & Maddock, 1953;Moody & Troutman, 2002).Along channel segments with significant long wavelength slopes (i.e., knickzones), channel width decreases and channel depth significantly increases, which is consistent with experimental and theoretical studies (Figures 9c, 9f, and 9h; Gardner, 1983;Whipple & Tucker, 1999).
Overall, this central channel has a sinuosity of 1.55, and is classified as meandering (Figure 10; Fisher et al., 2013;Leopold & Wolman, 1957).Meandering channels follow rules of self-similarity, where meander length is typically ∼11 times the channel width (Leopold & Wolman, 1960).The central channel is consistent with this prediction since channel width is 850 ± 150 m and meander length is 8.5 ± 0.5 km (Figure 9).Furthermore, the radius of curvature is ∼2.5 times channel width, consistent with expected values of 2-3 (Leopold & Wolman, 1960).
Vertical cross-sections located at significant meanders clearly illustrate channel asymmetry, consistent with observations of slope angle (Figures 7b and 10).Each cross-section exhibits what can be interpreted as a steep outer slope, an off-center thalweg, and a gently dipping point bar (Figures 10c-10e).Since the highest shear stresses are likely to be encountered along an outer concave bank (i.e., the cut bank), channel erosion can be expected to be enhanced which would generate a high slope angle and promote meander migration.At the same time, inward basal flow of the vortex is expected to sweep sedimentary particles toward the convex bank of a meander, adding to the low-angle point bar.Thus, in addition to fulfilling the criteria for a natural experiment, this  preliminary geomorphic investigation suggests that the drainage channels of Landscape R share many of the complex features described by examining present-day landscapes and laboratory experiments.

Spectral Analysis of Landscape R
There is longstanding debate concerning how fluvial channels acquire their longitudinal geometry (e.g., Dietrich et al., 2003;Howard et al., 1994;Rosenbloom & Anderson, 1994;Weissel & Seidl, 1998;Wolpert & Forte, 2021).This debate is centered upon the relative significance of regional uplift, which acts to elevate a channel, and erosional processes, which act to lower a channel.Roberts et al. (2019) showed that river profiles are dominated by long wavelength signals (i.e., ≳100 km) by calculating spectral power.Distance-averaged power spectra typically have a slope of 2 (i.e., red noise) at longer wavelengths and a slope of 1 (i.e., pink noise) at shorter wavelengths.They argue that these spectra imply long wavelength regional uplift is the primary external forcing that drives the macroscopic evolution of fluvial landscapes.The cross-over transition from red to pink noise at wavelengths ≲100 km suggests that complexity generated by lithologic changes, heterogeneous erodibility, and non-linear fluid flow instabilities arising from the Saint-Venant equations becomes significant.Conversely, emergent simplicity at long wavelengths suggests that landscape evolution can be quantified using a non-linear Langevin formulation with a noise driver which is referred to as the Kardar-Parisi-Zhang equation (Kardar et al., 1986).
Natural geomorphic experiments are an important means for testing these concepts since external forcing, duration and lithology are well constrained.Here, we spectrally analyze the drainage network extracted from Landscape R since it is considerably larger than the other landscapes and spectral power can be calculated over a significant range of wavelengths.Drainage extraction is carried out using a flow-routing algorithm (Tarboton, 1997).Power is calculated as a function of both distance and wavenumber using a wavelet transform (Roberts et al., 2019;Torrence & Compo, 1998).The wavelet transform, W x (s), is calculated discretely by summing wavelets over a range of horizontal scales, s, along a river profile where ψ is the mother wavelet, z x′ is the elevation at a given distance x′ along the profile, and δx = 200 m is the sampling interval (Torrence & Compo, 1998).A Morlet wavelet is assumed such that where the subscript zero refers to a wavelet that has not yet been scaled to s, ω 0 is a dimensionless frequency, and i 2 = 1.We use ω 0 = 6 but it has been shown by Roberts et al. (2019) that other values of ω 0 yield similar power spectra.The wavelet power spectrum is given by where the distance-average power spectrum, φ(s), is In order to compare wavelet and Fourier spectra, distance is converted into wavenumber (Torrence & Compo, 1998).Edge effects are minimized by mirroring river profiles about the z and x axis seven times (Roberts et al., 2019).
Spectral analysis of the eight principal drainage channels extracted from Landscape R is shown in Figure 11. Figure 11a demonstrates that whilst a continuous wavelet transform successfully recreates the river profile of the main channel, a good fit is also obtained when the 50 shortest wavelengths (i.e., <23.9 km) are omitted.This match demonstrates that spectral power of these river profiles resides at the longest wavelengths (Figure 11b).At short wavelengths, spectral power is localized at three prominent knickzones.The distance-average power spectrum, φ(s), for the main channel is presented alongside a moving average of the Fourier transform (Figure 11c).Both methods corroborate the fact that power resides at the longest wavelengths.Finally, the mean distance-averaged power spectrum for all eight channels is shown in Figure 11d.The slope of this normalized spectrum is 1 (i.e., pink noise).The results presented here are in close agreement with those of Roberts et al. (2019).Power is concentrated at the longest wavelengths and only extends to short wavelengths where sharp increases in slope occur.Roberts et al. ( 2019) find that at wavelengths of less than about 100 km, slopes of distance-averaged power spectra consistently change from 2 (i.e., red noise) to 1 (i.e., pink noise).Landscape R is too small to record wavelengths >100 km but clear spectral evidence for pink noise has significant implications.At the longest wavelengths, power spectra of longitudinal river profiles are expected to be dominated by red noise because the amplitude of topography scales with wavelength (Newman & Turcotte, 1990;Pelletier, 1999).The observed cross-over transition from red to pink noise at wavelengths of ∼100 km implies that self-similarity starts to be influenced by complex hydraulic and erosive processes that could be characterized by white, or possibly blue, noise (Birnir et al., 2001;Smith, Birnir, & Merchant, 1997;Smith, Merchant, & Birnir, 1997)

Development of Buried Landscapes
Landscape R exhibits the most extensive and best-resolved drainage network (Shaw-Champion et al., 2008).
Nevertheless, geomorphic analysis of landscapes S, T and U show that they have broadly similar characteristics .Landscape Q has not been analyzed in any detail because seismic reflection imaging has been compromised by extrusive and intrusive igneous features which degrade its spatial resolution.Time-to-depth and decompaction calculations are also subject to greater uncertainties.The drainage planform of Landscape S can be clearly picked out on a spatial gradient map (Figure 12b).This landscape is dominated by one large catchment whose main channel exhibits steeper incision along its northern edge (Figure 12c).There are six minor catchments.The spatial distribution of τ G demonstrates that the drainage planform behaves in accordance with the stream power formulation away from the edges of the mapped landscapes (Figure 12d).
Although Landscape T is located directly above Landscape S and has an almost identical geologic setting along the crest of a minor delta, its drainage planform is quite different (Figures 13a and 13b).First of all, it is compartmentalized into two separate basins each of which has four small catchments (Figure 12c).The northern basin is deeply incised with numerous elongated and eroded spurs which implies that it is approaching dynamic equilibrium.There is a steep, possibly eroding, cliff at the back of this basin bounding an undissected flat surface.The eastern basin is possibly younger since there is less evidence for isolated spurs.The two larger catchments have remarkably linear channels.The cliff at the back of this basin has a clearly defined sharp edge.
The spatial distribution of τ G suggests that both basins potentially record signals of similar duration (Figure 13d).
Landscape U is much more deeply incised than the older landscapes (Figures 14a and 14b).A trilete drainage pattern with three major catchments and three minor catchments is incised into an extensive surface of low relief which has an average elevation of ∼200 m (Figures 14a and 14c).Major catchments are dominated by dendritic drainage and flow into a steep-sided canyon that meanders in places.The spatial distribution of τ G is systematically organized (Figure 14d).This landscape generally has the steepest slopes which could imply a greater degree of disequilibrium.
All of these buried transient landscapes share important geomorphic characteristics-dendritic drainage patterns, steep asymmetric slopes on either side of any major channels, incised meanders, and long wavelength knickzones.From a tectonic perspective, each landscape should contain a recoverable record of regional uplift.Different studies have attempted to extract information about regional uplift histories from suites of longitudinal river profiles (e.g., Fox et al., 2014;Roberts & White, 2010;Rudge et al., 2015).The modified stream power formulation forms the basis of all of these studies which necessarily makes significant simplifying assumptions about rock erodibility and precipitation.Here, our intention is to exploit these stratigraphically constrained natural geomorphic experiments to investigate the quantitative relationship between tectonic forcing and landscape development.
Our approach is divided into two steps with the aim of demonstrating the particular importance of these buried landscapes as natural experiments.In the first step, complete drainage inventories of longitudinal river profiles are digitally extracted from each buried landscape.Even from a qualitative perspective, these suites of river profiles yield useful insights into the varying geomorphic expressions that result from tectonic forcing.A significant advantage of buried landscapes is that the geometry of both substrate and infill can be examined.In the second step, an inverse modeling algorithm is used to calculate uplift rate histories by minimizing the misfit between observed and calculated river profiles.Our goal is show that systematic differences in recovered uplift histories reflect geomorphic expression.In summary, a combination of lithologic and temporal constraints underline the geomorphic significance of this suite of landscapes as natural experiments.

Drainage Inventories
Extracting drainage networks from digital elevation models is a well-known operational procedure.Landscapes are first smoothed to remove anomalous spikes and pits.Secondly, drainage networks are identified using a standard flow-routing algorithm (Tarboton, 1997).Individual channels are identified by setting the minimum upstream drainage area to 10 4 m 2 .The extracted drainage channels and catchments are shown in Figures 7c, 12c, 13c, and 14c.Comprehensive suites of longitudinal river profiles from each landscape are presented in Figure 15.
In each case, the main river profile as well as about 10 tributary profiles are shown.Aside from minor artifacts, the recovered drainage planforms and profiles are well-resolved and appear to possess geomorphic realism.
All of the extracted river profiles appear disequilibrated since there are significant long wavelength slope excursions (i.e., knickzones; Figures 15b, 15f, 15j, and 15n).However, the degree and form of this disequilibration does vary from landscape to landscape.On Landscape R, nine river profiles display three spatially coherent knickzones at different horizontal ranges.Numerous minor and largely uncorrelated knickpoints are also evident.
The major knickzones are separated by low-relief surfaces, possibly peneplains, which occur at elevations of 250, 450 and 700 m (Figure 15n).All three knickzones have been cut into the lithologically homogeneous Lamba Formation, which implies that they were generated by regional base-level changes.Following ∼1 Myr of terrestrial exposure whose duration is stratigraphically constrained, Landscape R was rapidly infilled by the Flett and Balder formations (Figures 15o and 15p).
Nine river profiles were extracted from Landscape S.These profiles appear to be more disequilibrated since they have longer reaches that are more obviously convex upward (Figure 15j).Broader, but less obvious, knickzones occur between flatter reaches at elevations of 60 and 100 m.River profiles are cut into the lithologically uniform Stronsay Formation and the landscape is rapidly infilled by onlapping of the upper half of the same formation (Figures 15k and 15l).Reflectivity of the substrate does not correlate with loci of knickzones.Landscape T sits immediately above Landscape S but its drainage planform and the geometry of its 11 river profiles are quite different (Figures 15e and 15f).Each of the two catchments has a distinctive profile geometrythe northern catchment has five concave upward profiles with subtle knickzones between flatter reaches at elevations of 30 and 60 m while the eastern catchment has seven highly linear profiles with equally subtle knickzones between flatter reaches at elevations of 50 and 100 m.River profiles are cut into, and filled by, the largely lithologically featureless Stronsay Group (Figures 15g and 15h).
The youngest Landscape U appears to have the most disequilibrated river profiles-they are clearly convex upward (Figures 15a and 15b).Three knickzones occur between flatter reaches at elevations of 50, 80, and 150 m.Within the underlying substrate, there are significant changes in reflectivity but these changes are generally parallel to the river profiles (Figures 15c and 15d).

Regional Uplift Histories
Different methods are used to calculate regional uplift histories from fluvial landscapes.These methods are usually based upon a modified version of the stream power formulation given by Equations 1 and 2. Thus the geometry of a longitudinal river profile is controlled by regional uplift, U, and moderated by erosion, E, both of which can vary as a function of time, t, and distance, x.A simplified form of channel erosion can be employed where The values of v and m are determined from observational constraints.
A primary goal is to calculate U(x, t) which drives fluvial channelization and erosion.Different strategies are used to determine U as a function of time and space.One popular strategy exploits slope-area analysis by making simplifying assumptions about stream power.A second strategy determines U by inverse modeling.Here, we explore both of these strategies in order to shed light on the tectonic development of these particular buried landscapes.

Slope-Area Analysis
Numerous studies have quantitatively explored the links between uplift rate, channel steepness and upstream drainage area of a longitudinal river profile by imposing a steady state condition (i.e., ∂z/∂t = 0; Kirby & Whipple, 2012;Tarboton et al., 1989;Schoenbohm et al., 2004).This assumption means that Equations 1 and 2 can be rearranged so that By plotting observation of log(∂z/∂x) as a function of log(A), a linear regression yields a gradient of m/n and a y intercept of k s = (1/n) log(U/v) .In this way, U can be determined provided that the values of v and, say, m are known.
Here, we apply standard slope-area analysis to the main river profile from Landscape R since it is the longest and best resolved channel (Figure 16).This channel is discretized every 200 m and slope values are calculated using a spatial 9-point moving average along the channel.As is usually the practise, flat reaches of this profile, where ∂z/ ∂x = 0, are excluded.The results of slope-area analysis are shown in Figure 16b.As anticipated, calculated values of log(∂z/∂x) and of log(A) exhibit a broadly negative correlation.Linear regression of these values returns m/ n = 0.148 and U = 2.77 m Myr 1 , assuming that v = 10 which is the minimum recoverable value consistent with the duration of Landscape R. The recovered value of U is considerably smaller than that determined from stratigraphic constraints (Hartley et al., 2011).This value implies that Landscape R grows and decays over tens of millions of years which does not accord with stratigraphic constraints.Note that assuming a smaller value of v would lead to an even more reduced (and therefore unrealistic) value of U. If n = 1, the recovered value of m is 0.148.Alternatively, if m = 0.5, the recovered value of n is 3.4.
The reasons why slope-area analysis fails in this instance are two-fold.First, this approach is predicated upon the assumption of steady state which is unlikely to be correct for these natural experiments where stratigraphic constraints demonstrate that Landscape R underwent regional uplift, terrestrial exposure, and subsidence over a very short period of time (≤1 Myrs) during which fluvial equilibrium has demonstrably not been achieved.Secondly, slope-area analysis depends upon taking the derivative of discrete and noisy measurements and so is highly susceptible to the presence of even small uncertainties (Roberts et al., 2012).

Linearized Inverse Modeling
It has been shown that U can be recovered as a function of time and space under conditions where steady state is not assumed (i.e.∂z/∂t ≠ 0).Pritchard et al. (2009), Roberts and White (2010) are consistent with independent geologic observations.Their general strategy relies upon an integral solution which recovers U by minimizing the misfit between observed and calculated river profiles.Rudge et al. (2015) developed a linearized approach which uses the method of characteristic curves to solve Equation 1 for smooth histories of uplift rate.This methodology has been successfully applied to many modern fluvial landscapes (e.g., Conway-Jones et al., 2019;McNab et al., 2018;Tribaldos et al., 2017).The two appropriate integral solutions are and assuming boundary conditions x = x * , z = z * at t = 0, i.e. when landscape starts to form (12) and Uplift values, U, are discretized at spatial and temporal nodes with 1 km and 0.01 Myr intervals, respectively.The resultant elevation at a given position along the river profile, z, is given by where M is a matrix containing information required to calculate river profiles from a given uplift history (e.g., upstream drainage area, erosional parameters).Equation 14 can be inverted using a non-negative linear least squares approach, thus enabling U to be calculated from z.The matrix M is often under-determined which means that damping is required.A least squared approach is chosen such that the quantity is minimized where λ s and λ t are spatial and temporal smoothing parameters.The spatial smoothing matrix, S, is given by and the temporal smoothing matrix, T, is given by where s represents the two spatial dimensions.In this way, large numbers of river profiles can be simultaneous inverted to calculate a history of regional uplift as a function of time and space.

Parameter Calibration
A significant advantage of inverse modeling is the ability to calibrate erosional and smoothing parameters by repeatedly calculating the misfit between observed and calculated river profiles for different values of v (O' Malley et al., 2021;Rudge et al., 2015).Root-mean-squared (rms) misfit, χ 2 , is given by where J is the number of river profiles, I is the number of points along a given profile, K is the total number of points, z o and z c are the observed and calculated elevations, respectively, and σ = 10 m is elevation uncertainty, which equates to the resolution of the seismic reflection experiment.
A parameter sweep, in which m is systematically varied between 0.2 and 0.9 with an increment of 0.05, yielded a minimum misfit at 0.45 ≤ m ≤ 0.55 for all four buried landscapes (Figure 17a).This range of values is consistent with published estimates (Whipple & Tucker, 1999).The optimal value of v can be obtained from the duration of landscape exposure which is constrained by stratigraphic observations.The simplest approach exploits a scaling argument-the maximum value of τ G can be equated to the duration of sub-aerial exposure which yields 4 ≤ v ≤ 15.These estimates of v are greater than those obtained from modern fluvial landscapes probably because these buried transient landscapes are carved into unconsolidated (i.e., easily erodible) fine-grained sediments.Note that v is essentially a scaling parameter, which sets the pace of fluvial erosion.Our recovered range of values of v can be regarded as minima since some kinematic wave signals (i.e., knickzones) may have propagated all the way back to drainage headwaters with a concomitant loss of information.
Finally, the values of the smoothing parameters λ s and λ t are selected.Following Parker (1972), our goal is to seek the smoothest model which yields the lowest rms misfit by varying smoothing parameters λ s and λ t .Optimal values of λ s = 10 and λ t = 1 are obtained for χ 2 = 0.24-1.38 which occur at the elbow (i.e., change in slope) of the classic L-shaped curve (Figures 17b and 17c;Parker, 1994).

Calculated Uplift Histories
Inverse modeling of landscapes R-U recovers uplift rate histories which are presented in Figure 18.In each case, the misfit between observed and calculated river profiles has been minimized by varying uplift rate as a function of time and space (Figures 18a-18d).Residual misfits are generally close to or less than 1 and visual fits are good in that major knickzones are matched without overfitting minor knickpoints or any features that are smaller than a conservative spatial resolution of, say, 30 m.During the modeling process, initiation of regional uplift (i.e., the commencement of landscape emergence) is set at 0 Myrs.Model time increases toward the present day.In each case, geologic time can be obtained by consulting Figure 3 and Conway-Jones and White (2022).
The uplift rate history of Landscape R is recovered over a period of 2.1 Myrs, which is consistent with qualitative observations of major knickzone and low relief (i.e., probable peneplanation) surfaces (Figures 18e and 18i).An uplift rate of ∼0.2 mm yr 1 occurs between 1 and 1.4 Myrs.Between 1.5 and 1.8 Myrs, there is a second pulse of uplift with a rate of ∼1 mm yr 1 .A final pulse of uplift occurs between 2 and 2.1 Myrs which has a maximum rate of ∼4 mm yr 1 and is manifest as an inferred sea cliff located at the base of the buried landscape.This sea cliff probably represents significant and youthful base-level change.Cumulative uplift, which is obtained by integrating over the uplift rate history, reaches 625 m (Figure 18e).These results are consistent with those of Hartley et al. (2011).
Landscape S is half the extent of Landscape R but also records three phases of uplift at 0.2-0.35,0.39-0.55 and 0.55-0.58Myrs (Figure 18j).Each successive phase has a higher rate, culminating in a rate of 2 mm yr 1 .Cumulative uplift is 110 m (Figure 18f).Emergence of Landscape T, which has the smallest planform, is modeled over a period of 0.53 Myrs.Its drainage network also records three pulses of uplift at 0.1-0.28,0.31-0.4and 0.43-0.53Myrs (Figure 18k).The maximum recorded uplift rate is 0.8 mm yr 1 and cumulative uplift is 100 m (Figures 18g and 18k).There is evidence for spatial variation of U recovered from Landscape T, which reflects the fact that two discrete drainage basins were inverted (Figure 13c).Nevertheless, the uplift rate history is mostly consistent between these basins.Finally, inverse modeling of river profiles from Landscape U recovers three phases of uplift rate at 0.32-0.45,0.51-0.67,and 0.72-0.83Myrs.The most significant result concerns the peak uplift rate which reaches 2.6 mm yr 1 over a very short period of 0.1 Myrs (Figure 18l).Once again, the recovered spatial variation of U is negligible.All four buried landscapes are characterized by rapid uplift rates of 1-4 mm yr 1 and there is evidence that uplift rate histories occur in two or perhaps three phases.Stratigraphic constraints demonstrate that each of these landscapes is buried as rapidly as they are generated (Hartley et al., 2011).It is likely that this sequence of transient landscapes is generated by fluctuating behavior of the Icelandic mantle plume as postulated by Rudge et al. (2008) and by Conway-Jones and White (2022).Thus, these landscapes appear to indirectly record plume activity that is manifest by Neogene V-shaped ridges with the adjacent oceanic basins (Parnell-Turner et al., 2017; N. J. White & Lovell, 1997).In both cases, rapid uplift followed by equally rapid subsidence can be modeled by horizontal advective flow of a thermal anomaly beneath the lithospheric plate.Notwithstanding the overall commonality of landscape generation and decay, there are significant geomorphic differences from landscape to landscape which are manifest by the details of their drainage networks.These geomorphic differences may reflect minor changes in uplift rate history and underline the significance of this suite of landscapes as systematic natural experiments.

Conclusions and Future Work
Fluvial landscapes are complex multifactorial phenomena whose quantitative origins are much debated.Although it is now more widely appreciated that fluvial landscapes are externally forced by regional tectonic uplift, which varies as a function of space and time, the complexity of fluvial erosion at the hydrologic-lithologic interface makes recovery of this tectonic signal challenging.It is generally recognized that well-chosen natural experiments can play a useful role in testing landscape theory on sufficiently large space and time scales.Here, we identify, interpret and analyze a suite of buried transient landscapes, which fulfill the criteria for such natural experiments.Our results suggest that geomorphic response is primarily a function of the fluid dynamics of channel development rather than of lithology or precipitation which is an important insight.Despite the undoubted variation and complexity of this response, inverse modeling of river profiles extracted from four landscapes recovers uplift rate histories that are generally consistent with independent geologic observations about fluctuating activity of the Icelandic plume.Notwithstanding these simple insights, we recognize that the complex geomorphic details of these transient landscapes remain unexplained and are necessarily a topic of future research.

Appendix A
A suite of MP4 files are provided which present animated aerial flights through the five buried landscapes.The purpose of these animations is to illuminate, to enhance, and to clarify landscape descriptions provided in the main text.The following captions document the principal geomorphic features that are visible.In each case, elapsed times are quoted in seconds.Scale and orientation are provided by Figure 5.

A1. Landscape_Q_Fly.mp4
The oldest landscape has the poorest resolution due to a combination of depth and of predominantly volcanic infill.Nonetheless, the general outlines of the fluvial geomorphology of Landscape Q are clearly discernible.As the mouth of the main channel is approached, we can see a broad but incised valley which gently meanders upstream to the east.No significant sea cliff is evident.The sides of this valley are very steep, climbing onto a 200 m high, probably peneplained, surface characterized by numerous extensively eroded channels.There are excellent examples of fins and hoodoos, especially on the northern side.In the distance, a small number of tributaries emerging from an even higher surface are visible on both sides.At 28 s, this broad valley or channel dramatically narrows where a prominent knickpoint is visible, flanked by two spurs.Significantly larger spurs occur at elevation on either side of the undissected surface, especially on the southern side.At 32 s, the flight path veers northward to track the now much narrower channel which significantly meanders.By 37 s, we have emerged out onto a higher undissected surface.The flight ends with an eastward view across this peneplain toward a narrow incised channel which continues further into the distance.

A2. Landscape_R_Fly.mp4
This spectacular landscape is the most extensive and the best preserved.It also has the highest resolution, partly because the irregular and rugose surface is easier to map.The initial high aerial view dramatically reveals the meandering main channel and at least two extensive undissected surfaces.Its broad flood plain is flanked by steep escarpments that bound a 500 m high surface.At 18 s, there is an excellent view of the mouths of both the main channel and a subsidiary channel lying to the west.Both channels emerge onto a flat estuary which is bound by a steep sea cliff with an elevation of 200 m.This inferred sea cliff marks a broad inlet whose shape and orientation is consistent with the main channel.The flight path turns eastward at 24 s up the lower reaches of the main channel where a meander belt is clearly observed with a sinuosity index of ∼2.Vertical transects across the main channel demonstrate that individual meanders are incised.The outer, concave sides of these meanders have prominent river-cut cliffs (i.e., cut banks).At 29 s, a possible example of a meander cutoff is visible at the center of the view.Continuing upstream, the meandering loops become increasingly incised.Prominent and deeply incised tributaries occur on both valley flanks.At 43 s, there is an excellent view of a dramatic knickzone where the main channel narrows from a width of 5 km down to several hundred meters.To the southwest, a large deeply incised tributary is clearly visible.The flight path veers southwestward at 46 s up this incised tributary.In the distance, the main channel above the knickpoint broadens out into a wide valley.The final view to the southwest is of the highest undissected surface at an elevation of 600 m.

A3. Landscape_S_Fly.mp4
This landscape is less extensive and is characterized by two catchments, which differ considerably in terms of scale and shape.In the distance, a dramatic flat surface at 100 m is clearly visible and a continuous 50 m high sea cliff occurs in the foreground.At 4 s, the flight path turns northwestward in toward the smaller northern catchment.The clearest view of this catchment occurs at 14 s.The planform of the main channel curves from northeast toward the south.Its broad valley is not very incised but the many small tributaries that emerge from the high surface are much more incised.At 18 s, the flight path veers southeastward over a broad tilted surface where part of the drainage divide between the two catchments is visible.There is an excellent view of the incised mouth of the main channel at 20 s.The channel valley is narrow and flanked by steep to vertical escarpments that are 100 m high.Prominent eroded spurs, fins and hoodoos are clearly visible.At 24 s, the camera turns westward at which position the main channel narrows at a steep knickzone where there are numerous elongated spurs and fins.Above this knickpoint, the main channel occupies a broad flat valley bounded by moderately steep escarpments.
There is an excellent view of this valley at 26 s.The main channel narrows again at 28 s and then splits into a series of smaller tributary channels.The final view is of the highest, monotonous, undissected surface.

A4. Landscape_T_Fly.mp4
Landscape T resembles Landscape S to some extent in that the maximum relief is also about 100 m and there are two large north-facing catchment areas.The overall view shows a flat surface whose northern edge is dissected by incising drainage.At 9 s, the flight path turns toward the eastern catchment which is characterized by two largescale channels as well as other less well-resolved channels.One of these channels is remarkably linear with moderate incision.Toward the upper reaches of this channel it splits into a series of quasi-linear tributaries, one of which extents right up on top of the high surface.The northern large-scale channel has a more irregular morphology with many minor spurs and fins.Away from these large-scale channels, the spatial resolution is poorer but the drainage networks look more mature.There is an excellent view of both large-scale channels at 15 s.The flight path turns into the western catchment at 21 s.This catchment has a very different morphology as can be best seen at 24 s.Overall, the catchment is characterized by mature fluvial networks with concave upward channels.In the middle of this catchment, there are excellent examples of elongated fins that characterize the highly eroded interfluves.The back of the catchment is dramatically truncated by a steep convex escarpment which delineates the edge of the elevated surface.The catchment can be sub-divided into three sub-catchments, the western two of which are separated by a prominent ridge which is flown over at 33 s.All three sub-catchments have similar morphologies.

A5. Landscape_U_Fly.mp4
This landscape is considerably younger that the other four and it also has a very different fluvial morphology.The clearest overview occurs at 8 s where a single dramatically incised main channel is visible.This channel is characterized by a slot valley with several incised meanders.Both sides of this valley are oversteepened.A knickzone (possibly a sea cliff) is visible in the lowest reaches of the valley.Further back, the main channel splits into a series of three or four tributary channels, each of which is also deeply incised with steep valley walls.Isolated elongated spurs are visible in places.The northernmost tributary channel splits again into a host of minor channels that incise deeply into the edge of a dramatic undissected surface with an elevation of 200 m.At 23 s, there is a good view of a large-scale incised meander within the main channel.The edges of this meander are characterized by small-scale incision.Where the main channel splits, several prominent fins are visible.
Figures were generated using ArcPro, QGIS and GMT software tools.The Faroe-Shetland seismic MegaSurvey was generously provided by Petroleum Geo-Services (PGS) ASA to whom access enquires must be directed.Seismic interpretation and landscape extraction was carried out using the Petrel software package which is licensed by Schlumberger (SLB) Incorporated.Conversion from two-way travel time to depth and decompaction calculations were carried out using software written by Shaw-Champion et al. (2008).Borehole information is within the public domain and is freely available from the O&G authority website.Spectral analysis was carried out using software written by Roberts et al. (2019).Drainage inventories were extracted using the hydrology toolbox of ArcPro.Longitudinal river profiles were inverted using software written by Rudge et al. (2015).

Figure 2 .
Figure 2. Location of calibrated 3D seismic reflection surveys used to map buried landscapes.Map of region encompassing Faroe-Shetland trough.Open polygon = extent of Faroe-Shetland MegaSurvey generously provided by Petroleum Geoservices Ltd. (PGS); five labeled gray polygons = planforms of individual landscapes Q-U; labeled red lines = loci of 2D vertical slices from 3D survey (X-X′ shown in Figure 4; Y-Y′ shown in Figure 6); J and K are boreholes used for calibration; solid black circles = boreholes with checkshot surveys used to calibrate time-depth relationship for depth conversion; S = Shetland Islands.Inset at top left-hand side is schematic representation of 3D seismic volume showing vertical stacking of buried landscapes Q-U.Labeled red wiggly lines = individual mappable landscapes.

Figure 3 .
Figure 3. Paleocene to Early Eocene stratigraphic framework.At left-hand side, absolute time scale, nannoplankton (NP) zonation, dinoflagellate cyst (DC) zonation, and magnetic polarity chronology taken from Speijer et al. (2020).Buried landscapes P, Q, R, S, T, and U are shown in Figure 4. NS and F-S refer to North Sea and Faroe-Shetland stratigraphic templates taken from Mudge and Bujak (2001).Stratigraphic column for generalized lithologic log is adapted from Shaw-Champion et al. (2008).At right-hand side, timings of major bioevents are shown where LAD = last appearance datum and FAD = first appearance datum.Dinoflagellate cyst biostratigraphic framework is taken from Speijer et al. (2020).Note that Landscape U is younger than 50 Myrs as indicated by vertical arrow and that Landscape P, which crops out along coastline of West Greenland, is not described here.Redrawn from Conway-Jones and White (2022).

Figure 4 .
Figure 4. Vertical slice through seismic volume.(a) 2D vertical slice through 3D seismic reflection survey where Top Balder (54.8 Myrs) surface is flattened for clarity (see Figure 2 for location).Red/blue wiggles = reflections with positive/negative polarity; pairs of white arrow heads indicate loci of prominent channels; labeled drilling rig icons = projected loci of two boreholes used for stratigraphic correlation where J = 204/16 1 and K = 204/17 1.(b) Geologic interpretation of panel (a).Labeled colored polygons = Paleogene stratigraphic units; thick black lines = principal lithological boundaries and incised surfaces; thin black lines = deltaic foreset deposits that lie beneath upper pair of incised surfaces; encircled letters = four incised surfaces shown in Figure 5.

Figure 5 .
Figure 5. Four ancient buried landscapes.(a) Perspective view of buried Landscape Q. x-x′ shows orientation of vertical slice presented in Figure 4. (b)-(d) Same for landscapes R, S and T, respectively.

Figure 6 .
Figure 6.Visualization of Landscape U. (a) 2D vertical slice through 3D seismic reflection survey where topset surface of deltaic deposit is flattened for clarity (see Figure 2 for location).Red/blue wiggles = reflections with positive/negative polarity; pairs of white arrow heads indicate loci of prominent channels.(b) Geologic interpretation of panel (a).Blue/brown/ purple/blue polygons = mappable sub-units of Stronsay Formation; black lines = principal lithologic boundaries and incised surface; encircled letter = incised surface shown in panel (c).(c) Perspective view of Landscape U. Y-Y′ indicates orientation of vertical slice shown in panel (a).

Figure 7 .
Figure 7. Geomorphic analysis of Landscape R. (a) Topography of Landscape R. Black polygon = mappable planform of landscape; black rectangles = location of perspective slices shown in Figure 8.(b) Spatial gradient of Landscape R calculated from panel (a).(c) Drainage network extracted from Landscape R. Black lines = recovered fluvial channels; colored polygons = principal drainage catchments.(d) Calculated Gilbert Time, τ G .Small white patches denote lack of recoverability.Note that panels (c) and (d) do not exactly match since only clearest portions of drainage network are shown in panel (c).

Figure 8 .
Figure 8. 3D perspective slices across Landscape R. Four transects of central valley bisecting Landscape R. See Figure 7 for location.

Figure 9 .
Figure 9. Morphology of main channel extracted from Landscape R. (a) Detailed map showing main channel that bisects Landscape R. Black dotted lines with labels = loci of cross-sections shown in panels (c)-(h).(b) Spatial gradient of panel (a).(c)-(h) Six cross-sections of main channel whose loci are shown in panel (a).In each case, black line bounding brown polygon = elevation as function of distance across main channel; vertical black arrow = locus of main channel.Inset panel shows longitudinal river profile of main channel where: black line = elevation as function of distance from mouth; red/black dot = position of each cross-section at knickzone/ flat portion along river profile.

Figure 10 .
Figure 10.Main channel asymmetry across Landscape R. (a) Detailed map showing main channel that bisects Landscape R. Black dotted lines with labels = loci of cross-sections shown in panels (c)-(e).(b) Spatial gradient of panel (a).(c)-(e) Three cross-sections of main channel whose loci are shown in panel (a).In each case, black line bounding brown polygon = elevation as function of distance across main channel; black arrow points at steeper side of channel.Inset panel shows longitudinal river profile of main channel where: black line = elevation as function of distance from mouth; black dot = position of each cross-section along river profile.

Figure 11 .
Figure 11.Spectral analysis of river profiles extracted from Landscape R. (a) Longitudinal river profile of main channel extracted from Landscape R (see inset panel for location).Thick black line = observed river profile; white dots = inverse transform of wavelet power spectrum across all wavelengths; red line = inverse transform of wavelet power spectrum for wavelengths ≥23.9 km.(b) Morlet wavelet power spectrum for observed river profile shown in panel (a).Black line = wavelength of 23.9 km.(c) Comparison of Morlet wavelet and Fourier transform power spectra of observed river profile.Black line = distance-averaged Morlet wavelet power spectrum; gray line = five-point moving average for Fourier transform spectrum.(d) Morlet wavelet power spectrum that has been rectified and normalized with respect to power spectral slope of 2 (i.e., red noise).Black/gray lines = mean/extremal distance-averaged Morlet wavelet power spectra for all river channels shown in inset panel of panel (a); dotted pink lines = power spectral slope of 1 (i.e., pink noise).

Figure 12 .
Figure 12.Geomorphic analysis of Landscape S. (a) Topography of Landscape S. Black polygon = mappable planform of landscape.(b) Spatial gradient of Landscape S calculated from panel (a).(c) Drainage network extracted from Landscape S. Black lines = recovered fluvial channels; colored polygons = principal drainage catchments.(d) Calculated Gilbert Time, τ G .Small white patches denote lack of recoverability.

Figure 13 .
Figure 13.Geomorphic analysis of Landscape T. (a) Topography of Landscape T. Black polygon = mappable planform of landscape.(b) Spatial gradient of Landscape T calculated from panel (a).(c) Drainage network extracted from Landscape T. Black lines = recovered fluvial channels; colored polygons = principal drainage catchments.(d) Calculated Gilbert Time, τ G .Small white patches denote lack of recoverability.

Figure 14 .
Figure 14.Geomorphic analysis of Landscape U. (a) Topography of Landscape U. Black polygon = mappable planform of landscape.(b) Spatial gradient of Landscape U calculated from panel (a).(c) Drainage network extracted from Landscape U. Black lines = recovered fluvial channels; colored polygons = principal drainage catchments.(d) Calculated Gilbert Time, τ G .Small white patches denote lack of recoverability.

Figure 15 .
Figure 15.Longitudinal river profiles extracted from four mapped landscapes.(a) Drainage planform for Landscape U. Black polygon = extent of recoverable landscape; black lines = extracted river channels; red line = main river channel.(b) Principal river profiles extracted from this landscape.Black lines = major river profiles; red line = river profile shown in panels (c) and (d.) (c) Vertical slice through 3D seismic reflection survey along principal river profile.(d) Geologic interpretation of panel (c).Pale blue polygon = undivided Stronsay Formation; black lines = principal reflective surfaces interpreted to be intra-formational lithologic boundaries; horizontal arrows indicate stratal onlap.(e)-(p) Same for landscapes T through R. For panel (p), colored polygons = stratigraphic formations shown in Figure 3b.

Figure 16 .
Figure 16.Slope-area analysis of main channel extracted from Landscape R. (a) Extracted river profile.Black line = elevation plotted as function of distance; gray portions omitted from slope-area analysis; red line = drainage area as function of distance along river profile.(b) Measurements of slope as function of drainage area for river profile shown in panel (a).Circles with horizontal and vertical bars = measurements sampled every 200 m along river profile, assuming uncertainties of 20%; black line = best-fitting (r 2 = 0.068) linear relationship using parameter values given at bottom lefthand side.

Figure 17 .
Figure 17.Calibration of erosional parameters.(a) Normalized rms misfit as function of erosional parameter, m, for landscapes R-U.Solid/dashed/dotted/dot-dashed lines = misfit function for landscapes R/S/T/U; vertical gray bar = range of m values for which common minimum exists.(b) rms misfit and normalized spatial roughness as function of λ s for Landscape R. See right-hand edge for color scale.(c) rms misfit and normalized temporal roughness as function of λ t for Landscape R.

Figure 18 .
Figure 18.River profile inverse modeling.(a) Observed and calculated river profiles for Landscape R. Gray lines = observed profiles; red dots = calculated profiles.(b)-(d) Same for landscapes S-U.(e) Cumulative uplift history for Landscape R. Black line = calculated cumulative uplift; gray band = 2σ uncertainty.(f)-(h) Same for landscapes S-U.(i) Uplift rate history for Landscape R. Black line = calculated uplift rate; gray band = 2σ uncertainty; vertical arrows indicate discrete pulses of uplift.(j)-(l) Same for landscapes S-U.Note that model time starts at zero and increases from left to right (see Figure 3 and Conway-Jones and White (2022) for converted geologic time scale).
Wapenhans et al. (2021)Roberts et al., 2019;Wapenhans et al., 2021)ly homogeneous, we infer that at wavelengths shorter than 100 km the channelization process is affected by non-linear instabilities, shock behavior and hydraulic jumps.Significantly, landscape simulations generated by tectonic forcing yield power spectra that are completely dominated by red noise and do not show a cross-over transition to pink noise(O'Malley et al., 2021;Roberts et al., 2019;Wapenhans et al., 2021).Absence of a spectral transition reflects the fact that landscape simulations do not incorporate complex (i.e., realistic) fluid dynamics.A detailed discussion of the power spectral content of both observed and simulated fluvial drainage networks is provided byO'Malley et al. (2021)and byWapenhans et al. (2021).
Roberts et al. (2012)2012)developed and applied inverse modeling algorithms which demonstrated that it is possible to recover uplift rate histories that CONWAY-JONES AND WHITE