Tilting and Flexural Stresses in Basins Due to Glaciations—An Example from the Barents Sea

Many of the Earth’s sedimentary basins are affected by glaciations. Repeated glaciations over millions of years may have had a significant effect on the physical conditions in sedimentary basins and on basin structuring. This paper presents some of the major effects that ice sheets might have on sedimentary basins, and includes examples of quantifications of their significance. Among the most important effects are movements of the solid Earth caused by glacial loading and unloading, and the related flexural stresses. The driving factor of these movements is isostasy. Most of the production licenses on the Norwegian Continental Shelf are located inside the margin of the former Last Glacial Maximum (LGM) ice sheet. Isostatic modeling shows that sedimentary basins near the former ice margin can be tilted as much as 3 m/km which might significantly alter pathways of hydrocarbon migration. In an example from the SW Barents Sea we show that flexural stresses related to the isostatic uplift after LGM deglaciation can produce stress changes large enough to result in increased fracture-related permeability in the sedimentary basin, and lead to potential spillage of hydrocarbons out of potential reservoirs. The results demonstrate that future basin modeling should consider including the loading effect of glaciations when dealing with petroleum potential in former glaciated areas.


Introduction
It is commonly accepted that major changes in Earth's climate started in Gelasian (~2.6 million years ago), and initiated the growth of ice sheets in the Northern Hemisphere (e.g., [1]). Over the last 2.6 million years, there have probably been more than 30 glaciations with cycles of ice sheets advancing and retreating on time scales of 40,000 to 100,000 years. A glaciation cycle of 40,000 years dominated in the early period, while the 100,000 year cycle has dominated in the past million years or so [2]. The most extensive glaciations are traditionally referenced in the last million years (Late Quaternary; cf. [3]) and the covers on North America and Europe have been estimated to be around 2-4 km thick near the centers of maximum ice accumulation (cf. e.g., [4]), and tapered towards the ice sheet margins. Repeated episodes of growth and withdrawal of huge ice sheets, glacial erosion, and associated sediment deposition lead to geologically dramatic changes in surface loading. The lithosphere responds to these changes by subsiding under loading, and by uplifting when loading is removed. The driving mechanism for this is isostasy.
Differential vertical movement of the lithosphere related to glacial isostasy lead to repeated tilting of sedimentary formations and potential petroleum reservoirs therein, which may have greatly affected hydrocarbon migration pathways in the former glaciated areas [5]. In addition, the upward and downward bending of the lithosphere due to glaciation leads to flexural stresses [6] likely to affect faults and their permeability, which could add to the changes in hydrocarbon migration pathways.
In parts of the Barents Sea, there are clear indications of hydrocarbon spillage (e.g., [7]). Some of the hydrocarbon spillage is linked to isostatic movements due to cycles of Pliocene-Pleistocene ice sheet loading/unloading and glacial sediment redistribution [8][9][10][11][12]. It is suggested that the recent discovery of the giant Johan Sverdrup oil field in the Norwegian North Sea was oil charged during Quaternary, and that the area more than likely underwent tilting and possible leakage several times over the last one million years [13]. Other research, related to the North Sea, showed that influence of ice sheet loading does not result in a significant tilting [14], which may be the result of smaller gradients of the ice thickness. Anyway, detailed control on the glacial isostasy is important and so far an insufficiently utilized factor for identification of the remaining hydrocarbon resources in sedimentary basins formerly covered by ice sheets.
Faulting related to the last glaciation is reported from northern Scandinavia, where several fault with up to 30 m high fault scarps [15]. Examples include the 80 km long Stuoragurra fault in northern Norway [16] and the 150 km long Pärvie fault in northern Sweden [17]. There is also evidence of such faults in Denmark [18], Germany, Poland, and the United Kingdom ( [19][20][21][22][23]). Thus, such faults are found in the center of former ice sheets, as well as at their margins and beyond [24,25].
For the offshore areas, the situation is more complicated. However, the number of seismic events recorded in the Barents Sea in the last 20 years is significant ( [26,27]). For the Barents Sea, the impact of the glacial isostasy on triggering earthquakes has not yet been investigated. A better understanding of the origin of the seismicity will help to recognize which faults and related areas exhibit a risk of leakage of hydrocarbons from traps.
Stress changes induced by glaciations are related to two main effects: 1) direct load introduced by the weight of the ice on the surface and 2) flexural load caused by the subsidence and uplift of the lithosphere due to glacial isostasy [10,28]. Horizontal stress due to stress migration, which is depending on loading time and mantle viscosity, is an additional effect that is modeled by Steffen et al. [29]. Grollimund and Zoback [30] investigated, by a 2D approach, lithospheric flexure as an alternative mechanism to explain the local stress perturbations observed in the northern North Sea. They found that flexural stresses due to glaciations/deglaciation could explain the observed lateral stress variations.
Johnston [31] used the Mohr-Coulomb theory for the study of faulting and explained the lack of seismic activity in areas of large continental ice sheets. Turpeinen [32] use finite-element models of rheologically layered lithosphere with a thrust fault to investigate its slip evolution during the glacial/ postglacial period. Their modelling results suggest that the rate of thrusting decreases during ice sheet loading and strongly increases during deglaciation. Steffen et al. [33] introduced a new approach that is based on a Glacial Isostatic Adjustment (GIA) model. They combined the GIA model with a fault surface to investigate the fault slip and fault activation time during a glacial cycle. Their model allows estimation of the fault throw for former glaciated areas. One of their conclusions is that faults start to move close to the end of deglaciation.
Our approach is flexural isostasy modelling, and we are dealing with two effects of the glaciations which we believe future basin modelling has to consider:

