Histories of vertical lithospheric motions provide important clues about geodynamic processes. We present evidence of an ancient (c. 58–55 Ma) landscape that probably underwent rapid uplift and subsidence during the initiation of the Icelandic plume. Now buried beneath c. 0.4–0.8 km of rock in the North Bressay region in the North Sea, this landscape is located within a sedimentary basin on the margin of the North Atlantic Ocean. We use high-resolution 3D seismic reflection data to map this ancient surface. Correlation of stratigraphy with a survey in the Bressay region constrains the age and depositional environment. The landscape contains excellent evidence of meandering fluvial channels, some of which record avulsions, that terminate against a coastline to the east where deltaic landforms are identified. The landscape was depth-converted and decompacted to generate a digital elevation model from which river profiles were extracted. Their geometries indicate that the landscape was generated by three phases of uplift. This history of uplift and subsidence is analogous to similar-aged landscapes in the Judd area c. 400 km to the west and Bressay c. 30 km to the south, and appears to be another manifestation of lithospheric motions generated by the passage of warm thermal anomalies away from the Icelandic plume.

Histories of vertical lithospheric motions are important constraints for identifying processes responsible for driving Earth's surface evolution over geological timescales (e.g. shortening, extension, mantle convection). At present, the North Atlantic Ocean hosts a region of mantle convective upwelling centred on Iceland, resulting in 1–2 km of dynamic support and elevated continental margins at its fringes (e.g. White and McKenzie 1989; Rickers et al. 2013; Hoggard et al. 2016). The composition and age of basaltic magmatism, V-shaped ridges south of Iceland, and topographic, stratigraphic, tomographic and thermochronological observations have provided insights into the plume's evolution since its formation earlier than 60 Ma (e.g. Lawver and Müller 1994; Saunders et al. 1997; Japsen and Chalmers 2000; Jones et al. 2002; Mudge and Jones 2004; Storey et al. 2007; Shaw Champion et al. 2008; Parnell-Turner et al. 2014; Schoonman et al. 2017; Hardman et al. 2018). Nonetheless, we have a limited understanding of its planform through time or the influence of temporal fluctuations. One way to probe the evolution of the Icelandic plume at greater spatial and temporal resolutions is through the identification and study of buried terrestrial landscapes that were once uplifted by thermal anomalies sourced from the ancient plume (e.g. Underhill 2001; Shaw Champion et al. 2008; Hartley et al. 2011; Stucky de Quay et al. 2017; Conway-Jones and White 2022).

The Faroe–Shetland, Moray Firth and North Sea basins display stratigraphic records containing coarse Paleocene–Eocene terrestrial sandstones sandwiched between marine mudstones. Such observations indicate that these regions experienced short-lived (<3 myr) histories of marine–terrestrial–marine conditions close to the Paleocene–Eocene boundary (e.g. Mudge and Bujak 2001; Underhill 2001; Mudge and Jones 2004; Mackay et al. 2005; Shaw Champion et al. 2008; Hardman et al. 2018; Jolley et al. 2021; Conway-Jones and White 2022). There is continuing debate about the contributions of plate-tectonic processes (e.g. shortening and extension) and sub-plate support in generating vertical motions in these areas (see, e.g. Jolley et al. 2021). However, the amplitudes of uplift (up to c. 1 km), combined with its rapidity and, importantly, subsequent subsidence, indicate that these vertical motions cannot be attributed to glacio-eustasy, magmatic underplating or crustal tectonics (e.g. shortening or extension). Instead, they have been attributed to the passage of thermal, buoyant, anomalies in a low-viscosity channel beneath the lithosphere away from the centre of the Icelandic plume (e.g. Rudge et al. 2008; Shaw Champion et al. 2008).

