Numerical modeling of drying and consolidation of fine sediments and tailings

Session 10 - Desiccation Dewatering

NUMERICAL, iVlODELiWG OF DRYING AND CONSOLIDATIOW OF FINE

SEDIMENTS AND TAILINGS

Jochem van der Meulen\ Frits van Tol^ ^ Leon van Paassen\ Timo Heimovaara^ Delft University oflechnoiogy, Delft, Netherlands

^Deltares, Delft, Netherlands

A B S T R A C T

The extraction and processing of many mineral ores result in the generation of large volumes of fine-grained residue or tailings. These fine sediments are deposited as a slurry with very high water contents and lose water after deposition due to self-weight consolidation. When the surface is exposed to the atmosphere they dry due to evaporation. In thin lifts deposition both processes, consolidation and drying take place simultaneously. The understanding and a quantitative description of both processes is important for the prediction of the dewatering process of fine sediments and tailings. A numerical model, which is able to simulate one-dimensional vertical flow due to consolidation and drying by evaporation of fine sediments, is implemented in a IViATLAB code.

The proper functioning of the model was verified for several combinations of boundary conditions, including a closed boundary, a constant flux boundary, a constant pressure boundary for both the top or the bottom and deposition of a fresh iayer on top of a dried one. Subsequently the model was validated using data from experimental research on fine oil sand tailings. The model showed realistic outcomes for most boundary conditions and good correspondence with the observed behaviour in the drying and consolidation experiments. The paper describes the model, and the verification and validation simulations.

INTRODUCTION

Canada accommodates an enormous oil supply from oil sands. The extraction of oil from the sand creates a significant amount of slurry. The sludge, which consists of approximately 80 to 90% water, has been stored in temporarily basins. Dewatering by self-weight consolidation takes up decades. In order to accelerate the dewatering process, the use of so-called mud farming has been proposed. Mud farming also called sub-aerial drying represents the process where the sludge is

deposited in thin layers, in which every layer loses water mainly by evaporation. Thin lift deposition is a potential solution, and involves deposition in thin lifts that are subsequently allowed to dry and gain strength. When the layer has lost a sufficient amount of water, the next layer can be deposited on top resulting in consolidation of the layers below. To determine the thickness of the applied layers and the required time to dewater the deposited layers, a model that incorporates the drying and consolidation process is required.

Both the consolidation process as the drying process of the slurry involve large volume reduction. Large strain consolidation has been described by the finite strain consolidation (Gibson et al., (1967); Schiffmann et a!., (1988), to be able to model the self-weight consolidation of fine sediments. The drying process has been described by several authors mainly from the agriculture-related disciplines e.g. Reniersce (1983) who modeled the ripening of soils gained from the sea. More recently this process is described by, Wilson et al. (1994); Wilson et al., (1997), to model the atmospheric drying process of soils. The theory of the combined process of sub-areal drying, self-weight and overburden consolidation has been described by Kim et al, (1992).

This paper describes the implementation, verification and validation of this combined process in a MATLAB based computer code. For the validation simulations of experiments of the drying and consolidation process of fine oil sand tailings were used.

THE MODEL

In order to understand the calculation process of the model the basic theory as described by Kim et al, (1992) is repeated in this paragraph.

The model is based on the following correlations: • Water ratio and void ratio (shnnkage curve); • Water ratio and matrix potential (pF-curve);

a Water ratio and tiydrauiic conductivity (permeability curve).

The model basically is a flow model in which Darcy's law governs the process. After dividing a slurry layer in a finite number of sub layers the flux (i;) over all modelling layers with thickness ( A z ) can be calculated which gives the change in water ratio (Ai9) per unit of time ( A t ) . The loss of water per layer is equal to the change in water ratio per layer:

n = overburden component [kPa]

M

At

Av

[1] By using the shrinkage curve, see

Figure 1, the correlation between the water ratio and the void ratio the loss of water is translated to a loss of volume and in a 1D case to settlement.

Figure 1. The Correlation Between The Water Content And Void Ratio (Taken From Kim Et Al, 1992).

The fluxes between all modelling layers are being calculated with the use of Darcy's law;

V - -k-Vi

V i =

d(t)

d7

[2]

V = water flux [cm/day]

k = vertical permeability [cm/day]

Vt = gradiënt in water pressure [-]

dcj) = difference in water pressure [kPa] dz = thickness o f t h e modeling layer [cm]

