Isostasy as a tool to validate interpretations of regional geophysical datasets – application to the mid-Norwegian continental margin

The concept of isostasy has been applied to several regional studies worldwide, helping to estimate lithospheric strength and explain the shape and origin of topography (e.g. Watts 2001; Kaban et al. 2003). Local studies, however, which investigate the crustal structure in more detail by employing refraction and deep reflection seismic data as well as gravity modelling, often ignore the additional insight that could be gained from isostatic considerations. Positive exceptions to this are, for example, Watts et al. (1985), Zeyen & Fernàndez (1994), Maţenco et al. (1997), Hirsch et al. (2009) and Fullea et al. (2010), where isostasy has entered the model building and interpretation. Methods for joint analysis and inversion of gravity and isostasy have, for example, been developed by Zeyen & Fernàndez (1994), Stewart et al. (2000), Braitenberg et al. (2004) and Salem et al. (2014).

Only a few of the studies involving isostasy cover the mid-Norwegian margin (Kimbell et al. 2004; Ebbing 2007) and we show in this paper that simple isostatic considerations can also improve the understanding of the crustal structure.

Continental margins may be especially well suited for isostatic studies since the crust is thinned and its deeper structure and Moho are more easily imaged by seismic methods than, for example, those of cratonic lithosphere. Seismic mapping is especially extensive (and accessible) along hydrocarbon-bearing continental margins such as the mid-Norwegian one. Furthermore, thinned crust and lithosphere are generally weaker, implying that lateral density variations and Moho undulations will have a larger expression in the bathymetry or basin shape. Thermal weakening during rifting and breakup as well as rigidification during subsequent cooling additionally affect the flexure of the margin.

The enormous amount of geophysical and geological data that exists for the region of the NE Atlantic has been gathered, analysed and compiled in recent years into a comprehensive tectonostratigraphic atlas (the NAG-TEC Atlas, Hopper et al. 2014). Based on these data, a new seismic-derived Moho grid has been compiled (Funck et al. 2016) and subsequently further improved using gravity inversion (Haase et al. 2016). The latter study also constructed a first-order density model over the entire NE-Atlantic, which we utilize here to forward calculate an isostatically consistent bathymetry that reflects all internal loading of the model. The resulting differences from the true bathymetry enable us to detect consistencies and discrepancies between the datasets mapping the deep subsurface structure. Since gravity and isostatically balanced elevation are independent observables but both depend on the underlying density model, consistencies and discrepancies can point towards poorly understood or erroneously interpreted areas. The obtained elevation residuals are then discussed with respect to existing seismic profiles and structural interpretations along the margin. By considering residuals of shorter wavelengths only, we focus on the crustal structure, which is much better mapped and discussed than the underlying lithospheric mantle. Previous studies that investigated regional isostatic and gravity modelling in the NE Atlantic region (Kimbell et al. 2004; Reynisson et al. 2010) had less complete datasets at hand and focussed on somewhat different questions. This study merely tests an existing subsurface model and does not attempt to improve it.

The mid-Norwegian margin, like any passive rifted margin, consists of numerous rifted blocks and interjacent sedimentary basins and sub-basins (Fig. 1). The margin is subdivided into the three major segments of the Møre, Vøring and Lofoten margins. The large shelf areas of the North Sea and Barents Sea extend to the south and north, respectively. The structure of the mid-Norwegian margin reflects a long period of extension following the orogenic collapse of the Caledonian orogeny. Earliest extension during the post-orogenic collapse initiated in Devonian–Carboniferous times, forming large detachment faults and basins that are now partly exposed in onshore areas (e.g. Seranne & Seguret 1987; Osmundsen & Andersen 2001, and references therein) and possibly buried by younger sediments in the offshore domain. Carboniferous to early Mesozoic extension created the early, proximal rift basins (e.g. the Trøndelag Platform; Brekke 2000). Distal rift basins formed during Middle Jurassic to Early/middle-Cretaceous (e.g. the Lofoten, Vøring and Møre basins; Brekke 2000). Late Cretaceous to Paleocene extension shaped the outer ridges and sub-basin system defining the pre-breakup rift complex (Gernigon et al. 2003; Ren et al. 2003). Volcanic activity affected the distal margins in the Late Paleocene but was most significant later, during the breakup and volcanic margin formation in the early Eocene (c. 54 Ma, Talwani & Eldholm 1977; Skogseid & Eldholm 1987). This volcanic activity was part of the North Atlantic Igneous Province (Eldholm & Grue 1994; Meyer et al. 2007) and is reflected along the present-day mid-Norwegian outer margin by extensive regions with seaward-dipping reflectors and massive volcanic flows. Magmatic activity also led to intrusions into the thinned continental crust and underplating in the continent–ocean transition zone. Oceanic spreading along the Aegir and Mohn’s ridges led to the complete separation of the Greenland and Eurasian plates. Continued, contemporaneous Cenozoic rifting in East Greenland and associated breakup in late Oligocene (c. 25–24 Ma) led to the formation of a second spreading axis and to the separation of the Jan Mayen microcontinent from Greenland (Talwani & Eldholm 1977; Vogt et al. 1980; Nunns 1982; Gudlaugsson et al. 1988; Gaina et al. 2009; Gernigon et al. 2015). Today, the Kolbeinsey Ridge is the active spreading centre linking the Reykjanes Ridge in the south to the Mohn’s Ridge in the north.

