Estimating the Potential Differential Settlement of a Tailings Deposit Based on Consolidation Properties Heterogeneity

: Processing of extracted oil sands generates substantial volumes of tailings slurries. Due to the scale and inherent variability of the tailings properties, consolidation settlement is expected to occur at different rates and magnitudes across the tailings deposit. Estimating potential differential settlement of the consolidated deposit surface is an essential input for closure design. This paper presents a three-step methodology that generates multiple realizations of quasi-three-dimensional (3D) surfaces of the consolidated deposit based on the adjacent points. Each point is based on a stochastic one-dimensional (1D) large strain consolidation model developed with Monte Carlo techniques in GoldSim. The simulated surfaces provide early estimates of differential settlement based on the variability of consolidation properties expected in the tailings deposit. Comprehensive sensitivity analyses are performed for differently treated tailings material through 28 distinct scenarios to evaluate the sensitivity of the developed 1D and 3D models to consolidation input parameters over a 40-year time period. The analysis demonstrated that differential settlement is highly sensitive to tailings compressibility and hydraulic conductivity governed by the constitutive relationship parameters, and less sensitive to the solids content, speciﬁc gravity or thickness of a surcharge load. Tailings that underwent steady continuous settlement exhibited the largest degree of differential settlement.


Introduction
The Athabasca Oil Sands Region in Alberta, Canada, is home to the third-largest oil reserves after Venezuela and Saudi Arabia. As of 2016, Alberta's oil sands proven reserves were 165.4 billion barrels. Surface mining recovers 20% of this total volume and covers an area of about 4800 km 2 [1]. Fluid fine tailings (FFT) are a by-product of surface bitumen extraction activities, which have accumulated a total volume of 1.3 billion m 3 [2].
The government of Alberta views oil sands operations as temporary land use and requires oil sands operators to return disturbed land to equivalent land capability [3]. Reclamation efforts for such large-scale disturbances require the reconstruction of new landscape-scale ecosystems [4,5]. Revegetation is key to reclamation in a terrestrial landscape. Therefore, reclamation soil covers are designed to provide an adequate water supply for vegetation over dry summer periods [6,7]. However, most reclamation soil cover research does not explicitly consider the 3D nature of differential settlement within the context of large-scale landscapes, such as oil sand tailings deposits.
Understanding the dynamic nature in which reclamation takes place is fundamental in re-establishing ecosystem functionality. One of the constraints in landscape design is to limit the free water area to less than 10-20% of the reclaimed deposit [8,9]. However, differential settlement may result in shallow ponded water and wetlands or may deepen or expand the existing wet landscape features. There is a wealth of hydrological studies that quantify downslope movement of water within reclaimed upland slopes and temporal variability of soil water distribution [10][11][12]. However, little effort seems to be directed at the spatial variability of underlying tailings with respect to geotechnical stability and settlement. The large strain consolidation theory [13] is commonly used to model the nonlinear behavior of oil sand tailings consolidation. The spatial variability is of significance when reclamation covers overlay such highly compressible FFT. There is certainly a need to understand not only the settlement with time for these covers but also any potential 3D surface roughness or differential settlement, as it directly affects the contour of the final surface. In addition to reclamation and landscape designing, analyzing the consolidation of the treated tailings is also an essential input to predict the storage capacity of tailings storage facilities.
In any tailings storage facility, the three main stages involved in settlements of finegrained tailings material are flocculation, sedimentation and consolidation [14]. As the consolidation stage dominates the long-term settlement, the focus is typically on analyzing the consolidation stage [15]. A one-dimensional (1D) column consolidation is the basis for most of the current modeling practices and has been extrapolated to the whole deposits' consolidation behavior [16][17][18]. Another example of approximating the 3D consolidation based on a 1D large strain model is presented in [19], which studied the 3D analysis of a large strain consolidation model for tailings that are deposited in a mining pit. The goal was to estimate the filling time and tailings storage capacity. The authors performed Hydraulic Consolidation Test (HCT) to obtain the consolidation constitutive relationships, and generated approximate solutions of the 3D analysis. The proposed algorithm discretizes the pit volume in columns and layers to render upper and lower bound solutions based on a 1D large strain model.
The quasi-3D models have been developed to estimate the consolidation of tailings deposits, most of them based on a series of 1D columns of concentric cylinders or adjacent columns [20][21][22]. There exists a critical assumption in quasi-3D consolidation models; at the same elevation, the sides of the tailings storage facility face the same deformation as the tailings material faces in the center of the deposit. This assumption introduces estimation error to such 3D models, which is a function of the geometry of the tailings storage facility, boundary conditions, tailings material properties and the filling rate. For the tailings facilities with irregular impoundment geometry, non-homogeneous distribution of tailings and complex drainage paths, the error of quasi-3D consolidation models are magnified [23]. Due to these potentially large errors related to compressible foundation boundaries, the 3D large strain analyses are preferred over 1D large strain analysis for feasibility, and detailed level engineering work [24,25] developed an approach to eliminate the calculation errors caused by compressible boundaries of the tailings storage facilities. Reference [24] briefly discusses the common methods that are used in tailings consolidation modeling, with their relative accuracy and pitfalls, and specific recommendations and the applicability of each method based on the facility geometry, tailings properties and filling rates.
There are research works that address the multi-dimensional consolidation modeling without approximations based on 1D consolidation models. However, there are numerical challenges in multi-dimensional modeling of tailings consolidation. References [23,26] introduced a geotechnical numerical modeling approach to simulate the gradual tailings impoundment and its consolidation in multi-dimensional space based on the large strain consolidation theory. The proposed algorithm is based on the Norwegian Geotechnical Institute [NGI] modeling approach, and is developed on FLAC and FLAC3D, the commercially available specialized geotechnical modeling software. The authors validated the proposed algorithm using the benchmarks published in [27]. To model the gradual deposition of tailings, the slurry is divided into many continuous layers that are activated sequentially from the bottom. Once each new layer is activated, a large-strain consolidation stage is triggered and the consolidation time will be a function of discharge history and the layer volume. Published case studies were not found with statistical modeling capability to consider for heterogeneity and uncertainty of consolidating tailings deposits.
Reference [28] listed some of the effective techniques to measure the actual settlement of the field, including both remote sensing methods (photogrammetry, Light Detection and Ranging-Lidar) and on-tailings methods (hovercrafts and swamp machine). These post-deposition techniques provide a high accuracy of elevation data that is essential to generate a 3D model of tailings facilities for detailed design. The other group of techniques to generate a 3D model of the deposit surface includes different estimations of the potential surface profile using simulation techniques and numerical methods such as large strain consolidation theory [13]. These techniques work well as the tailings deposit can be fully characterized post deposition. However, successful closure planning begins during the early stages of the mine operations, before tailings are deposited or the mine operations break ground. Given the limited actual field deposit data described above that are available at the early stages of mine planning and operations, the impact of the variability of tailings properties within the deposit may not be fully appreciated. The merit of this paper is to generate a very first estimate of the surface roughness and differential settlements well before the tailings deposit is filled and construction of a cover begins.
The objective of this paper is to develop a framework for estimating the differential settlements of a consolidating tailings surface based on a series of adjacent consolidating 1D columns. This has been done through implementing a large strain consolidation model developed by [29,30]. The performance and utility of the developed model in this paper is demonstrated through several scenarios using a synthetic tailings deposit dataset. The sensitivity of numerical settlement predictions to solids contents, specific gravity and tailings constitutive relationships are also presented.
Different tailings types were investigated in this study, including composite tailings (CT), thickened tailings (TT), centrifuge cake and cyclone overflow. The presented framework is a tool to generate a first estimate of the probable surface roughness of a consolidating tailings deposit based on spatial heterogeneity in consolidation properties (including compressibility and hydraulic conductivity). As such, it is not aimed at estimating tailings storage volume. The developed model provides insight into the probable differential settlement behavior of a tailings deposit surface a priori, based on the anticipated range of variation in tailings properties.