The water pressure is built from 3 components (Philip, 1969), see Equation 3:

(^-y^-z + ip + D. (p = water pressure [kPa] z = gravitational component [m]

i/) = matric potential [kPa]

[3]

Yw specific weight of water [kN/m ]

In the original paper the dimensions of the components in equation 3 are given in cm (water pressure). The component of overburden is being determined by the following formulas:

P ( z ) = P ( 0 ) + j 7 , d z

Yb = _{1 + e}

[4]

[5]

[6] P(z) = total stress at depth z [kPa]

P(0) = external load [kPa]

Yb = bulk density [kN/m^]

dz = thickness of modelling layers [cm]

i9 = water ratio ( ^ ) e = void ratio ( ^ )

= specific weight of the grains [kN/m^] Kim et al, 1992 use the water retention curve as general relation between the water ratio and stress. Figure 2 presents a water retention curve (Van Genuchten, 1980), the stress strain relation that was found by Yao et al(2011) for thickened tailings, in consolidation experiments. For the saturated part they match quite well.

5 10"

t \

pF-cutve (Van Genuchten, 1980 (1)) Stress-strain relation linear (Yao et al., 2011 (2)) \

\

Water ratio, V, (1) - Void rato, V^/V^ (2) Figure 2. Correlation Between

Stress And Void Ratio.

Effective

The third relation is a relation between the permeability and the water ratio. As the soil becomes more compact, the permeability decreases. The permeability decreases even more when the pores become unsaturated. A similar relation is given in Figure 3.

waterraUo, V.,yV^

Figure 3. Permeability Which Has Been Used By Kim In Order To Model Ripening Of Marine Clay (Kim Et Al, 1991).

For one-dimensional flow the total flow equation is as written in Equation 7:

[7] For saturated compaction the decrease in water

fde equals the decrease in void ratio I — =

1 en—r = 0 , the general flow equation can be

simplified by equation 8.

[8] In case of large strain systems, usually one makes use of a material coordinate system (Smiles and Rosenthal, 1968; Philip, 1968). In these kind of systems the regular coordinate system ( d z ) will be written to a so called material coordinate system ( d m ) where the layer thickness will be kept at a constant value.

d z

dm = ——

1

dm = thickness of modeling layers in the material coordinate system [cm]

dz = thickness of the modeling layers in a

regular coordinate system [cm] The differential equation now becomes:

at

dv

dm [9]

Kim et al., 1992 reduce the general equation according one dimensional groundwater flow for non-steady state situations to Equation 10, in case without surface loading (P(0) = 0):

3J. = ^ \ K * ( ^ + SFl + SF2 • P)] [10] at dm L \dm dm/i With de

SFl = (l + e)-id + ys) —

. F 2 = ^ J > + r . ) d

K* = — , 5F1

In this equation m represents coordinate (positive upward)

the component of gravitational and overburden potential term and SF2 is the component of contribution by overburden potential arising from the difference in volume shrinkage (residual zone). In these differential equations only three parameters are unknown: e.ipenK. All three parameters can be written as a function of the water ratio (i9). Figure 1,2, and 3. These correlations depend on the type of soil and can be determined by lab testing.

BOUNDARY CONDITIONS

For both top and bottom of the model a boundary condition needs to be specified.

Top boundary condition

According to the top boundary, there are three different situations:

1. Under water consolidation (self-weight consolidation).

2. Drying by evaporation; 3. Precipitation.

Situation 1 (self-weigiit consolidation)

With regard to the first situation, the density of the top modelling layer does not change. There is no compaction, as there is no surface load or atmospheric suction. Zero compaction means no change in flux or matric potential xp =^ 0. For an already consolidated layer, deposed under water,