In the Judd region, located in the Faroe–Shetland Basin, an erosional unconformity has been mapped at the base of the coarse sandstone Flett Formation (Smallwood and Gill 2002; Shaw Champion et al. 2008). Depth conversion and decompaction of infilling stratigraphy and recorded lignite suggest that the unconformity was created by erosion of a terrestrial landscape with at least 900 m of relief (Shaw Champion et al. 2008; Hartley et al. 2011). Dinoflagellate and pollen chronostratigraphy indicate that this landscape formed by c. 55 Ma. Jolley et al. (2021) suggested that the unconformity might be a composite of two erosional surfaces that formed as a result of two phases of uplift between c. 58 and 56 Ma. Longitudinal profiles extracted from this surface have almost 1 km of relief and contain three large (>100 m of relief) knickzones (Hartley et al. 2011). Two younger and two older (c. 61.5–52.5 Ma) landscapes that each grew and were buried in a few million years have been mapped in the Judd area (Conway-Jones and White 2022). The history of vertical motions inferred from these observations and from similar-aged stratigraphy in surrounding regions has been interpreted as a record of thermal perturbations of the Icelandic plume (see, e.g. Conway-Jones and White 2022, their Section 2.5, for a recent summary).

To the east, in the Bressay region of the North Sea, a similar-aged (c. 58–55 Ma) terrestrial landscape, now buried by c. 1.5 km of rock, has been identified in well and seismic data (Fig. 1; Underhill 2001). The Bressay region is located on the East Shetland Platform of the northern North Sea. Buried extensional faults and backstripped well data indicate that the East Shetland Platform was stretched by a factor β = 1.3 ± 0.04 from the Late Jurassic for c. 60 myr (White 1990; Underhill 2001). It currently sits on the eastern fringe of the Icelandic plume (Rickers et al. 2013; Hoggard et al. 2016). Cuttings from wells that intersect the c. 58–55 Ma landscape contain angiosperm (flowering plant) debris, which, combined with mapped dendritic drainage patterns and coarse clastic material, indicates that it formed subaerially (Stucky de Quay et al. 2017). Similar to the Judd area, longitudinal river profiles extracted from this landscape contain three broad knickzones, which do not correlate with substrate, and instead indicate a staged uplift history. Dinoflagellate cysts in the marine stratigraphy below and above the landscape indicate that it formed in less than 3 myr and was rapidly drowned, again similar to the Judd area. These buried landscapes, scattered across fringes of the North Atlantic Ocean, thus provide an opportunity to better constrain Icelandic plume activity by providing multiple points of reference in space and time.

This study describes a hitherto unmapped Early Paleogene eroded surface c. 30 km north of the Bressay buried landscape (herein termed the North Bressay region; see Fig. 1). Similar to other sedimentary basins on the fringes of the North Atlantic Ocean, it contains stratigraphic evidence of rapid uplift, landscape evolution and subsidence during incipience of the Icelandic plume (e.g. Underhill 2001; Rudge et al. 2008; Shaw Champion et al. 2008; Hartley et al. 2011; Stucky de Quay et al. 2017; Jolley et al. 2021). In the following sections we first describe seismic data and examine stratigraphic, tectonic and geomorphological observations from the North Bressay area. We then extract, depth-convert and decompact a prominent erosional surface, to which we tentatively assign an age based on correlation with stratigraphy in the Bressay area to the south. Drainage networks are extracted from this landscape and longitudinal river profiles are inverted for a history of uplift. Finally, we discuss implications of the calculated uplift histories for evolution of the region during the initial stages of the Icelandic plume.

A high-resolution, 3D seismic survey (MC3D-ESP2015M) located directly north of the Bressay region was made available for the study by PGS (Petroleum Geo-Services) and Equinor. Dendritic and meandering drainage, unconformities and shallower (younger) polygonal faulting in North Bressay stratigraphy are clearly visible in slices through seismic variance (trace-to-trace variability; Fig. 2a and b). The meandering drainage networks at a depth of c. 600 ms (two-way travel time) coincide with the strongest reflector visible in the seismic cross-section shown in Figure 3a, where evidence for valley incision can also be observed (see red arrows in Fig. 3b).

To characterize the stratigraphy, well logs within the MC3D-ESP2015M seismic survey were also made available. However, their location to the east of a large unconformity (Figs 2c and 3) means that direct stratigraphic correlation between the wells and the eroded surface to the west is challenging. Fortunately, this survey abuts the northern edge of the PGS BBK survey, which was used to map a c. 58–55 Ma buried landscape in the Bressay region (Stucky de Quay et al. 2017). This previous study incorporated seven wells that intersect the eroded landscape and stratigraphy located towards the south. The well reports contained geophysical data (e.g. check-shots, gamma-ray logs) and biostratigraphic observations that were used to depth-convert seismic data, calibrate compaction parameters and quantify the age of stratigraphy and its uncertainties. The reader is referred to Stucky de Quay et al. (2017, and references therein) for detailed lithological and biostratigraphic descriptions and detailed age constraints. Similarities in acoustic impedance contrasts (e.g. seismic facies) between the MC3D-ESP2015M and BBK surveys give a basis for interpreting the stratigraphy in western North Bressay and for depth-converting and decompacting its stratigraphy using similar techniques.