A large number of seismic refraction profiles cover the mid-Norwegian margin but become more sparse to the north (Fig. 1). Some profiles are several decades old and large receiver spacing (10–20 km) often prevents the detailed mapping of smaller-scale features. Nevertheless, this extensive grid predominantly shaped our present-day understanding of the deep margin structures.

The 300 km wide Møre Margin is characterized by a narrow platform and a rapid thinning of the continental crust. The main area of crustal thinning is confined to a narrow zone with a width of only a few tens of kilometres (Olafsson et al. 1992; Raum 2000; Mjelde et al. 2008; Osmundsen & Ebbing 2008; Kvarven et al. 2014; Nirrengarten et al. 2014). West of this major necking zone lies the nearly 200 km wide and long Cretaceous Møre Basin, which reaches depths of more than 15 km. This basin is underlain by highly thinned crust (in places less than 10 km thick) and high-velocity, high-density lower crustal bodies (LCBs) have been mapped both outbound and inbound of the necking zone but not directly underneath it (Olafsson et al. 1992; Kvarven et al. 2014; Nirrengarten et al. 2014). The outer limit of the margin comprises the Møre Marginal High which is overlain by a series of volcanic flows. Whereas some authors interpret the basement of the Møre Marginal High as thick continental crust (e.g. Reynisson et al. 2010), others suggest a sediment succession of several kilometres thickness beneath the flood basalts (Nirrengarten et al. 2014).

The Jan Mayen corridor, which continues as the Jan Mayen Fracture Zone in the oceanic domain, offsets the much wider Vøring Margin (450 km) from the Møre Margin. The Vøring segment comprises the wide, proximal Trøndelag Platform, a domain of terraces and fault complexes, the 200 km wide Vøring Basin and the outer Vøring Marginal High (Brekke 2000). The crust in the Vøring segment thins gradually over a distance of more than 200 km from over 30 km under the platform domain in the east to a minimum of 8 km under the distal sag basin. Extensive LCBs underlie most of the outer margin and often constitute nearly 50% of the total crustal thickness (Mjelde et al. 2009a). P-wave velocities are generally between 7.2 and 7.5 km s⁻¹ within these bodies, but supposedly increase to significantly more than 8.0 km s⁻¹ in an area at the SE edge of the margin (Raum et al. 2006), although such values are not confirmed by other studies (Rouzo et al. 2006). Densities of these bodies are estimated to range between 3000 and 3300 kg m⁻³ (Mjelde et al. 2009a; Kvarven et al. 2014; Nirrengarten et al. 2014). Extensive lava flows and sills characterize the outer Vøring Basin and adjacent Vøring Marginal High, making it diffcult to accurately determine depth to both basement and Moho. Commonly used interpretations of the COB (continent–ocean boundary, here used as the outer limit of the continent–ocean transition zone) place it along the outer Vøring Margin, where thickened, probably underplated, crust is mapped; yet other interpretations exist (cf. Gernigon et al. 2004; Peron-Pinvidic et al. 2013; Funck et al. 2014; Eagles et al. 2015). Seaward of the COB of the outer Vøring Margin, a topographic high underlain by up to 15 km thick crust with velocities less than 7.2 km s⁻¹ constitutes the Vøring Spur (Breivik et al. 2008; Gernigon et al. 2009).

The Bivrost lineament marks the northern end of the Vøring Margin (Planke et al. 1991; Blystad et al. 1995; Olesen et al. 2002; Mjelde et al. 2003; Mjelde et al. 2005a). This lineament constitutes a distinct change in petrophysical properties (density, seismic velocity, susceptibility) in the offshore continental domain (Ebbing et al. 2006) but has no equivalent expression (e.g. a fracture zone) in the oceanic domain (Olesen et al. 2007). The northward located Lofoten Margin narrows from more than 250 km in the south to c. 50 km at the Senja Fracture Zone in the north. The Lofoten–Utrøst ridge system separates the inner Vestfjorden Basin from the outer Lofoten Basin. Crustal thickness undulates between 15 and 25 km below the Lofoten Archipelago, Lofoten Ridge and Utrøst Ridge and rapidly thins to less than 10 km under the Lofoten Basin (Avedik et al. 1984; Goldschmidt-Rokita et al. 1988; Mjelde et al. 1993). A pronounced shallowing of the Moho has been mapped under the Lofoten Ridge (Mjelde et al. 1993). In contrast, the Røst High to the west is associated with crustal thickening. The transition to the oceanic domain is relatively abrupt but seismic imaging is obstructed by extensive volcanic flows. No LCBs have been mapped here.