zero-flux-boundary-condition is used. In other terms; the flux over the bottom of the top modelling iayer is equal to the flux over the top of the top modelling layer. The delta flux over the top modelling layer is therefore equal to zero, see Figure 4.

Delta flux equals 0 and therefore no compaction occurs

Figure 4. In C a s e Of Underwater Consolidation, Delta Flux Of The Top Modelling Layer Needs To Be Equal To Zero.

Situation 2 (evaporation)

Evaporation can be modelled by making use of both a flux and a matric potential (Rijniersce, 1983). By setting the flux equal to the open water evaporation and setting the matric potential to a value corresponding to the relative humidity of the atmosphere, the outflow of water at the top boundary will never exceed the open water evaporation and also will never dry out more thqn the given matric potential corresponding to the relative humidity, see Equation 13.

[13]

Ri^atm = evaporation from the soil [cm/day]

E = maximum open water evaporation [cm/day]

i>atm = matric potential which corresponds to the atmospheric relative humidity [kPa]

rpt = matric potential which corresponds to the top

modelling layer [kPa]

= vertical permeability of the soil [cm/day]

dm = thickness o f t h e top modelling layer [cm]

Situation 3 (precipitation)

Precipitation can be modelled with the use of a potential or a flux.

Bottom boundary condition

According to the bottom boundary, there are three different situations:

1. The bottom boundary is closed off.

2. An open boundary which drains towards deeper layers without any constrain.

3. A boundary characterized by a fix seepage.

Situation 1 (zero flux)

The situation where no drainage or seepage occurs can be modelled by setting the bottom boundary to a flux equal to zero.

Situation 2 (drainage)

In order to model drainage, one needs to make use of a matric potential. Straight after deposing the sludge, water will drain at the bottom. This results into a water pressure equal to zero. This corresponds to a situation where the total pressure equals the effective pressure. The bulk weight of the upper lying sludge now rests on the grains at the bottom.

As long as the sludge remains saturated, the matric potential is equal to the effective stress. In order to set the right boundary condition this needs to be realized.

ttom ~ _1+e [14]

^bottom- matric potential at the bottom boundary

[kPa]

z = thickness of the soil layer [m]

Situation 3 (seepage)

In order to model seepage, a flux is used.

qbottom-^N = seepage [15]

Qbottom^N = flux over the bottom boundary [cm/dag]

VERIFICATION

In order to verify the model, multiple simulations have been made. The results have been verified by checking the mass balance. For the cases without evaporation (only self-weight consolidation) the expected effective stress has been compared to the calculated matric potential. The simulations are characterised by the boundary conditions in Figure 5.

CASE 2 CASE 3

1. Evaporation at the top boundary and outfiow at the bottom boundary.

2. Evaporation at the top boundary and zero flux at the bottom boundary.

3. Stagnant water at the top boundary and zero flux at the bottom boundary.

4. Similar to case 2, but a second fresh layer is being deposed on an partly dried layer.

Table 1 summarizes the used boundary conditions.

Figure 5. Verification C a s e s .

Table 1. Boundary Conditions for Verification C a s e s . Case Initial Water

ratio

Height Top'' Bottom

Flux Potential Flux Potential |V„A/,1 [cm] [cm/month] [cm] [cm/day] [cm]

1 5.8 170 5 -10 000

-

2 5.8 170 5 -10 000 0 ~

3 5.8 170 X 0.0 0

-4 5.8^> 250 5 -10 000

-

For the evaporation two boundary conditions are being used: a matric potential corresponding to the relative humidity and fhe maximum open wafer evaporation. The model works with fhe smallest of fhe two.

A fresh layer of sludge(170 cm) has been modelled on an already ripened layer (80 cm). The water ratio ofthe already ripened layer has been determined according case 1. It is assumed that the 2"''layer has been deposed after a period of 100 days.

C a s e 1

In case 1 the sludge is deposed on a dry permeable ground surface. It has been assumed that the sludge is able to drain without any resistance at the bottom. The calculation results are presented in Figures 6, 7, 8 and 9.