Methodology
A 1D large strain consolidation model based on Gibson's theory was developed by [29,30], which is the primary input to the quasi-3D tailings surface contour model in this paper. Zheng's model simulates the consolidation of a column of deposited tailings and generates a height interface of the column over time, using System Dynamics techniques in the GoldSim environment.
The rationale to include uncertainty in the modeling is to consider the variability of tailings properties in a large deposit. Sources of variability in oil sands tailings properties include bitumen content, mineral solids and clay content, which is evident both spatially and temporally in tailings deposits, as reported by [31]. The tailings consolidation and geotechnical behavior are significantly influenced by the amount and type of clay present in the tailings [32][33][34][35], while both sand and clay contents are variable with depth of the deposit. This variability results in considerable heterogeneity in type and amount of clay minerals and the particle size distribution (PSD) of the tailings profile [36][37][38], which consequently influences the compressibility and permeability of the consolidated deposit [20,39]. The developed methodology captures the inherent uncertainty of tailings consolidation through multiple replications of the model and sampling from tailings' properties distributions (i.e., solids content, specific gravity and large strain constitutive relationships).
To estimate the surface roughness of the consolidating FFT, three steps are followed in a sequence, as illustrated in Figure 1. This methodology generates multiple quasi-3D profiles of the probable final surface from the 1D column consolidation GoldSim model. In the first step, the 1D consolidation model is run for multiple realizations through changing the solids content, specific gravity and sampling for the parameters of the consolidation constitutive relationships from user-defined probability distribution functions. The height interface curves are generated at the end of step one. A settlement curve is randomly selected in the second step and assigned to 'Defined points' in a two-dimensional (2D) grid. The Defined points are the points' locations where the deposit profile is assumed, and the value of solids contents, specific gravity, fines contents and the settlement curves that govern the settlement rate over time are therefore assigned accordingly. In deposits that have been filled, the 'Defined points' are the locations for which the characteristics such as solid contents and specific gravity are known as a result of field samplings. In step three, the elevation of all the remaining points of the surface, i.e., the 'Intermediate points', are calculated based on the inverse distance weighted averaging of the Defined points. As a result of step three, the surface at each time step is estimated, and if the model is run for multiple realizations, multiple potential surfaces will be generated. The three steps are performed in two different, but related GoldSim models.
To estimate the surface roughness of the consolidating FFT, three steps are followed in a sequence, as illustrated in Figure 1. This methodology generates multiple quasi-3D profiles of the probable final surface from the 1D column consolidation GoldSim model. In the first step, the 1D consolidation model is run for multiple realizations through changing the solids content, specific gravity and sampling for the parameters of the consolidation constitutive relationships from user-defined probability distribution functions. The height interface curves are generated at the end of step one. A settlement curve is randomly selected in the second step and assigned to 'Defined points' in a twodimensional (2D) grid. The Defined points are the points' locations where the deposit profile is assumed, and the value of solids contents, specific gravity, fines contents and the settlement curves that govern the settlement rate over time are therefore assigned accordingly. In deposits that have been filled, the 'Defined points' are the locations for which the characteristics such as solid contents and specific gravity are known as a result of field samplings. In step three, the elevation of all the remaining points of the surface, i.e., the 'Intermediate points', are calculated based on the inverse distance weighted averaging of the Defined points. As a result of step three, the surface at each time step is estimated, and if the model is run for multiple realizations, multiple potential surfaces will be generated. The three steps are performed in two different, but related GoldSim models. The methodology is implemented through the following steps: the settlement curves generated by the 1D GoldSim model are stored in a data file. A second data file contains all the inputs required to run the quasi-3D model, including the X and Y coordinates of all the grid points in separate tables for the Defined and Intermediate points. The quasi-3D GoldSim model reads these two files, calculates the elevation of Intermediate points over time using inverse distance weighted (IDW) averaging and writes them automatically in the output text file. At the final step, a series of Matlab functions read the text results, analyze the differential settlement and visualize the quasi-3D surface profiles and calculate the surface roughness at each specified time step and replication.