Numerous geophysical surveys have been conducted over the past decades along the mid-Norwegian margin. A comprehensive tectonostratigraphic atlas for the entire NE Atlantic region has recently been compiled based on reflection and refraction seismic experiments, magnetic and gravity surveys, as well as borehole and heatflow measurements (the NAG-TEC atlas; Hopper et al. 2014). Datasets that did not enter the compilation include the vast amounts of proprietory data (2D/3D reflection seismic, borehole data) as well as a number of still unpublished aeromagnetic datasets (e.g. Gernigon et al. 2012b; Brönner et al. 2014).

The datasets relevant for the isostatic study presented here are the bathymetry, sedimentary and crustal horizons, as well as densities from a recent gravity inversion model (Haase et al. 2016). Their model covers the entire NE Atlantic and aims to refine the Moho horizon where it is not or only poorly covered by seismic profiles. It also includes the sub-Moho lithospheric thermal structure to first order. All datasets are described in detail in Haase et al. (2016) and are briefly summarized here for completion. The layers and respective density values of the model are listed in Table 1. The entire model has a lateral resolution of 10 × 10 km and extends to 250 km depth. A subset of this model is used for the isostatic calculations of the mid-Norwegian margin, covering an area of 1200 × 1550 km.

The newest, most detailed bathymetric grid consists of the IBCAO dataset (Jakobsson et al. 2012) and a more local compilation along the Norwegian margin (Olesen et al. 2010). For the gravity inversion and the current isostatic study, less resolution is required and the DTU10 elevation model is used, which is based on satellite altimetry (Andersen et al. 2010; Andersen 2010).

The top basement is derived from the total sediment sequence thickness, compiled from existing datasets and seismic refraction data compiled in the NAG-TEC atlas (Funck et al. 2014; Hopper this volume, in prep). Individual sedimentary layers are not differentiated, yet a vertical division into 2 km thick sublayers is introduced to be able to represent a density structure. The total sedimentary sequence also contains volcanic flows and in places, the top basement might merge with the top of the basalts rather than the top of the crystalline crust, especially around the continent–ocean transition. This should be kept in mind when interpreting isostatic compensation.

The Moho along the mid-Norwegian margin is constrained by numerous refraction seismic profiles (locations shown in Fig. 1). A refraction seismic-based grid for the Moho depth was constructed for the entire NE Atlantic by Funck et al. (2016) using kriging techniques supplemented by bathymetric and gravity data. This Moho horizon was further refined by gravity inversion (Haase et al. 2016), adhering to the seismic uncertainties along the seismic profiles. In the density model of the gravity inversion, the crystalline crust is divided into two layers, consistent with many seismic observations from thinned continental domains. Whereas the intra-crustal boundary is at constant 20 km depth in the onshore domains, it separates the offshore continental crust into equally thick upper and lower crustal layers (prior to inversion and adjustments of the Moho). This is a large simplification compared with the highly variable depth of this intra-crustal boundary as seen in seismic profiles, but such variations could not be included consistently across the entire study area owing to insuffcient seismic coverage. The high-velocity, high-density LCBs are also included according to the record of refraction seismic data (e.g. Mjelde et al. 2009a). A division into continental and oceanic domain is required, and the COB follows the interpretation of Funck et al. (2014), using magnetic lineations, tectonic reconstructions and the landward-most points of seismically defined oceanic crust. The lithosphere-asthenosphere boundary (LAB) is included as the 1300°C isotherm in the oceanic domain, based on the age of the crust by Gaina et al. (2016) and isotherm calculations after Sandwell (2001). In the continental domain, the LAB is also modelled as an isotherm, based on the definitions by Artemieva & Mooney (2001) and Artemieva (2006). The transition between the two regions has been smoothed to avoid long-wavelength effects (Haase et al. 2016). The LAB here represents the boundary between conductive and convective heat transfer, such that temperatures (and here also densities) are constant below the LAB. The sub-Moho lithosphere is subdivided along the 700 and 900°C isotherms with respective depths calculated according to Sandwell (2001, and references therein).

