Numerical simulation of the effect of laponite on saturated sandy ground in re-liquefaction events


Laponite is a nanomaterial that can transform the pore water into a viscous gel and continues to harden with time and can recover its strength after liquefaction. Initially, laponite suspension has very low viscosity, so it could be used in passive site remediation to improve the ground in highly populated areas. In this study, a constitutive model able to reproduce the behaviour of both granular and cohesive materials under different loading conditions is calibrated and validated using several laboratory experimental results. The calibrated parameters are then used to simulate the impact of a sequence of two earthquakes on a residential building with a shallow foundation on a loose saturated soil deposit containing sand that was treated with 1% laponite. This study attempts to provide insights on the level of improvement that laponite could provide as a liquefaction mitigation method in real applications and to show practitioners the potential of using this material for remediation purposes.


What is laponite?

Laponite is a non-toxic entirely synthetic layered silicate with a chemical composition equivalent to that of naturally occurring smectite clay minerals (similar shape and structure to Na-montmorillonite, but one order of magnitude smaller). Unlike natural clays, laponite has no impurities which make its suspension properties more controllable and stable. For instance, it is possible to prepare it with very low initial viscosity when injected and it has high yield shear stress with time. Compared with bentonite, laponite’s plasticity index has been reported to be 1200% (BYK Additives & Instruments 2014), twice that of bentonite. Its particle size is 5-10 times smaller than bentonite, and its initial viscosity and density are similar to water; these features favour laponite over bentonite for use in mitigating liquefaction. Besides, the use of highly plastic clays has been shown to increase the cyclic strength of sands (e.g. Ishihara 1993; Ishihara and Koseki 1989).

Laponite particles are disk-shaped, which makes its edges naturally positively charged, and its surface negatively charged. When laponite is in dry-powder form, its platelets arrange themselves electrostatically through sodium bonds (Figure 1a). However, when dispersed in water, the sodium ions are released, and the particles repel each other from their mutual positive charges. At this stage, the dispersion displays low viscosity and Newtonian behaviour (BYK Additives & Instruments 2014). When the equilibrium of the ions is reached, laponite particles start re-arranging themselves due to electroactive forces (the positive edge of some particles attracts the negative surface of other particles), and the formation of a “house of cards” structure starts taking place (Figure 1b). This change in internal structure is responsible for a change in the rheological behaviour of the fluid, i.e. a thixotropic gel. One of the key features of this behaviour relates with its yield shear stress, which increases with time. The gel does not flow unless a shear stress larger than this threshold is reached (); at this stage, the network breaks and the material flows, decreasing its viscosity (Figures 1c and 1d), but maintaining solid-like behaviour, which means that the elastic components are predominant over the viscous components. When the shear stress is suppressed, the network recovers and the resistance to flow increases, which is an indicator that laponite could recover after an earthquake and continue to harden with time.

a) Dispersed particle (BYK Additives & Instruments 2014) b) Scheme of house of cards formation c) Thixotropic behaviour (modified from Barnes 1997) d) Viscosity of a thixotropic material

Figure 1: Laponite particle dispersion and gel formation

This study

This study presents a numerical model using the finite element method to simulate the effects of two successive earthquakes on a loose saturated soil deposit containing sand with and without treatment of 1% laponite. The purpose of this numerical simulation is to provide an insight on the level of improvement that laponite could provide as an option to mitigate liquefaction after a series of earthquakes, and to show practitioners the potential of using this material for remediation purposes.

The application study is inspired by the 2010 Darfield and 2011 Christchurch Earthquakes, where the two earthquakes struck the same area 5 months apart (Figure 2). The main earthquake (Mw 7.1) occurred on 3 September 2010, 40 km west of Christchurch Central Business District (CBD) at 10 km depth. Liquefaction-induced damages were observed near Avon River (Figure 2). The aftershock (Mw 6.3) occurred on 21 February 2011, 10 km south of the CBD at a depth of 5 km. The aftershock, being shallower and closer to the CBD, caused even more damage and the extent of liquefaction was even larger (Figure 2).