1.
Possible glacially induced tilting of reservoirs. Here, we calculate the glacially induced tilting of the Norwegian Continental Shelf due to the last ice age, with particular focus on the SW Barents Sea.

2.
Assessment of the flexural stresses induced by the lithospheric flexure related to glacial isostasy and tilt, and their effect on faults during glaciations and interglacials, here illustrated by an example from the SW Barents Sea where the stress modelling is limited to the flexural stress effects of deglaciation after LGM. The stress effects of vertical ice load and related stress migration are not calculated here.

Glaciations
The last glacial period began about 110,000 years ago and ended 10,000 years ago. A model of the deglaciation history for our area is given in [34]. The glaciations that occurred during this glacial period covered many areas of the Northern Hemisphere. The extent of the Eurasian Ice Sheet during the Last Glacial Maximum (LGM) has been long debated in the literature. Researchers [35] argue for a highly asynchronous ice sheet with regional differences in the age of maximum ice-sheet positions. However, in this paper we use the LGM (~20,000 BP) reconstruction of [36], which suggests the SW Barents Sea was glaciated from 24,000 to 18,000 BP. In the following calculations, we assume the Earth reached isostatic equilibrium during this period.

Ice Thickness
Glacial isostasy is calculated based on the palaeo ice extent (from dated marginal moraines) and models of the ice thicknesses. Observational data increasingly constrain the extent of the Quaternary ice sheets. The thickness and volume of these ice sheets are much harder to reconstruct and generally need to be inferred from indirect evidence and modeling.
Several models of ice thickness of LGM and the late-glacial deglaciation are published; the loading scenarios are produced using different methods and sets of constraints. Some of them consider thermomechanical coupling (e.g., [37]). A number of models of the LGM ice thickness are presented in [38]. We follow the work of [36]; their LGM extent model is shown in Figure 1. period covered many areas of the Northern Hemisphere. The extent of the Eurasian Ice Sheet during the Last Glacial Maximum (LGM) has been long debated in the literature. Researchers [35] argue for a highly asynchronous ice sheet with regional differences in the age of maximum ice-sheet positions. However, in this paper we use the LGM (~20,000 BP) reconstruction of [36], which suggests the SW Barents Sea was glaciated from 24,000 to 18,000 BP. In the following calculations, we assume the Earth reached isostatic equilibrium during this period.

Ice Thickness
Glacial isostasy is calculated based on the palaeo ice extent (from dated marginal moraines) and models of the ice thicknesses. Observational data increasingly constrain the extent of the Quaternary ice sheets. The thickness and volume of these ice sheets are much harder to reconstruct and generally need to be inferred from indirect evidence and modeling.
Several models of ice thickness of LGM and the late-glacial deglaciation are published; the loading scenarios are produced using different methods and sets of constraints. Some of them consider thermomechanical coupling (e.g., [37]). A number of models of the LGM ice thickness are presented in [38]. We follow the work of [36]; their LGM extent model is shown in Figure 1.  [36]) and maximum thickness model for the Last Glacial Maximum (LGM), from [3]. The production licenses (from www.npd.no) in the Norwegian continental shelf are shown by grey polygons.