Sediment densities increase stepwise with burial depth mimicking an exponential trend. The respective densities for the 2 km thick intervals are 2200, 2270, 2340, 2650 and 2700 kg m⁻³, where the last value represents any sediments buried more than 8 km. These density values are typical for clastic sediments but do not consider the denser layers of basalts and volcanics that are present in parts of the study area. Crustal densities are 2750 and 2950 kg m⁻³ for the continental upper and lower crust, respectively. The density in the oceanic domain of 2850 kg m⁻³ is consistent with previous gravity models of the mid-Norwegian margin (e.g. Ebbing et al. 2006). The lower continental crust extends slightly under the oceanic crust and tapers over a distance of c. 40 km to ensure a gradual density change across the COB (representing a zone of mafically intruded crust). This is a crude approximation of the actual continent–ocean transition zone, of which the velocity structure has been seismically imaged in many places along the margin. The lower crustal bodies are assigned a constant density of 3100 kg m⁻³, although a range of densities has been suggested by previous gravity models (3000–3300 kg m⁻³, e.g. Mjelde et al. 2009a; Kvarven et al. 2014; Nirrengarten et al. 2014). Densities below the Moho consider the thermal regime to first order and decrease stepwise at the 700, 900 and 1300°C isotherm from 3300 kg m⁻³ at the Moho down to 3260, 3240 and 3200 kg m⁻³, respectively. The LAB of this model is thus only a thermal but not compositional boundary. The densities do not include pressure-related density increases through compaction or phase changes. Furthermore any lateral density changes are omitted which originate, for example, from the effects of the Iceland plume or from compositional variations.

We use the density model exactly as it has been presented by Haase et al. (2016). Because the density model was developed to match the regional gravity data of the entire NE Atlantic, but our study covers the smaller area of the mid-Norwegian margin, certain limitations and oversimplifications are expected and have to be kept in mind for the interpretation of the results. By applying isostatic considerations, we can point out areas where the density model does not fit. These might represent additional lateral density variation which had not been included in the subsurface model, located in sediments, crust or uppermost mantle. We do not attempt to improve the density model in this study but suggest a future joint inversion using gravity and isostatic calculations. Isostatic tests alone may not allow for new insights into the crustal structure of the margin but when evaluated in context of additional data – here for example regional seismic profiles – it can give clues as to the previously undetected variations in density, composition or thermal regime.

The gravity misfit for the regional density model is 22 mGal standard deviaton with the largest residuals remaining around the Greenland–Iceland–Faroe corridor, some oceanic ridges, as well as around the continent–ocean transitions. The long-wavelength residuals can be attributed to large-scale lateral variations in temperature and composition introduced by the Iceland plume and imaged, for example, by teleseismic tomography (Bijwaard & Spakman 1999; Weidle & Maupin 2008). Many of the smaller-wavelength residuals must be attributed to the crustal structure, which is only coarsely represented in this regional model. Deviations are thus to be expected around the very coarsely represented structures of the LCBs, the COB, the intra-crustal boundary and compositional density variation in the lithospheric mantle.

Long-wavelength residuals as introduced by an incomplete thermal model, the neglect of sublithospheric density variations or the forced smoothing of the LAB across the continent–ocean transition will affect the calculated long-wavelength bathymetric features but not the smaller-scale ones, which are subject of this study. We thereby acknowledge that the thermal regime assumed for the density model is only a first-order approximation. Additional modifications would lead to long-wavelength changes, which are not considered relevant for the smaller-scale structures discussed in this study.

In order to check the density model for isostatic equilibrium, we need to calculate the buoyancy of each vertical column of the model (each 10 × 10 km) and how it could be balanced by vertical movement, i.e. a modified bathymetry. We calculate the pressure at a compensation depth z_comp and its deviation from a reference pressure P_ref. We then calculate the flexural response of this load distribution with a given flexural strength of the lithosphere (effective elastic thickness) using the open source software tisc by Garcia-Castellanos (2002) and Garcia-Castellanos et al. (2003). The resulting deflection d_i thus represents the difference between the isostatically compensated bathymetry of the model and the true bathymetry, in the following sections referred to as elevation residual.

In the end-member case where lithosphere has no flexural strength and the system is locally isostatically balanced, each rock column is considered to be balanced independently of its neighbours and

where the index i refers to an individual rock column, P_i is the pressure at compensation depth, P_ref is a reference pressure, g is gravitational acceleration and Δρ = ρ_bottom−ρ_top is the density difference between material below the column (mantle) and above (water or air). The pressure at compensation depth is calculated by integrating the densities of each rock column from compensation depth to the surface

where h_i is the elevation. The deflection d_i is positive if P_i < P_ref and negative if P_i > P_ref. The new (forward calculated) isostatically balanced elevation h^bal is then given by

