Forecasting the long-term deterioration of a cut slope in high-plasticity clay using a numerical model

This paper details development of a numerical modelling approach that has been employed to forecast the long-term performance of a cut slope formed in high plasticity clay. It links hydrological and mechanical behaviour in a coupled saturated and unsaturated model. This is used to investigate the influence of combined dissipation of excavation-generated excess pore water pressures and seasonal weather-driven near-surface cyclic pore water pressures. Deterioration of slope performance is defined in terms of both slope deformations (i.e. service) and factor of safety against shear failure (i.e. safety). Uniquely, the modelling approach has been validated using 16 years of measured pore water pressure data from multiple locations in a London Clay cut slope. Slope deterioration was shown to be a function of both construction-induced pore water pressure dissipation and seasonal weather-driven pore water pressure cycles. These lead to both transient and permanent changes in factor of safety due to effective stress variation and mobilisation of post-peak strength reduction over time, respectively, ultimately causing shallow first-time progressive failure. It is demonstrated that this long-term (90 year) deterioration in slope performance is governed by the hydrological processes in the weathered near surface soil zone that forms following slope excavation. demonstrate a in the initial critical failure is and by post-excavation Bishop‟s effective stress framework (e.g. strain softening and hydraulic conductivity), and captures observed slope mechanisms and failure conditions. The model can be used as a tool to forecast long-term behaviour that leads to deterioration of slope performance for example due to climate change and the efficacy of engineering interventions in slowing deterioration. However, this study also demonstrates the need to include the development of the near-surface weathered zone, post excavation, so that the governing soil properties evolve in this zone over time. Work by Stirling et al. (2017, 2020) has described and started to quantify the magnitude and rate of material degradation mechanisms that result in these modified soil properties.