Correlation of seismic stratigraphy indicates that the prominent incised layer at c. 600 ms in North Bressay is of Early Cretaceous age (Cromer–Knoll–Shetland Group; Fig. 3). Overlying deposits are likely to be the Paleocene–Eocene aged Moray and Montrose groups (compare Bressay and Judd chronostratigraphy as given by Stucky de Quay et al. 2017, fig. 2). Well logs in the Bressay area suggest that the incised Cretaceous layer is composed of limestones and is overlain by interbedded sand and shales containing coarse, pebbly sandstone (see Stucky de Quay et al. 2017, fig. 4). Acoustic impedance contrasts within the Montrose and Moray groups are suggestive of interbedded sands, shales and conglomerates (Fig. 3). Above the Montrose and Moray groups, a succession of undifferentiated marine mudstones is mapped up to the seabed (Fig. 3b). Below Cretaceous stratigraphy in the Bressay area, the logs report sandstones, which are most probably of the Old Red Sandstone (Devonian) Group. Figure 4 shows a tentative interpretation of Paleocene–Eocene chronostratigraphy for North Bressay compared with the Bressay and Judd regions (after Shaw Champion et al. 2008; Stucky de Quay et al. 2017; Jolley et al. 2021). The fidelity of interpreted North Bressay stratigraphy could obviously be assessed and improved if data from wells penetrating the erosional surface become available.

Seismic mapping

The strong reflector at c. 600 ms in Figure 3a was mapped in two dimensions with a horizontal resolution of c. 10 m. Results using auto-tracking and manual mapping of in-lines and cross-lines yield similar results owing to the continuity of the mapped surface. The resulting mapped surface in two-way travel time is shown in Figure 5. It contains a dendritic drainage system in the NW comprising multiple tributaries. Unlike Bressay, these channels show evidence of meandering of the order of <1 km. Overlapping and interconnected meanders record channel avulsion events (e.g. southernmost channel in Fig. 5). The channels drain from higher elevations (c. 600 ms below sea-level; c. 400 ms below seabed) toward the east (>800 ms below sea-level). Incised valleys appear to stop along a palaeo-coastline, and localized deposits at the mouths of channels are most probably deltaic in origin (blue arrow in Fig. 5). These geomorphological features, combined with stratigraphic evidence for coarse fluvial sediments in overlying sediments, are suggestive of formation in a subaerial environment akin to the erosional landscape mapped at a similar stratigraphic level in the Bressay region. Preservation of meandering channels, interfluves, deltaic deposits and overlying sedimentary deposits suggests rapid burial and modest erosion during marine transgression. We note that quasi-linear depressions are observed between the mapped coastline and unconformity (Fig. 5, dashed and continuous lines). They do not obviously link with the meandering channels mapped upstream. We note that they sit between depositional highs, which we interpret as deltaic deposits. We are therefore cautious about interpreting them as downstream fluvial reaches and instead tentatively suggest that they might be marine channels.

North Bressay v. Bressay buried landscapes

The Bressay and North Bressay drainage systems are located c. 30 km apart. They share commonalities in their morphology including incised valleys sloping eastward toward palaeo-coastlines, and appear to contain similar stratigraphy, including coarse, pebbly channel infill. Stratigraphic correlations indicate that they are of a similar late Paleocene to early Eocene age (Fig. 4). The North Bressay landscape has a larger extent, greater resolution and more clearly defined valleys than the Bressay landscape. It contains a greater number of individual basins, channels and tributaries as well as evidence of coeval deposits at the valley mouths. Owing to its biostratigraphic records and abundance of available intersecting cores, data from the the Bressay area were used to depth-convert, decompact and constrain the evolution of the North Bressay landscape in the following section (Stucky de Quay et al. 2017).