The 1D Model: Large Strain Consolidation Theory
The 1D large strain consolidation model developed by [30] is the basis for the developed methodology of this paper. The 1D model implements a finite difference scheme [40,41] to solve the large strain model with Gibson theory [13] using System Dynamics (SD) techniques. It simulates the long-term soil water dynamics and self-weight consolidation process using the Monte Carlo technique in GoldSim environment [42]. The GoldSim model is a reproducible open-source and transparent model that enables the technical and non-technical stakeholders to understand the mechanism of the simulation model. Reference [29] explicitly models two types of uncertainties: the aleatory uncertainties that are due to inherent randomness of the parameters, and the epistemic uncertainties that arise because of the lack of knowledge about different aspects of the problem. The methodology is implemented through the following steps: the settlement curves generated by the 1D GoldSim model are stored in a data file. A second data file contains all the inputs required to run the quasi-3D model, including the X and Y coordinates of all the grid points in separate tables for the Defined and Intermediate points. The quasi-3D GoldSim model reads these two files, calculates the elevation of Intermediate points over time using inverse distance weighted (IDW) averaging and writes them automatically in the output text file. At the final step, a series of Matlab functions read the text results, analyze the differential settlement and visualize the quasi-3D surface profiles and calculate the surface roughness at each specified time step and replication.

The 1D Model: Large Strain Consolidation Theory
The 1D large strain consolidation model developed by [30] is the basis for the developed methodology of this paper. The 1D model implements a finite difference scheme [40,41] to solve the large strain model with Gibson theory [13] using System Dynamics (SD) techniques. It simulates the long-term soil water dynamics and self-weight consolidation process using the Monte Carlo technique in GoldSim environment [42]. The GoldSim model is a reproducible open-source and transparent model that enables the technical and non-technical stakeholders to understand the mechanism of the simulation model. Reference [29] explicitly models two types of uncertainties: the aleatory uncertainties that are due to inherent randomness of the parameters, and the epistemic uncertainties that arise because of the lack of knowledge about different aspects of the problem.
Equations (1) and (2) present the constitutive relationships in the 1D model [43] and govern the change in the height, void ratio and effective stress of each given column of tailings material over time: where e is the void ratio, σ is the effective stress, k is the hydraulic conductivity and A, B, C and D are the constitutive relationships parameters. These parameters are determined from curve fitting of experimental data and are unique to each tailings sample.

Inverse Distance Weighted (IDW) Averaging
The elevation of Intermediate points in the surface at any given time step is estimated based on the weighted average inverse distance calculation in GoldSim, as presented in Equation (3): where, d i is the Euclidean distance between point x (the point of interest) and point i in the XY plane, S is the number of all Defined points, m is the inverse distance order and Z x is the elevation of Intermediate points on the surface at any given time. The formulation simply implies that the elevation of Intermediate point x is the average of the elevation of all Defined points. The closer the Defined point to x, the greater impact it has on Z x . Parameter m quantifies the significance of closeness of Defined points to point x in determining Z x . Choosing higher values, e.g., m = 3 (cubed inverse distance) means that as the distance increases, the contribution of the farther points' elevation in the average is significantly lower than that of closer points, compared to the case where m = 1 (simple inverse distance). An example is illustrated in Figure 2 to show how parameter m affects the quasi-3D estimation result. (2) where is the void ratio, is the effective stress, is the hydraulic conductivity and A, B, C and D are the constitutive relationships parameters. These parameters are determined from curve fitting of experimental data and are unique to each tailings sample.

Inverse Distance Weighted (IDW) Averaging
The elevation of Intermediate points in the surface at any given time step is estimated based on the weighted average inverse distance calculation in GoldSim, as presented in Equation (3): where, is the Euclidean distance between point x (the point of interest) and point i in the XY plane, S is the number of all Defined points, m is the inverse distance order and is the elevation of Intermediate points on the surface at any given time. The formulation simply implies that the elevation of Intermediate point x is the average of the elevation of all Defined points. The closer the Defined point to x, the greater impact it has on . Parameter m quantifies the significance of closeness of Defined points to point x in determining . Choosing higher values, e.g., m = 3 (cubed inverse distance) means that as the distance increases, the contribution of the farther points' elevation in the average is significantly lower than that of closer points, compared to the case where m = 1 (simple inverse distance). An example is illustrated in Figure 2 to show how parameter m affects the quasi-3D estimation result.

Differential Settlement Analysis
Differential settlement is calculated to statistically analyze the distribution of the difference in settlement of points on a grid. This analysis shows the roughness and steepness of the slopes in the final surface, which is a critical input to surface water management of the reclamation cover.
Once the surface profile of the consolidated deposit is numerically modeled as indicated in Equation (3), the differential settlement can be calculated and analyzed through different measures that are available in the literature. Several formulations for surface roughness calculations are listed and reviewed by [44]. The same concepts can be implemented to measure the surface roughness. The differential settlement metrics that are used in this paper are listed in Table 1 and illustrated schematically in Figure 3.

Differential Settlement Analysis
Differential settlement is calculated to statistically analyze the distribution of the difference in settlement of points on a grid. This analysis shows the roughness and steepness of the slopes in the final surface, which is a critical input to surface water management of the reclamation cover.
Once the surface profile of the consolidated deposit is numerically modeled as indicated in Equation (3), the differential settlement can be calculated and analyzed through different measures that are available in the literature. Several formulations for surface roughness calculations are listed and reviewed by [44]. The same concepts can be implemented to measure the surface roughness. The differential settlement metrics that are used in this paper are listed in Table 1 and illustrated schematically in Figure 3.   Another simple measure of differential settlement is the 'Delta Z' that m calculates the difference in the height of any given point and its four adjacent point left, right, top and bottom. The distance, in the XY plane, between the given point an  Another simple measure of differential settlement is the 'Delta Z' that merely calculates the difference in the height of any given point and its four adjacent points at its left, right, top and bottom. The distance, in the XY plane, between the given point and the four adjacent points is constant and equal to the grid size. Hence, the distribution and the statistics of the Delta Z is very dependent on the grid resolution. The user can assign the grid resolution based on their site-specific understanding of anticipated deposit profiles.