Among the different types of structures, residential houses built on soft soils suffered more damage during the two earthquakes. Therefore, this numerical analysis is focused on the performance of a shallow foundation on a loose saturated sandy deposit subjected to a sequence of two ground motions (the Darfield/Christchurch earthquake sequence). The soil is represented by a multi-mechanism constitutive model capable of reproducing the behaviour of a wide range of soils under different loading conditions (drained, undrained, cyclic, monotonic). This constitutive model is used to represent the behaviour of the granular material used in this study (Mercer sand) and the addition of 1% laponite.

a) New Zealand map b) Satellite image showing epicentres of main earthquakes c) Distribution of liquefaction induced damage in Christchurch during the main earthquakes (modified from Orense et al. 2011)

Figure 2: Location of epicentres of main earthquakes and distribution of liquefaction around Christchurch


The model consists of a single-story structure on a shallow foundation on top of a 30 m thick layer of liquefiable soil deposit and underlain by bedrock (see Figure 3). The structure represents a one-story house, where the mass is distributed in the upper beam only, while the columns provide the strength and stiffness of the structure. The model is two-dimensional under plane strain condition.

At the bottom of the model, a 2.5 m thick elastic bedrock and a paraxial element was added to simulate the deformable unbounded bedrock and to prevent the reflection of seismic waves (Modaressi 1987). Only shear waves propagating vertically are considered. The computations are carried out in the time domain. The model was analysed in three steps: (1) Initialisation; (2) 2010 Earthquake; and (3) 2011 Earthquake. Figure 3 shows the boundary conditions used in each step. In the initialisation step, the model was consolidated under self-weight. In the next steps, the input motion was introduced through the paraxial element.

Initialisation step Earthquake loading (Steps 2 and 3)
a) Schematic of the model (not to scale) b) Boundary conditions

Figure 3: Model and boundary conditions

The soil deposit was divided into two layers, as depicted in Figure 4, allowing the possibility to use different constitutive models for the top 6 m soil and for the rest of the deposit for comparison purposes. In total, three models were analysed, in which Soils 1 and 2 (Figure 4) were represented by different combinations of pure sand and sand with 1% laponite, as summarised in Table 1. Model 1 is the control case, in which the whole soil deposit was modelled as pure sand (calibrated using the experimental results for Mercer sand).

The soil was modelled using the Hujeux constitutive model (Aubry et al. 1986; Hujeux 1985) implemented in the open source computational tool GEFDyn (Aubry et al. 1986). This is a multi-mechanism constitutive model that can reproduce the behaviour of granular and fine materials under static and dynamic loading. Table 2 shows a summary of the calibrated properties for each material.

Figure 4: Mesh for the different models

Table 1: Summary of models

Depth (m)

(top – bottom)

Layer ID Model 1:


Model 2


Model 3


0 – 6 Soil 1 Mercer sand 1% Laponite 1% Laponite
6 – 30 Soil 2 Mercer sand 1% Laponite Mercer sand

Table 2: Summary of the calibrated parameters used for the different materials

Pure Sand 1% Laponite Phase 1 1% Laponite Phase 2
Elastic Parameters Kref (MPa) 15.8 23 26.3
Gref (MPa) 9 13 15
ne 0.5 0.57 0.57
pref (MPa) 1 1 1
Loading function parameters φ’ (°) 33 19 19
ψ (°) 32 24.5 24.5
β 38 38 38
d 3.5 3.5 3.5
b 0.15 0.337 0.337
Hardening Parameters amon 0.02 0.048 0.048
acyc 0.001 0.001 0.001
cmon 0.00001 0.001 0.001
ccyc 0.004 0.004 0.004
αψ 0.4 0.5 0.5
Domain radii of behaviour rela 0.005 0.009 0.009
rhys 0.08 0.186 0.228
rmob 0.5 0.5 0.54
riso 0.01 0.001 0.001
xM 0.5 1.6 2