Glacial Isostasy and Tilting of Reservoirs
The Earth's response to deglaciation shows that the lithosphere acts as an elastic shell overlying a viscous mantle. If a load is applied to the elastic lithosphere, part of the applied load will be  [36]) and maximum thickness model for the Last Glacial Maximum (LGM), from [3]. The production licenses (from www.npd.no) in the Norwegian continental shelf are shown by grey polygons.

Glacial Isostasy and Tilting of Reservoirs
The Earth's response to deglaciation shows that the lithosphere acts as an elastic shell overlying a viscous mantle. If a load is applied to the elastic lithosphere, part of the applied load will be Geosciences 2019, 9, 474 4 of 15 supported by the elastic lithosphere, and part by buoyant forces of the mantle underneath, acting through the lithosphere.
In the literature, there is disagreement on the elastic thickness of the lithosphere; the elastic thickness of the lithosphere of Fennoscandia in GIA (Glacial Isostatic Adjustment) modeling varies from 30 to 160 km [39]. Our modeling of the post-glacial isostatic response, calibrated with the observations onshore Norway, gave best fit with an elastic thickness of 30 km; see [34]. Calibration with the Barents Sea gave similar results; see [40]. The method used in the isostatic calculations is described in [41] and in the Appendix A. The spatial resolution is 10 km.
There will be significant subsidence of potential sedimentary basins formerly covered by ice sheets. Sedimentary basins that are located in the periphery of the former ice sheet will be notably affected by tilting of the Earth's surface as the subsidence (or uplift after the melting) is gradually decreasing towards the periphery of the former ice sheet. Therefore, the largest tilts of the Earth's surface are found in the more peripheral areas of former ice sheets. This is clearly seen along the entire Norwegian coast and in the Barents Sea ( Figure 2). Sedimentary basins located in the peripheral areas in the SW Barents Sea could be tilted by up to 2.7 m/km due to the LGM ice sheet (Figure 2), and even more in the previous greater glaciations [3]. The ice sheets will act as seesaw during the total glaciation period (Quaternary); the Earth will subside under the ice sheets and rise again in interglacial periods. supported by the elastic lithosphere, and part by buoyant forces of the mantle underneath, acting through the lithosphere. In the literature, there is disagreement on the elastic thickness of the lithosphere; the elastic thickness of the lithosphere of Fennoscandia in GIA (Glacial Isostatic Adjustment) modeling varies from 30 to 160 km [39]. Our modeling of the post-glacial isostatic response, calibrated with the observations onshore Norway, gave best fit with an elastic thickness of 30 km; see [34]. Calibration with the Barents Sea gave similar results; see [40]. The method used in the isostatic calculations is described in [41] and in the Appendix. The spatial resolution is 10 km.
There will be significant subsidence of potential sedimentary basins formerly covered by ice sheets. Sedimentary basins that are located in the periphery of the former ice sheet will be notably affected by tilting of the Earth´s surface as the subsidence (or uplift after the melting) is gradually decreasing towards the periphery of the former ice sheet. Therefore, the largest tilts of the Earth´s surface are found in the more peripheral areas of former ice sheets. This is clearly seen along the entire m/km This process is, of course, important for migration modelling. Figure 3 shows how a simple structure behaves during one glacial cycle. It is assumed, for illustration purposes, that the structure is filled with hydrocarbons to the spill point at the onset of glaciation. When the lithosphere and the structure subside due to the glacial load, the reservoir is tilted hydrocarbons start to spill out of the structure. This goes on until the maximum tilting is reached ( Figure 3, upper right). The closure volume will be reduced. This process is, of course, important for migration modelling. Figure 3 shows how a simple structure behaves during one glacial cycle. It is assumed, for illustration purposes, that the structure is filled with hydrocarbons to the spill point at the onset of glaciation. When the lithosphere and the structure subside due to the glacial load, the reservoir is tilted hydrocarbons start to spill out of the structure. This goes on until the maximum tilting is reached ( Figure 3, upper right). The closure volume will be reduced.
When the lithosphere and structure are uplifted due to deglaciation, a new hydrocarbon-water contact (HCWC) will be established ( Figure 3, bottom). Paleo-HCWCs can theoretically be found at the spill point level for the positions before and after deglaciation. Residual oil can be found in the water zone down to the original spill point level. In a gas-filled reservoir with an oil leg originally, residual oil can also be found in the gas-filled part, limited by the uppermost HCWC.
The effect of the glacially induced hydrocarbon migration is determined by several factors: • The number of glacial events. Tilting will result from each isostatic event. A number of glacial events will create chaotic patterns of paleo contacts and residual oil which will be very difficult to interpret.

