Time series modelling: applications for groundwater control in urban tunnelling

Water ingress to tunnels may result in pore pressure drawdown and consolidation settlements in areas above tunnels founded on soft soil deposits, potentially causing damage to buildings and infrastructure. To limit pore pressure drawdown, requirements are set on water ingress to bedrock tunnels. To meet these requirements, pre-excavation grouting is often performed to reduce the hydraulic conductivity of the rock mass surrounding the tunnel. Real-time pore pressure monitoring may be used to document pore pressure drawdown during construction. However, the effect of tunnel water ingress can be difficult to distinguish from natural pore pressure fluctuations. This paper presents a tunnel case in Oslo, Norway, where time series modelling was applied to local pore pressure data using the transfer function model framework. The input to the models was daily meteorological data considering precipitation and evapotranspiration, and the output was simulated pore pressure levels with impulse response functions. The models were optimised with data from before tunnel excavation, and simulations were run during the tunnel excavation period. Simulated pore pressure levels were compared with observed pore pressure levels to assess tunnelling-induced drawdown. Model uncertainty ranges were used to produce upper, lower, and best estimates of the drawdown. The findings show that time series modelling with transfer function models may be used in tunnel projects to continuously assess the impact on the local groundwater environment, for better evaluation of the pre-grouting performance, and for quantifying both the temporary and long-term drawdown with increased accuracy.


Introduction
Urbanisation calls for increasing development of underground infrastructure services, such as roads, railroads, subway, sewage systems, and electric power.In urban areas founded on soft soil deposits, water ingress to underground construction projects can cause pore pressure drawdown and consolidation settlements.This challenge is present in many regions around of the world, including Scandinavia, Canada, USA, Mexico, and Southeast Asia, where risk of settlements needs to be thoroughly addressed to mitigate damage to existing buildings and infrastructure.Extensive settlements in connection to groundwater lowering during underground construction have been documented by many researchers (Olofsson 1994;Persson 1998;Karlsrud 2002;Yoo et al. 2012;Lopez-Fernandez et al. 2012;Garshol et al. 2012;Shen et al. 2014;Myrabø and Moss-Iversen 2014;Barton and Quadros 2019).
To mitigate settlements, pre-excavation grouting is performed to limit water ingress to hard rock tunnels.The process involves pumping a cement-water mixture with high pressure through drilled boreholes in front of the tunnel face to reduce the hydraulic conductivity of the rock mass (Gustafson and Stille 1996).If the pre-grouting is insufficient the water ingress will cause pore pressure drawdown in nearby aquifers.The pre-grouting efficiency is normally evaluated using various types of water ingress measurements performed in the tunnel.However, due to the uncertainties related to these measurements, new methods need to be explored to be able to assess the performance of the pre-grouting (Langford et al. 2022).This paper explores the potential of using real-time pore pressure data to assess the effects of tunnelling on pressure levels during construction.Monitored pore pressure drawdown in soil deposits overlying the tunnel could ideally be used to adjust the pre-grouting program as the tunnel progresses, to meet requirements on water ingress limits and to avoid settlements.However, pore pressure levels also fluctuate naturally, influenced by meteorological factors such as precipitation and evaporation.This makes it challenging to separate natural and un-natural fluctuations in pore pressure, and to assess the effects of tunnel excavation.
Time series analysis with transfer function models has been applied to simulate fluctuation of groundwater levels for different purposes (Hipel and McLeod 1994;Von Asmuth et al. 2002, 2008;Bakker et al. 2008;Von Asmuth 2012;Shapoori et al. 2015;Obergfell et al. 2019 among others).The use of transfer function models to simulate groundwater levels influenced by underground construction is less explored, limited to the recent work by Boström and Lindblom (2021).This paper further presents the use of transfer function models to simulate pore pressure levels during a recent tunnel project in Oslo, Norway.The aim is to demonstrate how the model simulations may be used in monitoring to assess the effect of the tunnel on the nearby groundwater levels during construction.