The sludge dewaters in both upward and downward directions in the first time steps (see Figure 9), resulting in a decrease of the water ratio simultaneously at the top and the bottom of the layer. Finally the clay has a water ratio equal to 1,5 and a matric potential of -1000 kPa. After a period of 4 years the initial layer of 170 cm is compressed to about 60 cm and the flux converges to zero. The vertical axes of Figures 6, 8 and 9 is given in calculation layers of 1 cm solids. To transfer to vertical axes in cm these values should be multiplied by 1+e, which gives initially 170 cm but changes with every calculation step in a different way over the depth.

.S -15 Figure 6. WATER RATIO t = 0 year t = 0.5 year t = 1 year I = 2 years t = 4 years 2 3 4 Water ratio, [V, A/ ]

Water Ratio Over The Height Of The Sample.

S E T T L E M E N T 160 140 120 u , 100 15 20 25 30 35 40 45 50 Time, [months] Figure 7. Settlement In Time.

MATRIC P O T E N T I A L c -10 '•E to o O ro -15 D) cr 2 Q. I = 0.5 year 1 = 0 year t = 1 year 1 = 2 years t = 4 years -1000 -800 -600 -400 -200 0

Matric Potential, [kPa]

Figure 8. Matric Potentiai Over The Height Of The Sample.

Case 2

Sludge is deposed on an impermeable layer. There is evaporation at the top boundary and no flow over the bottom boundary. The height of the sample is 170 cm with an initial water ratio of 5,8. The calculated matric potential and fluxes are presented in Figures 10, 11 and 12, The settlements are nearly equal to Figures 7 as the effect of bottom drainage on the settlements is relatively small.

However the difference between an open and closed bottom can clearly be seen in water ratio. Figures 6. In case 2 water evaporates and simultaneously the sludge is consolidates under it's own weight but the water ratio does not decrease as fast as in the case with the drained bottom. The sludge dewaters only upward direction. Figure 12.

F L U X E S t ~ 4 years t = 1 year ( • : = 0.6 year t = 2 years t = 0 year - 2 5 ¬ -0.1 0 0.05 0.1 Flux, [cm/dag] 0.15 0.2

Figure 9. _{Flux [Cm/Day] Over The Height Of} The Sample. W A T E R RATIO

I

Figure 10. Water Ratio Over The Height Of The Sample.

M A T R I C P O T E N T I A L W A T E R R A T I O I -10 !0 O O ra -15 ' C D c 2 O) -20 0) t = 4 years t = 0.5 year . t = O year t = 1 year t = 2 y e a r s -1000 -800 -600 ^ 0 0 -200

Matric Potential, [kPa]

-/ -- t = 20 years / /' ' t = 5 years / -t ~ 2 years / t = 1 year / .* /' / • t = 0 year '. / 3.6 4 4.5 6 Water ratio, [V N 1

F i g u r e n . Matric Potential Over The Height Of pigurelS. Water Ratio Over The Height Of The

The Sample. Sample.

!^ -10

F L U X E S

t = 4 years t = 1 year t = 0.5 year

- t = 2 years ,' 3 year • 1 • r • O 0.06 0.1 Flux, [cm/day]

Figure 12. Flux Over The Height Of The Sample.

S E T T L E I V I E N T

10 12 14

Time [years] Figure 14. Settlement In Time. C a s e 3

Case 3 is a situation of self-weight consolidation under water of a layer of 170 cm of slurry with the same initial water ratio of 5.8.

After a period of approximately 5 years (60 time steps) the flux is nearly equal to zero over the height o f t h e sample and the end of settlement has practically been reached. The matric potential after 20 years equals app. 4 kPa at the bottom, which is equal to the effective stress. As 12.21 kN/m^ a ; = (12,21-10) 1,7 = 3.75 kPa

MATRIC P O T E N T I A L t = 20 years t = 5years t = 2 years /• t = 1 year . . •'' 1 = 0 year -3 -2 -1

Matric Pctenlial, [l<Pa]

Figure 15. Matric Potential Over The Height Of The Sample. F L U X E S £2 f -5 £ -10 8 ro -15 CD

I

t

Figure 16. Flux Over The Height Of The Sample.

Case 4

Case 4 is represented by a situation where a fresh slurry layer is imposed on a drying layer (at t = 100 days). The top boundary is characterized by evaporation and the sample is able to drain at the bottom. The drying layer has been modelled according to case 1. The output has been generated at t = 100 days. Figures 17, 18 19 and 20 presents the results

W A T E R RATIO E -35 ^ 4 5 -t = 1 year f = 2 years t = 4 years t = 0 days , ' t = 1 month ^ 1 = 2 2.5 3 3.5 A 4.5 5 6.5 6 Water ratio, [V N ] _{^ w}

Figure 17. Water Ratio Over The Height Of The Sample.

S E T T L E M E N T

50

0 5 10 16 20 25 30 36 40 45

Time [months] Figure 18. Settlement In Time.

M A T R I C P O T E N T I A L F ° ra -5 Vi M -16 g -20 -8 -25 O ro -30 D) c 2 -35 D) ^ 0 t^ 5 t = 4 years t = 1 year t = 2 years t = 0 days t = 1 month OJ -50 -1O00 -800 .600 ^ 0 0 -200 0