Notwithstanding these similarities, we note that are there subtle but important differences between the two landscapes and datasets. First, the Bressay landscape sits below prominent gas columns that partially obscure deeper stratigraphy; they are far less prominent or are absent atop the North Bressay landscape. Second, the MC3D-ESP2015M survey does not extend far enough to the west for river heads to be mapped. Despite these challenges, the proximity of the nearby Bressay region provides an opportunity to constrain the relief of the North Bressay landscape as it formed and make a tentative assessment of its uplift and erosional history.

Depth conversion and decompaction

Five wells in the Bressay area had check-shot data available, which range from zero to c. 1200 ms (c. 1400 m) below the seabed. These data were used to constrain the following optimal relationship between two-way travel time (t, milliseconds) and depth (z, metres):
Equation (1) fits the available check-shot data with r2 = 0.99 (see Stucky de Quay et al. 2017, fig. 9a). One standard deviation was used to assess uncertainties in calculated depths.
The depth-converted landscape was corrected for post-depositional compaction using a simple widely used parameterization of porosity (ϕ),
where ϕo is initial porosity (e.g. at the seabed), λ is compaction wavelength and z is depth. Values of compaction parameters (ϕo, λ) were estimated by fitting check-shot data using an empirical relationship that includes depth, velocity and compaction parameters, which yielded optimal (best-fitting) parameters of ϕo = 0.77 and λ = 2.4 km (Wyllie et al. 1956; Sclater and Christie 1980; Christensen 1982; Czarnota et al. 2013). It should be noted that, although a global minimum is identified, there is a trade-off between ϕo and λ, which affects assessment of uncertainties in compaction parameters (see Stucky de Quay et al. 2017, fig. 9b and c). Decompaction was performed assuming that the solid grain fraction is conserved; that is, 1ϕ(z)dz is the same between compacted and decompacted columns. Depositional thickness, z3, is given by solving the integrals to yield
where z2 – z1 is the compacted thickness of the layer and z3 is its depositional thickness. It should be noted that z3 appears on both sides of the equation, which was solved by Newton–Raphson iteration; we set z3 = z2 – z1 in the first iteration, and values of z3 tend to converge to 3 decimal places. by 3–5 iterations (see, e.g. Stucky de Quay et al. 2017, for details). The depth-converted and decompacted eroded surface is shown in Figure 6a. Uncertainties in relief (c. ±20 m) vary as a function of the velocity model and compaction parameterization. The total relief of the landscape is c. 500 m, and contains evidence of rivers draining east towards a palaeo-coastline.

Drainage patterns

To recreate the landscape's hydrology at the time of subaerial exposure, Esri flow routing algorithms were used to extract drainage networks from the mapped landscape (Fig. 6a). First, flow directions were calculated using the steepest descent from each cell in the digital elevation model of the landscape (Tarboton 1997). Subsequently, flow accumulation to each cell was measured and flow lengths were calculated. Drainage networks were extracted from the six main basins and a subset of 107 longitudinal river profiles (elevation as a function of distance) was generated (Fig. 6b–g). The extracted drainage patterns highlight the presence of eastward-flowing, dendritic and meandering networks atop the eroded landscape. The mapped rivers are up to c. 13 km long and have a maximum relief of c. 380 m. Stucky de Quay et al. (2017) showed that uncertainties in the velocity model and compaction parameters stretched or compressed calculated river profiles (vertically) by ±25 m and that the location and relief of knickzones are largely unaffected by these uncertainties. We therefore extract longitudinal river profiles from the digital elevation model produced using the optimal velocity model and compaction parameters.

The river profiles shown in Figure 6b–g are all broadly convex-upward and most contain large knickzones with relief up to the order of 100 m. The geometries of these profiles indicate that the landscape was eroding in response to changes in its history of uplift before it was rapidly drowned and buried by sediment (see, e.g. Howard et al. 1994; Rosenbloom and Anderson 1994; Whipple and Tucker 1999; Pritchard et al. 2009). The mapped surface is essentially a snapshot of a short-lived (less than a few million years) evolving landscape. In the following section we invert the four longest complete longitudinal river profiles for a history of uplift of the North Bressay landscape.