Project area and data
In Oslo, Norway, large parts of the city centre are built on man-made fill overlying deposits of soft, marine clay (Bjerrum 1967).The clay deposits in the project area have a thickness of maximum 15 m.The underlying bedrock consists largely of Cambro-Silurian shale and limestone cut by Permian igneous dikes and sills.It is common to find a layer of moraine between the bedrock surface and the overlying clay.The soil deposits can be grouped into an upper aquifer in the topsoil (manmade fill material), and a lower aquifer in the moraine layer just above bedrock.The two aquifers are separated by the low permeable clay layer, acting as an aquitard.The natural pore pressure fluctuation is governed by factors such as amount of precipitation, evaporation, amount of surface run-off, snow melt, ground freezing, local topography, the presence of drainage systems, type of surface cover and soil type, etc.The lower aquifer is confined, and pore pressure levels are sensitive to discharge and recharge.Under these conditions, water ingress to tunnels and caverns in bedrock may lower the pore pressure in the confined aquifer.If the pore pressure reduction is sustained, a consolidation process will propagate in the clay, leading to ground settlements and potentially cause damage to buildings and structures (Karlsrud 2002).The magnitude of the settlements depends on the amount of pore pressure drawdown, the thickness of the clay layer, and the consolidation parameters of the clay deposit.
The project area, located in downtown Oslo, is schematically illustrated in Fig. 1.A bedrock tunnel was constructed, and pre-excavation grouting was performed continuously around the tunnel face ahead of excavation.Two piezometers, PZ1 and PZ2, were installed to monitor the pore pressures in the confined aquifer at the bedrock surface.Infiltration Fig. 1 Schematic illustrational showing bedrock tunnelling with pre-grouting in typical Oslo ground conditions.Piezometers are installed in the lower, confined aquifer at bedrock level, and artificial infiltration is conducted in bedrock wells wells IW1 and IW2 were installed into bedrock in settlementsensitive areas to counteract water ingress to the tunnel.IW1 is installed close to the deep depression where PZ1 is located, and IW2 is installed in the vicinity of PZ2.
Tunnel excavation was started on January 1, 2020 and the whole project was finished around May 1, 2022.Pore pressures were monitored in PZ1 and PZ2 from the end of March 2017 until May 1, 2022 (Table 1).The piezometers were connected to remote sensing units, transferring pore pressure data in real time to a local server.The pore pressure data had hourly logging, but the data was resampled to daily average values.Piezometer PZ2 has some missing data in 2020, caused by instrument failure and replacement.The observed (measured) pore pressure time series are displayed in Fig. 2a, b; installation data is given in Table 1.
The pore pressure levels were assumed to be influenced by precipitation and evaporation.Meteorological data was retrieved from the closest meteorological station in Oslo, located about 3-4 km from the piezometers.The recorded precipitation is displayed in Fig. 2c.A snow storage model was implemented to correct for the effect of snowfall during the winter months.The model corrects for snowfall by storing measured precipitation during days with negative temperature, releasing the stored precipitation when the temperature eventually rises above 0°.The corrected precipitation (net precipitation) equals measured precipitation on days with positive temperatures plus release of stored precipitation (snow) as meltwater.Daily evaporation was estimated using the Penman-Monteith equation for reference evapotranspiration from grass-covered surfaces (Penman 1948).The calculations include daily air temperature, wind speed, relative humidity, air pressure, solar radiation, site elevation, site latitude, and day of the year; results are shown in Fig. 2d.
The tunnel construction was split into two different parts, sites A and B, where water ingress was measured separately at the two sites (Fig. 2e).Underground construction at site A was started in January 2020, while construction at site B was started in May 2021.For both sites, tunnel water ingress was measured at specific times, typically during holidays and on Sundays, as the ingress measurements require a period of no activity on site (no use of water on site).A continuous time series of water ingress from the two sites was obtained by interpolation between the individual measurements.There are uncertainties related to the measured amount of water ingress at both sites due to difficulties in performing measurements consistently throughout the construction period.The measurements taken at site A are considered more uncertain than measurements taken at site B.
Water infiltration flux (litre/minute) was measured continuously in infiltration wells IW1 and IW2 (Fig. 2f, g).Infiltration well IW1 was installed in the beginning of 2020, while well IW2 was installed in 2021.Infiltration tests were conducted in 2020 and 2021, seen as short spikes in the data series.The infiltration flux was logged every hour, but the data was resampled to daily average values.
The pore pressure time series of 2020-2022 indicate natural seasonal pore pressure fluctuation of about 1 m in PZ1 and PZ2 prior to construction.It is also seen that the pore pressure in PZ1 responded well to the infiltration test in IW1 during March 2020.Data from PZ2 was not available at this time.A gradual decline in pore pressure was observed in PZ1 during spring 2020, when there was also an increase in water ingress rate from about 0 to 30 L/ min at site A. The pore pressure decline coincided with a drought period in Oslo, and the cause of the pore pressure decline was thus uncertain.The water ingress rate at site A remained constant at 30 L/min during the rest of the project period.In spring of 2021, infiltration tests were conducted in wells IW1 and IW2.A clear response was seen in pore pressure levels of both PZ1 and PZ2.In August 2021, large amounts of water ingress were recorded at the site B tunnel.To counteract the drainage to the tunnel site, water infiltration was increased in both wells to compensate for pore pressure drawdown.During the rest of 2021 fluctuating pore pressure was seen in both PZ1 and PZ2 as water ingress and water infiltration proceeded simultaneously.In 2022, water ingress to A and B was reduced as a result of grouting efforts.Water infiltration was ended during the first 2 months of 2022.In May 2022, pore pressure levels at both PZ1 and PZ2 were about 1-1.5 m lover than levels recorded before start of tunnelling.

