Comprehensive risk assessment of groundwater drawdown induced subsidence

We present a method for risk assessment of groundwater drawdown induced land subsidence when planning for sub-surface infrastructure. Since groundwater drawdown and related subsidence can occur at large distances from the points of inflow, the large spatial extent often implies heterogeneous geological conditions that cannot be described in complete detail. This calls for estimation of uncertainties in all components of the cause-effect chain with probabilistic methods. In this study, we couple four probabilistic methods into a comprehensive model for economic risk quantification: a geostatistical soil-stratification model, an inverse calibrated groundwater model, an elasto-plastic subsidence model, and a model describing the resulting damages and costs on individual buildings and constructions. Groundwater head measurements, hydraulic tests, statistical analyses of stratification and soil properties and an inventory of buildings are inputs to the models. In the coupled method, different design alternatives for risk reduction measures are evaluated. Integration of probabilities and damage costs result in an economic risk estimate for each alternative. Compared with the risk for a reference alternative, the best prior alternative is identified as the alternative with the highest expected net benefit. The results include spatial probabilistic risk estimates for each alternative where areas with significant risk are distinguished from low-risk areas. The efficiency and usefulness of this modelling approach as a tool for communication to stakeholders, decision support for prioritization of risk reducing measures, and identification of the need for further investigations and monitoring are demonstrated with a case study of a planned railway tunnel in Varberg, Sweden.