If isostatic adjustment raises a column from below sea-level (h⁰ < 0) to above sea-level (h^bal > 0), the equations need to be modified:

where ρ_a and ρ_w are the densities of the asthenosphere (sublithospheric mantle) and water, respectively. Air density is considered to be negligibly small.

If flexural isostasy is assumed, the load ΔP that causes a deflection of the lithosphere, must be convoluted with a flexural filter. This filter effectively suppresses short wavelengths and depends on the flexural rigidity D (Turcotte & Schubert 1982; Watts 2001). In the spectral domain, it can be written as

where Y(k) and L(k) are the Fourier transforms of the deflection and the load, respectively. In our case, the load is given by the pressure distribution P_i at the base of the compensation depth.

Φ_e(k) represents the flexural response function:

with D being the flexural rigidity of an elastic plate. This parameter can also be expressed through the effective elastic thickness T_eff, the Young’s modulus E and the Poisson ratio ν:

Values of T_eff range from approximately 5–100 km (e.g. Tesauro et al. 2013). These values are determined from natural loading experiments (removal of ice loads, large lakes, seamounts) or spectral correlation of gravity and topographic data. For oceanic domains, the effective elastic thickness roughly corresponds to the depth of the 450°C isotherm (Watts 2001). For continental domains T_eff is usually much smaller than the thickness of the lithosphere measured with seismic methods (or even the crust in some places). These different values of seismic (or thermal) and effective elastic thickness of the lithosphere arise because the lithosphere is not a homogenous elastic layer. Instead, it consists of layers of multiple lithologies with a net response that resembles that of a single, yet thinner, elastic layer. We test a range of effective elastic thicknesses from 5 to 20 km.

The compensation depth z_comp is chosen to be 250 km, which is far below the lithosphere–asthenosphere boundary (usually less than 100 km in oceanic and offshore domains) and therefore well within the region of viscoplastic deformation.

The reference pressure is commonly calculated from a reference column. This is often chosen to be a column of oceanic lithosphere with well-known density distribution or a pure asthenospheric column representative of a mid-oceanic ridge. In this study, the density structure of the lithospheric mantle is highly simplified and thereby not directly comparable with the absolute values. We thus use the average pressure at compensation depth across the study area as the reference pressure. This reference pressure is representative for large parts of the oceanic domains and the resulting elevation adjustment will here accordingly be zero. This approach also leads to fairly equally distributed positive and negative elevation adjustments along the margin and is therefore well suited to investigate the local density structures in the study area. A different reference pressure would lead to an overall shift in the calculated bathymetry and thus in the residuals.

When comparing a forward calculated dataset with an observed one, it is important to consider data content of similar wavelength. The existing bathymetric data has thus been regridded to a 10 × 10 km grid in order to be comparable with the forward calculated bathymetry (Fig. 2a).

Figure 2b and c show the calculated bathymetry with T_eff = 5 km and the corresponding residuals, respectively. The overall elevation is well reproduced, yet residuals of several hundreds of metres remain in the continent–ocean transition area around the Jan Mayen microcontinent, the SW Barents Sea and the mid-Norwegian margin. Calculations with higher effective elastic thickness of T_eff = 20 km (Fig. 2d) highlight the more regional residuals of longer wavelength than in the case with T_eff = 5 km. Calculations based on the pre-inversion density model with a purely seismically derived Moho (not shown here) display the same patterns of residuals, yet much higher magnitudes. This shows that, although isostasy was not considered in the gravity inversion, the modifications are in line with isostatic balancing. Both gravity and elevation are thus sensitive to the same density structure.

Applying a constant effective elastic thickness for the entire study area is a strong simplification. Lateral variations must exist for oceanic and continental domains, for highly stretched and very thick crust. Furthermore, a temporal change of the rigidity (weak when young and hot, strong when old and cold) is to be expected. This would lead to different compensation of the crystalline crustal structure and the subsequent basin fill, as explained, for example, by Kimbell et al. (2004). In this study, where only the regional trends are analysed, a constant effective elastic thickness is appropriate. But further, local investigations need to take potential spatial and temporal variations into account.

The compensation depth z_comp to be used in isostatic modelling is not a well-defined value. Here, we take the view that the convective asthenospheric mantle provides no isostatic support and thus choose a compensation depth below the LAB at 250 km. Since no density variations are expected below this depth, an even larger depth for z_comp would not change the result. Results of calculations with shallower compensation depths are shown in Figure 3.