The transfer function model framework
Time series modelling and simulation of groundwater levels may be done using different methodologies, as described by Hipel and McLeod (1994).In this study, time series modelling of the pore pressure levels of the confined aquifer was done using a transfer-function modelling framework.The transfer-function model was first described by Box and Jenkins (1970), being a type of model that describes "the dynamic relationship between two time series."In the case of hydrogeology, the time-dependent relationship between observed pore pressure and external stress series is used to make simulations.
A transfer function model may be described in the following way (Von Asmuth et al. 2002): For continuous time series, the contribution of stress p on the pore pressure at time t, h p (t) , is found by calculating the integral of the product of the stress time series p(t) and the impulse response function for that specific stress, p from 0 to time t (Von Asmuth et al. 2002): When using natural aquifer recharge as input stress, recharge models are used to convert precipitation and evaporation into a single stress series.A linear recharge model was used in this study due to its simplicity and ease of interpretation.The linear recharge model assumes that the relationship between precipitation, evaporation, and aquifer recharge is linear.In this configuration, the daily (natural) aquifer recharge rch is calculated as precipitation p minus evaporation e, where the evaporation is multiplied by an evaporation factor f (Von Asmuth et al. 2008;Collenteur et al. 2019): The impulse response function p describes the timeresponse of pore pressure to an external stress.For the purpose of this study, the time-response of pore pressure from natural recharge was simulated with an exponential function.The exponential function is suitable for systems where response to an impulse is near instantaneous.The shape of the exponential function is determined by two parameters, A and a: Autocorrelation of the model residuals is a common challenge related to time series modelling and can make it problematic to assess uncertainty and make statistical inference (Bakker and Schaars 2019).A noise model may be used on the residual time series to overcome this challenge, converting the model residuals into a series which more closely resembles white noise.Adding a noise model ( 1) was abandoned in this study due to inconsistency in model calibration and parameter optimisation.The effect may be caused by the high frequency of the data (daily values), combined with a relatively short time span of the data, but the actual cause of the inconsistency is still unclear.
The time series models were fitted to observed pore pressure time series using a nonlinear least squares solver from SciPy (Virtanen et al. 2020).The least squares solver optimises the model parameters to minimise the cost function.The cost function, in our case, is the sum of squared residuals from the calibration dataset.The open-source python package Pastas was used as the modelling framework.The background and basis of the Pastas package are described in Collenteur et al. (2019).

Results
A time series model was generated with natural recharge as a single input stress series.The model was calibrated on PZ1 and PZ2 data in the period before tunnelling (2017-2019).The optimised model parameters were used to make simulations of natural, undisturbed pore pressure levels during the following 2 years of tunnelling (2020-2022).