The materials were calibrated using the results of the cyclic simple shear tests, in which laponite-treated specimens were subjected to a series of two consecutive shear loadings applied until the samples reached liquefaction (defined as double amplitude shear strain ). The calibrated parameters were validated using the results of the other element testing, such as consolidation, drained and undrained simple shear tests; further details are provided by Tobar (2019). Figure 5 shows the results of the calibration for cyclic simple shear tests, in terms of the shear strain curves versus the number of cycles to liquefaction for the simulated and the experimental (in black) results at different levels of cyclic stress ratios . All the tests were performed at an initial vertical effective stress of 100 kPa.

(a) Sand (b) Laponite Phase 1 (c) Laponite Phase 2
















Figure 5: Calibration of the parameters based on CSST

(a) 2010 Earthquake (b) 2011 Earthquake

Figure 6: Records of earthquake motion

The selected motions were those recorded in one of the strong motion stations in Christchurch (RHSC station in Figure 2). As the earthquake motions were recorded at the surface of the soil deposit, it was necessary to de-convolute them to obtain the equivalent motions at a depth of 30 m. For this purpose, a frequency domain analysis with equivalent linear assumption and visco-elasticity was performed, while the degradation of shear modulus and damping were also considered. Additionally, each component was zero-padded to increase the length of the earthquake record (see Figure 6), in order to allow the excess pore water pressure to drain after the earthquakes.



Figure 7 shows the evolution of the vertical displacement 10 cm below the foundation, in which the time scale has been split into two to show the co-seismic settlement in more detailed scale, and the settlements after the earthquake have been presented in a larger scale to show the overall behaviour.

During the initialisation phase (Figure 7a), the soil started with an initial stress state very similar to the stress conditions due to self-weight. The construction of the structure started at 4000 s, and the three models experienced similar settlements due to the weight of the structure, which is between 6.1 − 6.6 cm. However, the models with laponite (1Lap and 1Lap6m) experienced the highest settlement at Dv = 6.6 cm.

Figure 7: Vertical displacement under the foundation

Regarding co-seismic settlements, the soil deposits experienced more deformation during the 2010 Earthquake than the 2011 Earthquake. That is, the vertical displacement during the effective duration of the first earthquake was between 7 − 11 cm, while during the second earthquake, the deformation was almost negligible because the sand particles had already densified considerably following the previous quake, and had little room for more deformation, so the co-seismic settlement varied between 0.1 − 0.4 cm. Most of the deformation occurred after the motions have finished due to the dissipation of the excess pore water pressure, which led to a total settlement of the soil deposit between 23 − 29 cm at the end of the simulation.

Excess pore water pressure ratio

Figures 8a and 8b show the vertical profiles of the excess pore water pressure ratio under the foundation for the 2010 and 2011 Earthquakes, respectively. Moreover, Figures 8c and 8d show the vertical profiles of the excess pore water pressure ratio in the free field (24 m away from the centre of the structure) for the two earthquakes. In these figures, the lines represent the envelope of all the excess of pore water pressure ratios for each model. In both earthquakes and locations (i.e. under the foundation and in free field), the case with pure sand developed the highest level of excess pore water pressure ratio, reaching a maximum value close to ru = 1 in free field during the 2010 Earthquake. On the other hand, the case with 1% laponite developed a maximum ru = 0.73 in the same conditions.

Under the foundation Free field
  1. 2010
  1. 2011
  1. 2010
  1. 2011

Figure 8: Pore water pressure ratio profiles under the foundation and free field

The shape of the ru distribution is similar among the models. In the free field, the peak value occurs near the surface, while under the foundation, it occurs at a depth of about m. Below this depth, the excess pore water pressure is minor, being less than 0.15 at a depth of m, except for the case with 1% laponite where the increase in pore water pressure is distributed in a wider zone, affecting deeper layers. Additionally, from Figures 8a and 8c, it is possible to observe that at the end of the motions, the excess pore water pressure had dissipated considerably in all the models before the 2011 earthquake started.

The development of the excess pore water pressure during 2011 Earthquake was less compared with that during the 2010 Earthquake, but a similar trend was obtained. Thus, the case with pure sand exhibited higher excess pore water pressure, while the case with 1% Laponite in the whole soil deposit performed better.

Table 3: Summary of the results of the simulation

