The impact of evaporation induced cracks and precipitation on temporal slope stability

The stability of a dike is influenced strongly by its water content, via both changing its weight and strength. While safety calculations using both analytical and numerical methods are well studied, the impact of surface boundaries exposed to natural conditions is rarely considered, nor is the fact that this surface is covered in vegetation and is susceptible to cracking. This paper presents a numerical study of stability of dikes covered with grass, subject to meteorological forcing and crack formation due to drying conditions. Building on a previous study and adding the impact of cracking, a crop model and a Finite Element Method (FEM) model are integrated together using an optimisation method to ensure mass balance and consistency. The crop model, used to simulate vegetation growth and infiltration/evaporation in response to meteorological forcing, is modified to consider preferential flow due to cracking. The FEM model, used to simulate the dike stability and hydro-mechanical behaviour, has the material properties modified to simulate the impact of cracks. Results simulating a ten-year period indicate a strong impact of cracking on the factor of safety. The vegetation was found to be responsive to both crack presence and an increase in the amount of cracks, which suggests that monitoring vegetation could be a useful tool to identify cracked (vulnerable to cracking) locations along dikes.


Introduction
More than 50% of the Netherlands is below the rivers or the sea level [1]. To prevent these areas from flooding, dikes are built along river banks and the coast. Inland waterways are protected from flooding the land by so-called secondary dikes of which there are about 14 000 km of in the Netherlands. Though dike failures are seldom, their consequences can be substantial. For example, the peat dike in Wilnis failed in August 2003 when a weight reduction of the dike, due to drought, led to a horizontal shear failure [1]. Around 300 houses were flooded as a result, and material damages amounted to around EUR 10 million [2].
The majority of secondary dikes have a vegetated surface which significantly affects the water infiltration and evaporation into and out of the dike [3,4]. During extended periods of drought, significant water loss can occur which can significantly impact the mechanical stability of the dike. For example, peat levees that become dehydrated can shrink and lose a significant portion of their self-weight, a situation that can predispose them to instability (due to uplift or low shear strength due to low confining pressure) or overtopping and ultimately lead to a breach [5]. Dikes which have other swelling and shrinking soils (soils which change in volume in response to a change in water content) will not lose as high a proportion of mass, but may suffer from cracking. The amount of volume change depends on the amount and type of clay minerals and water content change. In addition, as with all soils, the stress history affects the volume change behaviour. In dikes, shrinkage due to extreme drying may result in the occurrence of shrinkage cracks, which can weaken the soil structure and provide favorable conditions for rain infiltration. As a result, the overall shear strength of the soil and the factor of safety (FoS) of the slope can drop significantly [6,7]. Additionally, rainfall infiltration into a dike body through surface fractures will occur faster and will increase the weight (and overturning forces) of the dike and reduce shear strength derived from soil suction. Assessing the impact of cracks on the infiltration at the soil surface and subsequent redistribution of water within the soil is important to characterise hydro-mechanical behaviour. Both processes are different compared to non-shrinking soil, for example, due to changes in surface runoff and preferential flow in the cracks [8].
Most of the current hydro-mechanical models for slope stability analysis are based on the continuum modeling approach, as explicitly simulating the cracking process and preferential flow is difficult at a structure scale. However, some numerical analyses have been conducted to study smaller scale soil cracking using the finite element method (FEM) (e.g. [9][10][11]), discrete element method (e.g. [12,13]) and universal distinct element codes (e.g. [14]). These studies focused on the initiation, development and pattern of cracks, as well as the mechanical properties of the cracked soil [15].
In the previous study [16], a crop growth model was used to simulate the growth of grass cover on a dike surface, and this was integrated to a FEM model to quantify the influence of vegetation on the FoS. It was shown that vegetation is strongly coupled to the moisture available in a dike, and particularly in the root zone [16]. Meteorological aspects also govern this value, therefore a seasonal change is also seen. Some studies have highlighted the value of using the condition of the vegetation as an indicator of subsurface conditions. For instance, Hasan et al. [17] concluded that grass growing in areas with cracks and fractures were stressed during winter and early spring due to a lack of moisture compared to grass in adjacent areas.
In [16], the FEM model is a fully coupled hydro-mechanical model. The crop model, however, utilises meteorological data (e.g. precipitation and radiation) to simulate the root zone water balance and vegetation growth and does not include other coupled proceeses, such as the temperature distribution. The two models are one way coupled together, i.e. the output of the crop model is used as an input into the FEM model. Other models, e.g. [18] include tightly coupled Thermo-Hydraulic surface boundary conditions and model in a tightly coupled manner, although they do not consider the effect of vegetation on the boundary fluxes. This current work follows the one-direction coupled approach.
One of the limitations of the previous study [16] was that it did not account for cracks. The objective of this study is to further develop the approach of [16] to account for both the development of cracks and their impact, i.e. preferential flow and the reduction in strength of the soil. As in the previous study, a crop model and a geotechnical model have been integrated together, resulting in a quantification of the soilvegetation-atmosphere (SVA) interaction and the temporal FoS. Both the crop model and geotechnical model have been modified to allow for consideration of cracks and cracking. A case study is then provided, with real atmospheric data and an idealised dike, to demonstrate the performance of the model.