Simulating natural pore pressure levels
A time series model, from now on called the "basic" model, was calibrated to observed data from PZ1 and PZ2 from piezometer installation to the start of tunnel construction (2017-2019).The modelled and observed pore pressures for PZ1 and PZ2 are shown in Fig. 3.
The coefficient of determination (R 2 ) and root mean square error (RMSE) of the model fit are 0.79 and 0.11 for PZ1, and 0.79 and 0.09 for PZ2 (Table 2).The results are reasonably good; however, the models overestimate and underestimate the pressure levels during the most prevalent peaks and throughs.The lowest throughs might be caused by freezing of the upper soil layers during the coldest winter months.The underestimation during peak periods may be caused by large spatial variations during heavy rainfall and/ or snow melting events, being a cause for potential large discrepancy between measured precipitation at the meteorological station compared to actual rainfall at the site of the piezometers.
The block and step response curves of the recharge stress for the two models are shown in Fig. 4. The block response curve describes the time-response of the pore pressure from an infiltration impulse of 1 (the impulse response function).The step response curve describes the time-response of the pore pressure from a continuous recharge input of 1.We can see that the block responses for the two piezometers are very similar, both in amplitude ◂ and shape.This indicates that the pore pressure response and decline after recharge events are similar in the two locations.The step response is, however, different for the two piezometers.The impact of continuous recharge is higher for piezometer PZ1 than for PZ2, indicating that the pore pressure of PZ1 might be more sensitive to changes in long-term recharge rates than at PZ2.The basic models contain four optimised parameters, the A and a parameter of the exponential impulse response function, the evaporation factor f, and the base level d.The parameter optimal values and standard errors from the least-squared optimisation are shown in Tables 3 and 4. The standard errors of the optimised parameters are low (< 5%), indicating that the fitted values are reasonably consistent throughout the calibration period.The base level values (d) are in line with the observed values of PZ1 and PZ2 and considered reasonable.

Identifying tunnelling-induced drawdown
The optimised parameters of the calibration period were used together with daily measured meteorological parameters to  generate simulations of pore pressure during the tunnelling period of 2020-2022.The results are seen in Figs. 5 and 6, with a vertical black dotted line indicating start of the tunnelling period.
The simulated pore pressures of 2020-2022 are based on calibration on data before the tunnelling period.Assuming that the observed pore pressure behaves naturally, the simulated levels of 2020-2022 are estimations of natural pore pressure that would have occurred without any influence from the nearby tunnel project.Based on these assumptions, the difference between the observed pore pressure and the simulated pore pressure is a direct estimate of artificial pore pressure drawdown during the tunnelling period.
As seen in Fig. 5, the difference between observed and simulated pore pressure for PZ1 increases during spring-summer 2020, and then decrease somewhat during the last half of 2020.Heavy water ingress combined with increased water infiltration causes a somewhat chaotic picture for PZ1 and PZ2 in 2021 (Figs. 5 and 6).However, as water infiltration is shut off during winter 2022, the difference between the time series increases to maximum levels towards the end of the project period.The estimated tunnelinduced drawdown in May, 2022 amounts to approximately 1 m for PZ1 and 0.65 m for PZ2.However, this is just a snapshot from May, 2022, and the final, permanent drawdown may be calculated as the average difference between simulations and observations over the following years.

Visualising uncertainty
The model simulations might be run together with the observed time series for direct estimations of both temporary and permanent drawdown.However, the method does not describe uncertainty in the estimations or the model fit.Results from the calibration period show that the time series models sometimes over-and underestimate the pore pressure levels throughout the 2-year calibration period.The following methodology is proposed to visualise this uncertainty.
In Fig. 7, the basic model output was plotted against the pore pressure observations of the calibration period of 2017-2019 for PZ1 and PZ2.The scatter is an indication of model uncertainty during the calibration period.Upper and   lower bounds were set based on the assumption of linearity.
The upper and lower bounds were added to the model simulations of 2020-2022 (Fig. 8).The uncertainty ranges for the PZ1 and PZ2 model simulations is ± 0.28 m and ± 0.20 m, respectively.The plot may be used directly for increased control during pore pressure monitoring, as observations falling outside the grey area point to un-natural levels not previously seen in the calibration period (assuming a generalised model, not overfit to the calibration dataset).
Figure 9 shows a more detailed view of observed and simulated pore pressure of PZ1 in 2020.During spring-summer of 2020, the observed pore pressure fell after an initial test of infiltration well IW1.This also coincided with a period of drought in Oslo, and a decrease in pore pressure was thus expected.However, as seen in Fig. 9, the observed decline was larger than the simulated pore pressure decline, the difference increased through April.The observed pore pressure continued to fall, eventually reaching levels below the uncertainty range of the model simulations.This would have been a clear indication that the pore pressure decline was caused by increased tunnel water ingress.A pore pressure drawdown like that of March-May of 2020 is difficult to discover just by visually monitoring the pore pressure time series.This example demonstrates how model simulations may be used by engineers on site for early detection of drawdown.Early identification of drawdown is useful, as the pre-grouting design may be adjusted to match the current geological conditions.Hypothetically, the total pore pressure drawdown at the end of the project may have been reduced if appropriate adjustments to the pre-grouting program had been done in time during March-May 2020.
Comparing pore pressure observations and model simulations may also be used to estimate the impact on the nearby groundwater environment.Upper, lower, and best estimates of tunnelling-induced drawdown may be calculated by taking the difference between observed pore pressure and the upper, lower, and best estimates of the simulated pore pressure (uncertainty range, Fig. 9).These metrics may be used for daily monitoring purposes during the construction phase to identify drawdown.Alarm thresholds could be set on these metrics to signal the need for intervention and pregrouting adjustment.If the pore pressure monitoring and model simulations are extended beyond the end of the project period, upper, lower, and best estimates of permanent tunnelling drawdown may be calculated following the same methodology.