Model 1: Sand Model 2: 1Lap Model 3: 1Lap6m
(2010EQ) (cm) 22 20 21
(2011EQ) (cm) 0.97 0.58 0.64
Total Settlement (cm) 29 27 28
Max (FF, 2010EQ) (-) 1 0.6 0.7
Max (FF, 2011EQ) (-) 0.62 0.24 0.20
PGA (FF, 2010EQ) (g) 0.23 0.16 0.22
PGA (FF, 2011EQ) (g) 0.37 0.32 0.28
PGD (FF, 2010EQ) (cm) 46 43 45
PGD (FF,2011EQ) (cm) 30 26 29

Table 3 shows a summary with the maximum values obtained for the different models analysed: (i) pure sand; (ii) sand with 1% laponite; and (iii) sand with 1% laponite only in the top 6 m. In the table, the notation FF stands for free field; and 2010EQ and 2011EQ denote the 2010 and 2011 Earthquakes, respectively. Moreover, for each variable, the cells corresponding to the best performing model are highlighted in grey. In general, the control case (pure sand) performed worse than those models with 1% laponite. The model with 1% laponite throughout the depth had better results in more categories, but it has the excess of pore water pressure distributed over a broader zone, reaching deeper than the rest of the models.


This study presented a numerical simulation inspired by the 2010- 2011 Canterbury earthquake sequence that caused great damage due to liquefaction, especially in residential houses. For this reason, a model of a structure with a shallow foundation on liquefiable ground was analysed. The soil deposit was represented using a multi-mechanism constitutive model that was calibrated using experimental data.

The results of the simulation show the potential of laponite as a possible mitigation option on large-scale applications. Laponite could be injected into the groundwater and, based on the results presented, it would be more beneficial to cover as deep as possible, because laponite could reduce the peak excess pore water pressure which could also lead into a smaller shear deformation.

It is important to note that the results obtained in the numerical simulation are a direct consequence of the quality of the calibration parameters selected to represent each of the materials. The set of parameters used in this study were considered acceptable because they allowed reproducing the macro-response of the materials under both cyclic and monotonic stress paths. However, the viscous behaviour of the laponite suspension or its capacity to harden with time or to recover after the removal of shear loading was not modelled. To date, there is no constitutive model able to capture its behaviour and its interaction with the pore matrix.


Aubry, D., Chouvet, D., Modaressi, A. and Modaressi, H. 1986. “GEFDYN : Logiciel d’analyse de comportement mécanique des sols par eléments finis avec prise en compte du couplage sol-eau-air.” Manuel Scientifique, Ecole Centrale Paris, LMSS-Mat.

Barnes, H. A. 1997. “Thixotropy – a review.” Journal of Non-Newtonian Fluid Mechanics 70(97):1–33.

BYK Additives & Instruments. 2014. “Laponite: performance additives.” Technical Information B-RI 2 1.

Hujeux, JC. 1985. “A constitutive law for the cyclic loading of the soils.” Génie Parasismique, edited by V. Davidovici and Al. Presses ENPC, 287–302.

Ishihara, K. 1993. “Liquefaction and flow failure during earthquakes.” Geotechnique 43(3):351–415.

Ishihara, K. and Koseki, J. 1989. “Cyclic shear strength of fines-containing sands.” Proc, 12th Int. Conference on Soil Mechanics and Foundation Engineering. Vol. 27. Rio de Janeiro, 101–106.

Modaressi, H. 1987. “Numerical modelling of wave propagation in an elastic porous media.” PhD Thesis, Ecole Centrale de Paris.

Orense, R. P., Kiyota, T., Yamada, S., Cubrinovski, M., Hosono, Y., Okamura, M., and Yasuda, S. 2011. “Comparison of liquefaction features observed during the 2010 and 2011 Canterbury earthquakes.” Seismological Research Letters 82(6):905–18.

Tobar, G.P. 2019. “Mitigation of soil liquefaction with two sustainable materials: Biochar and Laponite”, PhD Thesis, University of Auckland.

Tags : #Laponite#re-liquefaction

G.P. Tobar, R.P. Orense
NZGS Symposium>21st NZGS Symposium
Conference, Proceedings

Leave a Reply