Method
The idealised regional dike is used in this study is shown in Fig. 1. A permanent grass cover is on the top surface of the dike, which has a fixed root zone depth, in this example of 40 cm (red area in Fig. 1). The dike is of limited height (2 m) and is 41 m wide. The dike is used to retain fixed water levels, of 0.5 m above the ground level on the left hand side and a relatively high ground water level (0.5 m below the ground surface) on the right hand side.
The soil surface is a boundary with energy fluxes and water fluxes forced by meteorological conditions, the bottom boundary is assumed to be impermeable and the left and right boundaries are controlled by the water levels in the groundwater or water body adjacent to the dike. The vadose zone is of limited thickness and consists of the root zone and a portion of the upper dike material.
The strategy for the modelling is: (i) A 1D crop model is used to model the water balance in the root zone, including the impact of the climate and vegetation. (ii) A 1D geotechnical model is used to optimise the hydraulic material properties of the root zone so that the drainage from the bottom of the root zone matches in both models. (iii) The hydro-mechanical behaviour of the dike is simulated in a 2D geotechnical model, giving the temporal FoS, pore water pressure (and water content) and deformation. The description of the model and workflow are described in detail in our previous studies [16,19]. The overview of this process is given in the flowchart in Fig. 2. The individual sub-models which have been used in this study, i.e. LINGRA and PLAXIS, are considered validated ( [20,21]).
The modelling strategy and sub-models have been updated from the previous work [16] so that the crop model predicts when the root zone cracks, tracks the amount of cracks and calculates the drainage through both the cracks and the soil; the 1D geotechnical model is optimised based on total drainage from the crop model at different stages for each cracking event; and the soil shear strength is dependent on the amount of fractures in the root zone. It is assumed that cracking only occurs in the root zone with a constant thickness over the full depth of the root zone.
The outputs from the model (shown in the flow chart) are: drainage from the root zone (from the crop model (D L ) and 1D geotechnical model (D P )); the leaf area index (LAI): a measure of the amount of vegetation; crack volume and area in the root zone (V A , crack crack ); boundary net flux (Q net ) which is "effective precipitation" (precipitation minus interception) minus soil evaporation and transpiration (collectively referred to as evapotranspiration); displacement and Factor of Safety (FoS) which are outputs from the 2D geotechnical model. The workflow of integrating these two models is controlled via Python.
In the following sections the updates for simulating cracks in both the crop model and the geotechnical model are described.

Crop model
The crop model LINGRA (LINTUL GRAssland) [20] is used to solve the water balance in the root zone of a grass cover. The main components of interest are the water balance and leaf growth. More information on the model can be found in [22,23]. In this research, LINGRA has been adopted to incorporate the development of cracks and the impact of cracks in the water balance in the root zone, and consequentially the leaf growth.