•
The geometry of the receiving area. If a saddle area is opened, the fluids can migrate into higher lying traps or be lost to the surface. A surrounding flat area, on the other hand, can become part of the trap during the subsidence, and fluid can migrate back into the trap under uplift.

•
Additional potential factors that determine the migration effect are: i) the geometry of the deformation, ii) the geometry of the hydrocarbon filled structures, iii) the initial fill of the trap, and iv) the spill point location and orientation of the structure.
The isostatic response on the Pleistocene sediment redistribution and ice loading has been evaluated for the Bjørnøyrenna Fault Complex in the Barents Sea [11]. The isostatic impact on hydrocarbon trap capacity changes and hydrocarbon maximum spillage was assessed in that study and the conclusion was that tilting together with gas volume expansion might have been responsible for some part of the hydrocarbon loss during the Cenozoic.
This process is, of course, important for migration modelling. Figure 3 shows how a simple structure behaves during one glacial cycle. It is assumed, for illustration purposes, that the structure is filled with hydrocarbons to the spill point at the onset of glaciation. When the lithosphere and the structure subside due to the glacial load, the reservoir is tilted hydrocarbons start to spill out of the structure. This goes on until the maximum tilting is reached (Figure 3, upper right). The closure volume will be reduced. Illustration of tilting during one glacial cycle (redrawn from [5]). Upper left shows the situation prior to tilting; upper right during tilting. Lower panel shows the situation after deglaciation. "Spill point" is the structurally lowest point in a hydrocarbon trap that can retain hydrocarbons. Once a trap has been filled to its spill point, further storage or retention of hydrocarbons will not occur. The illustration assumes the trap filled to the spill point. HCWC is the hydrocarbon-water contact.

Flexure-Related Stress Effects
As mentioned above, the flexural stresses that occur during glaciations will be related to the isostatic response of the lithosphere, which is a time-dependent process. The lithosphere is subsiding under the weight of ice, and uplifted back to original state during unloading. During this process, the stresses in the elastic part of the lithosphere will change.
The estimation of stress regime in the lithosphere due to flexural stress during glaciations depends on the timing of crustal stress equilibrium. Two elastic models have been proposed for the estimation of stress changes by flexural stresses during glaciations. Stephansson [42] assumes that the lithosphere is in equilibrium before the ice load is introduced (Figure 4, left). The resulting stress state during ice loading then becomes compressive underneath the ice sheet, and extensional outside the margin of the ice sheet (forebulge area). During ice withdrawal, the horizontal stresses slowly diminish until the uplift ceases and a new equilibrium is reached. Stein et al. [43] present a model where the lithosphere is assumed to have reached equilibrium after some time of ice loading (Figure 4, right). Thus, during melting and uplift the area originally covered by ice experiences extensional stresses, and the area outside the former margin of the ice sheet becomes subject to compressive stresses. These stress changes decrease with depth and are inverted in the lower half of the lithosphere.
If we assume that the duration of ice loading and unloading is long enough for the lithosphere to reach equilibrium during both glaciation and deglaciation, then both models are applicable during one glacial cycle. The average duration of glaciations and interglacials during Quaternary has been 40,000 to 100,000 years. The last ice age ended around 10,000 years ago, and at present the crust is nearly in isostatic equilibrium. The assumption of crustal equilibrium during both glaciations and interglacials is therefore not unrealistic, as the relaxation time for the Fennoscandian uplift is less than 4000 years [34,44]. Geosciences 2019, 9,

