Numerical modelling of desiccation cracking

INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS IN GEOMECHANICS Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 Published online 11 March 2010 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/nag.894 Numerical modelling of desiccation cracking Aruna L. Amarasiri, Jayantha K. Kodikara∗, †, ‡ and Susanga Costa Department of Civil Engineering, Monash University, Australia SUMMARY The ability to model and predict the formation of desiccation cracks is potentially beneficial in many applications such as clay liner design, earth dam construction, and crop science, etc. However, most studies have focused on statistical analysis of crack patterns and qualitative study of contributing factors to crack development rather than prediction. Because it is exceedingly difficult to capture the nonlinear processes during desiccation in analytical modelling, most such models handle crack formation without considering variation of material properties with time, and are unattractive to use in realistic modelling. The data obtained from laboratory experiments on clay soil desiccating in moulds were used as a basis to develop a more refined model of desiccation cracking. In this study, the properties, such as matric suction, stiffness and tensile strength of soil, and base adhesion, could be expressed approximately as functions of moisture content. The initial conditions and the development of suction due to desiccation and the varying material properties were inputted to UDEC, a distinct element code, using its internal programming language FISH. The model was able to capture some essential physical aspects of crack evolution in soil contained in moulds with varying lengths, heights, and materials of construction. Extension of this methodology is potentially beneficial not only for modelling desiccation cracking in clay, but also in other systems with evolving material properties such as concrete structures and road pavements. Copyright q 2010 John Wiley & Sons, Ltd. Received 22 October 2008; Revised 8 September 2009; Accepted 31 December 2009 KEY WORDS: desiccation cracking; numerical modelling; moisture content 1. INTRODUCTION Climate change is affecting the world in many ways. It has made many regions in the world drier, and some regions wetter than before [1]. Though wetting and drying due to seasonal weather changes can be a driving force for desiccation cracking, climate change could potentially exacerbate the harmful effects. There have been reports that increase in the occurrences of breakages in pipe buried in reactive soil and climate change could be a contributing factor [2]. Soil shrinks as it loses moisture, and magnitude of shrinkage can be very high for reactive soils. Clay liners, which are commonly used to contain hazardous and solid waste, may also crack upon drying, rendering them ineffective as a hydraulic barrier [3]. With these possibilities in mind, it is appropriate that tools to model the formation of desiccation cracks be developed and validated. The heterogeneity and nonlinearity of soils, as well as the fact that it is not a material made under factory conditions like many other engineering materials has been a factor in making modelling its exact behaviour difficult. The capability to model desiccation cracking has lagged behind most other geo-engineering fields ∗ Correspondence to: Jayantha K. Kodikara, Department of Civil Engineering, Monash University, Australia. E-mail: Jayantha.Kodikara@eng.monash.edu.au ‡ Associate Professor. † Contract/grant sponsor: Australian Research Council Project; contract/number: ID DP 0773861 Copyright q 2010 John Wiley & Sons, Ltd. 83 NUMERICAL MODELLING OF DESICCATION CRACKING of study. There are many nonlinear processes involved in desiccation. As the soil dries, the formation of menisci causes very high negative water pressures (suction) to develop in the soil. The stiffness of the soil, which is characterized by parameters such as Young’s modulus and shear modulus, increases with increasing suction. The resistance to cracking, i.e. the tensile strength or fracture toughness of the soil also changes upon drying. Finally, the adhesion at interfaces, which is essential in providing the restraint for desiccation cracks to form, changes with moisture content [4, 5]. The four phenomena mentioned above sometimes cause change of several orders in magnitude in material properties, making it virtually impossible for rigorous analytical solutions to be formed. Therefore, it is not totally surprising that, after many decades of research, published literature worldwide still focus on qualitative, semi-quantitative, and statistical analysis of crack patterns [6–9]. Some studies analyse cracking without taking into account the changes in material properties. In the research reported herein, the distinct element code UDEC was used, along with the embedded programming language FISH to successfully model desiccation while capturing the nonlinear behaviour and considerable changes in properties while desiccation was in progress. The laboratory desiccation tests reported by Nahlawi and Kodikara [10] were used as a basis for validation and calibration of the model. The same methodology can be used to study field applications involving drying soil, as well as systems with time-varying properties in fields as diverse as paint, concrete, and pavement industries. 2. LABORATORY EXPERIMENTS 2.1. Soil and desiccation test procedure The basaltic clay used in the laboratory tests described above was obtained from a natural deposit in Werribee, a western suburb of Melbourne, Australia. X-ray diffraction revealed that the main phases present in the soil were Ca-Smectite (42%), quartz (30%), feldspars (8%), illite (10%), and kaolinite (10%) [11]. The soil for the above experimentation was obtained at depths ranging from 1.5 to 2 m below the ground surface. The Werribee clay is considered to be a reactive soil with high shrink swell potential. The desiccation tests were carried out in an environmental chamber within which the temperature and relative humidity both could be controlled. Tests were conducted in perspex and metal moulds with rectangular cross section. The details of the desiccation tests are shown in Table I. The laboratory procedure in carrying out the desiccation tests will be described hereafter [10]. The soil used for making the moulds was prepared from material passing the 425 m sieve using the procedure given in AS 1289.1 [12]. The moisture content of the soil was brought close to the liquid limit of 127 % and then the soil was left covered for 12 h for the moisture content within the sample to stabilize. The specimens were prepared by placing and moulding soil using a spatula. The sides of the moulds were greased to prevent adhesion of the clay, whereas the base was not Table I. Details of desiccation tests on slurry clay [10]. Test no.∗ 0 1 2 3 4 5 Mould dimensions (L × W × D) (mm) Number of moulds used Dry density of soil† (kg/m3 ) Test duration (h) Initial water content (%) Relative humidity (% RH) Mould type 250×25×12.5 250×25×5 250×25×5 600×50×10 600×25×30 600×25×20 4 2 2 2 2 4 — — — 540 590 620 191.5 48.25 29 147.75 272 49 127 127 127 127 134 127 40 20 20 25 25 25 Perspex Perspex Perspex Metal Metal Metal ∗ Tests 0–3 were carried out at 18◦ C temperature and Tests 4 and 5 were carried out at 20◦ C. † The dry density was computed at the start of the test. ‘—’ indicates that dry density was not recorded, but it is considered to be similar to other tests owing to the same test procedure used. Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag 84 A. L. AMARASIRI, J. K. KODIKARA AND S. COSTA greased, and there could be adhesion between the clay and the mould there. This adhesion provided the restraint necessary for the formation of desiccation cracks. Replicate specimens of clay drying in moulds were prepared for the purpose of taking moisture content measurements at intervals so that the moisture content in the moulds where cracking is being observed could be estimated. The slurry clay specimens had an initial dry density ranging between 540 and 620 kg/m3 , as determined using the gravimetric method. 2.2. Soil drying and crack evolution Upon drying, predominantly parallel cracks developed in the moulds, as was expected because of the use of long rectangular moulds for the preparing of specimens. Cracks could be classified as primary, secondary, and tertiary. Primary cracks are the first to develop as the soil dries, and cause the breaking up of the intact material into separate blocks. Often, a single primary crack appears in the centre of the mould breaking the clay into two blocks. Secondary cracks appear subdividing blocks formed by primary cracks, while sometimes tertiary cracks appear that subdivides blocks formed by secondary cracking. Infrequently, cracks will start at one edge but not extend fully to the other side of the mould. In general, the cracks were well distributed over the entire length of the specimen with some variation in crack spacing. Data from the observations made about the study of cracking are shown in Table II. A typical plot of water content measurements against time is shown in Figure 1. They show the average values obtained from the top and bottom layers. The water content of the lower layers Table II. Number of cracks and dimensions of blocks formed [10]. Number of moulds used Mould dimensions (L × W × D) (mm) Number of cracks/mould number Average number of cracks Mean spacing (mm) Final thickness of the soil (mm) 0 1 4 2 250×25×12.5 250×25×5 4 (average) 6 8 1, 2 4 7 35.8 26.6 8.5 2.5 2 2 250×25×5 8.5 19.5 3 3 2 600×50×10 11 35.2 6.5 4 2 600×25×30 5 68 20 5 4 600×25×20 9 8 1, 2 11 , 11 1 2 5 1 5, 7, 6, 4 1 2 3 4 67.7 13 Test no. 5.5 Figure 1. Moisture content variation with time for a Test 3 specimen. Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag 85 NUMERICAL MODELLING OF DESICCATION CRACKING Table III. Time to first crack, cracking water content and desiccation coefficient. Initial soil thickness (mm) Test Number Time to first crack (h) Cracking water content (%) Desiccation coefficient K (h−1 ) 30 4 10 20 5 8 12.5 0 11 10 3 6.5 5 2 5 5 1 3.5 wt 110.42 wb 129.22 wt 102.41 wb 109.24 wt 95.33 wb 101.64 wt 107.4 wb 109.83 wt 94.77 wb 96.69 wt 96.48 wb 115.9 0.022 0.021 0.037 0.036 0.040 0.038 0.057 0.055 0.126 0.124 0.171 0.170 was slightly higher than that of the upper layers for most of the duration of the tests. It is thus apparent that drying took place from the top of the specimen. Some curling of blocks formed was observed. This phenomenon can be attributed to differential drying between the layers. According to the water content measurements, it took about 29-272 h for the desiccation process to reach equilibrium, depending on the thickness of the specimen and the drying environment. It was observed that specimens with greater thickness took a longer time to develop cracks, as can be seen from data shown in Table III. The difference between the top and bottom water contents was generally higher for thicker specimens. However, it was observed that the moisture contents of all the specimens reached an equilibrium moisture content of about 10%. 3. OVERVIEW OF DISTINCT ELEMENT METHOD The software UDEC was used to model the desiccation cracking. UDEC is a distinct element programme by ITASCA, which has been widely used in the industry, primarily for rock engineering [13]. UDEC models a continuum using overlapping constant strain finite difference triangles [14]. It also uses explicit time marching and, therefore, it applies Newton’s equations of motion on the lumped masses at the grid points of the zones over small time steps rather than form a global stiffness matrix. Damping is introduced to the system to allow it to come into static equilibrium. Distinct element programmes can model the breaking apart of material as well as the formation of new contacts. Because rock materials are usually of high strength, its behaviour in bulk is greatly affected by weaker planes, i.e. joints. As such, UDEC is designed to address very efficiently modelling requirements in jointed material, and can be used effectively to analyse the opening of cracks. However, one of the limitations of UDEC is that it can accommodate fluid flow only through joints, and not through media. The limited capability of UDEC to handle fluid flow has limited its use outside analysis of rock. However, if the equilibrium moisture contents throughout the medium are known, or as in the case of the modelling reported here, the time progression of the moisture content in the regions of the medium is known, UDEC is an attractive computer programme for the analysis of desiccation cracking. However, common issues encountered when carrying out explicit analyses such as allowing sufficient time steps for the system to come into equilibrium under each suction increment and keeping each suction increment small enough not to ‘shock’ the system had to be addressed in this numerical analysis. Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag 86 A. L. AMARASIRI, J. K. KODIKARA AND S. COSTA 4. APPLICATION OF UDEC TO MODEL LABORATORY EXPERIMENTS 4.1. Model characteristics and parameters To analyse desiccation cracking using UDEC, representative material properties have to be inputted to the model. The important properties are tensile strength, shear modulus, bulk modulus, and suction of the soil and shear strength and shear stiffness of the soil–base interface. Kodikara et al. [15] and Kodikara and Choi [4] reported findings from laboratory experiments on Werribee clay that are appropriate for use in formulating the material properties. However, it should be noted that all the material properties need to be input as functions of the moisture content in order to capture their variation during the course of desiccation. The programming language FISH, which is embedded in the UDEC code, allows access to all the material properties and facilitates capturing of the variation of properties with moisture condition. The soil block material was taken to be nonlinear elastic with the material properties’ functions of the moisture content at that point of time in desiccation. This is acceptable so long as reversal of strains is not predominant, as in the continuous drying. UDEC requires the input of the bulk modulus K and shear modulus G as the material stiffness properties, rather than the elastic modulus E, and Poisson’s ratio . The joints assigned in modelling Test 1 is shown in Figure 2. The base was considered to be 2 mm thick. The interface between the soil and the base, and the vertical joints in the soil that have been introduced as locations of potential cracks can be seen in this figure. The soil has been split into two layers, each with its own desiccation coefficient k, by a horizontal joint. It should be noted that although this is the generic mesh used in the studies, the dependency of the simulated results on the mesh was studied as a separate item. 4.2. Simulation of drying, suction, and strength development The desiccation progression was input to the program using the internal programming language, FISH. The definition adopted to mathematically describe the desiccation rate was w −wr = (wi −wr )e−kt (1) dw = −k(w −wr ) dt (2) and therefore, Figure 2. Plot of joints assigned in modelling Test 1. Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag NUMERICAL MODELLING OF DESICCATION CRACKING 87 Figure 3. Soil water characteristic curve for the Werribee clay. where wi is the initial water content, wr is the residual water content, and k is the desiccation coefficient with dimensions of [T ]−1 . The values of k which best fits the moisture loss curves with time as given by Nahlawi and Kodikara [10] were used in the modelling. The primary mechanism associated with desiccation is the development of soil suction, which leads to increase in stiffness and strength of fine-grained soils. In general, suction increases as water content decreases and the variation of suction with moisture content is given by the soil water characteristic curve. The experimental points for the soil water characteristic curve for the Werribee clay starting from a slurry has been given by Kodikara et al. [15], and its theoretical approximation used in the simulation is given in Figure 3. The suction  (kPa), in the soil at a given moisture content may be approximated by the following exponential formula:  = exp((136.5−w)/9.92) (3) as shown in Figure 3. As the suction on the soil increases due to drying, its stiffness increases. Kodikara and Choi’s [4] recommendation for the Werribee clay, H = 240.95 (4) where H (kPa) is a modulus which with respect to suction was used. Note that an exponent of unity would be more suitable for Equation (4), especially for materials close to saturation, but the exact regression developed by Kodikara and Choi is used here. Next, Young’s modulus, E, was estimated as recommended by Fredlund and Rahardjo [16] by using the equation E = H (1−2) (5) where  is Poisson’s ratio of the soil. Poisson’s ratio of 0.42 was used since the soils remain close to saturation. Once Young’s modulus has been estimated for a given moisture content, the bulk modulus K and shear modulus, G, follow from the well-known relationships of elasticity. During desiccation, suction stresses tend to cause shrinkage, and restraints may then produce tensile stresses in the soil causing cracking space. In the laboratory tests modelled in the work described herein, the only restraint is provided by the base to the soil immediately above it. Nahlawi and Kodikara [10] reported test results of tensile strength tests carried out on the Werribee clay at various moisture contents (Figure 4). Based on these tests, they recommended the following empirical formula for the tensile strength T (kPa) of the Werribee clay at a given moisture content, w(%) T = 149, 194w −2.3744 Copyright q 2010 John Wiley & Sons, Ltd. (6) Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag 88 A. L. AMARASIRI, J. K. KODIKARA AND S. COSTA Figure 4. Tensile strength data (adapted from Nahlawi [11]). Figure 5. Interface shear strength data. 4.3. Soil/mould interface properties Another parameter that poses significant impact on desiccation cracking is the adhesive strength of the soil–base interface. Researchers have reported that if the base is greased, as in the free shrinkage tests, there may not be any desiccation cracks at all [15, 17]. Nahlawi [11] carried out tests using a shear box to determine the shear strength and the shear stiffness of the base and the Werribee clay at various moisture contents, and the results of these test are shown in Figure 5. The base adhesive strength C (kPa) was then proposed using a Bezier equation as C = 1 (1−u)4 +2 (7) u = (w −5)/3 (8) where The parameters 1 , 2 , and 3 were taken as 200, 4 kPa and 125%, respectively. During the parametric study carried out described fully later on in this paper, it was found that the above formula may overestimate the adhesive strength of the soil–base interface based on the number of cracks that open up during the course of desiccation. It should be noted that the shear strength tests were carried out at relatively rapid rates (peak shear stress reached within about 2 min), whereas Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag 89 NUMERICAL MODELLING OF DESICCATION CRACKING Table IV. Experimental interface properties. Water content (%) Shear stress at failure (kPa) Shear stiffness (103 kPa/m) 17.54 4.4 3.2 2 3.98 6.88 12.31 1.85 27.5 21.9 17.4 12.6 shrinkage due to desiccation happens much slower, with first cracks appearing up to 2 h after initiation of drying. Furthermore, certainly during the initial stages of drying, the soil is close to its liquid limit, and viscous effects may be present. 1 , 2 , and 3 were changed to 125, 3.4 kPa, and 129% for all tests in Perspex moulds and 100, 2.6 kPa, and 129% for all tests in steel moulds, which gave reasonably consistent results with the experiments. Besides the shear strength or cohesion of the interface, another parameter that was thought could affect the behaviour of cracking is the shear stiffness of the soil–base interface. A relatively high shear stiffness would lead to the shear stresses being developed with relatively smaller shear movement between the base and the soil, hence producing stiffer restraint to shrinking soil and vice versa. Additionally, UDEC being a programme that uses the explicit method of solution, the normal and shear stiffness of joints are used in calculating the critical time step. The critical time step has to be small compared with the natural time period corresponding to the largest stiffness to mass ratio of the system. This condition is automatically fulfilled by UDEC when it calculates the time step over which it allows each lumped mass to move under the out of balance forces that may be present in the system. Nahlawi [11] obtained the shear stiffness values at corresponding moisture contents as shown in Table IV, and were used as approximate input to the program. 4.4. Soil/soil interface properties During the parametric study, it was found that the normal stiffness of the soil joints affects the cracking behaviour. This phenomenon is understandable, because the tensile stresses developed at a joint in the model is equal to the multiplication of the normal stiffness of the joint and the opening between the two faces of the joint. It was found that if the normal stiffness (jkn) of the joint was set at the upper limit recommended by the UDEC users’ guide, jkn10(K +4G/3)z, where z is the smallest dimension of all zones on either side of the joint, all joints contained within the middle 10th of the specimen opened up simultaneously when test no. 3 with 600 mm length and 10 mm height soil was modelled. This unnatural modelling result (not found while replicating any of the other tests) is obtained because there is a long flat peak in the tensile stress distribution near the top of the soil layer due to the soil layer being long and thin, as shown in Figure 6. This issue will not be that important for thicker layers due to more less steep tensile stress distribution. Furthermore, in reality the flaws giving low tensile strength within the soil will be activated in the region where stresses are relatively uniform. During the parametric study, the normal stiffness of soil was varied, and it was found that when it was set at one 60th of the upper limit given by the UDEC user’s guide, acceptable crack patterns could be obtained in test no. 3. This value was then adopted for normal stiffness in soil for all tests. The shear stiffness of the joints between soil blocks, which is not considered a significant parameter because cracks open in tension than in shear, had to be prescribed thereafter. It was found that if they are set at the values recommended by the UDEC manual, the execution of the programme was very slow. This inefficiency is due to the wide disparity in joint stiffness values of the interface and cracks within the soil. A compromise where the stiffness values within the soil were set at one-tenth of the upper limit recommended by the UDEC manual led to better solution times. UDEC offers the use of many joint and block constitutive models. The joint model used at the soil–base interface was Mohr–Coulomb failure criteria with residual properties. An alternative joint constitutive model that could have been used was Mohr–Coulomb with perfectly plastic behaviour after the shear strength of the joint has been reached. However, it was though that there should Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag 90 A. L. AMARASIRI, J. K. KODIKARA AND S. COSTA Figure 6. Tensile stress variation along the top of the desiccating soil for Test 3. be some reduction in the shear stresses at the interface after slip. This behaviour is particularly important for the interface between the base and the soil. Therefore, a residual value of the cohesion after shear slip must be specified which would be a fraction of the intact value. As soil dries, it becomes more brittle, and for soil of relatively low moisture content, say close to the plastic limit, the residual cohesion can reasonably be taken to be a small value. On the other hand, when the soil is more moist, perhaps close to the liquid limit, the residual strength can intuitively thought to be a more significant portion of the intact strength that when the soil was drier and the model would approach that of the elastic-ideal plastic Mohr–Coulomb model. The ratio between the residual and intact shear strength, R, was taken to be 1.0 for soil wetter than the liquid limit, and 0.0 for soil drier than the plastic limit. Thereafter, the variation of R when the moisture content was between the liquid limit and plastic limit was tied to the liquidity index LI where w −wp LI = (9) wl −wp and wl and wp are the liquid limit and the plastic limit (%) of the soil, respectively. The variation of R was modelled as a cos function, R = 0.5−0.5 cos(L I ) (10) 2 cos−1 (L I )  (11) and a cos inverse function R = 1.0− Based on consistency between tests, and time progression of cracks, etc., the cos function was selected as the one which gave the most consistent results. 5. MODELLING RESULTS AND DISCUSSION OF RESULTS 5.1. Replication of laboratory tests A comparison of the experimental test data with the modelling results in terms of number of cracks opening up for each test, the residual height of the soil, the total area of cracks generated, and cracking water content is shown in Table V. Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag 91 NUMERICAL MODELLING OF DESICCATION CRACKING Table V. Summary of results. Test no. Item No. of cracks Residual soil height (mm) Cracking moisture content (%)∗ Average crack opening width (mm) Data type 0 1 2 3 4 5 Experimental Modelling Experimental Modelling Experimental Modelling Experimental Modelling 4 3 8.5 6.5 98.5 98.4 — 22.0 7 10 2.5 2.6 106.2 111.3 5.75 8.65 8.5 8 3 2.7 95.2 111.2 4.75 10.9 11 11 6.5 5.2 108.6 113.3 14.7 19.3 5 3 20 14.8 119.8 111.0 32.2 60.7 5.5 5 13 11.8 105.8 111.4 24.7 30.0 ∗ Average of top and bottom moisture content measurements. When examining the number of cracks observed experimentally and from numerical analysis, it can be seen that there is a fair match, though the number of cracks predicted numerically is smaller for the tests where the thickness to length ratio is highest (Tests 0 and 4). In the opinion of the authors, it may be possible for the methodology outlined to be further refined by better determination of interface properties, more precise moisture content-material properties’ correlations, etc. to obtain a closer match in crack numbers obtained experimentally and predicted numerically. The residual soil height obtained by numerical analysis seems to follow what was obtained experimentally though the residual soil height has been somewhat lower for a few tests. It could be that the stiffness of the soil, which was inputted to the model by specifying Young’s modulus as a function of the moisture content may be underestimated. Conversely, it may be that the suction at the lower moisture contents is somewhat overestimated by the model. Such deviations are not surprising considering that the model attempted to capture variations in suction and stiffness over four orders of magnitude each. The average crack width opening predicted by the model follows the trends of the crack widths observed experimentally though the model seems to currently overestimate crack width openings. As in the case of the residual height of the soil, a lower suction or higher elastic modulus at the residual moisture content is indicated. There is scope for future improvements in the model in this regard. The cracking moisture content is of significance in real-world applications such as clay liners, where cracking may signify the start of diminishing of its functionality. The modelling versus experimental cracking moisture content is shown in Figure 7. It can be seen that the model has been successful in approximately predicting the onset of cracking because all the points are close to the 1:1 line plotted on the figure. The differences between experimental and modelling cracking water content are not excessive when it is considered that the modelling captured the moisture content variation from the start of the test, 127% (134% for Test 4) to the residual moisture content of 10%. A comparison between experimental and modelling development of the number of cracks versus moisture content is shown in Figures 8 and 9 for Tests 3 and 5, respectively. Examining Figure 8, it can be seen that the model correctly predicts that the number of cracks opening up would be 11 and that crack evolution will stabilize at a moisture content of approximately 60%. However, the model seems to predict that the cracks opening up at higher moisture contents appear somewhat earlier than experimentally observed. It can be seen in Figure 9 that the model has predicted the number of cracks and the moisture content at which the crack evolution stabilizes for Test 5 as well. The evolution of number of cracks with moisture content in the model follows the trends observed experimentally. Some curling is predicted by the model with some lift off at the edges and is shown in Figure 10. Significant curling was reported in the laboratory experiments by Nahlawi and Kodikara [10] and a theoretical basis for curling has been given by Kodikara et al. [15]. Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag 92 A. L. AMARASIRI, J. K. KODIKARA AND S. COSTA Figure 7. Moisture contents of initiation of cracking. Figure 8. Crack evolution in Test 3. Figure 9. Crack evolution in Test 5. Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag NUMERICAL MODELLING OF DESICCATION CRACKING 93 Figure 10. Curling of soil blocks in Test 4. Figure 11. Crack progression for Test 2. On the basis of the numerical simulations, the progression of breaking the intact soil into individual blocks is shown schematically in Figure 11 for Test 2 and Figure 12 for Test 5. The model shows the formation of primary, secondary, and tertiary cracks as found in desiccation tests of long moulds, as discussed by Nahlawi and Kodikara [10]. It should be noted that absolute symmetry is not maintained. Indeed particle image velocimetry (PIV) studies have shown that soil desiccating in long moulds tend to move more towards either the left or the right while the block is still intact, i.e. there are sometimes asymmetrical processes involved. One reason could be that once the shear stresses at a side in a desiccating block reaches the shear strength and slippage Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag 94 A. L. AMARASIRI, J. K. KODIKARA AND S. COSTA Figure 12. Crack progression for Test 5. takes place, further slippage may be biased towards the side that has slipped rather than the intact side as desiccation progresses. Therefore, considering the complexity of the actual behaviour of these desiccation tests, the proposed modelling results can be considered to capture some of the essential physics of the problem. 6. SENSITIVITY ANALYSIS A sensitivity analysis of three key input parameters, i.e. number of possible cracks in the soil block, shear strength at the soil–base interface, and tensile strength of the soil were performed. Experiment Test 1 was selected for this study. 6.1. Number of potential cracks The number of vertical joints provided in the soil block used to obtain the results shown in Table V for Test 1 was 49. Analyses were carried out by changing the number of potential cracks to 15, 19, and 69 while leaving all other parameters constant. The number of cracks opening up were 11, 10, 10, and 10 for number of potential cracks 15, 19, 49, and 69, respectively. It was noted while manipulating the number of potential cracks for other runs that the number of cracks opening up would change with increasing joint number and then reach a plateau. Sufficient joints were provided in modelling carried out reported herein to ensure that less than 30% of the joints opened up to form cracks. It is desirable to place a large number of joints as practically possible to allow the soil to crack at points of maximum tensile stress. 6.2. Shear strength at the soil–base interface The shear strength at the soil–base interface for Test 1 was given by Equation (7) with 1 , 2 , and 3 taken as 125, 3.4 kPa and 129%, respectively. Tests were done keeping all other parameters constant and changing the shear strength to 50, 75, 125, and 150% of this value. The number of cracks opening up with the shear strength at 50, 75, 100, 125, and 150% of the value given by Equation (7) were 1, 8, 10, 11, and 12, respectively. As expected, the number of cracks decreased with decreasing shear strength. Thus, it is apparent that the number of cracks opening up is sensitive to the shear strength, with reduction in shear strength from the value used in the model sharply reducing the number of cracks, and the increase in shear strength from the value moderately increasing the number of cracks. Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag NUMERICAL MODELLING OF DESICCATION CRACKING 95 6.3. Tensile strength To ascertain the influence of tensile strength inputted to the number of cracks opening up, modelling was carried out keeping all parameters stationary, but changing the tensile strength to 50, 75, 100%, and 125 and 150% of the value given by Equation (6). The number of cracks opening up were, respectively, 19, 10, 10, 8, 8. First, it is apparent that an increase in tensile strength causes a reduction in the number of cracks opening up as expected. The reduction in shear strength at the interface by 50% caused a reduction in the number of cracks opening up by an order of magnitude, whereas the changes in the tensile strength had caused a change in the same number within a factor of two, which indicates that the model may be more sensitive to shear strength at the soil–base interface than tensile strength. 7. CONCLUSIONS Quantitative prediction of desiccation cracking has not been entirely successful in the past because of the difficulty in capturing the manyfold variation of parameters upon drying. Most studies have dealt with qualitative and statistical analysis of desiccation cracking, or reported numerical analyses where the great changes happen in soil upon drying has not been captured. This paper reports numerical modelling of desiccation cracking of a clay where the required material properties for modelling have been obtained by laboratory experimentation. Laboratory desiccation tests carried out in long rectangular moulds were then modelled using UDEC with the changes in properties incorporated with the aid of the embedded programming language FISH. Attempts were made to replicate six laboratory tests with varying mould construction materials, dimensions, and drying environments with numerical modelling. The modelling was successful in capturing the essence of desiccation crack evolution in terms of number of cracks formed, residual soil height, width of cracks, and crack initiation moisture content, etc., perhaps the first of this kind reported in literature. A sensitivity analysis showed that the number of cracks was not very sensitive to the number of joints provided to the model, sensitive to the tensile strength inputted to the soil, and could be highly sensitive to the shear strength at the soil–base interface within the limits of the changes in the parameters. Herein, the location of potential cracks must be known a priori, but this is not a major limitation because of the ease with which UDEC allows cracks to be included in a block of material in orientations, locations, and spacings of the user’s choice. This methodology can potentially be extended to predict desiccation cracking in earth dams and clay liners, etc., and cracking in other materials with properties varying with time and temperature such as concrete and asphalt. Though this study was carried out in two dimensions, the methodology can be extended to three-dimensional analyses though the geometry and spatial distribution of the cracks may become more irregular and mixed mode fracture may also contribute to the propagation of some of the cracks. The choice of suitable distribution of pre-existing joints and with the ‘right’ joint properties may still require a certain amount of work. ACKNOWLEDGEMENTS This paper reports research carried out with funding from the Australian Research Council Project ID DP 0773861. The kind advice of Dr Mike Coulthard is gratefully acknowledged. REFERENCES 1. Gore A. An Inconvenient Truth: The Planetary Emergency of Global Warming and What We Can Do About It. Rodale Press Inc.: Pennsylvania, U.S.A., 2006. 2. Gould S, Moglia M, Davis P, Burn S. Determining the Economic Lifetime of a Pipeline. AWA Ozwater 2007 Convention & Exhibition. AWA, Sydney Convention & Exhibition Centre, 2007. 3. Konrad JM, Ayad R. An idealized framework for the analysis of cohesive soils undergoing desiccation. Canadian Geotechnical Journal 1997; 34:477–488. Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag 96 A. L. AMARASIRI, J. K. KODIKARA AND S. COSTA 4. Kodikara JK, Choi X. A simplified analytical model for desiccation cracking of clay layers in laboratory tests. In Proceedings of UNSAT 2006 Conference. Miller GA, Zapata CE, Houston SL, Fredlund DG (eds). ASCE Geotechnical Special Publication, Unsaturated Soils, Carefree, AZ, vol. 2, 2006; 2558–2567. 5. Corte A, Higashi A. Experimental Research on Desiccation Cracking in Soil. Research Report, U.S. Army Snow Ice and Permafrost Research Establishment, IL, U.S.A., 1964. 6. Chertkov VY. Modelling cracking stages of saturated soils as they dry and shrink. European Journal of Soil Science 2002; 53:105–118. 7. Morris PH, Graham J, Williams DJ. Cracking in drying soils. Canadian Geotechnical Journal 1992; 29:263–277. 8. Ayad R, Konrad JM, Soulie M. Desiccation of a sensitive clay: application of the model CRACK. Canadian Geotechnical Journal 1997; 34:943–951. 9. Peron H, Hu LB, Laloui L, Hueckel T. Mechanisms of desiccation cracking of soil: validation. In Numerical Models in Geomechanics, Pande G, Pietrusczczac S (eds). Taylor & Francis Group: London, 2007. 10. Nahlawi H, Kodikara JK. Laboratory experiments on desiccation cracking of thin soil layers. Geotechnical and Geological Engineering 2006; 24:1641–1664. 11. Nahlawi H. Behaviour of a reactive soil during desiccation. M. Eng Thesis, Department of Civil Engineering. Monash University, Clayton, Victoria, Australia, 2004. 12. AS 1289.1 Preparation of Disturbed Soil Samples for Testing. Standards Australia, Sydney, 1991. 13. ITASCA. UDEC Users’ Guide. ITASCA: Minnesota, U.S.A., 2004. 14. ITASCA, UDEC Theory and Background. ITASCA: Minnesota, U.S.A., 2004. 15. Kodikara JK, Nahlawi H, Bouazza A. Modelling of curling in desiccating clay. Canadian Geotechnical Journal 2004; 41:560–566. 16. Fredlund DG, Rahardjo H. Soil Mechanics for Unsaturated Soils. Wiley: New York, 1993. 17. Hu LB, Hueckel T, Peron H, Laloui L. Desiccation shrinkage of unconstrained soil in the saturated phase. Proceedings of the First European Conference on Unsaturated Soils (E-UNSAT 2008), Durham, U.K., 2–4 July 2008. Copyright q 2010 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2011; 35:82–96 DOI: 10.1002/nag