Simulating crack development in the crop model
It is assumed that the cracking occurs in the upper layers of the soil, and therefore the soil is under very low confining pressure and any potential volume loss results in either cracking or subsidence, rather than volume expansion of the soil matrix, which would cause tension. This assumption is not valid at greater soil depths and therefore is limited to the root zone. Also, it is assumed here that there are no initial cracks or macropores due to worm and root holes.
The intact soil is considered to be composed of solid material and pores (see Fig. 3). As a soil shrinks, the solid particles stay the same size, move and rearrange so that the void space is reduced and the soil shrinks ( Fig. 3(a)). Soil shrinkage can occur in both the vertical and horizontal direction. Vertical shrinkage generally leads to soil surface subsidence and horizontal shrinkage results in cracks, as shown in the right hand side image of Fig. 3(a).
The volume fractions are shown in Fig. 3(b). In an intact soil, the total volume is made up of pores and solids. When a soil shrinks due to moisture loss (see Fig. 3(b) right hand side), the solid volume fraction (V solid ) remains the same in reference to the original volume and a portion of the pores reduces in volume, i.e. shrinks. Note, that it would also be possible to consider volume fractions in respect to the soil matrix, or the original volume soil matrix including the fractures and subsidence. The crack volume fraction with respect to the original volume is calculated from overall and vertical shrinkage as: where V V ,   subsidence volume fraction respectively, as shown in Fig. 3(b). It is also useful to calculate the volume fraction of the cracks in reference to the soil matrix and cracks volume, as this is also equal to the area fraction of fractures on the surface of the dike and therefore can be used to calculate flow through the fractures. This can be calculated as: The relation between the proportion of soil shrinkage that results in cracking and subsidence can be governed via: where r s is a dimensionless geometry factor which determines the partition of total volume change over change in layer thickness and change in crack volume [24]. In case of three-dimensional isotropic shrinkage, r s = 3. When cracking dominates subsidence r s > 3, when subsidence dominates cracking 1 < r s < 3. In case of subsidence only, r s = 1.
The shrinkage volume fraction, V shrinkage , is equal to the fraction loss of the pore volume, which can also be calculated from the void ratio ( = e V V / pores solid ) and the volumetric water content at saturation:

Shrinkage curve
The matrix shrinkage is a function of volumetric moisture content and material shrinkage behaviour. A shrinkage curve or shrinkage characteristic curve describes the relationship between soil volume and soil moisture content [8]. Initially, a very loose saturated soil may shrink and remain almost fully saturated. As the amount of water in the soil reduces, the soil will typically shrink less in proportion to the amount of water reduction, resulting in de-saturation.
Starting from a completely saturated shrinking soil and drying Bronswijk, [25] identified the following stages: (1) structural shrinkage; (2) normal shrinkage; (3) residual shrinkage, and (4) zero shrinkage. According to Bronswijk [25], the first stage would only occur in wellstructured soils. In this stage, the macropores empty without   considerable change in aggregate volume, and air enters the macropores. Many forms of shrinkage curve exist. A convenient one relates the void ratio to the moisture ratio ( = μ V V / water solid ), i.e. both the amount of water and voids are related to the solid soil volume as a reference. The moisture ratio can also be calculated from the more well recognised volumetric moisture content ( = θ V V / water soil ), by dividing by + e (1 ). In this work, the approach for the shrinkage curve of Kim [26,27] is followed: where e 0 is the void ratio at = μ β 0, K and γ K are dimensionless fitting parameters and μ sat is the moisture ratio at saturation.
It is assumed that cracking occurs when the soil shrinks and it does not recover or seal. Therefore, the shrinkage curve is only used for drying, and the void ratio is only updated when the moisture ratio is lower than that which the soil has previously experienced.

Water balance in cracked matrix
LINGRA uses a tipping bucket approach to solve water balance in the root zone [20]. Hence, ponding in the cracks is disregarded. A fraction of precipitation infiltrates into the crack matrix (I P crack , ) and goes directly to the lower layers, and the rest infiltrates into the soil matrix (I P matrix , ) or runs off, if the flow capacity is exceeded. The daily changes in the amount of water stored in the soil (matrix) is therefore: where In is leaf interception, ET matrix is the evapotranspiration, D matrix is the drainage rate and Rn matrix is the runoff, all from the soil matrix. All quantities are in m − day 1 (due to the model being 1D). These are quantified using the approach of LINGRA [20].
The proportion flowing into the cracks is assumed to be equal to the surface area of the cracks, which in turn is equal to the volume fracture, as the cracks are assumed to be only vertical.
where A matrix is the surface area of the soil which is equal to − A 1 crack . When the amount of water in the soil matrix reaches the field capacity (θ fc ), which is the maximum water storage capacity of the root zone, defined as the volumetric water content at a soil moisture suction of 10 kPa or pF 2.0 [28], the excess water drains from the lower boundary of the root zone. The total amount of drainage (D L ) to the layer below the root zone, as shown in Fig. 4, includes drainage from both the soil matrix (D matrix ) and cracks (D crack ), calculated as: Due to cracking, the field capacity changes. Here, as the computational domain remains the same size, the field capacity relates to the total (original) volume.
The field capacity for the cracked soil matrix (θ fc matrix ( ) ) is calculated assuming that the degree of saturation ( = Sr V V / water pores ) at field capacity for the soil matrix (Fig. 3) remains constant in the intact and cracked soil: where Sr intact and Sr matrix are the degree of saturation in the intact and cracked soil, respectively. The field capacity is related to the degree of saturation by: (10) and therefore, where the subscripts intact ( ) and cracked ( ) related to the intact and cracked soil, respectively.