Introduction
Sub-surface projects are generally constructed in materials formed and impacted by complex geological and anthropogenic processes (Lundman 2011), creating highly uncertain and variable ground conditions.These conditions can result in a wide range of risks including groundwater drawdown induced subsidence initiated by groundwater leakage into the sub-surface constructions.Since groundwater drawdown can affect large areas (km 2 ), see e.g.Burbey (2002), Huang et al. (2012), the number of affected buildings and constructions can be substantial.There are many examples where groundwater drawdown induced land subsidence have led to severe consequences, e.g.Shanghai (Xue et al. 2005), Mexico City (Ortega-Guerrero et al. 1999), Bangkok (Phien-wej et al. 2006), Las Vegas (Burbey 2002), Los Angeles (Bryan et al. 2018) and the Scandinavian cities Stockholm, Gothenburg and Oslo (Karlsrud 1999;Olofsson 1994).
The cause-effect chain from groundwater leakage, subsidence and resulting damage and cost is an interaction between several processes.The chain (Fig. 1) starts with groundwater extraction via leakage of groundwater into a sub-surface construction in bedrock (1a) or soil (1b).The Electronic supplementary material The online version of this article (https://doi.org/10.1007/s00477-018-01647-x)contains supplementary material, which is available to authorized users.
& Jonas Sundell jonas.sundell@trafikverket.se;jonas.sundell@chalmers.seprocess proceeds with reduction of groundwater heads (2), reduction of pore pressure in clay or other compressible deposit (3) and ground subsidence (4).The response and sensitivity of the constructions affected by the subsidence determines the extent of damages (5).As a final point, the cost associated with the damage determines the economic consequences (6) (Sundell 2018;Sundell et al. 2017).
When planning future sub-surface constructions, there are different alternatives for risk reducing measures such as improved sealing to avoid leakage and infiltration of water to maintain groundwater levels.The decision problem is to identify the best of these options given a decision criterion, which in this paper is the highest net benefit (NB) by means of cost-benefit analysis (CBA).To evaluate the highest NB, the risk for subsidence damage as a function of probability of occurrence and damage costs need to be estimated for each alternative.Since the cause-effect chain is complex and processes occur in highly heterogeneous systems, accounting for uncertainties in all stages within the cause-effect chain are necessary for risk assessment.
As the information varies in character between each part of the chain, different methods for uncertainty estimation need to be combined for a holistic assessment of risk.In urban areas, large numbers of borehole logs from previous construction projects describing soil stratigraphy are often easily obtainable (de Rienzo et al. 2008;Marache et al. 2009;Velasco et al. 2013).This geologic information is necessary for an understanding of groundwater flow and the subsidence process in different materials.With high spatial density of the boreholes, the geostatistical method Kriging (Matheron 1963) can be used for estimation of spatial uncertainty of the stratigraphy, see e.g.Asa et al. (2012), Bourgine et al. (2006) or Thierry et al. (2009).Samples and tests describing hydrogeological properties are on the other hand often sparse, making the Kriging method unfeasible.Instead, inverse calibrated stochastic simulations (Burrows and Doherty 2015;Carrera et al. 2005;Doherty 2003;Li and Zhang 2018;Siade et al. 2017;Tonkin and Doherty 2009) can be used for an estimation of groundwater flow conditions.Uncertainties in spatially sparse compression properties can be evaluated by detrending the samples and estimate probability density functions (PDFs) (Sundell et al. 2017).Information on building characteristics such as construction and foundation type can (in Sweden) be found at the municipality's archive and their resilience can be classed from damage schemes, see e.g.Skempton and Macdonald (1956), Bjerrum (1963), Rankin (1988) and Son and Cording (2005).Finally, the associated damage cost can be valued from an uncertainty estimation of historical records of subsidence in respective damage category.
The aim of this study is to present a comprehensive risk assessment procedure that facilitates cost-benefit analysis of risk reducing measures and considers the entire chain of events, from the initiating groundwater drawdown to the resulting consequences.This procedure considers uncertainties in all parts of the chain and is adapted to the large scale of the problem and different information types.Despite a comprehensive approach, the uncertainty estimation is limited to parameter classification whereas model uncertainties are not considered, see e.g.Refsgaard et al.
Fig. 1 The cause-effect chain of groundwater drawdown induced subsidence damages.The pink area illustrates bedrock, green: coarse grained material, yellow: soft clay and grey coarsegrained filling material.The hatched line at 1a illustrates a fracture zone in the bedrock (Sundell 2018;Sundell et al. 2017) (2007).The procedure is exemplified with a case study of a planned railway tunnel in Varberg, Sweden.

Method
This paper summarizes a complex modelling procedure and workflow, presented in Fig. 2. In a first step soil stratification and bedrock levels are modelled from borehole data (Fig. 2a).The second model (Fig. 2b) is a finitedifference groundwater model based on the geometry given by the first step together with observations of e.g.groundwater heads, precipitation and evaluated hydraulic conductivities.The effect of groundwater drawdown and changes in water balance of different design alternatives is then modelled in each of these plausible model solutions resulting in different drawdown scenarios (Sect.2.2).In the third model (Fig. 2c), changes in pore pressure and effective stress are calculated as well as subsidence from the difference in groundwater head between the calibrated models and each of the corresponding drawdown solutions (Sects.2.3 and 2.4).Finally, the resulting subsidence simulation is combined with a cost function to estimate the economic risk for subsidence damages for the different design alternatives (A 0 , A 1 , A 2 in Fig. 2D), see Sect.2.5.1.The best prior alternative can then be identified as the alternative with the highest expected net benefit compared with a reference alternative, see Sect.2.5.2.

Bedrock and soil stratification model
The soil stratification and bedrock level model describes stratigraphy between boreholes and surface mapping.The modelling procedure has previously been presented in Sundell et al. (2016) and Sundell et al. (2017) but a summary with minor modifications of the method follows here.In the case-study, the soil stratification is simplified to three continuous layers: filling or coarse-grained material (topmost), clay (mainly silty or sandy) and coarse-grained material (glacial till and glacio-fluvial deposits) above crystalline bedrock (Fig. 3a).Since the boreholes contain different types of information (Fig. 3b), the model follows a stepwise procedure to consider all available information and dependencies between layers.
In the first step, a variogram is modelled from boreholes with bedrock level information (Fig. 3c).Kriging of bedrock levels results in a grid with average values and standard deviations.With this result, possible deviations in bedrock levels are estimated (red lines in Fig. 3d).In the same manner, standard deviations and mean values are calculated from boreholes containing both bedrock levels and lowest level without reaching bedrock (Fig. 3e).The latter dataset contains both boreholes with a deep lowest level close to the bedrock surface and boreholes with a shallow lowest level at larger vertical distance from the bedrock surface.In order to consider only deep and not shallow levels, the same quantile is always simulated for Fig. 3d and e (e.g. the median in D is compared with the median in E) and the lowest level gives the bedrock level in the simulation (Fig. 3f).
From boreholes with soil stratification information, the parameter z pa is defined as the proportion of the topmost layer of the total soil thickness.This proportion is then transformed to a normal distribution from the probability integral of the standardized normal distribution N(0,1).Similarly, z pb is defined as the proportion of clay and the two bottommost layers (clay/[clay ?coarse grained]) and transformed to normality analogous z pa .These transformations are necessary for independent simulations between soil layers and previous bedrock level simulation.From the transformations, variograms are modelled (Fig. 3g and j) and Kriging interpolations produce average values and standard deviations of respective parameter.Based on this result, possible deviations of respective material proportion are simulated (Fig. 3h and K).The thickness of the topmost layer is then given by multiplying the simulated proportion with the total soil thickness (Fig. 3i).The total soil thickness is calculated by subtracting the surface level with the modelled bedrock level.From the total soil thickness, the thickness of the topmost layer and z pb the thickness of the clay layer is calculated (Fig. 3l) together with the levels of respective layer (Fig. 3n). Figure 3m illustrates max and min values of clay thickness from combined simulation results.
In each simulation sequence the total vertical stress (r v ) is estimated based on the density of respective material (Fig. 3o).The resolution in the case study is chosen to 20 9 20 m, which is considered sufficient to cover soil heterogeneities and individual risk objects for the area of the case study.At each of these grid points a vector with 0.1 m vertical resolution is created for calculation of stresses in the soil, r v .

Groundwater model
Since many different model solutions can be acceptable and consistent with available observations, this part forms an ill-posed problem with more independent than dependent parameters (Beven 2006).A randomized inverse calibration results in several plausible model solutions with different water balance, calibrated heads and parameter fields.In this model, the numerical groundwater model MODFLOW (Harbaugh 2005) with the NWT solver (Niswonger et al. 2011) is combined with the PEST (Parameter ESTimation code) sub space Monte Carlo (SSMC) (Tonkin and Doherty 2009) technique.The overall process follows principles and steps suggested by e.g.Reilly (2001), Freeze et al. (1990) and LeGrand and Rose ´n (2000), including: definition of project goal, collection of data, definitions of conceptual and numerical models, model parameterization, calibration with available observations, simulation model of the problem at hand and finally, a design suggestion based on the model results.From the previous soil and bedrock stratification model, the numerical model is discretized with different materials and layers in a three-dimensional grid, using mean layer thicknesses.The hydraulic transmissivity (T = K * d), determines flow in an individual layer.T can be adjusted either with respect to the layer thickness (d) or with the hydraulic conductivity (K), see e.g.Fetter (2001).With the choice of only including the most likely stratification Except for K of the different materials, the model is parameterized with recharge (RCH).Within the different materials, significant heterogeneity can be present.To model this heterogeneity, fields of material properties are assigned from a 2D scatter point set (pilot-points or PP) representing different positions within each material (Doherty 2003).From the expected variability of each material property, prior parameter estimates are assigned to the PP.If a specific parameter is considered reliable at a certain location (for example an evaluated K from a reliable pump test), a fixed value can be assigned to this PP.
The steady state groundwater model is calibrated against available observations of e.g.groundwater heads with PEST SSMC resulting in several plausible solutions, which is documented in a wide range of literature (Burrows and Doherty 2015;Doherty 2011;Fienen et al. 2013;Hou et al. 2015;Rossi et al. 2014;Woodward et al. 2016).After the SSMC procedure, solutions that are not able to fulfill a calibration criterion despite the re-estimation are removed.These include deviations in water balance greater than 10% between in-and outflow, together with differences between modelled and observed heads greater than 1 m in any observation well.In the remaining n number of solutions, the effect of changed drainage conditions from different design alternatives are modelled in each solution (see right part of Fig. 4).The design alternatives (A i , with index i = 1…m, A 0 is the reference alternative and m is the number of alternatives) include the sub-surface construction and safety measures to reduce the effect of the leakage.

Change in pore water pressure and effective stress
The subsidence model is initiated by calculating change in pore pressure and effective stress between each of the calibrated groundwater model solutions (simulation 1-n in the left part of Fig. 4) and the corresponding simulation with a design alternative (right part of Fig. 4).Although it would be possible to model changes in pore-pressure in the clay layer directly in the MODFLOW model e.g. with the SUB package (Galloway and Burbey 2011), this model does not take into account variability in clay thickness from the probabilistic stratification model.After a soil stratification profile and r v has been simulated (see Sect. 2.1) groundwater heads in layer 1 and 3 (above and below the clay layer) are simulated from the accepted calibrations with a uniform distribution, U(1,n).From these, the pore pressure profile (u) is calculated as a linear approximation between the pressure heads at the bottom and the top of the clay layer (Fig. 4A), see e.g.Persson (2007) and Zeitoun and Wakshal (2013).If layer 1 is dry and the pressure head in layer 3 is below the upper edge of the clay, u is calculated hydrostatically with the head in layer 3. The effective stress (r 0 v0 ) (which governs the deformation of saturated granular medium containing water within its voids) is calculated from r vu (Fig. 4B).The change in u (Fig. 4C) and r 0 v (Fig. 4D) in a simulation with a design alternative is calculated in the same way as in the calibrated solution.As with r v , u and r 0 v are simulated at each of the vertical grid points and at every 0.1-m interval.Note that in the example in Fig. 4 there is no change in head in the topmost layer but only in layer 3.
The above calculation of u and r 0 v considers the steady state conditions at the end of the consolidation process.Since the groundwater model is steady-state, it is conservatively assumed that the drawdown in the coarse-grained material in layer 1 and 3 is instantaneous.The consolidation process in the clay in layer 2 is calculated from Terzaghi's one dimensional consolidation theory (Terzaghi 1943) with a solution based on Fourier series as described by Taylor (1948).In this solution, a consolidation coefficient (c v ) and a time factor (T v ) calculates a consolidation coefficient (U z ) for every 0.1 m increment.

Compression parameters and subsidence model
To calculate subsidence, we use a simple one-dimensional elasto-plastic model by Larsson and Sa ¨llfors (1986) (similar to the tangent modulus approach by Janbu (1967)).Soil compressibility is evaluated from diagrams where vertical strain (e v ) is plotted against effective vertical stress (r 0 v ).The diagrams are based on constant rate of strain (CRS) oedometer tested clay samples.The CRS-test is described in different standards (ASTM International 2012; European Committee for Standardization 1997b; Swedish Standard Institute 1991; The European Union Per Regulation 305/ 2011 2004) but is internationally less common than incremental loading oedometer tests.The benefits of CRStests include generation of a continuous stress-strain curve and a shorter test period.Disadvantages include inability to evaluate creep (secondary or delayed consolidation with no change in stress state) (Larsson and Sa ¨llfors 1986) and dependence between the measured response and the applied strain rate (Pu and Fox 2016).Although the method is demonstrated with CRS-test results evaluated with a modulus approach, the principle is applicable on results omitted from incremental loading tests with the more recognized compression index method, see e.g.(Fang 2013).The method has previously been applied for probabilistic simulations of subsidence in Sundell et al. (2017) but a summary of the method as applied in this study, its parameters, statistical analysis (data transformation, regression analysis, ANOVA, T test) as well as details on the simulations are presented in Appendix A.
From the previous described simulations of changes in r 0 v and PDFs describing the relationship between different compression parameters (r 0 c r 0 L , M L , M 0 and M', see Appendix A), subsidence is simulated at each vertical vector and 0.1 m interval on the horizontal grid.The entire simulation sequence is repeated for all accepted groundwater model runs in each design alternative using 1000 random draws from the bedrock and soil stratigraphic model.These computations result in a distribution of subsidence magnitudes where the combined results of all grid points are used for spatial mapping of subsidence risks in terms of magnitude and probability.

Risk expressed as a function of damage and cost
The extent of a damage depends on the subsidence magnitude together with the response of the risk object.In this study, risk objects are limited to buildings although other constructions could be damaged (e.g.roads and pipes).The construction of the buildings and their foundation, together with any historical damages, determine the building's response to subsidence.In most damage schemes, the occurrence of cracks is correlated with movements such as absolute and differential subsidence, angular distortion and lateral strain, see e.g.Skempton and Macdonald (1956), Bjerrum (1963), Rankin (1988), Cooper (2008) and Son and Cording (2005).Minor cracks are associated with esthetical damages with small restoration costs whereas larger cracks can affect the function or the structure of a building with substantial damage and cost.The transition Fig. 4 Calculation of pore pressure and effective stress from calibrated groundwater models and models with drawdowns between these damage extents can be described as a continuous function where the economic risk of subsidence, R s , is given by a combination of the economic cost, C s , and the probability of a damage, f s : R s is first calculated for each building.The sum of R s for all buildings gives the total economic risk for each design alternative.To be in accordance with recommendations in damage schemes, the continuous function is simplified to a stepwise function, Fig. 5.In this simplification three broad categories: (1) esthetical, (2) function and (3) stability are applied according to Driscoll (1995).

Subsidence and damage
Detailed calculations of subsidence movements and building response, see e.g.Giardina et al. (2015) and Schuster et al. (2009) are beyond the scope of this paper since the presented method aims to give an overview of the risk at an early stage.More advanced models of building response could be useful at a later stage, when following up high-risk objects identified by the approach presented in this study.In the applied subsidence model, only vertical movements (total subsidence) are considered.If the subsidence is evenly distributed over a building's floor area, large total subsidence does not necessarily mean a large damage.Nevertheless, a large total subsidence does often result in a movement that causes damages.The lower limit of a harmful absolute subsidence is defined to 10 mm in Rankin (1988) and Son and Cording (2005), 20 mm in Chiocchio et al. (1997), 50 mm in Eurocode (European Committee for Standardization 1997a), 20 mm for the City-tunnel project in Malmo ¨, Sweden (case M 487-04, Land and Environmental Court of Appeal 2004).For damages in the more severe categories, Chiocchio et al. (1997) defines 30 mm as an upper limit for restauration of esthetical damages and 100 mm as a lower limit for rehabilitation measures related to functional damages.Rankin (1988) and Son and Cording (2005) defines the limit for possible structural damage to 50 mm.The lower limit for severe damages where evacuation is needed is set to 150 mm in Chiocchio et al. (1997).In Rankin (1988) the lower limit for expected structural damage is defined to 75 mm.From these numbers, lower limits for respective damage category are estimated: esthetical; 10 mm, functional; 30 mm, structural; 75 mm.From these definitions and the simulation of subsidence, f s is calculated.

Damage costs
Damage costs can both be direct (costs for repairing the damage) and indirect (e.g.delay costs of the project, loss of goodwill, reduced market value of the damaged building and inconveniences for tenants).For the case-study, only direct costs are considered.Since insurance policies in general do not cover subsidence in Sweden, no rigorous database exists.Instead, the valuation is based on a few recent legal cases in Sweden.Not all cases are related to groundwater-induced subsidence but the damages and compensations are relevant for the valuation, see Table 1.
As additional information to Table 1, the gross charge for new production was 33,500 SEK/m 2 for single family houses and 45,500 SEK/m 2 for apartment buildings in Sweden year 2016 (Statistiska Centralbyra ˚n 2017).From these numbers, the average cost per m 2 gross floor area is estimated to 400, 14,000 and 38,000 for each respective damage category 1-3.Since the cases are few, no differentiation between foundation-and construction types is possible.Another consequence of the situation with few cases is that PDFs cannot be defined directly.Instead, a highest reasonable damage cost for respective category is estimated.This highest value is assumed to represent the 95th percentile in respective category.For class 1, 2 and 3 this value is estimated to 1000, 30,000 and 60,000 SEK respectively.The percentiles were defined based on the information given in Table 2 and expert judgment in the project team.From these numbers, a log-normal PDF describes C s for each damage class, LN(l,r): 1; LN(5.99, 0.557), 2; LN(9.55, 0.463), 3; LN(10.55,0.277).A log- normal PDF was chosen to assure non-negative costs and to represent larger uncertainties in the right tail (Hansen and Singleton 1983).
The definitions of f s and C s , calculates R s for all risk objects in each design alternative.With one alternative defined as the null alternative, the benefit of an alternative is given by: Each alternative has an investment costs, c i (the cost relative A 0 ).The net benefit (NB) of an alternative is given by: To evaluate a situation where there is a stream of costs and benefits over time and where the benefits and costs do not appear simultaneously, Eq. 3 needs to be adjusted to include time-discounting for calculation of U i as a Net Present Value (NPV), se e.g.Boardman et al. (2011), So ¨derqvist et al. (2015).Discounting is necessary in situations of e.g.longer consolidation times where subsidence damages (and costs) propagate over time.

Varberg case study
The method is applied for a part of a planned railway shaft and tunnel in Varberg, Sweden, see Fig. 6.The bedrock geology is dominated by gneissic granite and charnockite, both of magmatic origin (1700 and 1 400 million years).The gneissic granite has become gneissic, veined and folded due to several periods of deformation and metamorphism.During the Svecofennian orogeny (1 400 million years), the charnockite was subjected to granulitefacies (Lundqvist and Kero 2008).Where the two bedrock types intersect (location of pump test in Fig. 6), several fracture zones are observed (Sundlo ¨f et al. 2016).In general, hydraulic conductivity decreases with depth below surface since weathering effects declines on deeper levels together with reduced fracture apertures since rock stress increases at greater depths (Gustafson 2012).As a result, the bedrock is assumed to be more conductive at its more weathered topmost part and along fracture zones.
As mentioned in the method, the soil stratigraphy is simplified to three continuous layers.The bottom layer is dominated by a sandy glacial till and sandy glaciofluvial deposits.Above this layer, glacio-marine clay is present in the low altitude areas close to the shoreline.Laminates of silt and rests of shells have been observed in some clay samples (see supplementary material SM 3.2).Postglacial coarse grained gravely beach deposits formed by redeposition of glacial deposits dominates the topmost layer (Pa ˚sse 1990).In the west parts close to the shore, filling material is deposited on top of the gravely deposits.
Figure 6 shows observations of bedrock outcrops, boreholes with bedrock levels and clay thickness, groundwater observation wells, pump tests (Martinsson et al. 2016) and piston samples of clay (Hurtig et al. 2016;Trafikverket 2016).The borehole and groundwater head data is collected from investigations made for the planned tunnel, an inventory of the municipal archive in October 2017 and well data, openly available from the Swedish  is modelled from 44 boreholes where the entire clay layer is penetrated, 19 boreholes where it cannot be confirmed if the entire clay layer has been penetrated and 12 boreholes where it can be confirmed that clay is not present for the entire soil profile.The variograms (from a larger dataset within Varberg than the model area, see SM1.1) for z pa are represented by an exponential function whereas z pb is represented by a spherical function.The modelled variograms are fitted to the experimental variograms with least square to minimize the quadratic sum of the difference between the these, see Cressie (1985).Both variograms are modelled with a nugget effect because of sample uncertainties.Opposite to the bedrock levels, no anisotropy is revealed for the transformed proportions.From the variograms, borehole data and observations of bedrock outcrops, grids with expected values and standard deviations are modelled using kriging.From these grids, uncertainties in soil stratigraphy is modelled according to the procedure presented in Sect.2.1.The result of the model for clay thickness is presented in Fig. 7.As can be seen in the figure, clay thickness is greater in the low-altitude western parts close to the shoreline.At locations with bedrock outcrops, clay is nonexistent in all simulations.

Groundwater model
Groundwater flows from higher areas in the east to lower areas in the west.Along the corridor of the planned railway, the groundwater heads are about 1 MASL.Low groundwater gradients are present between this corridor and the sea.Groundwater recharge is expected to the higher altitude areas where clay is less likely to be present (nr1 in Fig. 8).From this area, most of the groundwater flows in the coarse-grained materials and in the fractured topmost part of the bedrock (nr2) and further to the sea (nr4).Less groundwater flows in the fracture zones (nr3) before it reaches the sea (nr4).The model area partly covers the sea which level cannot be affected by the planned groundwater drawdown.
From the bedrock and soil stratification model and the location of fracture zones a geometrical model for groundwater modelling is defined, see SM2.2.Since continuous soil layers are modelled also where these are not present, layers thinner than 0.3 m are removed.The bedrock is divided into three materials; fracture zones, a 2 m thick fractured topmost bedrock and less fractured bedrock on deeper levels.The grid for the numerical model consists of 7 layers with a general 30 9 30 m cell resolution, which is refined at the location of the railway.

Probabilistic model calibration
Within each material, PP are defined every 100th meter with prior parameter estimates presented in Table 2.A large interval for the filling material is motivated by the variability of fine grained and coarse-grained fractions at different locations.The hydrogeological investigations (Sundlo ¨f et al. 2016) for the coarse grained material and the bedrock demonstrate a high variability of evaluated K values even at short distances with the exception of the marked pump test in Fig. 6.At this location, K is evaluated to 2 9 10 -5 m/s and a fixed PP with this value is defined for the coarse-grained material (see SM2.2).Because of the presence of clay and discharge of groundwater in the western parts, a RCH area is defined for a part of the total model area (shaded area in SM2.2 nr1).The net precipitation in the area is about 350 mm/year (1.1 9 10 -8 m/s).The large range for RCH (Table 2) is motivated by significant spatial variation of recharge because of surface runoff, variability in infiltration capacity and leakage from water and sewage pipes.Along the shoreline, the sea level is represented by a constant head boundary (CHD) in layers 1 and 2. Groundwater head observations from 38 observation wells in different layers (see Fig. 6 and SM2.1) are used to calibrate a steady state groundwater model.Additional groundwater wells are present in the area but each of these are closely located to an observation well and have observations in the same magnitude.A steady-state groundwater model is motivated by the short and not always overlapping time periods of the available observations.Nevertheless, a few wells have longer observation series (a few years) where the head variation is about 1 m.To represent this variation, the calibration target is defined to 1 m between observed and modelled heads.In addition to this target, the deviation between out-and inflow should not be greater than 10%.
Following the description in Sect.2.2 a groundwater model is first manually calibrated, then recalibrated with PEST following a SSMC calibration with 1000 runs.In total 731 of these models fulfilled the calibration criteria.The water balance for these models varies within the RCH area between 24-68 mm/year with a median of 40 mm/ year, see SM2.5.Max-and min values of calibrated K and RCH fields (interpolated from PP) are presented in SM2.3 and SM2.4.These posterior parameter fields reveal a significant heterogeneity in the different materials and wide parameter ranges for some locations and parameters, which demonstrates that several parameter combinations can calibrate the model equally well.A clear effect of the SSMC calibration is observed in layer 3 where the variability in K is lower close to the observation wells compared with locations at larger distances from the wells.

Modelling of design alternatives
From the calibrated models, the effect of three different design alternatives (A 0 , A 1 , A 2 ) for sealing of the tunnel and shaft is evaluated: A0.Without sealing of fracture zones and soil layers in shaft.A1.Sealing of shaft and fracture zones to K = 1 9 10 -7 m/s.A2.Sealing of and fracture zones to K = 1 9 10 -8 m/s and a draining layer below the sealed shaft.
The drain layer in A 2 is necessary to avoid damming of the groundwater flow.The railway is located between -4 and -15 MASL from the northern to the southern part, which defines the drain levels in all alternatives.The resulting drawdowns in all alternatives are presented in Fig. 9 with large drawdowns in A 0 and less drawdowns in A 1 and A 2 .The water balance in A 0 (SM2.6)results in drainage between 60 and 200 mm/year and an large leakage from the CHD boundary representing the sea.In A 1 , the inflow is reduced to 30-45 mm/year, which also reduces the leakage from the sea.In A 2 , the sealing reduces the inflow to 15-18 mm/year, which eliminates the leakage from the sea.Coarse grained material 4 9 10 -7 0.02 Fractured topmost bedrock 5 9 10 -8 5 9 10 -5 Less fractured bedrock 1 9 10 -10 1 9 10 -8 Fracture zones 5 9 10 -10 1 9 10 -4 RCH 9 9 10 -13 6 9 10 -8 Compression parameters from 13 locations (46 levels) have been defined from Hurtig et al. (2016) and Trafikverket (2016) and plotted against depth in Fig. 10 (see also SM3, not all samples contain all information).As can be seen in Fig. 10, all parameters in the top row reveal a weak trend with higher values against depth.

Data processing and statistical analysis of parameters
Following the procedure in Sect.2.3, u is calculated from groundwater observations together with r 0 v and OCR.At 6 m depth in 14T2046G, OCR is calculated to \ 1 as a result of disturbances during sampling, lab evaluation, an underestimated u and/or an overestimated r v and therefore excluded for further analysis.Dependencies and homoscedastic errors are accounted for through ln(OCR -1), ln(r 0 L /r 0 c -1), ln(M0/ML), 10log(k) and ln(Density), see Table 3.In addition to these, r 0 L and M L reveal a strong linear dependency (R 2 = 0.85), which introduces ln(M L /r 0 L ).M' is kept in its original form since it is independent of other parameters and can be described with a normal distribution (Sundell et al. 2017).Following these transformations, two parameters, ln(OCR -1) and ln(r 0 L /r 0 c -1) reveal a vertical trend along depth that cannot be neglected (R 2 [ 0.05) since it partly describes parameter uncertainties.The other parameters have no significant trend, which means that the regression result is ignored in the next step.
Three different group divisions are tested for significant differences on the 0.05 level: clay type (1), sample disturbance (2) and area (3).The first division is grouped into samples that only contain silty clay (siCl) and all other samples (mostly siCl with sand and shells), the second into three groups of sampling quality based on volume change at reconsolidation relative water content (Larsson et al. 2007;Lunne et al. 1997) and the third in four areas along the railway, see SM3.In the third division, area 3 is within the model area characterized by a built urban environment.Area 1 and 4 are unbuilt greenfield sites and area 2 along the present railway track.Since the first division contains only two groups, it has been tested with a T-test instead of ANOVA.Parameters with significant and insignificant Fig. 9 5, 50 and 95th percentiles of groundwater drawdown for design alternatives A 0 , A 1 , A 2 .Note that increased groundwater levels also can occur as a result of barrier effects of the planned tunnel differences are marked with ''Y'' (rejected H0) and ''N'' (H0 cannot be rejected) respectively in Table 3.Although significant differences are present for divisions between clay types, these distinctions cannot be observed from boreholes where clay is observed to build the soil stratification model (Sect.2.1).Therefore, the observed differences between soil types in the clay samples are not considered in the proceeding model.The observed differences between groups of sampling quality are highly correlated with the groups of soil types.Here, higher quality is related to siCl samples, while lower quality is exhibited for samples containing sand and shells.To consider uncertainties under the current state of knowledge where the sampling quality is dependent on the clay type (which cannot be distinguished in the stratification model), the observed differences are disregarded in the modelling stage.For the area division, significant differences are observed between most groups, most likely because of different load history between the areas.The post hoc test reveals that results from areas 1 and 2 can be excluded from area 3 (model area) for four parameters (Table 3).With these areas excluded, the remaining samples are tested for vertical trends (last column in Table 3), exhibiting a significant trend for ln(M 0 /M L ).

Definition of PDFs
From the result of the transformation and the statistical analysis, PDFs for simulation of subsidence are defined, Table 4.For the cases with a significant vertical trend, the PDFs describe the residual of the regression line.The regression line is described with y = a * x ?b where y is the parameter value, x depth below surface and a and b the regression coefficients presented in Table 4.
The resulting regression lines and confidence intervals of future observations are presented in Fig. 11 and the bottom row in Fig. 10.

Simulation of subsidence
From the soil stratification model, the groundwater model, the PDFs of the parameters in Table 4 and the simulation scheme in Appendix A, subsidence is simulated at each 20 9 20 m grid point (the result of the groundwater model is rescaled to this grid size), see Fig. 12. Simulations are carried out for 6 months and final subsidence (only final subsidence is presented in Fig. 12) for each of the three design alternatives/groundwater drawdown scenarios.As can be seen in Fig. 12, larger subsidence magnitudes are modelled at locations with thicker clay layers (Fig. 7) and larger groundwater drawdown magnitudes (Fig. 9).

Risk model
From the subsidence calculations and the PDFs for C s in 2.5.2,R s is calculated for each building and design alternative, see Fig. 13.In this calculation, buildings that are planned to be demolished or are founded directly on bedrock are excluded from the analysis (marked as not sensitive for subsidence in different figures).No existing damages on buildings have been observed within the area.As can be seen in Fig. 13, R s for individual buildings are correlated to locations with larger subsidence magnitudes.Since the cost functions are expressed per square meter (see Sect. 2.5.2),larger subsidence on larger buildings also results in larger R s .
In Fig. 14, the economic cost, C s , and the probability of a damage, f s for any building in each design alternative is expressed.As can be seen in Fig. 14, larger costs are more common in A 0 than in A 1 and A 2 .Since no or very small subsidence is simulated for most buildings, failure costs that approximate zero are the most common.
When R s is aggregated for all buildings in each design alternative, the total R s and B i after 6 months and final subsidence are calculated, Table 5.Relative the reference alternative A 0 , the benefit (B i ) is substantial in both A 1 and

123
A 2 .With a positive NB as a criterion to invest in an alternative, the large B i values result in generous margins for investment (c i ) in both A 1 and A 2 .Despite the large benefits, R s are still high in A 1 and A 2 , meaning that additional safety measures can be motivated depending on their investment costs.Such measures include additional sealing and infiltration to maintain groundwater levels.These measures should be prioritized to areas with high R s .
Fig. 11 Transformed parameters with regression lines for all samples (grey) and samples that cannot be excluded from area 3. The hatched lines show 90% prediction interval of samples that cannot be excluded from area 3 If the investment costs for A 1 and A 2 are the same, A 2 is identified as the best prior alternative since it has the highest net benefit.

Discussion
The presented method is comprehensive since all parts of the leakage-subsidence-damage-cost chain are considered.Despite this, simplifications have been made by only considering parameter uncertainties in chosen models.
One simplification is that loads from buildings (with or without basements) and considerations of different foundation types are omitted in the calculations.The effective stress can thus be expected to be greater at locations of buildings, which means that here subsidence can be  underestimated.On the other hand, if a building is founded on piles, the magnitude of subsidence can be overestimated.In addition, none of the soil samples are taken from below buildings, which could mean that the load history and compression parameters differ from the available samples.Despite these constraints, together with the fact that creep (secondary consolidation) and other movements than vertical are not considered, the order of magnitude of the subsidence calculations are expected to be reasonable for an initial risk assessment.Regardless of simplifications, the calculated risk costs (R s ) for individual buildings are reasonable when compared with the recent legal cases in Sect.2.5.2.Nevertheless, uncertainties in subsidence magnitude (see Fig. 12), building response and resulting damage cost are high; additional investigations with detailed models that consider different foundation types and building loads are recommended for buildings with high R s .Additional investigations can be prioritized with means of Value of Information Analysis (VOIA), where the costs for collecting new information are compared with the expected benefits of reduced risk of making an erroneous decision relative to a reference alternative, see e.g.Sundell et al. (2019) and Zetterlund et al. (2011).Based on the result of the VOIA, additional investigations, the best prior alternative or additional safety measures are prioritized.

Conclusions
This article presents a novel and comprehensive risk assessment procedure that facilitates cost-benefit analysis of risk reducing measures.The procedure is comprehensive since the entire chain of events, from the initiating groundwater drawdown to the resulting consequences (damage costs) are considered together with parameter uncertainties in all parts.The procedure facilitates a spatial CBA of risk expressed by a combination of the economic cost and the probability of damage for multiple design alternatives.The risk cost presented on maps, can be an important decision-support tool regarding planning of additional investigations, monitoring and safety measures.
The case-study results indicate high economic benefits, relative a reference alternative, for two design alternatives that include additional sealing.Despite high benefits, i.e. large reductions of economic risks, individual buildings with high residual risks are present in the alternatives.For these highlighted buildings, detailed investigations of building response are recommended.Based on the outcome of these investigations, implementation of additional safety measures should be considered.
Future research is recommended to evaluate not only parameter uncertainties but also uncertainties on different assumptions, conceptual and numerical models.Inclusion of refined models of building load and subsidence response are also suggested.In addition, the cost function should be refined by including more damage cost appraisals as well as different building types.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative commons.org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Detailed information on subsidence model
Different parameters are evaluated from the over-and normal consolidated part of the CRS plot.The over consolidated part is described by a recompression or swelling line (RCL) and the normal consolidated part at higher stress by a normal compression line (NCL) where the shape of these curves relates to the stress history of the soil.The preconsolidation pressure (r 0 c ) describes the maximum historical effective pressure.If the change in r 0 v results in stresses below r 0 c , the process is considered linear elastic.Instead of evaluating swelling-and compression indexes (as in the compression index method), different modulus of compressibility (M = dr 0 v /de v ) are evaluated.The evaluation of the modulus at the RCL part of the CRS-curve (M 0 ) often results in an underestimation (Fang 2013;Olsson 2010).As recommended by Sa ¨llfors (2001), the evaluated M 0 is multiplied with a factor 3. Despite this simplification, subsidence magnitudes at the RCL part are generally small.The modulus at the NCL part, the second constant modulus (M L ) is evaluated instead of a compression index as in the compression index method.In addition to the RCL and NCL states, the CRS test reveal a third state at higher stress levels above r 0 L where the simplification of constant modulus is no longer valid.Above r 0 L , the modulus number M' is evaluated as DM/Dr 0 (Larsson and Sa ¨llfors 1986).To evaluate time-dependent subsidence k (hydraulic conductivity of clay, necessary for calculation of c v ), is evaluated from the CRS-test at stress-levels equal to r 0 c .Depending on the state of r 0 v in relation with r 0 c and r 0 L different equations for calculation of subsidence are used (Eq I-III in Table 6).Eq IV integrates the result from the different 0.1 m segments of the vertical vectors.

Data processing and statistical analysis of parameters
The compression parameters r 0 c r 0 L , M L , M 0 and M' are evaluated at different locations and depths.To parameterize each vector of the 20x20 m horizontal grid from a relatively few samples, it is necessary to investigate how representative the samples are for the investigated area and if any trends and dependencies in the data is present.This investigation follows a procedure presented in Sundell (2016).The whole process is exemplified for r 0 c and r 0 v with data from the case study in Fig. 15.In a first step, dependencies between the parameters on the CRS-curve and r 0 v are investigated (Fig. 15a and b).The relationship r 0 c /r 0 vo defines the over consolidation ratio (OCR), Fig. 15c.Since r 0 c represents the highest historical effective stress, r 0 v0 cannot be higher than this point which conditions OCR [ 1.Because of the stress-strain relationship on the CRS-curve, the conditions r 0 L /r 0 c [ 1 and M 0 /M L [ 1 are introduced.In addition, dependencies between r 0 L and M L are investigated with regression analysis since these are related on the stress-strain curve.To avoid ratios \ 1 in subsequent simulations, OCR and r 0 L /r 0 c are subtracted by 1.In the next step, the subtracted values are ln-transformed to assure homoscedastic errors (equal variance at each vertical interval), Fig. 15d.Presence of vertical trends are investigated with linear regression (Fig. 15e), which results in normally distributed residuals (Fig. 15f).Normality is tested with residual-and normal score plots together with a Kolmogorov-Smirnov (KS) test, see e.g.Johnson et al. (2013).When the coefficient of determination (R 2 ) is close to zero, no vertical trends are assumed.In all other cases, R 2 partly describes parameter uncertainties.The conditioning of OCR means that values close to 1 are sensitive for small variations in both r 0 c (sample and lab disturbances) and r 0 0 (over/underestimated u and r).As a result, an outlier with OCR \ 1 is rejected (Fig. 15c).Because of the ln-transformations, these variations also result in a larger spread in the upper tail of the residuals (Fig. 15f).For the time dependent calculations, k is transformed to normality with the 10-logarithm.
Although sampling is spatially scarce, it is necessary to investigate if a sample is relevant for the investigated area.By dividing the samples in relevant groups based on e.g.location and evaluated sample disturbance, differences between these groups can be investigated with ANOVA, see e.g.Marx and Larsen (2006).Since ANVOA requires normally distributed data with equal variances, the result of the previous transformation is used.Equal means among groups defines the null hypothesis (H0), which is rejected on a 0.05 significance level.For tests with rejected H0, the Bonferroni method (Dunn 1961) is used as a post hoc test to compare differences between means.If the post hoc test reveals significant differences between groups (at 0.05 level), these groups are used to define probability density functions (PDFs) of parameters for subsequent simulations.If H0 cannot be rejected, difference between groups cannot be distinguished and all samples are used for definition of PDFs for the area of interest.

Simulation of subsidence
From the previous described simulations of changes in r 0 v and PDFs for the different compression parameters, subsidence is simulated at each vertical vector on the horizontal grid.This process is described in Fig. 16 with data from the Varberg case study.The simulation is initiated with a bedrock level (Fig. 16b) and upper and lower limits of the clay layer (Fig. 16c).From r 0 v before groundwater drawdown (Fig. 16d) and a PDF describing the residuals of ln(OCR -1), r 0 c is simulated along the soil profile (Fig. 16e).The shaded area describes the 90-percent prediction interval for observations of OCR.The red hatched line illustrates a simulation of OCR along the soil profile from the exemplified iteration (red dot) of the PDF for ln(OCR -1).The previous result together with a PDF describing the residuals of ln(r 0 L /r 0 c -1), simulates r 0 L (Fig. 16f).In the same manner, M L is simulated from a PDF describing ln(M L /r 0 L ) and the previous r 0 L (Fig. 16g).
Table 6 Equations for calculation of subsidence according to Larsson and Sa ¨llfors (1986) Case Equation Eq I r 0 v ?Dr \ r 0 In the case study, Du = Dr since it is only groundwater drawdown that causes a change in load.Eq IV integrates the total subsidence for each segment Stochastic Environmental Research and Risk Assessment 123 M 0 is simulated analogously from a PDF of ln(M 0 /M L ) and M L (Fig. 16g).M' (Fig. 16i) and 10log(k) (not illustrated in the figure) are simulated independent of the other parameters and of depth.
The simulated values of r 0 c and r 0 L are compared with the change in r 0 v as a result of groundwater drawdown.From this comparison, the adequate equation I-III in Table 6 is selected for calculation of subsidence both for final consolidation, Du(t max ) and after a certain time Du(t) at each 0.1-meter interval.In Fig. 16j, Eq I is used for the top part where r 0 v ?Du \ r 0 c (RCL-part), meanwhile Eq II is used for the bottom of the soil profile where r 0 c \ r 0 v ?Du \ r 0 L (NCL-part).The total subsidence in each grid point is approximated with the trapezoidal rule (Eq.IV).The whole simulation sequence is repeated for all accepted groundwater model runs with 1000 random draws from the bedrock and soil stratigraphic model.The simulation is repeated at each grid point for each groundwater drawdown scenario (originating from different design alternatives) (Fig. 16k).These computations results in a distribution of subsidence magnitudes (Fig. 16l).

Fig. 2
Fig. 2 A coupled model for probabilistic risk modelling of groundwater drawdown induced subsidence.The first part (a) models soil stratification and bedrock levels, the second part (b) groundwater drawdown, the third (c) subsidence and the final part (d) risk in terms of economic cost (C s ) and the probability of damage to a certain degree (f s ) for different design alternatives

Fig. 5
Fig. 5 Economic risk (R s ) expressed as a combination of economic cost (C s ) and the probability of damage to a certain degree (f s ) in a continuous (left) and a stepwise function (right)

Fig. 6
Fig.6Model area (red square), available data and location of profile in Fig.8

Fig. 10
Fig. 10 Parameter values with depth below ground surface with average values for all samples (grey line in bottom row).The average value of samples that cannot be excluded from area 3 is indicated with

Fig. 13
Fig. 13 Economic risk, R s , per building and design alternative A 0 , A 1 , A 2 for 6 months and final subsidence.10SEK & €1

Fig. 15
Fig. 15 Data processing to control dependencies, vertical trends and homoscedastic errors exemplified for r 0 v and r 0

Table 1
Classed cases with damage costs, 10 SEK approximates €1 Smith (1999)and soil stratification modelBedrock levels within the model area is modelled from 230 boreholes confirming bedrock level, 175 to coarse grained soil and 64 without confirmation of neither coarse grained soil nor bedrock.Experimental and modelled spherical variograms, see e.g.LloydSmith (1999), for bedrock levels in two horizontal directions is presented in SM1.1.The variogram demonstrates an anisotropy with the lowest values in direction -119 (NNE) and the highest values in -29 (WNW), which is in accordance with the main direction of the general tectonic structure pattern, indicated by the soil filled valleys in NNE.Similarly, clay thickness

Table 3
Transformed parameters and result of statistical analysis

Table 5
Economic risk, R s , per design alternative A 0 , A 1 , A