Discussion
As demonstrated, time series modelling with transfer function models may be used as a tool in urban tunnel projects to separate natural and un-natural pore pressure fluctuation, enabling increased accuracy when quantifying temporary and long-term drawdown.
The transfer function modelling framework is favourable in this use case due to its simplicity and physical basis.The framework is not a black box, as the mathematical model is based on known physical relationships.Additionally, the limited number of parameters compared to other model types (i.e.neural networks) makes it less susceptible to overfitting.The length and quality of the input data impact the overall quality and reliability of the results.In a typical tunnelling project in Norway, the piezometers are installed 1-3 years prior to start of the tunnel construction period.The amount of data is thus often limited.Overfitting is likely if the length of the calibration period is short.The risk of overfitting must be evaluated in each case.In this study, the results from the fitted models (number of model parameters, low variance in optimised values) indicate well-generalised models.Longer data series are, however, still preferred.
When applying time series modelling on data from real projects, the main "drivers", i.e. the main input stresses must be included in the model.If not, the model will not be able to successfully simulate the system, resulting in a bad model fit.A bad fit might also be caused by a complex relationship between pore pressure levels and the external stresses, not being suitable for implementation into the proposed model setup.In some areas, the dynamics between measured precipitation, evaporation, and aquifer recharge is not well described with a linear recharge model.In this case, more complex nonlinear recharge models might be needed to adequately simulate aquifer recharge (Peterson and Western 2014).The time series models may be expanded to include data such as water ingress measurements and water infiltration rates.Simulations of pore pressure response from different combinations of water ingress and water infiltration may be done after an initial calibration, as demonstrated by Boström and Lindblom (2021).Additionally, the pore pressure response to climatic changes may be explored by adjusting meteorological input levels, for instance increasing or decreasing the amount of rainfall or temperatures over a given period.Such simulations come with its own uncertainties; however, the possibilities should be explored further.
Time series modelling with transfer function models may be used for other purposes than demonstrated here, including estimation of annual aquifer recharge rates (Obergfell 2019;Collenteur et al. 2021), extension of piezometer data series by model simulation (forward or backwards), filling missing data, to name a few.Time series modelling might also be useful in combination with well infiltration and pumping tests.It might be possible to use observed data during well testing to make simulations of pore pressure drawdown in an aquifer during tunnel design.From outside the world of tunnelling, using transfer function models may also prove useful in the field of geohazards, using high frequency forecasts of pore pressure levels in early warning systems of landslide hazard, or in the field of engineering for simulations of groundwater levels as input to geotechnical design.These concepts should be explored further.

Fig. 2
Fig. 2 Data used in the study.a Observed pore pressures PZ1, b observed pore pressures PZ2, c daily net precipitation, d reference evapotranspiration, e measured tunnel water ingress at site A and B, f water infiltration flux in well IW1, and g Water infiltration flux in well IW2

Fig. 4
Fig. 4 Block and step response curves from recharge for the piezometers PZ1 and PZ2

Fig. 5
Fig. 5 Top: PZ1 basic model output.Observations shown with a black line.Model simulations shown with a red dotted line.Black vertical line marks start of the tunnelling period

Fig. 7
Fig. 7 Basic model output vs. pore pressure observations during the calibration period of 2017-2019

Fig. 8
Fig. 8 Basic model simulations during the tunnel period of 2020-2022 with uncertainty range in grey

Table 1
Installation data for piezometers PZ1 and PZ2 is the number of external stresses, h i (t) is the simulated pore pressure at time t based on input stress i, d is the base pore pressure level, and r(t) is the model residual at time t.The model residual is the difference between observed and simulated values (the model error).

Table 3
Optimised model parameters with standard errors for PZ1 basic model