Geotechnical model
The commercial finite element code, PLAXIS 2D, was used in this study. The workflow is controlled via the PLAXIS Python interface. Further information on hydromechanical and safety analysis of this model can be found in [21].
Cracks, in addition to providing preferential flow channels which increase the soil permeability, also decrease the soil strength [29]. Furthermore, cracks can form a part of the critical surface and therefore can ultimately influence the stability. In this study, the impact of cracks is considered as a bulk effect via a change in the shear strength parameters. This allows for a relatively simple method of quantifying the impact of cracks, without predicting complex crack patterns, orientations or very local changes in the critical failure surface. Moreover, cracks are considered to only extend over the root zone. The calculated crack volume from the crop model is used in this geotechnical model to update the mechanical strength of the cracked root zone. The shear strength parameters, cohesion (c) and friction angle (ϕ) have been reduced according to the crack volume, using as a first approximation, a linear relationship. The values of c and ϕ for an intact soil and the minimum value of c and ϕ for a maximum V crack are input by the user.

Case study
The example dike, shown in Fig. 1, has been investigated. The dike is considered to be covered with perennial ryegrass which has the majority (85%) of its root system in the shallow soil layer of 0-40 cm below soil surface [23], therefore the root depth is considered fixed at 40 cm.
The numerical experiment was performed for ten years with climatic data, from 2009 to 2019, used to obtain a time series of FoS to investigate the influence of meteorological conditions on the soil shrinkage and cracking behaviour and, consequently, the dike safety. Furthermore, the impact of cracking on vegetation growth was investigated.

Input data
The meteorological data were obtained from the Royal Netherlands Meteorological Institute (KNMI) station at Schiphol Airport (Amsterdam) (52°′ 19 04°′ 47 OL). In Figs. 5(a)-(d) precipitation (P), average air temperature, wind speed and solar radiation for the 10 years simulation period is shown. These data are used as inputs for the crop model.
The key material parameters used by the crop model for both the soil and vegetation are listed in Table 1. In the soil parameters, the intact (volumetric) water content field capacity (θ fc intact ( ) ) is the Fig. 4. Updated root zone water balance modified from [16] for cracked soil matrix to include water flux in cracks. maximum water storage capacity of the root zone. The water content at the wilting point (θ wp ) is the limit of water content, below which plant water uptake ceases and plants start to wilt. Below the critical water content (θ cr ), transpiration is reduced due to water stress. Both θ wp and θ cr are assumed not to change for intact and cracked soil. Finally, drainage is limited by the maximum drainage rate (DRATE) of the subsoil.
For the vegetation (perennial ryegrass) the Specific Leaf Area (SLA; leaf area/ leaf mass) determines how much new leaf area to deploy for each unit of biomass produced. The critical leaf area (LAICR) is the value beyond which death due to self-shading occurs [30]. The values used here are typical for Dutch soil conditions and typical grass cover, based on reported values by [23].
Typical parameters required for the shrinkage curve (Eq. 5) were obtained from literature ( [27,31] . In Fig. 6 the calculated shrinkage curve for this study is shown by the solid blue line; the two dashed lines are selected measured shrinkage characteristics of clay and peaty soils in the Netherlands, as described by [32,33], respectively; and the dotted green line is the saturated line. In this study, the isotropic shrinkage has been considered, so r s = 3. As shown in Fig. 1 the example dike consists of the root zone and the soil of the dike body. Constitutive and hydraulic input parameters for those parts of the dike are listed in Table 2 for the intact soil matrix. The values are based on the default soil properties from the PLAXIS library for the root zone (silt clay) and for the dike body (organic clay), except for the hydraulic values of the root zone which are obtained from the optimisation code. The initial parameters optimised for intact soil are shown here. Since the cracking occurs only in the root zone, the dike body parameters do not change as A crack increases. The shear strength parameters (c and ϕ) and the hydraulic parameters change for the root   ) was extracted (10 %). Then shear strength values were picked to ensure that the model had a FoS > 1, so that it would continue to run. The value of shear strength parameters were selected arbitrarily for demonstration purposes and for more realistic analyses they should be calibrated for real cases.