Identification of Glacially Flexured Areas
To identify if a certain study area has been affected by glacially-induced flexure, we need information of the glaciation history. The lithospheric flexure is calculated based on the glacial isostasy ( Figure 5), which is derived from the ice extent and thickness, as well as the properties and thickness of the lithosphere and mantle.
Here, we will provide numerical models of flexure-related stress effects in the Hoop Fault Complex, located in the SW Barents Sea (Figure 6). The modelling will be limited to the stress effects related to the isostatic response during the deglaciation after LGM. For identifying the areas of SW Barents Sea that are affected by flexural stresses, we use the glacial isostasy of the LGM ice sheet calculated in Section 2 as input.
Based on the resulting glacial isostasy, we can identify the areas of maximum bending of the crust, and the hinge zones with minimum flexure. The derivation results of a function are illustrated in Figure 5. The first function illustrates the glacial isostasy on a vertical section through the lithosphere near the ice margin. The gradient (or tilt) can be derived from the first derivative of the glacial isostasy, the curvature or flexure of the lithosphere can be found by the second derivative.

Identification of Glacially Flexured Areas
To identify if a certain study area has been affected by glacially-induced flexure, we need information of the glaciation history. The lithospheric flexure is calculated based on the glacial isostasy ( Figure 5), which is derived from the ice extent and thickness, as well as the properties and thickness of the lithosphere and mantle.

Model 1 [42]
Model 2 [43]  The resulting stresses are large during loading and diminish during unloading. The resulting stress state is compressional beneath the ice sheet and extensional in the hinge zones. Right: Flexural stress model ( [43]) assuming crustal stress equilibrium after a period of ice loading. From [28].

Identification of Glacially Flexured Areas
To identify if a certain study area has been affected by glacially-induced flexure, we need information of the glaciation history. The lithospheric flexure is calculated based on the glacial isostasy ( Figure 5), which is derived from the ice extent and thickness, as well as the properties and thickness of the lithosphere and mantle.
Here, we will provide numerical models of flexure-related stress effects in the Hoop Fault Complex, located in the SW Barents Sea (Figure 6). The modelling will be limited to the stress effects related to the isostatic response during the deglaciation after LGM. For identifying the areas of SW Barents Sea that are affected by flexural stresses, we use the glacial isostasy of the LGM ice sheet calculated in Section 2 as input.
Based on the resulting glacial isostasy, we can identify the areas of maximum bending of the crust, and the hinge zones with minimum flexure. The derivation results of a function are illustrated in Figure 5. The first function illustrates the glacial isostasy on a vertical section through the lithosphere near the ice margin. The gradient (or tilt) can be derived from the first derivative of the glacial isostasy, the curvature or flexure of the lithosphere can be found by the second derivative.  Here, we will provide numerical models of flexure-related stress effects in the Hoop Fault Complex, located in the SW Barents Sea ( Figure 6). The modelling will be limited to the stress effects related to the isostatic response during the deglaciation after LGM. For identifying the areas of SW Barents Sea that are affected by flexural stresses, we use the glacial isostasy of the LGM ice sheet calculated in Section 2 as input.
Based on the resulting glacial isostasy, we can identify the areas of maximum bending of the crust, and the hinge zones with minimum flexure. The derivation results of a function are illustrated in Figure 5. The first function illustrates the glacial isostasy on a vertical section through the lithosphere near the ice margin. The gradient (or tilt) can be derived from the first derivative of the glacial isostasy, the curvature or flexure of the lithosphere can be found by the second derivative.
The second derivative of the glacial isostasy of the LGM ice for the SW Barents Sea is shown in Figure 6. According to this, it is clear that the SW Barents Sea was greatly affected by lithospheric flexure during the deglaciation after LGM. There is a pronounced hinge zone (white area) extending from west of Bjørnøya over the Stappen High and into Loppa High. From here it extends eastwards over the Samson Dome onto the Finnmark Platform. Areas in the hinge zone are less likely to be affected by the flexural stress after the LGM deglaciation.  Figure 2). Yellow, orange, and green areas have high curvature and are likely to be influenced by stress effects during and after the glaciation. (cf. Figure 5). The black square marks the area included in the stress modeling, and the blue arrows show the orientation of maximum flexure-related extension after LGM deglaciation within that area.
The second derivative of the glacial isostasy of the LGM ice for the SW Barents Sea is shown in Figure 6. According to this, it is clear that the SW Barents Sea was greatly affected by lithospheric flexure during the deglaciation after LGM. There is a pronounced hinge zone (white area) extending from west of Bjørnøya over the Stappen High and into Loppa High. From here it extends eastwards over the Samson Dome onto the Finnmark Platform. Areas in the hinge zone are less likely to be affected by the flexural stress after the LGM deglaciation.