Pritchard et al. (2009) and Roberts and White (2010) showed that longitudinal river profiles could be successfully inverted for a history of uplift rate. Following Stucky de Quay et al. (2017), we use a simplified version of the stream power model to invert rivers draining the North Bressay landscape for uplift rate, U, as a function of time, t,
where erosion rate, ∂z/∂t, is parameterized using the advective term on the right-hand side of equation (4). This term controls the velocity of kinematic erosional waves as they propagate upstream (e.g. Rosenbloom and Anderson 1994; Whipple and Tucker 1999). In this parameterization, erosional velocities, and ultimately the duration of landscape response to uplift, depend on upstream drainage area, A, which is measured from the extracted landscape, and the values of the erosional parameters, v and m. z is elevation and x is distance along the river.

The lack of direct biostratigrapic control means that the age of the channelized North Bressay surface cannot be constrained with the precision of the Bressay and Judd areas. However, the overlying sedimentary rock appears to be of a similar age to sediments that infilled the Bressay landscape (Late Paleocene–Early Eocene). Because there are no other heavily channelized surfaces in either region, and considering their proximity (c. 30 km), presence of infilling coarse sandstones and comparable palaeo-coastlines, we consider it likely that they they formed in a similar way (subaerially) around the same time. The biostratigraphic records discussed above indicate that the duration of exposure is <3 myr. Joint inverse modelling of the Bressay drainage inventory indicates that the best-fitting value of erosional parameter m = 0.5. Calibrated values of v in the Judd and Bressay regions are 4.75 ± 0.1 Ma−1 and 1.6 ± 0.1 Ma−1, respectively (Hartley et al. 2011; Stucky de Quay et al. 2017). We invert North Bressay river profiles for histories of uplift using this range of erosional parameter values.

The inverse problem is solved by seeking the minimum value of the objective function, H, which incorporates data (root mean squared, r.m.s.) misfit, model roughness and a positivity penalty function, P,
where N is the total number of measurements of elevation and distance in the drainage inventory. Data variance σi is set to 20 m, which is a conservative estimate of the vertical resolution of the seismic data. In practice, the model seeks the smoothest uplift rate history that yields the lowest residual misfit between observed and theoretical river profiles (Parker 1994; Roberts and White 2010). U and U are the first and second derivatives of calculated uplift, which is parameterized by selecting M discrete values (see Roberts and White 2010). Consistent with Stucky de Quay et al. (2017), we set smoothing parameters Wn = 0.1. The positivity penalty function P = cosh(U) – 1 if U < 0 and P = 0 if U ≥ 0. The model uses Powell's algorithm, a conjugate gradient method, to invert for uplift rate histories, which yields stable solutions (Roberts and White 2010; Hartley et al. 2011; Stucky de Quay et al. 2017). In the starting model z(x) = 0 and U(t) = 0. Two rivers from the two largest basins were jointly inverted for their uplift rate histories (Fig. 6a).

Residual misfit between observed (black curves) and best-fitting theoretical (red curves) river profiles is low (r.m.s. ∼ 1; Fig. 7a and c). Changing the value of erosional parameter v scales the timing of calculated uplift. As expected, the relatively low value of v = 1.6 Ma–1 (calibrated in the Bressay region) results in a calculated uplift history that is about three times longer-lived than when the faster velocity, based on calibration from the Judd area, v = 4.75 Ma–1, is used (compare the black and grey curves in Fig. 7b and d). Given the proximity of the study area to Bressay, we favour the slower advective velocity v, which we note yields broadly similar uplift histories to those calculated using rivers in the Bressay area (compare Stucky de Quay et al. 2017, fig. 12b). However, we acknowledge that both parameterizations yield uplift histories <3 myr in duration and that no direct biostratigraphic control means that evidence to discriminate between the two uplift histories is lacking.

Seismic reflection data in the North Bressay region of the North Sea contain evidence of an ancient fluvial landscape now buried beneath c. 0.4–0.8 km of marine deposits. The quality of the survey permits extraction of a high-resolution (of the order of 10 m) digital elevation model of the landscape. Correlation of stratigraphy observed within this survey and an adjacent survey in the Bressay region indicates that the landscape formed close to the Paleocene–Eocene boundary (c. 55 Ma) within a few million years (<3 myr). Mapped meandering rivers show evidence of avulsion in the southern sector of the survey. The rivers appear to terminate in the east against what we have mapped as an undulating coastline. At least two deltaic landforms abut the mouths of rivers mapped in the northern section of the survey. The fluvial landscape extended beyond the western limit of the survey (see Fig. 5).