Model Validation
The performance of a proposed model can be validated through comparing the model results against the real field data, or the results of available software packages that are already validated through different benchmarks. The 1D model component of the proposed consolidation model is validated by [29,30] through three sets of scenarios, based on the reported cases by [27] as the well-accepted benchmarks. The benchmark cases have been widely used in the literature for evaluation of numerical large strain consolidation models, and are the results of nine different models that predicted the consolidation behavior of four waste clay disposal scenarios [23]. In [29,30], the authors simulated the tailings consolidation under different permeability and compressibility conditions for three types of tailings, i.e., untreated mature fine tailings (MFT), thickened tailings (TT) and phosphate tailings (PT). The authors considered a wide range of solid contents and initial deposit height to investigate the robustness of the model. The primary performance indicators used for validation include settlement over time, effective stress and void ratio profile. The commercial software FSCA is used to compare the key performance indicators for thickened tailings, based on the large strain consolidation theory.
The interface height for the validation cases are illustrated in Figure 4. The 1D model component provides acceptable settlement predictions, labeled TMSim and TMSim-consol in Figure 4, when benchmarked against commercial software and laboratory datasets. Readers are directed to [29,30] for comprehensive details about the validation of the 1D model. datasets. Readers are directed to [29,30] for comprehensive details about the validation of the 1D model. The quasi-3D model is tested for special conditions to verify its numerical results. As expected, it generates zero differential settlement (a flat surface) when all of the defined points follow a unique settlement curve. However, there are no available data about the actual surface of a deposit during the active life of a tailings storage facility or pos deposition after the pond is filled. A review of the publicly available tailings reports published by oil sands operators clearly states that during the active life cycle of a tailings storage facilities, because the deposits are constantly filled with the tailings, the operators do not track the consolidation: "Settlement and consolidation of fluid tailings are no tracked in active ponds, due to the complexity and inaccuracy inherent in measuring and calculating these values." [45]. Therefore, there are no valid data available to be compared against the results of the model. On the other hand, the developed model here generates a pseudo-3D estimation of the potential surface of the consolidated pond once the active deposition is over and the treated tailings has been consolidated. The model is no intended to predict the final surface profiles or for stage-volume predictions, but rather to inform the closure landform design process.
In the literature, the 2D and 3D models are rarely validated independently from 1D models [46]. Two of the reasons from the literature are: (1) the 2D effects become significant when the impoundment is relatively deep (with a width-to-height ratio in the order of 5 or less, associating with side drainages), whereas the majority of the tailings storage facilities in the oil sands region have much higher width-to-height ratios [47], and (2) while a 1D model satisfies most of the requirements [39], there are limited detailed published data to compare 3D models. Therefore, the developed model utilizes the previously validated 1D consolidation model to infer probable surface roughness or differential settlement statistics across a deposit under different scenarios due to uncertainties.