Model Setup
The stress modelling is done in 2D using the state of plane stress, and carried out using the commercial software Comsol Multiphysics 5.2 and the add-on structural mechanics module. Comsol Multiphysics is a FEA (Finite Element Analysis) software, where the problem domain is divided by a mesh of triangular elements in which the calculations are done. The set-up of a stress model in Comsol Multiphysics includes defining model geometry, then adding material properties, loads and boundary conditions. The model geometry includes faults geometries based on a published continental shelf map of Norway [45] with scale 1:2 mill, interpreted at the base of Top Jurassic. The model includes an area of approximately 300 × 300 km in the SW Barents Sea that was covered by ice during LGM. All fault geometries are manually traced in a CAD tool, and then imported into Comsol Multiphysics.
Every geometry unit of the model is assigned material properties, that is, Young´s modulus and Poisson´s ratio for the lithologies included. In the models presented here, the host rock is modelled as soft shale with a Young´s modulus of 5 GPa and a Poisson´s ratio of 0.25, which is within the typical values for shales [46]. Fault zones normally consist of numerous fractures and smaller faults which lower the general stiffness, or Young´s modulus, of the host rock [47]. The faults are therefore given a Young´s modulus of 0.1 GPa and Poisson´s ratio of 0.25.  Figure 2). Yellow, orange, and green areas have high curvature and are likely to be influenced by stress effects during and after the glaciation. (cf. Figure 5). The black square marks the area included in the stress modeling, and the blue arrows show the orientation of maximum flexure-related extension after LGM deglaciation within that area.

Model Setup
The stress modelling is done in 2D using the state of plane stress, and carried out using the commercial software Comsol Multiphysics 5.2 and the add-on structural mechanics module. Comsol Multiphysics is a FEA (Finite Element Analysis) software, where the problem domain is divided by a mesh of triangular elements in which the calculations are done. The set-up of a stress model in Comsol Multiphysics includes defining model geometry, then adding material properties, loads and boundary conditions. The model geometry includes faults geometries based on a published continental shelf map of Norway [45] with scale 1:2 mill, interpreted at the base of Top Jurassic. The model includes an area of approximately 300 × 300 km in the SW Barents Sea that was covered by ice during LGM. All fault geometries are manually traced in a CAD tool, and then imported into Comsol Multiphysics.
Every geometry unit of the model is assigned material properties, that is, Young's modulus and Poisson's ratio for the lithologies included. In the models presented here, the host rock is modelled as soft shale with a Young's modulus of 5 GPa and a Poisson's ratio of 0.25, which is within the typical values for shales [46]. Fault zones normally consist of numerous fractures and smaller faults which lower the general stiffness, or Young's modulus, of the host rock [47]. The faults are therefore given a Young's modulus of 0.1 GPa and Poisson's ratio of 0.25.
Next, boundary conditions are assigned to the model. As the model represents the base of a geological unit at depth, the sides of the model cannot deform freely. Two opposite boundaries of the model are therefore fixed. The other two boundaries are assigned a load representing the flexure-related extension originating from the isostatic response during the deglaciation after LGM. The orientation of the flexure-related extension is perpendicular to the curvature trends, with an orientation~SW-NE. (Figure 6). In a study of the post-glacial lithospheric flexure in the North Sea, analytical and numerical calculation of glacial unloading (assuming isostatic equilibrium before the onset of unloading) of a 1 km thick ice sheet and a lithosphere with elastic thickness of 30 km resulted in bending stresses in the order of 20 MPa [30]. We have used 5 MPa, which is a more conservative magnitude for this effect. The model does not take into consideration possible tectonic background stress at the time of glacial unloading. This modelling may therefore be considered as a study of the effect of stress changes due to the glacially induced flexure.