cycles leading to strength deterioration and progressive failure means that highplasticity clay slopes can remain stable during one wet event but fail a number of years later due to a similar magnitude event (Take and Bolton, 2011;Postill et al., 2019). This study uses numerical modelling to show how deterioration of slope performance, as defined using both slope deformations and factor of safety against shear failure respectively, is driven by a combination of post-excavation pore water pressure dissipation and long-term cycling of weather driven near surface pore pressures that produce reductions in soil strength (Stirling et al., , 2020. Elia et al. (2017) illustrate the difficulty of modelling soil-vegetation-atmospheredriven long-term cyclic pore water pressures and highlight the limited data sets available to validate model behaviour. In addition, Elia et al. (2017) discuss the significance of including coupled hydro-mechanical behaviour in order to capture long-term soil strength deterioration due to shrink-swell stress cycles that ultimately lead to progressive failure. To date, numerical simulations of long-term clay slope behaviour have had one or more of the following characteristics; non-softening behaviour that can only analyse transient changes in effective stress (Oh and Lu, 2015;Tsiampousi et al., 2017); use has been made of simple wetting and drying boundary conditions to explore stress cycles; or where progressive failure is considered, this is undertaken within a saturated model (Kovacevic et al., 2001; 3 an infrastructure embankment in high-plasticity clay. Accumulated seasonal rainfall infiltration has also been shown to have an effect on slope stability in a number of modelling studies (Rahardjo et al., 2001;Alonso et al., 2003;Calvello et al., 2008;Zhang et al., 2014;Cotecchia et al., 2019). The importance of near surface behaviour, particularly for the modelling of pore pressures, is emphasised in Rahardjo et al. (2013). The methodology derived by Rouainia et al. (2009) was later used to investigate the stability of a natural slope in glacial till (Davies et al., 2014) and applied to a cut slope. This cut slope model also began to introduce deterioration into the approach adopted here by demonstrating the reduction in residual factor on the ultimate failure surface over time (Rouainia et al., 2020). The importance of the inclusion of the air phase and the effect on model results is further demonstrated by Cho (2016). Postill et al. (2019), also using a coupled two-phase flow and mechanical formulation, showed that the mechanism of ratcheting due to seasonal pore water pressure cycles, as observed in centrifuge experiments, could be successfully modelled numerically. That work demonstrated the effect of seasonal stress cycles on plastic strain accumulation and progressive failure due to prolonged wetting as well as continued seasonal stress and included a validation of both the pore water pressure cycles and the mechanical response.
A significant further advance is made in the current paper using a two-phase flow numerical procedure (after Postill et al., 2019) to study long-term deterioration of a highway cut slope in high plasticity clay near Newbury, UK. Importantly, it also responds to the review carried out by Elia et al. (2017) by validating the numerical model outputs using a 16-year time series of field measurements from the slope (Smethurst et al., 2012). The primary mechanism of deterioration studied is strain-

Numerical Modelling Approach
The numerical analyses reported in this paper used the finite difference software FLAC two-phase flow (FLAC-TP) (Itasca, 2016), which has been modified to include user-specified algorithms described in this section and in Postill et al. (2019). FLAC-TP is formulated for the simulation of wetting and non-wetting fluids flowing within a continuum model, allowing hydro-mechanical coupling. The approach adopted uses a displacement-pore pressure ( ) formulation and assumes the fluids (i.e. water and air) within the soil are homogeneous and idealised (i.e. fluids are immiscible). The software and approach were used and discussed by Postill et al. (2019) for the modelling of the seasonal ratcheting mechanism validated against centrifuge experiments and the London Clay local-strain softening was validated in Rouainia et al. (2020). In brief, inter-particle stresses are described by Bishop"s generalised effective stress, , as follows (Bishop, 1959): Equation (1. where , is the effective stress parameter. In this model, the degree of saturation, , is used as an approximation for . It has been recognised in the literature that can be used to quantify the Bishop's effective stress parameter at high degrees of saturation, such as those found in the typical range of pore water suctions predicted in this work (Khalili et al., 2004;Nuth and Laloui, 2008).
where and are fitting parameters.
To account for the reduction in hydraulic conductivity ( ) as the soil becomes unsaturated, relative hydraulic conductivity of the water ( ) and air ( ) phases are calculated as per Equations 3 and 4.
where and are the fluid and gas dynamic viscosities, respectively.
In this work, slope performance is defined in terms of factor of safety (FoS) against ultimate limit state failure and is classified as being unacceptable once failure has occurred (FoS 1). The effective stress changes due to dissipation of construction induced pore pressure and seasonal wetting and drying cycles would occur regardless of soil type. However, for progressive failure, which is the critical mechanism driving deterioration of the soil shear strength in high plasticity clays, a strain softening constitutive model is required. The adopted approach could also be applied to lower plasticity materials that undergo similar stress cycles, however the magnitude of strength reduction will be lower due to the increased residual strength in lower plasticity soils.
However, analyses using strain-softening models can be compromised by mesh dependency Ba ant and Pijaudier-Cabot, 1988;Vermeer and Brinkgreve, 1994;Di Prisco and Imposimato, 2003;Galavi and Schweiger, 2010). To minimise this in the numerical solution, a nonlocal regularisation technique was adopted, where nonlocal plastic strains ( ) are obtained by averaging local plastic strains ( ), calculated using the Mohr-Coulomb constitutive model (Summersgill et al., 2017b;Postill et al., 2019). in FLAC is derived as follows (Itasca, 2016): In which ( ) are the principle plastic strain increments.

8
Plastic strain averaging is calculated using Equation (6): Where, ∑ ′ = weighted volume. The weighting function ( ) adopted in this work is that proposed by Galavi and Schweiger (2010).
where is the distance between the calculation point and the neighbouring stress point and is the internal length which controls the distance over which neighbouring strains effect the strain calculation.
This weighting function was shown to be the least mesh-dependent of a number considered by Summersgill et al. (2017b) for analysis of clay cut slopes.

Material Properties
Newbury cutting is formed within London Clay; and strain-softening behaviour from previous modelling studies in this material has informed the parameters adopted in the modelling (Potts et al., 1997;Kovacevic et al., 2007Kovacevic et al., , 2013Ellis and O"Brien, 2007;Rouainia et al., 2009;Summersgill et al., 2018) as well as laboratory and field data (Bromhead, 1978;Bromhead and Dixon, 1986).
The elastic behaviour makes use of a mean effective stress dependent model, where Young"s modulus ( ) is derived from Bishop"s generalised mean effective stress ( ) where is the reference Young"s modulus at zero mean effective stress and is the minimum allowable Young"s modulus.
The stiffness model was calibrated against swelling tests undertaken on London Clay samples  taken at a depth similar to the cutting height modelled in this work (8.9 m). The swelling response from this and other published work on the modelling of excavations in London Clay are also included by way of comparison (Potts et al., 1997;Kovacevic et al., 2007). Figure 2 shows the ability of FLAC to replicate the swelling behaviour shown in Potts et al (1997), however when this is compared to the more recent data of Hight et al. (2007), it tends to overestimate the magnitude of swelling at vertical stresses between approximately 20 and 100 kPa. As such a minimum bulk and shear modulus were specified based on the values adopted by Jurečič et al. (2013). The resulting swelling profile more closely matches the field data. Figure 2: Calibration of the stiffness model against Oedometer swelling tests (Apted, 1977;Hight et al., 2007) and comparison with previously adopted stiffness models (Potts et al., 1997;Kovacevic et al., 2007;Rouainia et al., 2020).
Strength is modelled using a Mohr-Coulomb strength envelope that combines internal friction angle ( ) and apparent cohesion ( ). Shear strength/shear strain behaviour is incorporated via definition of peak, post peak and large displacement residual strength properties and the respective strains required for mobilisation.
It should be noted that in this work, the post peak strength adopted is equal to the average strength at failure from the back analysis of slope failures in the field (Bromhead, 1978;Bromhead and Dixon, 1986). This is often described as residual strength by other authors Potts et al., 1997;Kovacevic et al., 2007;Ellis and O"Brien, 2007;Summersgill et al., 2018). In this work the adopted residual strength is derived from ring shear testing of London Clay at large displacements (see for example Bishop et al., 1971;Lupini et al., 1981) and includes a total loss of cohesive strength and a further reduction in friction angle as used in Rouainia et al. (2020) where the local strain softening model was also validated for London Clay.
As described in section Numerical Modelling Approach3 the non-local strain softening model adopted here introduces an additional parameter (the internal length; ) that influences the softening rate once plastic strains occur. As such it was necessary to calibrate the model to ensure the adopted softening criterion were appropriate to the adopted mesh (as the non-local model will reduce, but not eliminate, mesh dependency). The calibration was based on the methodology softening model. It can be seen that in the local softening models, the mesh size has a large effect on the softening rate with the rate decreasing as mesh size increases.
The non-local approach reduces, but does not eliminate, the mesh dependency of the model results especially in coupled hydro-mechanical models. The non-local results for elements of 0.5 and 0.75 m bracket the accepted local softening rate with the 0.5 m model being more conservative. As such, 0.5 m elements were adopted in this work.   (Potts et al., 1997;Kovacevic et al., 2007Kovacevic et al., , 2013Ellis and O"Brien, 2007;after Rouainia et al., 2009after Rouainia et al., , 2020Jurečič et al., 2013;Summersgill et al., 2017aSummersgill et al., , 2018  The presence of vegetation cover on the slope influences soil strength in two ways. The first effect is hydrological, due to the suctions generated by the uptake of water by plant roots, which causes changes in the effective stress, and therefore, soil The second effect is due to the direct influence of roots on the mechanical strength of the soil. Whilst it is known root reinforcement has an influence on near-surface soil strength, there is not an accepted consensus in the literature on how this can best be implemented in numerical models (Cohen et al., 2009). Where parameters are suggested (for example as an increase in cohesion) the potential range of values is large and can therefore be the dominating factor that has significant impacts on slope response. Therefore, the mechanical component of soil shear strength enhancement due to vegetation has been omitted from the simulations. As such the modelling results can be considered conservative.
The hydrogeological properties used to represent Newbury cutting are shown in demonstrates the changes in SWRC due to weather-induced soil deterioration in both the laboratory and field. To reflect this, an intact and a "weathered" London Clay SWRC were adopted in the analysis, see Figure 5c. The SWRCs used for intact and "weathered" London Clay were adopted from Briggs et al. (2013) and Croney (1977).
The region over which the near-surface saturated hydraulic conductivity and J o u r n a l P r e -p r o o f Journal Pre-proof weathered soil water retention properties are applied in the slope model are shown schematically within Figure 6. In the numerical model, prior to excavation, the saturated hydraulic conductivity is set to the profile shown in Figure 5a with the reference datum being the original ground level. The adopted profile is informed by laboratory and field permeability data for the London Clay (Garga, 1970;Dixon and Bromhead, 1999;Dixon et al., 2019). Once excavation was complete, the increased saturated hydraulic conductivity profile is added to the near-surface of the model with a revised datum representing the excavated ground surface.
The model employed accounts directly, and continuously, for hydraulic conductivity change due to desaturation; but does not incorporate time-dependent changes due to the development of desiccation cracking (Anderson et al., 1982) or the process of progressive SWRC and saturated hydraulic conductivity change resulting from material deterioration mechanisms (i.e. weathering). There is currently little information in the literature on how soil parameters evolve in a developing weathered zone, and so deriving and validating a hydraulic deterioration model that could be employed in the analysis would not be possible at present. However, the effects of these processes were included by assuming they had occurred by the end of excavation (i.e. the newly excavated slope has a fully developed weathered nearsurface zone). While there are simplifications in the modelling regarding these complex processes, the work presented advances the current state of modelling long-term behaviour through the inclusion of unsaturated behaviour, deteriorated near-surface hydrogeological properties and validation against long time series monitored data.

J o u r n a l P r e -p r o o f
Journal Pre-proof

Initial Conditions and Cutting Excavation
In-situ stresses were initialised assuming a coefficient of earth pressure at rest of , and a phreatic surface 1 m below the pre-excavation model ground surface (see Figure 6a). This in-situ stress is equal to the values used for prior modelling of natural London Clay (Potts et al., 1997;Tsiampousi et al., 2017;Rouainia et al., 2020) and from available in-situ measurements . The adopted value best captures the swelling response measured by Hight et al. (2007). See The initial groundwater table is set equal to the UK winter average at 1 m below the surface (Vaughan, 1994) matching that adopted in previous numerical modelling work (Potts et al., 1997;Rouainia et al., 2020). These stresses were then used to derive the initial bulk and shear stiffness distributions ( 5000 kPa) as per Kovacevic et al. (2013). The cut slope excavation was modelled in eight 1 m increments with swelling allowed to occur between each layer for a period of 9 days (region excavated is shown in Figure 6a). The total excavation time was 72 days, the same excavation rate as modelled for a cut slope in London Clay by Potts et al. (1997). During the post excavation swelling phase and during weather driven cycling of pore water pressures, the stiffness distributions varied using the relationships presented in Table 1.

Boundary Conditions
The lateral boundaries were fixed to prevent horizontal movements, and the base was fixed to prevent both horizontal and vertical displacements. The lateral and horizontal boundaries were made impermeable. Surface boundary conditions were applied as daily discharges, obtained through a soil water balance calculation (Blight, 2003) utilising the weather data obtained during monitoring at Newbury cutting (Smethurst et al., 2012), with regional weather data used to infill missing data (i.e. prior to the start of monitoring in 2003). The water balance calculation is shown conceptually within Figure 6b and given by: where is the soil moisture deficit, is potential evapotranspiration, is rainfall and is runoff.
From observations, Smethurst et al. (2012) showed run-off occurred only when the soil was at, or very close, to field capacity (i.e. = 0 mm). Therefore, within the water balance calculation, run-off was assumed to be zero until the soil was at field capacity.
In order to calculate a daily flux boundary condition (m/s), the daily change in the soil water balance was calculated. This flux was applied to the numerical analysis and solved for a one-day time increment, the discharge value updated for the next day and the model solved for another day, and so on. Daily flux increments were calculated as shown in Equation (8).
When the model surface became fully saturated, the surface boundary condition was switched from a discharge boundary to a fixed pore water pressure of 0 kPa until a boundary condition removing water was applied (i.e. when the surface pore water pressure equalled 0 kPa, the boundary condition changed to a fixed pore water pressure of 0 kPa until a discharge condition that removed water was stipulated, after which discharge boundary conditions were reinstated on the boundary). This is a similar approach to that implemented by Smith (2003). This algorithm was applied to each surface node individually so that a mix of surface pore water pressure fixity J o u r n a l P r e -p r o o f 20 and/or discharge could occur at a given time increment depending on the local degree of saturation/available water storage.

Past Weather Data
For the period 2003 to 2019, data from the onsite weather station for Newbury cutting was used. This included rainfall, run-off and reference potential evapotranspiration ( ) calculated using the FAO Penman-Monteith equation (Allen et al., 1998). See Equation (11.
Where is the slope of the vapour pressure curve, is the net solar radiation, is the surface heat flux, is the psychrometric constant, is the air temperature, and are the saturation and actual vapour pressures and is the wind speed at 2 m height.
To infill gaps in the weather time series, regional weather data and the Blaney-Criddle (1950) method to derive reference potential evapotranspiration was used.
The sources of regional weather data are summarised in Table 3 Table 4.
Where, is a calibration factor.      Table 4: Site-specific calibration factors for Blaney-Criddle (1950) calculation (after Doorenbos and Pruitt, 1977) To convert reference potential evapotranspiration into actual evapotranspiration for the soil water balance calculation, the effect of vegetation (using a crop factor, = 1 for grass) and soil water retention behaviour were accounted for using the method outlined by Smethurst et al. (2006). This takes into account the readily available and total available water ( and ) in the soil as a function of the stress and wilting points of the vegetation in relation to the soil water retention properties. As suggested by Smethurst et al. (2006), = 144 mm and = 58 mm. Using this data, actual evapotranspiration ( ) can be calculated in line with Equation (13) (Clarke et al., 1998;Smethurst et al., 2006).

} (13)
A comparison of annual cumulative rainfall and potential evapotranspiration for regional and site-specific weather data are presented in Figure 9a and b, with the derived profile shown in Figure 9c. Figure 9 shows that regional and site specific weather is comparable and for this exercise is adequate to model the history The boundary conditions for the analysis were generated from the as described in Equation (9).

Future Weather Data
To forecast future performance of Newbury cutting and gain an understanding of long-term behaviour, a future weather sequence was required. Stochastic synthetic daily meteorological data was generated using past seasonal variation in weather patterns using the statistical down-scaling approach described by Wilby et al. (2014).
Baseline regional weather (see Table 3) for the period 1961 to 1990 (see Figure 10) was used to calibrate the synthetic weather generated. Baseline weather was used to obtain monthly data with cumulative probability distributions created to establish the range and probability of occurrence of temperature and precipitation on any given day of a month. The probability of wet days in each month was also obtained for the baseline period.
The synthetic weather data was then used with the Blaney-Criddle (1950) method and the steps described previously to obtain the and boundary conditions for the analysis. Figure 11 shows the for the past and future synthetic weather sequence, where the seasonal cycles of for the future synthetic weather are comparable to the past weather.

Results of Modelling Newbury Cutting
The following section presents the results of the numerical analysis of Newbury cutting using the methodology, material properties and boundary conditions described. Initially, the period covering the field monitored data is considered (2003 to 2019), then the forecast longer-term behaviour of Newbury cutting is presented.
The model geometry and mesh used for the analysis is shown in Figure 6 along with the locations of the field-monitored and modelled pore water pressures.

Hydrogeological Response
For the period in which pore water pressure records are available, the near-surface hydrogeological behaviour of the numerical analysis and monitored data were compared for an exemplar location with the results presented in Figure  In addition to comparing modelled and monitored pore water pressures across the slope, pore water pressures below the toe of the slope were considered specifically in order to demonstrate the capabilities of the model to capture post-excavation dissipation of excess pore water pressures. This is shown in Figure 13   However, pore water pressure dissipation towards steady state is still occurring at depth (i.e. 5.0 m and 10.0 m) 100 years post-excavation. This is broadly consistent with the times described in Chandler and Skempton (1974) for London Clay cut slopes. Figure 13 shows the inter-related nature of long-term behaviour and annual cycles of wetting and drying. Extreme wet and dry periods can be seen to both cause transient fluctuations in pore water pressures and affect the rate of long-term dissipation. However, the overall trend is a dissipation of post-excavation suppressed pore water pressures, which was captured by the numerical model.

Mechanical Response
The modelled mechanical behaviour as a response to the seasonal pore water pressure cycles and post-excavation dissipation of excess pore water pressures from construction to the end of 2019 is plotted in Figure 14. Figure 14a illustrates the modelled pore water pressure record at piezometer D at 1.0 m depth, with mid-slope ground surface displacements shown in Figure 14b and toe surface displacements in Figure 14c. For the time period plotted, seasonal shrink-swell displacements give a net upwards and outwards swelling of the slope at the toe and at mid-slope. This initial swelling can be associated with post-excavation pore water pressure dissipation and stress relief, the majority of which at these shallow depths occurs in the first 5 to 7 years (see Figure 13a and Figure 14a).  Figure 13a and Figure 14a). This is largely a function of the non-uniform effective stress change and generation of negative pore pressures during excavation which suppress pore pressures to a decreasing extent with increasing distance up the slope.

Forecast of Long-Term Behaviour
The This is a form of the strength reduction approach (see for example Griffiths and Lane, 1999  As such, the recorded FoS represent a series of discrete records of the slope stability at individual points in time. Full details of the approach are given in (Itasca, 2016) to 20 years corresponding to the time taken for excess pore water pressure dissipation at shallow depth (i.e. 2.5 m to 5.0 m) below the toe of the slope in the numerical analysis. At greater depth (10 m) this increases significantly (40 to 50 years) and is consistent with the findings of Vaughan & Walbancke (1973) for similar excavations in low hydraulic conductivity materials such as London Clay.
The 5  37 purely a result of effective stress changes due to dissipation of construction induced negative pore water pressures. Conversely, the additional reduction in FoS for the softening model must then be attributable to the additional reduction in strength caused by strain-softening. This assumption is used to compare the proportion of total deterioration that can be attributed to pore water pressure dissipation (derived from the non-softening deterioration curve) with that which can be attributed to the strain softening behaviour (the difference between the total deterioration in the softening model and the non-softening model). See Figure 15f. This shows that the initial rapid dissipation of pore pressure is largely responsible for the initial decrease in FoS during the first 10 years following excavation which is also reflected in Figure   15e where the FoS for the non-softening and softening models is identical. From this point, the relative importance of strain softening begins to progressively increase over time.

Discussion
An approach has been presented for modelling long-term deterioration of highplasticity clay cut slopes with seasonal stress cycles driven by daily boundary conditions.
Hydrogeological behaviour of the numerical analysis has been validated against transient pore water pressure cycles measured in the field over a 16-year period.
Mechanical behaviour has also been investigated although it has not been possible to validate this against field data so instead laboratory tests have been replicated.
Furthermore, the ability of the slope modelling approach adopted here to replicate displacements driven by surface water flux boundary conditions has been previously Seasonal pore water pressure driven stress cycles are controlled predominantly by the near-surface hydraulic conductivity and the surface-atmosphere interaction, whereas the long-term dissipation of the construction induced pore water pressures is dominated by the deeper, lower hydraulic conductivity material, as shown by Vaughan and Walbancke (1973). These distinct pore water pressure changes are J o u r n a l P r e -p r o o f Journal Pre-proof shown to interact, whereby there is superposition of the short timescale near-surface changes onto the pore water pressure response at depth. This causes suppression of the rate of pore water pressure recovery during drying periods and an acceleration during wet periods. The wetting periods correlate with an increase in the magnitude of annual displacements due to the transient reduction in effective stress and hence shear strength (see for example Lade, 2010). This in turn leads to plastic strain development and depending on the magnitude of this strain, softening to occur.
Therefore, while the decrease in FoS due to wet periods may be transient due to effective stress change, this in turn can cause permanent reductions in shear strength and hence FoS due to softening (Take and Bolton, 2011;Postill et al., 2019;Rouainia et al., 2020). This critical long-term behaviour has been captured within the presented modelling approach.
The model results demonstrate a change in potential failure mechanism through the life of the slope, whereby the initial critical failure surface is deep seated and rotational, with the primary deterioration being caused by stress relief and post-excavation pore water pressure recovery. With increasing numbers of seasonal cycles, stability of the near-surface zone, which experiences the largest annual stress cycles, begins to deteriorate due to softening and there is The field evidence from Newbury cutting indicates that a zone of elevated nearsurface saturated hydraulic conductivity developed rapidly after construction (within 20 years) along the full length of the slope . This hydrological deterioration has a significant impact on the surface water balance and pore water pressure cycling within the slope. In order to replicate this observation, the modelling included an elevated saturated hydraulic conductivity and adjustments to the soil water retention properties of the material in this weathered zone. This allowed the numerical model to replicate measured pore water pressure cycles in the nearsurface zone with good accuracy. However, this was included in the model immediately after construction, whereas in practice this weathered near-surface zone would develop over time due to vegetation growth and cycles of wetting and drying.
The hydrological parameters will change concurrently with mechanical parameters (Stirling et al., , 2020. In future work, improvement of the numerical modelling approach may be possible by employing a constitutive model that is able to replicate this behaviour and more fully capture the coupled mechanical and hydrological aspects of the near-surface weathering process. At present, such a model does not exist.
The modelling also demonstrates the importance of prolonged "wet" periods, defined as a period of limited summer suctions that are dissipated entirely during a subsequent wet winter. This combination leads to elevated pore water pressures for a significant proportion of the year. The magnitude and duration of the associated effective stress cycles cause a change in the mechanical slope behaviour with increased displacements during both the swelling and ratcheting phases. Larger seasonal displacements in wet periods increase the generation of plastic shear J o u r n a l P r e -p r o o f Journal Pre-proof strains, leading to increased strain-softening, redistribution of stresses and progressive failure. This process increases the rate of stability deterioration and hence reduces the long-term performance of the slope.
The average size of annual displacement increments was also shown to increase with time, even though the annual pore water pressure cycle size (see Figure 15a), and the average annual pore water pressures (Figure 17d to h) are very similar once the post-excavation equilibration has completed. This change in annual displacement is primarily a function of strength change rather than changes in amplitude in effective stress cycles. This behaviour along with the difference between the strain-softening and perfectly plastic Mohr-Coulomb models (Figure 15e and f) indicates that the long-term reduction in FoS (i.e. slope performance) is a result of a permanent deterioration in soil strength due to material softening rather than solely a function of transient changes in pore water pressures. This can be seen in the FoS vs time plot (Figure 15e), with the FoS reducing equally for both the strain-softening and perfectly plastic models for the first 5 years of the life of the slope, followed by an increase in the rate of deterioration past this age, which is characteristic of progressive failure. However the majority of the deterioration in FoS for ¾ of the slope"s life is due to effective stress changes caused by negative pore pressure dissipation and annual pore pressure cycling along with a small but gradually increasing contribution from strength reduction due to strain-softening.
Following on from this period where the behaviour is dominated by pore pressure change, is an increase in the rate of deterioration which occurs once the near surface reduction in strength leads to a change in the critical failure mechanism.
When this occurs, the near surface seasonal ratcheting mechanism becomes J o u r n a l P r e -p r o o f Journal Pre-proof dominant over deep-seated shear strains due to dissipation of excavation excess pore water pressures. The seasonal ratcheting causes progressive reduction in strength in the near surface ultimately leading to a shallower (3.0 m depth, compared to 7.5 m for the initial critical depth <25 years after excavation) first-time failure within the slope which is of a type often observed in grass-covered clay highway slopes (Briggs et al., 2017). Field monitoring and modelling work (Smethurst et al., 2015;Tsiampousi et al., 2017) has demonstrated that the type of vegetation present can have a positive impact on the stability of slopes in high plasticity clays (while having a negative impact on their serviceability), and conversely vegetation removal can have a negative impact or even be shown to trigger failure. In this work, the vegetation is assumed to be a grass through the life of the slope. In reality this is a simplifying assumption and it is likely that shrubs and in turn trees will develop over time. As the presence of trees has been shown to contribute to stability, this assumption of grass cover throughout is conservative in terms of overall stability but may have implications for the size of annual pore water pressure cycles and further near surface deterioration.
The modelling approach presented captures the primary factors influencing longterm strength deterioration in cutting slopes formed in high plasticity clays (i.e. strainsoftening with large reduction in strength from peak to residual and progressive failure). It has been shown that the near-surface weathered zone on the slope plays a fundamental role in long-term behaviour by influencing the in-and out-flow of water at this boundary. However, it is currently not possible to model weathering processes that lead to near-surface hydrogeological deterioration following excavation of the slope and how this develops over time. In order to further develop understanding of J o u r n a l P r e -p r o o f Journal Pre-proof Elsevier Production Office Reference: ENGEO 105912 Article reference: ENGEO_2020_347 49 slope deterioration behaviour, new constitutive models are required that allow the progressive development of this near-surface zone to be replicated, including the linked changes in saturated hydraulic conductivity, void ratio, water retention behaviour, stiffness and shear strength (Stirling et al., , 2020. The effect of changing vegetation over time and the impact of mechanical root reinforcement will also influence stability and the development of a combined root hydrological and mechanical model would be a significant advance on current capability.

Conclusions
The paper details a numerical modelling approach that allows forecast of long-term deterioration in the performance of cut slopes formed in high plasticity clays subject to seasonal weather sequences. The following are the main findings:  Uniquely, hydrological behaviour of the model has been validated using 16 years of field measurements made at multiple locations in a slope. Daily weather inputs have been used to replicate seasonal cycles of pore water pressures.
 Coupled hydrological-mechanical behaviour using an approach validated previously, has been used to forecast long-term deterioration of slope performance as measured by deformation and factor of safety against ultimate limit state failure. It is shown that a seasonal near-surface ratcheting mechanism leads to development of a shallow softened zone (i.e. with reduced shear strength).
 Seasonal pore water pressure cycles in the slope are regulated by surfaceatmospheric interactions coupled with elevated hydraulic conductivity of the near-surface weathered zone that forms post excavation.
 These cycles of pore water pressure, and hence effective stresses, lead to strain softening, strain accumulation and progressive reduction in soil strength in a zone that propagates from the toe of the slope with number of cycles.
 Prolonged periods of elevated pore water pressures (i.e. during wet winters) leads to increased magnitude of slope displacements.
 Cycling of near surface pore water pressures also reduces the rate of dissipation of excavation generated negative excess pore water pressures at depth during dry periods and increase the dissipation rate in wet periods.
 Over a period of decades, the temporal change in pore water pressure regime in the slope alters the dominant potential failure mechanism from deep seated during the initial response to excavation stress changes to a shallow mechanism driven by weather slope surface interactions.
 Modelling outputs show that long-term reduction in factor of safety against shear failure is a function of permanent deterioration of soil shear strength and not just a direct response to seasonal cycles of pore water pressure. This explains observations that a slope can fail under wet winter conditions even though it has remained stable in the past when subject to the same conditions.
This study advances the numerical modelling capability for analysing long-term behaviour of clay cut slopes. It incorporates the governing soil behaviour in a J o u r n a l P r e -p r o o f Journal Pre-proof