Matric Potential, [kPa]

Figure 19. Matric Potential Over The Height Of The Sample.

F L U X E S = 2 years t = 4 years I = 1 year ^ t = O d a y s -0.1 -0.05 0 0.05 0.1 0.15 0.2 Flux, [ c m / d a y ]

Figure 20. Flux Over The Height Of The Sample.

Figure 17 shows that the water ratio at the top of the bottom layer increases in the first time steps. The water from the imposed layer drains into the bottom layer. At the bottom of the lower layer the water ratio decreases; caused by the overburden of the imposed layer. Due to the high suction pressure in the top of the present layer, the evaporation becomes less. When the new slurry layer is being imposed, the matric potential at the top of the lower layer immediately becomes less negative due to the infiltration form the new layer.

Figures 21 up to 24 present the calculation results. The comparison of the settlements is shown in Figure 25.

The dots in Figure 25 represent the measured results and the lines represent the calculated results. The measured and calculated profile are very similar. The calculated settlement profile which makes use of the pF-curve and permeabilty curve come to a lesser degree of similarity than the fitted profile. The water ratio over the depth of the profile has been measured at the end of the test at t=100 days. In Figure 26 the measured and calculated waterratio at t=100 days over the depth ofthe sample are plotted

W A T E R RATIO • -0.6 -1 . - 1 5 -2 •2-5 -3 ' -3.5 I -4 -i.s • -5 Water ratio, [V^,

VALIDATION

Besides verification, the model needed to be validated. Not all cases from the previous chapter have been physically modelled. Weijermars (2011) performed consolidation and sedimentation tests on thickened tailings from the Albian Sands Muskeg River Mine in Fort McMurray. In order to monitor the self-weight consolidation, a constant water head has been applied. The boundary condition are presented in table 2

Table 2. Boundary Conditions for Validation

C a s e W a t e r ratio H e i g h t T o p " Flux Potential B o t t o m Flux Potential E n d

I V J V . l [cml Icm/month] [cm] fcnVdayl [cm] [days]

V a l . 6.3 34 - 0.0 0 - 100

^' For the evaporation tv;o boundary conditions are being used: a matric potential corresponding to the relative humidity and the m a x i m u m allowable open w a t e r evaporation. T h e model works with the smaller of the two.

Figure 21. Water Ratio Over The Height Of The Sample.

S E T T L E M E N T

0 10 20 30 40 50 60 70 80 90 100

Time [days] Figure 22. Settlement In Time

M A T R I C P O T E N T I A L > . -1.6 -i 1 - / ^ ) 1 I = 100 days -- _{t = 60 days} • -

_/

Matric Potential, [kPa]

Figure 23. Matric Potential Over The Height Of The Sample. F L U X E S >- -IS -2.6 £ -4.6 -0.1 1 = 0 days 0 0.1 0.2 0.3 Flux, [cm/dag] 0.4

Figure 24. Flux Over The Height Of The Sample.