A compensation depth of 100 km (Fig. 3a) yields a similar distribution of the local elevation residuals along the mid-Norwegian margin. However, the long-wavelength residuals, in particular around the mid-oceanic ridges, are much higher. This reflects that the thermal structure of the oceanic domain and the passive margin comprises major density variations at 100 km depth around the spreading ridges, but relatively uniform densities deep under the continental margin. Kimbell et al. (2004) suggest that a compensation depth of 90 km is suffcient, yet they employ a slightly different thermal structure with maximum temperatures of 1100°C and no significant lateral density variations at greater depths.

Figure 3b shows the elevation residual if isostatic compensation is assumed to occur at Moho depth. This is here implemented by a constant density of the underlying mantle of 3300 kg m⁻³. The residuals clearly show the influence of the thermal structure on elevation. Omitting the lithospheric mantle in isostatic studies can thus give undesired long-wavelength residuals. This is even more the case if compositional density changes are present within the lithospheric mantle. These are also not included in the density model used here. The calculations furthermore do not include the dynamic pressure exerted by flowing mantle (strict definition of dynamic topography), merely the thermally induced density differences in the lithosphere and asthenosphere.

We now focus on the short-wavelength residuals (T_eff = 5 km) displayed on the mid-Norwegian margin where a multitude of refraction seismic profiles are available that constrain the subsurface model and allow us to interpret the elevation residuals. Four areas of high residuals are highlighted (Fig. 4a) and will be discussed in the following. These ‘unbalanced’ areas are also marked on the map of gravity residuals (Fig. 4b), which result from the inversion by Haase et al. (2016).

Area 1 shows a large mass deficit in an area where seismic data coverage is high and subsurface geometries should therefore be well determined. The gravity response of the subsurface model equally indicates a mass deficit in the density model here (Fig. 4b). As Haase et al. (2016) point out, the densities assumed for the sedimentary layers are representative for clastic sediments but do not consider the contribution of sills and lava flows. These flood basalts and sills are massive along the mid-Norwegian margin (location shown as green transparent area on Fig. 5) and additionally hamper seismic imaging on the outer margin. Thus, seismic penetration is poor and the structures below the volcanics – in particular the thickness of the basalts themselves and of the underlying sediments as well as the Moho depth – are poorly resolved. These uncertainties could justify significant modifications to the subsurface model, which could yield a better fit to gravity data and isostatically compensated bathymetry. The current mass deficit in the model could, for example, be compensated by thinner sediments (shallower basement) or a shallower Moho.

In addition to the volcanics near the surface, magmatic activity has probably affected the lower crust as well. The seismic velocities around Area 1 are increased with respect to the more landward, stretched continental crust (Funck et al. 2014). These high velocities probably go hand in hand with higher densities, but regions of transitional crust are not represented in the current density model. An inversion of the gravity residual for crustal densities also indicates a need for higher densities in this area (Haase et al. 2016).

An alternative explanation of the regional residuals can be found in the LCBs (high-velocity, high-density lower crustal bodies), which in large parts of the mid-Norwegian margin correlate with the extent of the volcanic flows (Reynisson et al. 2010). Figure 5 illustrates that their extent fairly well matches the extent of Area 1, where residuals for both bathymetry and gravity are high. Calculations with an LCB density of 3150 kg m⁻³ instead of 3100 kg m⁻³ show only a slight improvement (10–20% of the gravity and elevation residual). The LCB density increase alone is thus not suffcient to explain the residuals. Higher densities for the sedimentary section in this area strongly improves the gravity residuals but has little effect on the elevation residuals.

Seaward of Area 1 exists a small area of positive elevation and gravity residuals, which is crossed by several refraction seismic profiles (Mjelde et al. 2001, 2005b; Mjelde et al. 2007; Breivik et al. 2008). An exemplary refraction seismic profile is shown in Figure 6 together with the crustal structure of the density model. The elevation and gravity residuals correlate with the positive bathymetry outboard of the mapped COB (Figs 1 & 2a). This indicates that the density structure across the COB may not be as simple as included in the present subsurface model and perhaps even that a further seaward located COB might be more appropriate (Gernigon et al. 2012a; Peron-Pinvidic et al. 2013; Eagles et al. 2015).

An overall mass excess is seen in an elongated, coast-parallel region, which is crossed by a few refraction seismic profiles (Fig. 5). The gravity residual indicates a slight mass deficit for Areas 3a and 3c and a generally good fit for the central Area 3b. Since elevation and gravity residuals are derived from the same density model, these seemingly contradicting results must be explained by the different methods. Isostatic adjustment should be less sensitive to small wavelength features as their respective loads can be supported by a suffciently rigid lithospheric plate. Yet, even with a stronger lithosphere (T_eff = 20 km), the calculated elevation is too low (Fig. 2d) and a significant residual remains. This indicates that the excess mass in the density model is a more regional feature. Because the gravity signal is more sensitive to shallower depths and the isostatic compensation does not depend on the vertical distribution of mass within a model column, the mass excess revealed in Areas 3a and 3c probably stems from greater depth.