Inverse modelling of the two largest river systems extracted from the landscape indicates that fluvial geometries are best explained by uplift occurring in three distinct stages (I, II and III); this is strikingly analogous to uplift histories calculated by inverting river profiles from similar-aged (now buried) landscapes c. 30 km to the south (Bressay) and c. 400 km to the west (Judd; Fig. 8). Rapid uplift and subsidence of these landscapes indicates that neither glacio-eustasy (which is too low-amplitude), underplating (which generates uplift but no or very little subsidence) nor crustal tectonics (e.g. shortening followed by extension) played a governing role in generating the landscape. Instead, simple isostatic calculations indicate that these landscapes could have been generated by the passage of a warm thermal anomaly (c. 130–50°C excess temperature) away from the Icelandic plume in a low-viscosity channel beneath the plate (Rudge et al. 2008).

The proximity of North Bressay to Bressay, and their stratigraphic similarities, indicates that uplift may also have been generated by the passage of a warm thermal anomaly beneath the plate in c. 1–3 myr. If we assume that isostasy prevails and that compressibility is negligible, the average excess temperature in a low-viscosity channel, T¯, required to generate the observed uplift, u, is given by
where h is the thickness of the anomalously warm layer, To is background temperature and α = 3.28 × 10–5 K–1 is thermal expansivity. If the background temperature of the plume head is 1400°C and the thickness of the hot layer is 200 km, an excess temperature of 44°C is required to generate the c. 300 m of uplift predicted by inverting river profiles. It should be noted that this calculation is insensitive to assumed background temperature. The calculated excess temperature is within the range of values estimated for similar (c. 61.5–52.5 Ma), now buried, landscapes in the Judd area (16–120°C; Conway-Jones and White 2022).

We tentatively suggest that the buried landscape mapped in this study is another record of rapid transient vertical motions generated by mantle convection (Fig. 9). Mapping of ancient buried landscapes is a rich source of geomorphological information including constraints on erosion rates and planform responses to rapid uplift and subsidence. The discovery and analysis of spatially distributed buried landscapes through seismic mapping holds great potential for providing quantitative constraints on the timing and extent of the Icelandic plume through time.

In this study we present a newly identified and previously unmapped ancient buried landscape located in the North Bressay region of the northern North Sea. The landscape was identified by mapping seismic reflectors on a 3D survey provided by PGS and Equinor. It contains excellent evidence of well-preserved, incised meandering channels, avulsion and a palaeo-coastline with deposition of sediments in deltas. Correlation of stratigraphy in the region with an existing survey to the south indicates that the eroded surface was terrestrial and formed in a few million years between c. 58 and 55 Ma. Longitudinal river profiles extracted from the depth-converted and decompacted landscape contain convex-upward reaches, which indicate that the landscape was youthful when it was subsequently buried by marine sediment during subsidence. Inverse modelling of the longitudinal profiles indicates that uplift of the landscape occurred in three stages, similar to uplift histories calculated for similar-aged landscapes in the Judd area c. 400 km to the west and Bressay c. 30 km to the south. The uplift and subsidence history of the North Bressay landscape is consistent with the notion that excess thermal anomalies travelling beneath the plate away from the Icelandic plume generated rapid Paleocene–Eocene lithospheric uplift and subsidence.

Seismic and well data were made available by PGS (survey MC3D-ESP2015M) and Equinor, and we thank them for their help. We thank R. Bell, C. Jackson, G. Hampson, L. Lonergan, T. Weight and N. White for their help. N. Schofield, T. Sømme and two anonymous reviewers are thanked for their insightful reviews.

GSQ: conceptualization (equal), data curation (equal), formal analysis (equal), investigation (equal), methodology (equal), software (equal), writing – review & editing (equal); GGR: conceptualization (equal), formal analysis (equal), funding acquisition (lead), investigation (equal), methodology (equal), project administration (equal), resources (lead), software (equal), supervision (lead), validation (lead), writing – original draft (lead), writing – review & editing (equal)

G.S.Q. was supported by the National Environmental Research Council (NERC) Grantham Institute SSCP DTP (grant number NE/L002515/1).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Requests for access to seismic data should be made to PGS (Petroleum Geo-Services). Data, input files and source code for the inverse models are available from the authors.

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)