Both settlement profiles correspond reasonably well. The difference in both profiles can be explained by the fact that the measured profile height at t=100 was slightly less than the calculated profile height at t= 100. Therefore the measured water ratio is a somewhat less than the calculated water ratio. A possible explanation is that the chosen permeability for the calculation was too low, and also the relation was simplified (linear at log-scale). It can also be seen that the water ratio at the top and bottom of the sample correspond exactly, this means the water retention curve curve (pF-curve) was estimated quite well.

S E T T L E M E N T

Elapsed Time [days] • Measured settlement, w=211%

Matlab results, w/=211%

Figure 25. Calculated And Measured Settlement Of Under Water Consolidation Of Thickened Tailings. W A T E R R A T I O (t=100) • y / • / water ratio (-) • M e a s u r e d , W = 2 1 1 % • M a t l a b r e s u l t a a t , w = 2 1 1 %

Figure 26. The Measured And Calculated Water Ratio.

CONCLUSION

A numerical model, which is able to simulate one-dimensional vertical flow due to consolidation and drying by evaporation of fine sediments, was successfully implemented in a MATLAB code. The model calculates per time step the settlement, water ratio, the matrix potential and the flux over the depth of the sample. By a number of verification calculations it was shown that model gives consistent results. Validations using the results of laboratory tests could be simulated rather well.

During the development of the model some interesting findings were made: / the pF-curve can be interpreted in the same way as Terzaghi's stress-strain relation. For the fully saturated part the matric potential equals the effective stress. //This makes it possible to apply a matric potential at the bottom boundary of the sample, which equals the effective stress (case 1 and 3). This same applies for the upper boundary condition. If there is no load or evaporation, the matric potential at the top boundary can be set at zero. /// By the application of an overburden component, the pore water will flow upwards as a consequence of self-weight consolidation.

With the model it was possible to simulate successfully subsequent deposited layers of slurry which is promising for the prediction of the drying process in thin lifts. For the time being the model uses only one stress strain relation for drying and rewetting. This short coming can be avoided by using different retention curves for drying and rewetting.

R E F E R E N C E S

Genuchten, Van, M. Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil. Sci. Soc. A m . J., 44:892¬ 898.

Gibson, R. E., G. L. England & M.J.L. Hussey. 1967. The theory of one-dimensional consolidation Kim, D. J., J. Diels & J. Feyen. 1992. Water movement associated with overburden potential in a shrinkage marine clay soil. J. Hydrol. 133:179¬ 200.

Philip, J. R., 1968. Kinetics of sorption and volume change in clay-colloid pastes. Aust. J. Soil Res., 6:249-267.

Philip, J. R., 1969. Hydrostatics and hydrodynamics of swelling soils. Water Resour.Res. 5: 1070-1077.

Rijniersce, K., 1983, A simulation model for physical soil ripening in the IJsselmeerpolders. Lelystad, The Netherlands.

Schiffman, R. L., Vick, S. G. & Gibson, R. E., 1988. Behaviour and properties of hydraulic fills. Proc, Hydraulic Fill Structures.

Smiles, D. E. & M. J. Rosenthal. 1968. The movement of water in swelling materials. Aust. J. Soil Res., 6: 237-248.

Weijermars, E., 2011. Sedimentation of oil-sand tailings. Bachelor Thesis: TU Delft, Geosciences. Wilson, G. W., D. G. Fredlund & S. L. Barbour. 1994. Coupled soil athmosphere modeling for soil evaporation. Canadian Geotechnical Journal, 31(2): 151-161.

Wilson, G. W., D. G. Fredlund & S. L. Barbour. 1997. The effect of soil suction on evaporative fluxes from soil surfaces. Canadian Geotechnical Joumal, 34(1):145-157. Doi:10.1139/cgj-34-1-145. Yao, Y., L. A. van Paassen, A. F. van ToI, H. J. Everts & A. Mulder. 2010. Experimental research on accelerated consolidation using filter jackets in oil sands tailings. Proceedings of Second

International Oil Sands Tailings Conference, Edmonton, Alberta, Canada.