Results
In the resulting models (Figure 7), the stress is calculated based on the assumption of a linear elastic material. Rocks normally behave elastic up to 1%-2% strain at low temperature and pressure [46,48,49]. The elastic behavior of rocks is limited by their strength, that is, the maximum applied compressive or tensile stress the rock can withstand before failure occurs. The tensile strength of most rocks is in the range 0.5 to 6 MPa, most commonly 2-3 MPa [46,50]. This means that crustal tensile stress exceeding these values is not very common, as it is likely to result in tensile failure of the rock. The shear strength of rocks is twice their tensile strength, as follows from the Griffith failure criterion [46,51]. Thus, in the models, the tensile stress is truncated at 6 MPa, and the shear stress at 12 MPa. Areas concentrating stress above these values are likely to fail in tension or shear. However, due to uncertainties related to the material properties for the involved sedimentary units and the flexure-related load, the stress magnitudes in the models below should be read in terms of relative values. The stress is likely to be released by failure in the high magnitude areas (marked in yellow and red), and less likely to lead to failure in low magnitude areas (marked in blue).
The first model (Figure 7, top) shows the resulting maximum principal tensile stress (σ 3 ) when the area is subject to extension during LGM deglaciation. Tensile stress can be described as stretching stress, and can lead to opening up of new extension fractures (when the tensile strength of the rock is exceeded) or opening up of preexisting fractures. The SW-NE extension causes large areas between the faults to concentrate high magnitude tensile stress, whereas other parts of the host rock are located in a "stress shadow". Here, the host rock surrounding SW-NE trending faults accumulate high magnitude tensile stresses, whereas the fault zones themselves show minimal effect. This indicates that due to the LGM deglaciation, the host rock surrounding SW-NE trending faults is likely to experience increased local fracture-related permeability, which may lead potential hydrocarbons towards and along the faults. Focusing on the exploration wells, the result suggests that wells 7324/2-1 and 7324/8-2 are placed within areas of high fracture-related permeability during the LGM deglaciation. Potential hydrocarbon may have migrated along faults out of a potential reservoir during the deglaciation. According to the results, the other wells are located in areas of low tensile stress, and therefore less likely to have been affected by flexural stresses during and after LGM.
The von Mises shear stress results (Figure 7, bottom) show similar tendency. Most of the wells are located in areas of low magnitude shear stress, whereas wells 7324/2-1 and 7324/8-2 are located in areas of higher shear stress magnitudes. Well data for exploration wells drilled in the area show that all wells included here are oil or gas discoveries, except for well 7324/2-1 and 7324/8-2 which are dry with shows, indicating former presence of hydrocarbons.

Glaciation Model
Our glaciation model is not intended to include time-varying ice sheet geometries. However, LGM was the ice sheet configuration that was latest influencing the SW Barents Sea. The SW Barents Sea was glaciated around 24,000 BP and starting to melt away at 18,000 BP [36]. The LGM ice sheet is thus a reasonable ice sheet configuration for modelling the flexural stress effects in SW Barents Sea.

Earth Model
In the isostatic calculations we have assumed an elastic lithosphere thickness of 30 km (flexural rigidity 5 × 10 23 Nm). As mentioned above this Earth model has been shown to match other Fennoscandian post-glacial data, but most of the published GIA Earth models of Fennoscandia have a thick elastic lithosphere, between 80 and 160 km [34,39]. However, a thick lithosphere is not a viable option for data on tilted palaeo shorelines of coastal southwestern Norway and not for the Barents Sea [34,52]. However, a recent study of the Tapes transgression along the Norwegian coast indicates a thicker effective elastic lithosphere thickness of 70 km (flexural rigidity 3 × 10 24 Nm) in Northern Norway [53].
The glacial isostatic tilt caused by a 70 km thick elastic lithosphere is shown in Figure 8. Sedimentary basins located in SW Barents Sea could be tilted by up to 2.0 m/km due to the LGM ice sheet, which is about 70% of the tilt with lithospheric thickness of 30 km (cf. Figure 2). The flexural stress due to a 70 km thick lithosphere is thus somewhat lower than for a thinner lithosphere, but the stresses will be located basically in the same area. We believe, however, that a model with thin elastic lithosphere of 30 km better describes the Earth rheology of the Barents Sea than a model with thicker elastic lithosphere (see [52,54]).