Results
In Fig. 7, temporal results from the crop model for soil which is able to crack is compared with soil which cannot crack, i.e. the model presented in [16]. As shown in Fig. 7(a), the cracks cannot seal during wet periods, but only increase in conditions drier than previously encountered. This assumption ensures that the worst-case scenario has been considered. In Fig. 7(a) the crack area increases from spring 2009 and gradually increases from 0 to 6 % by May. Wet days from May -August 2009 ensure that the crack amount remains constant until the next dry period in June -July 2010 during which cracks increase to 7.5 %. During summer 2011, cracks grow again (9.3 %) and the soil experiences the next drier condition in the summer 2018, when cracks again grow (10%). As the crack expands only in drier conditions than have been previously encountered, the time between cracking events gets longer as the analysis progresses.
In Fig. 7(b), temporal variations in LAI are shown. The LAI is highest in spring and summer, since reduced solar radiation limits growth in the autumn and winter. Higher LAI values in the summer lead to higher evapotranspiration, and hence a reduction in the amount of water flux into the dike (Fig. 7(c)). The mowing events on 15 June and 15 August were imposed in the crop model based on the mowing schedule for regional dikes in the Netherlands [16], which can be seen by the sudden decrease in LAI. The vegetation growth can be seen to be influenced by the presence of cracks, due to a portion of the precipitation draining directly to the lower soil layer and a reduction in the maximum stored water (seen in Fig. 7(d)). In the case of a cracked soil, the LAI is lower or equal at almost all times than the case without cracks.
Positive Q net values occur in wet periods when precipitation exceeds evapotranspiration demand. When there are cracks, the boundary net flux (Q net ) is seen to be higher than for the un-cracked soil (Fig. 7(c)). In the cracked soil, the combined effect of drier root zone and lower LAI cause lower evapotranspiration and leaf interception.
In Fig. 7(d) it is shown that the water content in the root zone (θ rz ) decreases during the summer due to high levels of evapotranspiration. On the other hand, during wet periods with a consistently high Q net , the root zone reaches the field capacity, and extra water drains to the lower layer. In the cracked matrix, θ rz is lower as some rainfall passes directly through the cracks and does not enter the soil root zone. Additionally, the field capacity is reduced, therefore the maximum amount of water stored is reduced. In summer 2010, the crack area increases due to the dry period which influences the vegetation growth, i.e. after first mowing in June 2010, due to the very low water content in the root zone ( Fig. 7(d)), vegetation cannot recover in the growing season. The same situation happens in May -June 2011 and June -August 2018, when the average LAI is very low over almost 3 months in both cracked and un-cracked cases.
Drainage (D L ) occurs when there is a positive (downward) Q net and θ rz reaches the field capacity ( Fig. 7(e)). This can generally be seen in the winter months. For example, spikes in D L in August 2010 and 2017 correspond to the heavy rainfall and therefore high Q net . Infiltration of precipitation through cracks and the reduction in the field capacity in cracked soil, which both increase with A crack rising, causes D L to increase.
The results of the optimisation procedure for the hydraulic parameters of the geotechnical model are shown in Table 3. It is seen that, in general, as the crack area increases, the hydraulic conductivity increases. In addition, different soil water retention curves (SWRC) for  These are shown in Fig. 8, in which Sr is the degree of saturation, and h is the suction height above the phreatic level. In general, more cracks are associated with a drier root zone. In Fig. 9 the geotechnical model outputs are illustrated. Fig. 9(a) shows the crack percentage (from the crop model) for convenience. Fig. 9(b) and (c) show the pore water pressure (pwp) at points A (in the root zone) and B (in the dike body), shown in Fig. 1, respectively. Positive values are for compressive pressures and negative values indicate suction. As expected, high levels of drainage (August 2010 and 2017), or long periods of cumulative drainage (winter 2009-2010), lead to higher pwp in the root zone and dike body. In both locations (A and B) pwp is higher in case of cracked soil as more water reaches the soil system via the higher Q net and more D L in the cracked soil. As A crack increases, pwp rises and decreases more slowly.
In Fig. 9(d) and (e) the magnitude of displacement of points A and B is shown. The displacement rises following large Q net and recovers between precipitation events. A slight accumulation of displacement over time is observed, due to plastic displacement. The displacement of the points increase as crack grows during time, which depicts the effect of shrinkage behaviour in the root zone (where cracks exist) and more drainage into the dike body. By increasing the crack area the difference in displacement between the crack and un-crack condition gets more noticeable.
The temporal variation of FoS is shown in Fig. 9(f) from 2009 to 2019. Safety in the dike is responsive to the climate and vegetation condition. The safety of the dike in the cracked condition is significantly lower than the case without cracks under the combined effect of more infiltration into the dike (Fig. 7)) and lower shear strength induced by modified cohesion and friction angle. The maximum crack area leads to the minimum shear strength parameters (Table 3), thereby generally lower FoS. During the simulated period, results from Figs. 7 and 9 suggest that more cracks lead to a lower amount of vegetation (LAI decreases) and a lower amount of stored water in the root zone. In general, this leads to lower safety in the dike. As found in the previous studies [16,19], heavy rainfall events cause a dramatic decrease in the safety. Therefore, it is the combination of cracking, due to drought, which reduces the strength and general level of safety, and heavy rainfall events which significantly lower the safety temporarily.