The velocity structure of the refraction seismic profile crossing Area 3a is shown in Figure 7 together with the crustal structure of the density model. The continental shelf is underlain by thinned continental crust that forms several basement highs, neighboured by sedimentary basins. Whereas the Røst High and Marmæle Spur are underlain by thickened crust, the Lofoten Ridge is located above a Moho high. This scenario is in disagreement with local isostasy. Similar to these results, Reynisson et al. (2010) also calculated that for achieving complete isostatic balancing, the respective Moho needed to be several kilometres deeper than the seismically imaged Moho. A deeper Moho however, would not be compatible with the gravity data but could be compensated, for example, by a different composition of the overlying basement high (Reynisson et al. 2010). It was originally proposed that the Lofoten Ridge is a large-scale horst structure (Mjelde et al. 1993) and that the isostatic imbalance is too small to be locally compensated. Recent studies (yet still based on older seismic data) interpret a high-velocity body in this area (Mjelde et al. 2013). The reflector originally interpreted as the Moho would then be the top of a high-density body with mantle-type velocities. While the refraction seismic data did not detect lateral variations of subcrustal velocities, a deep mantle reflector c. 10 km below the detected Moho suggests that some deeper structures are still present below the Lofoten Ridge (Mjelde et al. 1993). Reflection seismic data show yet a different Moho horizon with an up to 4 km shallower Moho east of the Lofoten Ridge (Fig. 7, Mjelde et al. 1993). This shows that seismic methods do not unambiguously map the deep structures here, and differences in determining seismic velocities may play a role.

A similar scenario can be brought forward for Area 3c (Fig. 8). The refraction seismic line of Mjelde et al. (2008) shows again that a relative basement high is underlain by a Moho high, but the gravity inversion modified this Moho high to a generally landward-dipping horizon. Further outboard, the crust is drastically thinned and comprises a high-velocity lower layer, which pinches out towards the shelf. Seismic investigations from the last few years (Kvarven et al. 2014; Nirrengarten et al. 2014) and older, expanding spread profiles and sonobuoy data from neighbouring areas (Olafsson et al. 1992) describe an additional high-velocity lens in the area showing the mass excess (Fig. 8b, c). Slightly higher subcrustal velocities were also detected in the refraction seismic profile entering the current subsurface model (Fig. 8a, Mjelde et al. 2008), but only much more recently interpreted by the same authors as lower crustal eclogitic bodies (Mjelde et al. 2013). This high-velocity, high-density lower crustal lens is thought to be a remnant of old eclogized crust which was once part of the Caledonian orogen but broke up and transformed during rifting (Nirrengarten et al. 2014). While these new interpretations did not enter the subsurface starting model for the gravity inversion, the respective results nevertheless show a deeper Moho in this area (Fig. 8a).

For the central Area 3b, no indications for LCBs have been brought forward (Breivik et al. 2011). Here, the Trøndelag Platform is underlain by relatively thick continental crust (more than 20 km) and the simple structure assumed in the gravity model is probably not adequate to capture the effects of internal layering and inhomogeneities on the gravity and isostasy signal. A better mapping of the margin structure and more detailed modelling is needed here.

While the suggested Moho modifications for Area 3 (additional LCB, lower Moho) would be sufficient to reduce the gravity residual, a positive elevation residual still remains that points to excess mass (missing buoyancy). A local, further lowering of the Moho could remove this mismatch, yet the densities of this lower crustal structure remain unresolved. Whereas isostatic considerations require low crustal densities, gravity modelling by Kvarven et al. (2014) suggests densities as high as 3300 kg m⁻³. On the other hand, Nirrengarten et al. (2014) suggest a density of 3100 kg m⁻³, closer to the values used in our study. A separate explanation could be lateral, compositional density variations in the uppermost mantle. Indications for such variations are seen in the refraction seismic data (see Fig. 8a, b) and discussed, for example, by Kimbell et al. (2004) for an area NW of Scotland.

Area 4 exhibits a strong mass excess, which is seen in both the elevation and the gravity residuals (Fig. 4a, b). This region is situated in the corner between the outer Møre and Vøring margins in the outer part of the Jan Mayen corridor. Reynisson et al. (2010) mapped an equal mass excess (isostatic Moho shallower than refraction seismic Moho) slightly landward, on the outer Møre Margin. The area under investigation is considered to extend across the COB, covering both oceanic and transitional crust. It is also highly affected by the shearing and differential opening along the Jan Mayen Fracture Zone (Talwani & Eldholm 1977; Skogseid & Eldholm 1987; Gernigon et al. 2015). Refraction seismic profiles that partially cover this area are shown in Figure 9. The crust is generally thin (c. 7 km) and no LCB has been mapped. For the area under discussion here, it is also plausible that the crustal densities are not accurately represented in the model, namely that the actual densities are lower than those assumed. This could be explained by a different, farther seaward located course of the COB or a different, more abrupt, density transition across the COB.