Glaciation Model
Our glaciation model is not intended to include time-varying ice sheet geometries. However, LGM was the ice sheet configuration that was latest influencing the SW Barents Sea. The SW Barents Sea was glaciated around 24,000 BP and starting to melt away at 18,000 BP [35]. The LGM ice sheet is thus a reasonable ice sheet configuration for modelling the flexural stress effects in SW Barents Sea.

Earth Model
In the isostatic calculations we have assumed an elastic lithosphere thickness of 30 km (flexural rigidity 5 × 10 23 Nm). As mentioned above this Earth model has been shown to match other Fennoscandian post-glacial data, but most of the published GIA Earth models of Fennoscandia have a thick elastic lithosphere, between 80 and 160 km [33,38]. However, a thick lithosphere is not a viable option for data on tilted palaeo shorelines of coastal southwestern Norway and not for the Barents Sea [33,52]. However, a recent study of the Tapes transgression along the Norwegian coast indicates a thicker effective elastic lithosphere thickness of 70 km (flexural rigidity 3 × 10 24 Nm) in Northern Norway [53].
The glacial isostatic tilt caused by a 70 km thick elastic lithosphere is shown in Figure 8. Sedimentary basins located in SW Barents Sea could be tilted by up to 2.0 m/km due to the LGM ice sheet, which is about 70% of the tilt with lithospheric thickness of 30 km (cf. Figure 2). The flexural stress due to a 70 km thick lithosphere is thus somewhat lower than for a thinner lithosphere, but the stresses will be located basically in the same area. We believe, however, that a model with thin elastic lithosphere of 30 km better describes the Earth rheology of the Barents Sea than a model with thicker elastic lithosphere (see [52,54]).

Stress Models
The way the stress distributes in a material is depending on existing fractures and weaknesses and their orientation to the added load. In the stress model examples we have shown here, a

Stress Models
The way the stress distributes in a material is depending on existing fractures and weaknesses and their orientation to the added load. In the stress model examples we have shown here, a published regional interpreted fault map was used as basis. A more detailed fault interpretation will be needed for a more local study of the stress effects of glacially induced flexure.

Implications for Basin Modeling
Repeated ice loading during the Quaternary could have significant impact on sedimentary basins and add a degree of difficulty to petroleum exploration. The isostatic effects of ice sheets, (in addition to possible glacial erosion/deposition) could give significant tilting of sedimentary basins in the peripheral areas of former glaciated area. The tilting of the reservoirs at the Norwegian Continental Shelf could be up to almost 3 m/km and could alter hydrocarbon migration pathways. The ice sheets will act as seesaw during glacial/interglacial period, while the isostatic response of erosion and sedimentation will grow iteratively during the glaciations, depending on the stability of the depocenters.
Differential vertical movement of the lithosphere related to glacial isostasy leads to repeated tilting of sedimentary formations and potential petroleum reservoirs therein, which may have greatly affected hydrocarbon migration pathways in the former glaciated area [3]. In addition, the upward and downward bending of the lithosphere lead to flexural stresses likely to affect faults and their permeability, which could add to changes in hydrocarbon migration pathways. The flexural stresses are adding to background stresses in the area.
The map of the flexural effect based on glacial isostasy shows that the SW Barents Sea is affected by flexural stresses due to the last ice age. Stress models suggest that this effect may have caused high magnitude stresses to accumulate and therefore lead to an increased fracture-related permeability and escape of hydrocarbons in certain areas.
Most of the production licenses on the Norwegian Continental Shelf are located inside the former Last Glacial Maximum ice sheet (Figure 1). The results of our study demonstrate the importance of including the effect of glaciations in basin modeling studies of areas of petroleum potential. Adding to this are the temperature effects of glaciations and the effects of glacial erosion on temperature and stress discussed in more detail in [3]. Acknowledgments: The authors would like to thank four anonymous reviewers for constructive and fruitful feedback, which improved the paper.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Glacial Isostatic Model
The ultimate isostatic compensation achieved by any harmonic component is given by [44,55]: where F(k) = transformed ice load, ρ = density of the upper mantle, and g = gravity.