Using vegetation as an indicator for dike health
In the current study it is shown that vegetation responds to the presence of cracks, which influences the available water in the root   zone and therefore makes more cracking likely. Consistent with visual observations from dike inspectors in the Netherlands, in summer 2018 the water content in the root zone was extremely low and the LAI was low for an extended period of time (Figs. 7(b) and (d)). Visually it was seen that a substantial amount of the grass cover died and took several months to recover, this is shown in the simulation (Fig. 7(b)). In the analysis, it was also shown that A crack increased in August 2018, Fig. 7(a), after about 7 years of no increase. This summer was the warmest summer during the simulation period. After the mowing events in 2018, the vegetation is seen to be able to partially recover in the simulations without cracks, whereas it cannot in the cracked root zone (Fig. 7(b)). Therefore, it is proposed that vegetation growth could be used as an indicator of crack presence, also indicated by [17]. The cumulative precipitation, root zone saturation and LAI for 2017 and 2018 are plotted in Figs. 10(a)-(c), respectively. The former year considered as a 'wet' year and the latter one is the driest year in the 10year simulation. Before the first mowing, the amount of precipitation and consequentially the available water in the root zone was similar. However, in the summer, growing season, the precipitation in 2018 was less than 2017 which led to a dry condition and consequently the crack volume grew in July 2018 (red line in Fig. 10). The difference in vegetation growth is significant and seen by the difference between LAI in the following months. In first days of September, the rainfall event for both years is almost the same and in both case the water content reaches the field capacity in the root zone ( Fig. 10(b)), but due to the larger crack area in 2018, the vegetation cannot recover as much as in September 2017. This indicates that the LAI could also be used directly to identify cracked dikes which need maintenance. However, this does not seem to occur consistently throughout all years (see Fig. 7).
To further investigate the use of vegetation as an indicator in more detail,two differential LAI values are shown in Fig. 11. The first differential, i.e. the velocity or growth rate, is (the value is not divided by the time window), where t is the current day and x is the time window (Fig. 11(a)). The absolute of the second differential, i.e. the absolute value of vegetation growth acceleration or rate of change of growth rate (again, note that the value is not divided by the time window), with a time window of 15 days, is shown in Fig. 11 [ ]|) are identical for the cracked and un-cracked cases. As the crack area increases over time, the differential LAI time series (Fig. 11) can be categorized into periods where: (1) the two lines are virtually indistinguishable; (2) the cracked simulation exhibits a lower acceleration; or (3) the cracked simulation exhibits more variability in both cases of velocity and acceleration growth. These categories occur at different times of year and under different degrees of water stress, as highlighted in Fig. 11 Fig. 11) in the following periods: January -August 2009, January -April 2010, January 2012 -January 2013, January -August 2016. This is observed to be when either the crack area is starting to grow (January -August 2009) or there is a moderate amount of LAI and the root zone water content is reasonably high. In particular, in the whole year 2012 and the June -July 2016, the water content of the root zone remained higher than other years, and it is high enough for the vegetation to grow even over cracked areas. These are seen to be in situations with a low water content in the root zone, and low values of LAI, which implies dry periods. In particular, between mowing times, there seems to be significant differences between cracked and un-cracked simulation results in dry years (2011,2014,2015,2018). In dry periods, less water is available in the root zone in the cracked area than in the uncracked area (Fig. 7(d)). Therefore, vegetation cannot regrow after mowing over the cracked area as much as it regrows over the healthy areas. This low variability mainly occurs in the summer, and depending on the extent of the drought extends through the following year. However, this does not happen in the wet years (2012, 2013, 2016) as explained in the previous paragraph. At the times there is generally a moderate amount of LAI and the root zone water content is relatively quickly increasing or decreasing due to heavy rainfall after a dry period and couple of days in a row with the negative Q net , respectively (Fig. 7(c)). This is generally observed in the spring and autumn periods, when the energy for vegetation growth is limited and LAI variation is mainly responsive to the θ rz variation.
This suggests that monitoring to investigate the dike should be timed accordingly. In periods of moderate LAI and precipitation (generally winter), almost no differences are likely to be observed. In the summer periods there are more significant differences in dry years and in the spring and autumn much more variability is seen in the cracked soils in almost all years. Therefore, monitoring prior to the summer (when more cracks may occur) or prior to the winter (when the lowest safety is seen) is advisable.
Results in Figs. 10 and 11, support the argument that vegetation could be used as an indicator to distinguish whether a dike is significantly cracked. By observing anomalies in vegetation, further more targeted investigation can be planned. In addition, as discussed in the results section, displacement can be used as a proxy for both saturation (short term changes) and for accumulation of cracks (long term changes), although long term changes may also indicate subsidence or other processes. Displacement and vegetation indices can be obtained from Earth Observation (EO) data. Interferometric Synthetic Aperture Radar (InSAR) [34,35] is a technique that maps millimeter-scale deformations of the Earth' surface from satellite images. This could be used to monitor both short and long term changes, depending on the satellite overpass frequency [36]. Vegetation indices can be measured from both optical and radar images with fine resolution provided by satellites [e.g. 37,17]. While no absolute value of vegetation indices can be predictive of cracking, anomaly detection identify vulnerable areas that warrant further investigation.