Isostatic calculations are an easy way to provide additional support for seismic interpretation and gravity modelling. Especially where new concepts meet older data interpretations, the consideration of isostatic elevation residuals may provide additional insight. This is, for example, the case for the eclogitic bodies along the proximal mid-Norwegian margin that have been discussed in recent publications and whose presence is now also confirmed by isostatic calculations. Even where no other interpretations are available, elevation residuals may point out areas that require renewed investigation through more data or renewed data analysis. This is, for example, the case for the outer Vøring Marginal High (Area 2) or the complex area SW of the Vøring Margin (Area 4).

Elevation residuals can stem from inadequate isostatic calculations (e.g. inappropriate effective elastic thickness T_eff or compensation depth) or an inadequate subsurface density model. If the erroneous density structure is relatively shallow, it will also affect the gravity residuals (see e.g. Area 1). If it is deep, it will probably remain undetected by the gravity calculations. An erroneous subsurface model can result from wrongly assigned densities, wrong interpretation of refraction seismic data and poor quality or even lack of data, which in turn may lead to erroneous data interpretations. The latter seems to be the case for the deep LCBs along the mid-Norwegian margin (Areas 3a and c).

The deep structure, in particular the mantle densities, and the choice of compensation depth have a large effect on the long-wavelength elevation residual. The choice of the reference column, however, determines only the absolute, not the relative, values. In areas where thermal anomalies create lateral density changes, these need to be taken into account. Isostatic compensation is usually not reached at the Moho level.

The mantle densities of the current subsurface model decrease with depth, following the thermal trend. They do not include pressure-dependent density increase with depth as is the case in the natural system. This pressure dependency is a nearly uniform density increase without significant lateral variations. It therefore does not affect the gravity residuals nor the relative elevation residuals. Lateral density variations in the uppermost mantle, however, which could originate from different degrees of melting and depletion, can have a major effect.

Isostatic calculations comprise a simple, additional tool to validate interpretations of deep seismic data and gravity models. They can either be applied separately, as in this study, or included in a joint data analysis or joint inversion with gravity. Isostasy is only sensitive to regional structures, but, different from other methods, is equally sensitive to shallow and deep structures.

With this simple study we hope to have demonstrated that interpretation of a single refraction seismic profile can give misleading results, where high uncertainties arise from resolving deep structures (e.g. Moho, LCB and mantle reflections) or sub-basalt features. Different seismic methods (here expanding spread profiles, reflection seismic) and additional geophysical datasets (gravity, bathymetry) should always be considered for the interpretation.

This work formed part of the NAG-TEC project and we acknowledge the support of the industry sponsors (in alphabetical order): Bayerngas Norge AS; BP Exploration Operating Company Limited, Bundesanstalt für Geowissenschaften und Rohstoffe (BGR); Chevron East Greenland Exploration A/S; ConocoPhillips Skandinavia AS; DEA Norge AS; Det norske oljeselskap ASA; DONG E&P A/S; E.ON Norge AS; ExxonMobil Exploration and Production Norway AS; Japan Oil, Gas and Metals National Corporation (JOGMEC); Maersk Oil; Nalcor Energy Oil and Gas Inc.; Nexen Energy ULC, Norwegian Energy Company ASA (Noreco); Repsol Exploration Norge AS; Statoil (UK) Limited, and Wintershall Holding GmbH. We furthermore would like to thank the reviewers for their critical evaluation of the manuscript and Laurent Gernigon for further improvements of the text.

This work formed part of the NAG-TEC project and we acknowledge the support of the industry sponsors (in alphabetical order): Bayerngas Norge AS; BP Exploration Operating Company Limited, Bundesanstalt für Geowissenschaften und Rohstoffe (BGR); Chevron East Greenland Exploration A/S; ConocoPhillips Skandinavia AS; DEA Norge AS; Det norske oljeselskap ASA; DONG E&P A/S; E.ON Norge AS; ExxonMobil Exploration and Production Norway AS; Japan Oil, Gas and Metals National Corporation (JOGMEC); Maersk Oil; Nalcor Energy Oil and Gas Inc.; Nexen Energy ULC, Norwegian Energy Company ASA (Noreco); Repsol Exploration Norge AS; Statoil (UK) Limited, and Wintershall Holding GmbH. We furthermore would like to thank the reviewers for their critical evaluation of the manuscript and Laurent Gernigon for further improvements of the text.