The Case Study
The authors reviewed a number of recent fluid tailings management reports from some of the largest Alberta oil sands operators [45,[48][49][50][51][52] that are publicly available and published in accordance with Directive 085 [1]. These reports include tables containing annual estimations of volumetric values of different slurry components such as fines, sand and water that will be generated in oil sands processing and pumped into the tailings storage facilities. The estimated values will be updated on an annual basis once the actua oil sands production and processing happens. The reports also include satellite maps cross-sectional and plan views of the tailings storage facilities with the sampling point The quasi-3D model is tested for special conditions to verify its numerical results. As expected, it generates zero differential settlement (a flat surface) when all of the defined points follow a unique settlement curve. However, there are no available data about the actual surface of a deposit during the active life of a tailings storage facility or post deposition after the pond is filled. A review of the publicly available tailings reports published by oil sands operators clearly states that during the active life cycle of a tailings storage facilities, because the deposits are constantly filled with the tailings, the operators do not track the consolidation: "Settlement and consolidation of fluid tailings are not tracked in active ponds, due to the complexity and inaccuracy inherent in measuring and calculating these values." [45]. Therefore, there are no valid data available to be compared against the results of the model. On the other hand, the developed model here generates a pseudo-3D estimation of the potential surface of the consolidated pond once the active deposition is over and the treated tailings has been consolidated. The model is not intended to predict the final surface profiles or for stage-volume predictions, but rather to inform the closure landform design process.
In the literature, the 2D and 3D models are rarely validated independently from 1D models [46]. Two of the reasons from the literature are: (1) the 2D effects become significant when the impoundment is relatively deep (with a width-to-height ratio in the order of 5 or less, associating with side drainages), whereas the majority of the tailings storage facilities in the oil sands region have much higher width-to-height ratios [47], and (2) while a 1D model satisfies most of the requirements [39], there are limited detailed published data to compare 3D models. Therefore, the developed model utilizes the previously validated 1D consolidation model to infer probable surface roughness or differential settlement statistics across a deposit under different scenarios due to uncertainties.

The Case Study
The authors reviewed a number of recent fluid tailings management reports from some of the largest Alberta oil sands operators [45,[48][49][50][51][52] that are publicly available and published in accordance with Directive 085 [1]. These reports include tables containing annual estimations of volumetric values of different slurry components such as fines, sand and water that will be generated in oil sands processing and pumped into the tailings storage facilities. The estimated values will be updated on an annual basis once the actual oil sands production and processing happens. The reports also include satellite maps, crosssectional and plan views of the tailings storage facilities with the sampling point coordinates marked on the maps, all in PDF format. The spatial data, including the sampling data such as the solid contents, specific gravity and the parameters required for consolidation modeling at marked coordinates, cannot be retrieved from the published reports and maps. Therefore, in the absence of publicly available data on the characterization of full-scale tailings deposits, a case study was developed as an analogue for model development that represents a potential real tailings deposit.
The hypothetical case includes the essential characteristics of a tailings deposit, such as the deposit size, the grid resolution, the coordinates of the grid points, the initial elevation of the deposit surface and the parameters used in the 1D consolidation model. An example of the hypothesized grid is illustrated in Figure 5. Larger grid cells (lower resolutions) can be considered for more homogeneous deposits. A 1 km 2 deposit size was chosen to reduce model runtime but could easily be scaled to larger deposit sizes. Given the very large lateral extents of the tailings ponds (several km 2 ) compared to the depth (50-100 m), boundary effects are expected to be minimal and drainage will predominantly be vertical. In this example, there are a total of 625 grid cells and 25 Defined points in a 1 km by 1 km square deposit. coordinates marked on the maps, all in PDF format. The spatial data, including the sampling data such as the solid contents, specific gravity and the parameters required for consolidation modeling at marked coordinates, cannot be retrieved from the published reports and maps. Therefore, in the absence of publicly available data on the characterization of full-scale tailings deposits, a case study was developed as an analogue for model development that represents a potential real tailings deposit. The hypothetical case includes the essential characteristics of a tailings deposit, such as the deposit size, the grid resolution, the coordinates of the grid points, the initial elevation of the deposit surface and the parameters used in the 1D consolidation model. An example of the hypothesized grid is illustrated in Figure 5. Larger grid cells (lower resolutions) can be considered for more homogeneous deposits. A 1 km 2 deposit size was chosen to reduce model runtime but could easily be scaled to larger deposit sizes. Given the very large lateral extents of the tailings ponds (several km 2 ) compared to the depth (50-100 m), boundary effects are expected to be minimal and drainage will predominantly be vertical. In this example, there are a total of 625 grid cells and 25 Defined points in a 1 km by 1 km square deposit.

Scenarios
The first step of the methodology is to run the 1D model and generate the required settlement curves for the quasi-3D surface estimation runs. The 1D model is run for a total of 67 cases, and the results are investigated to select the scenarios for quasi-3D runs. This resulted in selecting 28 cases for more comparison on the effect of different solids contents and specific gravity (9 cases), surcharge (6 cases) and different material types reflecting a range of A, B, C and D parameters (13 cases). The parameters for these cases are presented in Table 2. Among these cases, A1 and N1 are tailings derived from a centrifuge dewatering process, F1 and F2 are cyclone overflow, C1 and C2 are composite tailings, and the rest are flocculated and thickened tailings from inline or conventional tank thickeners. The centrifuged tailings are representative of tailings with higher fines (~90% fines < 44 μm) and clay contents while the composite tailings are much coarser (~20 to 25% fines < 44 μm) [53].

Scenarios
The first step of the methodology is to run the 1D model and generate the required settlement curves for the quasi-3D surface estimation runs. The 1D model is run for a total of 67 cases, and the results are investigated to select the scenarios for quasi-3D runs. This resulted in selecting 28 cases for more comparison on the effect of different solids contents and specific gravity (9 cases), surcharge (6 cases) and different material types reflecting a range of A, B, C and D parameters (13 cases). The parameters for these cases are presented in Table 2. Among these cases, A1 and N1 are tailings derived from a centrifuge dewatering process, F1 and F2 are cyclone overflow, C1 and C2 are composite tailings, and the rest are flocculated and thickened tailings from inline or conventional tank thickeners. The centrifuged tailings are representative of tailings with higher fines (~90% fines < 44 µm) and clay contents while the composite tailings are much coarser (~20 to 25% fines < 44 µm) [53]. For any of these cases, the stochastic 1D model is run for five replications of nested Monte Carlo simulation for 100 years of simulation time. The stochastic A, B, C and D parameters for the inner 1D models, as originally developed by [29], are sampled from BETA_PERT distributions within 10% deviation around the outer A, B, C and D parameters. This deviation envelopes the expected variation of model parameters and is user-adjustable for different cases based on the statistical distribution of the material properties. The inner and outer 'A' parameters are shown in Figure 6. The three groups of settlement curves, each associated with one of the scenarios, are presented in Figures 7-9.
According to Figure 7, as the solids content increases from 50% to 70%, the column of treated tailings settles less over time due to lower void ratios associated with higher solids contents. Figure 7 also shows the settlement is larger for higher specific gravity regardless of the solids content. Figure 8 shows that having a surcharge load will slightly increase the settlement of the column. In both of the two cases (colored in blue and red), the solid lines on top show the situation where there is no surcharge load, the broken lines at the middle represent the cases with two-meter surcharge cap on top of the deposit and the broken lines at the bottom show the cases with a 4 m surcharge cap. However, comparing the two studied scenarios shows that the addition of a 2 m cap does not significantly change the settlement rate over time. That is because there is a very small difference between the settlement of the material with 2 m versus 4 m caps for both solids contents of 60% and 70%. Figure 9 shows how various material types with different constitutive relationship parameters will settle over time. The various range of tailings material in this figure show the applicability of the developed 1D and quasi-3D consolidation models on comparable cases. The scenarios presented in Figure 9 span a large range of tailings material, treated with different technologies and in different stages, including the centrifuge cake, thickened tailings, cyclone overflow and composite tailings. The solids contents and specific gravity of these cases are provided in Table 2. Based on the settlement curves, different materials show one of the five settlement behaviors: the first group includes scenario A1, the centrifuge cake, showing very little settlement over time. The second group of scenarios consists of the thickened tailings, including cases G1 and R1, showing a continuous settlement over time. The third group contains cases F1 and F2, the cyclone overflow, showing a very steep curve in the very early periods, then a flat curve without any more settlements over the rest of the years. The fourth group includes cases C1 and C2, which are composite tailings with a relatively rapid settlement up to year 10, and then continues with a flat curve without any considerable settlement for the rest of years. The fifth group, again thickened tailings material, contains S1 and S11 with a rapid settlement profile with time up to 20 years at which time it plateaus. Figure 10 shows the range in settlement for a particular tailings case with stochastic constitutive relationship parameters that are modeled in GoldSim. The shaded range shows the extent of settlement between the upper and lower bounds, each representing a set of consolidation constitutive relations parameters within the 10% range around the A, B, C and D parameters. The maximum of 5.4 m difference in settlement range in year 40 shows how the settlement is sensitive to the variability of the parameters. Figure 11 presents the constitutive relationship curves to understand the settlement rate of these scenarios as the relations of effective stress, void ratio and hydraulic conductivity for the same scenarios presented in Figure 9. The impact of PSD is captured in these constitutive relationships. Coarser tailings would tend to have higher hydraulic conductivities and lower compressibilities, whereas high fines and clay content tailings would exhibit lower hydraulic conductivities and greater compressibility [53]. For any of these cases, the stochastic 1D model is run for five replications of nested Monte Carlo simulation for 100 years of simulation time. The stochastic A, B, C and D parameters for the inner 1D models, as originally developed by [29], are sampled from BETA_PERT distributions within 10% deviation around the outer A, B, C and D parameters. This deviation envelopes the expected variation of model parameters and is user-adjustable for different cases based on the statistical distribution of the material properties. The inner and outer 'A' parameters are shown in Figure 6. The three groups of settlement curves, each associated with one of the scenarios, are presented in Figures 7-9.      According to Figure 7, as the solids content increases from 50% to 70%, the column of treated tailings settles less over time due to lower void ratios associated with higher solids contents. Figure 7 also shows the settlement is larger for higher specific gravity    According to Figure 7, as the solids content increases from 50% to 70%, the column of treated tailings settles less over time due to lower void ratios associated with higher solids contents. Figure 7 also shows the settlement is larger for higher specific gravity regardless of the solids content. Figure 8 shows that having a surcharge load will slightly increase the settlement of the column. In both of the two cases (colored in blue and red), the solid lines on top show the situation where there is no surcharge load, the broken lines at the middle represent the cases with two-meter surcharge cap on top of the deposit and the broken lines at the bottom show the cases with a 4 m surcharge cap. However, comparing the two studied scenarios shows that the addition of a 2 m cap does not significantly change the settlement rate over time. That is because there is a very small difference between the settlement of the material with 2 m versus 4 m caps for both solids contents of 60% and 70%. Figure 9 shows how various material types with different constitutive relationship parameters will settle over time. The various range of tailings material in this figure show the applicability of the developed 1D and quasi-3D consolidation models on comparable cases. The scenarios presented in Figure 9 span a large range of tailings material, treated with different technologies and in different stages, including the centrifuge cake, thickened tailings, cyclone overflow and composite tailings. The solids contents and specific gravity of these cases are provided in Table 2. Based on the settlement curves, different materials show one of the five settlement behaviors: the first group includes scenario A1, the centrifuge cake, showing very little settlement over time. The second group of scenarios consists of the thickened tailings, including cases G1 and R1, showing a continuous settlement over time. The third group contains cases F1 and F2, the cyclone overflow, showing a very steep curve in the very early periods, then a flat curve without any more settlements over the rest of the years. The fourth group includes cases C1 and C2, which are composite tailings with a relatively rapid settlement up to year 10, and then continues with a flat curve without any considerable settlement for the rest of years. The fifth group, again thickened tailings material, contains S1 and S11 with a rapid settlement profile with time up to 20 years at which time it plateaus. Figure  10 shows the range in settlement for a particular tailings case with stochastic constitutive relationship parameters that are modeled in GoldSim. The shaded range shows the extent of settlement between the upper and lower bounds, each representing a set of consolidation constitutive relations parameters within the 10% range around the A, B, C and D parameters. The maximum of 5.4 m difference in settlement range in year 40 shows how the settlement is sensitive to the variability of the parameters. Figure 11 presents the constitutive relationship curves to understand the settlement rate of these scenarios as the relations of effective stress, void ratio and hydraulic conductivity for the same scenarios presented in Figure 9. The impact of PSD is captured in these constitutive relationships. Coarser tailings would tend to have higher hydraulic conductivities and lower compressibilities, whereas high fines and clay content tailings would exhibit lower hydraulic conductivities and greater compressibility [53].

Sensitivity Analysis-1D Consolidation Model
The presented scenarios in the previous section shows the sensitivity of the 1D consolidation model and the settlement to some of the selected parameters, such as the surcharge load, constitutive relation parameters, solid contents and the specific gravity. However, there are more parameters involved in 1D and multidimensional consolidation. Reference [54] performed a comprehensive parametric sensitivity analysis to examine if, and to what extent, the consolidation results are sensitive to each of the material property functions and design variables such as compressibility constants, hydraulic conductivity constants, rate of rise, specific gravity, years of deposition, surcharge load application, and initial solids content.
The hydraulic conductivity and the compressibility functions have a combined effect on the consolidation behavior. To capture this combined effect, a full consolidation model must be run that simultaneously accounts for physical characteristics, such as specific gravity and initial solids content, pore water pressure and the stresses associated with tailings deposition. Reference [54] modeled a full consolidation model using the 1D finite strain program FSConsol © . The authors run over 700 models in total, each for a 1000-year simulation time, and the following three results of each model run were recorded for further analysis: (1) height after 1000 years, as an indicator of volume reduction; (2) time to reach 90% consolidation, as an indicator of reclamation timeframes; and (3) solids content after 1000 years, as an indicator of shear strength. Among these three recorded results, the statistical analysis regarding the second measure, i.e., the time to reach 90% consolidation, is presented in Figure 12, after the original work in [54].

Sensitivity Analysis-1D Consolidation Model
The presented scenarios in the previous section shows the sensitivity of the 1D consolidation model and the settlement to some of the selected parameters, such as the surcharge load, constitutive relation parameters, solid contents and the specific gravity. However, there are more parameters involved in 1D and multidimensional consolidation. Reference [54] performed a comprehensive parametric sensitivity analysis to examine if, and to what extent, the consolidation results are sensitive to each of the material property functions and design variables such as compressibility constants, hydraulic conductivity constants, rate of rise, specific gravity, years of deposition, surcharge load application, and initial solids content.
The hydraulic conductivity and the compressibility functions have a combined effect on the consolidation behavior. To capture this combined effect, a full consolidation model must be run that simultaneously accounts for physical characteristics, such as specific gravity and initial solids content, pore water pressure and the stresses associated with tailings deposition. Reference [54] modeled a full consolidation model using the 1D finite strain program FSConsol © . The authors run over 700 models in total, each for a 1000-year simulation time, and the following three results of each model run were recorded for further analysis: (1) height after 1000 years, as an indicator of volume reduction; (2) time to reach 90% consolidation, as an indicator of reclamation timeframes; and (3) solids content after 1000 years, as an indicator of shear strength. Among these three recorded results, the statistical analysis regarding the second measure, i.e., the time to reach 90% consolidation, is presented in Figure 12, after the original work in [54].
Appl. Sci. 2022, 12, 6206 13 of 19 Figure 12. Years to 90% consolidation-after [54]. Figure 12 shows the effect of the model parameters on the range of time to reach 90% consolidation. The first two rows show that changing all of the constitutive relation parameters ends up in a wide range in results. The overall depth of deposit is determined by the Rate of Rise (RoR) and years of depositions (Years) together, and hence, have a significant effect on consolidation times too. The C and D parameters are very influential, Figure 12. Years to 90% consolidation-after [54]. Figure 12 shows the effect of the model parameters on the range of time to reach 90% consolidation. The first two rows show that changing all of the constitutive relation parameters ends up in a wide range in results. The overall depth of deposit is determined by the Rate of Rise (RoR) and years of depositions (Years) together, and hence, have a significant effect on consolidation times too. The C and D parameters are very influential, and over the studied cases, C has a more significant effect on consolidation than D. The parameters A and B are less influential, to some extent because the 90% consolidation endpoint accounts partially for the compressibility curve. The figure shows that the initial solids content has a relatively strong effect, while the surcharge loading has a very small contribution in consolidation time. The authors conclude that the hydraulic conductivity function is very important in determining the consolidation time, compressibility and hydraulic conductivity inputs have a complex influence on the coefficients of consolidation, and the relationships between rate of raise and the extent of initial consolidation is another important factor in consolidation [54].

Quasi-3D Settlement Results
The final step of the methodology in Figure 1 is to generate the quasi-3D surfaces through sampling the settlement curves from the 1D results, generating multiple realizations of the surface profile and analyzing the differential settlement. Among the listed 1D scenarios in Table 2, eleven distinct cases are selected to represent the ranges of different material types, solids contents and specific gravity. The settlement curves of the chosen cases are randomly assigned to the Defined points of the hypothesized grid. The quasi-3D scenarios that are run for 40 years of simulation time and 10 replications are listed in Table 3.
The simulated surface topography is a quasi-3D surface that is generated from replicating the 1D column consolidation model on the 25 defined points, and calculating the simulated Z values at all other undefined points across the deposit surface (1 km by 1 km) based on an extrapolation of at the defined points. The quasi-3D model does not consider the lateral influences among the columns or between columns and facility side walls, as should be considered in real 3D consolidation modeling. This modeling approach has been replicated in the literature, for example in [18,19], and is a valid approach when the goal is to have an early estimation of the potential surface or expected differential settlement in the deposit for long-term decisions where high accuracy of the results is not intended. It is assumed that all 25 columns at the defined points have the same bottom elevation. With this limiting assumption, the simulated quasi-3D surface can be considered as an upper bound for the probable consolidated surface. Additionally, the model assumes the lateral extents are sufficiently large, and hence, the deformation of the side walls of the facility are not considered in consolidation modeling. The real tailings ponds are several kilometers wide with dams 50-100 m high, mainly with one-dimensional drainage. If the pond configurations are different from the assumptions of the presented model, i.e., for ponds with highly variable base topography, relatively small ponds with binding side walls or under valley infill scenarios, the model would not be representative.    The simulated surface topography is a quasi-3D surface that is generated from replicating the 1D column consolidation model on the 25 defined points, and calculating the simulated Z values at all other undefined points across the deposit surface (1 km by 1 km) based on an extrapolation of at the defined points. The quasi-3D model does not consider the lateral influences among the columns or between columns and facility side walls, as should be considered in real 3D consolidation modeling. This modeling approach has been replicated in the literature, for example in [18,19], and is a valid approach when the goal is to have an early estimation of the potential surface or expected differential settlement in the deposit for long-term decisions where high accuracy of the results is not intended. It is assumed that all 25 columns at the defined points have the same bottom elevation. With this limiting assumption, the simulated quasi-3D surface can be considered as an upper bound for the probable consolidated surface. Additionally, the model assumes the lateral extents are sufficiently large, and hence, the deformation of the side walls of the facility are not considered in consolidation modeling. The real tailings ponds are several kilometers wide with dams 50-100 m high, mainly with onedimensional drainage. If the pond configurations are different from the assumptions of the presented model, i.e., for ponds with highly variable base topography, relatively small ponds with binding side walls or under valley infill scenarios, the model would not be representative.
The analytical results of the eleven quasi-3D scenarios are compared in Table 4. The reported values show the mean of the metrics over ten replications for each experiment at year 40. Among the listed metrics, the average elevation and the maximum settlement are governed by the 1D consolidation model, while the maximum slope, the standard deviation of elevations around the average height and the roughness show the The analytical results of the eleven quasi-3D scenarios are compared in Table 4. The reported values show the mean of the metrics over ten replications for each experiment at year 40. Among the listed metrics, the average elevation and the maximum settlement are governed by the 1D consolidation model, while the maximum slope, the standard deviation of elevations around the average height and the roughness show the unevenness of the surface and are more affected by the quasi-3D model. Other metrics, such as the maximum valley depth and peak height and the distances between peaks and valleys, can also be viewed as roughness measures. However, the main metrics that are more representative of the material settlement and surface roughness are the maximum settlement (meters) and the maximum gradient slope (%). The results of the three solids contents scenarios are compared in Figure 14A. As expected, the case S1 with 50% solids contents settles the most due to the higher volume of water that it contains and expels over time, compared to the cases S2 and S3. Figure 14B shows the same pattern of maximum settlement values for specific gravity experiments. In this case, the material with the highest Gs (S7) shows the highest settlement because it contains heavier contents (less bitumen), which tends to settle quicker than the other two cases S1 and S4. However, the maximum slope does not show a direct relation to the solids contents or the Gs. The results of tailings material scenarios can help in explaining the maximum slope results. The results of the three solids contents scenarios are compared in Figure 14A. As expected, the case S1 with 50% solids contents settles the most due to the higher volume of water that it contains and expels over time, compared to the cases S2 and S3. Figure 14B shows the same pattern of maximum settlement values for specific gravity experiments. In this case, the material with the highest Gs (S7) shows the highest settlement because it contains heavier contents (less bitumen), which tends to settle quicker than the other two cases S1 and S4. However, the maximum slope does not show a direct relation to the solids contents or the Gs. The results of tailings material scenarios can help in explaining the maximum slope results. According to Figure 15, the type of tailings material affects both settlement of the thickened tailings over time and the roughness of the final deposit surface. In fact, different constitutive relationship parameters (the A, B, C and D parameters in Equations (1) and (2)) will directly determine the settlement behavior of the material and indirectly the roughness of the final surface of the consolidated deposit. According to Figure 15, the type of tailings material affects both settlement of the thickened tailings over time and the roughness of the final deposit surface. In fact, different constitutive relationship parameters (the A, B, C and D parameters in Equations (1) and (2)) will directly determine the settlement behavior of the material and indirectly the roughness of the final surface of the consolidated deposit.
Among the numerical metrics of the surface differential settlements, slope analysis is an important consideration for closure landform decision making. Assessing the magnitude of slope grades offers insight into the potential future peaks, valleys, surface water flow and possibility for seasonal small and large pools of water. The visual analysis of the quasi-3D scenario results shows that the roughness of the final surface is more affected by the type of tailings materials, rather than the solids contents or the specific gravity. In order to analyze the slope distribution statistically, all the combinations of the points on the final surface are paired, and the slopes of the straight lines between all the pairs are calculated. Then, the slope values are sorted based on the straight distance between the points and grouped in 100-m classes. Figure 16 shows a sample slope distribution as box plots for one of the replications of case R1. The box plots show that the slopes are relatively steeper in short distances (smaller than 400 m). The median of slope values is 1% for a 100-m distance, with a maximum slope of 5.1%, which is also reported in Table 4 and Figure 15. The median of slope values is reduced to around 0.5% for 500 m distances, then to 0.25% in 800 m and remains the same for farther distances up to 1500 m. This means that the large slope values, i.e., the outlier red plus signs around the box plots, are seen in smaller distances around the points. As the distance increases, the steep slopes are replaced with slight slopes. This general pattern is seen in all cases. However, the difference between the more uneven surfaces (such as case R1) and the smoother surfaces (like G1) is the different values of median and the maximum slopes. Among the numerical metrics of the surface differential settlements, slope analysis is an important consideration for closure landform decision making. Assessing the magnitude of slope grades offers insight into the potential future peaks, valleys, surface water flow and possibility for seasonal small and large pools of water. The visual analysis of the quasi-3D scenario results shows that the roughness of the final surface is more affected by the type of tailings materials, rather than the solids contents or the specific gravity. In order to analyze the slope distribution statistically, all the combinations of the points on the final surface are paired, and the slopes of the straight lines between all the pairs are calculated. Then, the slope values are sorted based on the straight distance between the points and grouped in 100-m classes. Figure 16 shows a sample slope distribution as box plots for one of the replications of case R1. The box plots show that the slopes are relatively steeper in short distances (smaller than 400 m). The median of slope values is 1% for a 100-m distance, with a maximum slope of 5.1%, which is also reported in Table 4 and Figure 15. The median of slope values is reduced to around 0.5% for 500 m distances, then to 0.25% in 800 m and remains the same for farther distances up to 1500 m. This means that the large slope values, i.e., the outlier red plus signs around the box plots, are seen in smaller distances around the points. As the distance increases, the steep slopes are replaced with slight slopes. This general pattern is seen in all cases. However, the difference between the more uneven surfaces (such as case R1) and the smoother surfaces (like G1) is the different values of median and the maximum slopes.

Conclusions
Integration of closure planning and design through the full life cycle of a mine is a key component of a successful tailings management plan. An early understanding of the range of potential tailings deposit surface roughness that may be encountered is an important input to landscape design for mine closure planning. The main contribution of an important consideration for closure landform decision making. Assessing the magnitude of slope grades offers insight into the potential future peaks, valleys, surface water flow and possibility for seasonal small and large pools of water. The visual analysis of the quasi-3D scenario results shows that the roughness of the final surface is more affected by the type of tailings materials, rather than the solids contents or the specific gravity. In order to analyze the slope distribution statistically, all the combinations of the points on the final surface are paired, and the slopes of the straight lines between all the pairs are calculated. Then, the slope values are sorted based on the straight distance between the points and grouped in 100-m classes. Figure 16 shows a sample slope distribution as box plots for one of the replications of case R1. The box plots show that the slopes are relatively steeper in short distances (smaller than 400 m). The median of slope values is 1% for a 100-m distance, with a maximum slope of 5.1%, which is also reported in Table 4 and Figure 15. The median of slope values is reduced to around 0.5% for 500 m distances, then to 0.25% in 800 m and remains the same for farther distances up to 1500 m. This means that the large slope values, i.e., the outlier red plus signs around the box plots, are seen in smaller distances around the points. As the distance increases, the steep slopes are replaced with slight slopes. This general pattern is seen in all cases. However, the difference between the more uneven surfaces (such as case R1) and the smoother surfaces (like G1) is the different values of median and the maximum slopes.

Conclusions
Integration of closure planning and design through the full life cycle of a mine is a key component of a successful tailings management plan. An early understanding of the range of potential tailings deposit surface roughness that may be encountered is an important input to landscape design for mine closure planning. The main contribution of

Conclusions
Integration of closure planning and design through the full life cycle of a mine is a key component of a successful tailings management plan. An early understanding of the range of potential tailings deposit surface roughness that may be encountered is an important input to landscape design for mine closure planning. The main contribution of this paper is the development of a combined 1D and quasi-3D simulation framework to evaluate the surface profile variability of a tailings deposit along with various sensitivity analyses. The 1D and quasi-3D models are developed in GoldSim as a reproducible open-source model that makes them understandable for technical and non-technical stakeholders, rather than typical black-box approaches of commercial software tools. The developed models are implemented under different tailings material type scenarios (thickened tailings, centrifuge cake, composite tailings and cyclone overflow) on a hypothesized grid set.
The numerical results demonstrated tailings that undergo continuous settlement throughout the simulation time (thickened tailings, R1 and G1) will exhibit the largest degree of differential settlement when considering elevation differences between peaks and valleys, greatest maximum slope and potential roughness. Tailings that undergo rapid settlement during time of interest (thickened tailings S1, S11 or composite tailings C1) or very little settlement (centrifuge cake, A1) exhibit much less differential settlement across the deposit. Differential settlements were less sensitive to the tailings properties, such as the specific gravity, which is a function of site geology, and the solids contents for the tailings evaluated (S1).
The developed framework to investigate differential settlement provides insight for landform designers so that designs can be prepared early in the mine life to accommodate the differential settlement or plans developed for ongoing maintenance as needed. These evaluations can provide insight on the impacts to surface drainage features and the potential for unplanned ponding as a result of heterogeneity of the underlying tailings deposit. While the model provides an assessment of the magnitude of differential settlement, it cannot predict the exact locations of flow paths, areas of high and low relief or the gradients of the future surface. Ongoing monitoring and post deposition measurements of the tailings deposit should be conducted and integrated into the closure design process throughout the life of the tailings deposit to update and inform final detailed closure landform design.