Perspectives on validation
Validation of the modelling methods and results is needed, but has not been possible during this study. The following aspects could be investigated: (i) cracking of a vegetated embankment surface, (ii) the consequential additional inflow, (iii) impact on the bulk shear strength, and (iv) the overall influence on the stability. A number of these aspects are scale dependent and are therefore difficult to observe in laboratory experiments; for example, the cracking is influenced by the vegetation rooting depth, the soil grain size and the root zone properties, which themselves are governed by the atmospheric conditions. The bulk reduction in shear strength properties, while convenient for numerical analysis, is difficult to validate as it relies on knowing the failure surface size, orientation and interaction with individual cracks. The influence on the overall stability of cracks could be validated via either scale model tests in the laboratory, or via full scale failure tests [e.g. 38]. The qualitative behaviour is well supported by, albeit limited, literature. One important field test was the BIONICS research embankment [39] which provided a full scale test where the hydro-mechanical behaviour was monitored, although it was not brought to failure. The additional inflow into cracked vegetated embankments by the use of Electrical Resistivity Tomography was shown by Stirling et al. [40] and the fractures in the same embankment were limited approximately to < 400 mm [41], approximately the same depth as the root zone [40]. It is clear that further experimental validation is needed.

Conclusion
The integrated model framework composed of a crop model and a geotechnical model including the impact of cracking was used to illustrate the sensitivity of the factor of safety to root zone soil moisture and vegetation cover in an idealized dike. This extended a previous study where cracking was omitted. Here, simple modifications or parameterisation was included in both sub-models to account for the formation of cracks. This provides a means to account for the preferential flow into the dike that is associated with cracks in the cover layer and a reduction in shear strength. This represents a step forward to understand soil-vegetation-atmosphere interactions in grass-covered dikes. Simulations with the new integrated model were used to compare vegetation growth and safety under intact and cracked soil conditions. To the authors' knowledge, neither comprehensive field or laboratory data are available to validate this numerical research and it is suggested for the future studies. It is shown that the presence of cracks affects the water flux into the dike and the shear strength, both of which impact the factor of safety (FoS). The history of the precipitation, root zone water content and LAI have an impact on crack propagation. Therefore, vegetation condition (Leaf area index or comparable indicator) and root zone water content could be useful as indicators to detect cracked areas along a dike and also of increasing crack volume at an early stage. Results suggest that monitoring in the spring or autumn may provide the most reliable and useful results.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.