Modelling and constructing engineered soil for post-mining landforms to create optimal ecological outcomes: application of a computer based Landscape Evolution Model (SSSPAM)

SSSPAM weathering module theoretical framework

The weathering module of SSSPAM is comprised of two main components: (1) the weathering geometry, which governs the grading of daughter particles; and (2) the rate of weathering for various soil layers, which controls the rate of weathering of the parent material. The weathering rate of each soil layer generally (but not necessarily) varies with subsurface depth. SSSPAM focuses solely on physical weathering to degrade both the surface armour layer and the subsurface soil layers, ensuring mass conservation as a weathering assessment demonstrated this to be the key mechanism (Welivitiya and Hancock 2022, 2024). Processes like chemical dissolution are not currently accounted for, although the model’s design allows for the inclusion of chemical weathering if necessary.

Weathering framework of SSSPAM

In SSSPAM, the soil grading at a given time and depth is described by the state vector, g_. Each entry, g_i, in this vector represents the proportion of material within the grading size range, i. Transformations in soil grading resulting from various geomorphological processes, including weathering, are modelled using transition matrices.

SSSPAM’s weathering mechanism model was employed to calculate the weathering rates and processes for different spoil materials collected from the field. This was achieved by calibrating the module parameters to align them with experimental observations (Welivitiya and Hancock 2022, 2024). The weathering process of SSSPAM is described by using a mathematical framework known as the weathering transition matrix. The matrix entries, B, reflect changes in soil grading size classes, based on the fracture mechanism that is associated with the weathering process. This approach enables SSSPAM to simulate a range of particle fracture behaviours, representing different weathering mechanisms.

Throughout the weathering process, larger parent particles break down into smaller daughter particles, causing each grading size class to lose some of its mass to smaller classes. The daughter products may be distributed into one or more smaller grading classes, based on the particle size range produced when parent particles break down. The quantity of material transferred to each smaller class is determined by the size distribution of the grading classes, the fracture mechanism, and the size properties of the daughter particles.

In SSSPAM’s weathering module, soil grading transformations due to weathering are represented by a matrix equation:

(1)gw_=[I+(W_Δt)B]gr_

where gr_ is the soil grading vector before the weathering process, gw_, is the soil grading vector after the weathering process, and I and W_ are the identity matrix and the weathering rate vector, respectively. The weathering rate vector indicates that different sized particles may undergo varying degrees of weathering, with each element of W_ and W_i corresponding to the weathering rate for a particular grading size class. Δt is to SSSPAM’s time step duration when weathering takes place, and B is the weathering transition matrix.

To generalise the fracture behaviours, consider a parent particle with a diameter of d that breaks into one larger daughter particle with a diameter of d₁ and n − 1 smaller daughter particles with diameters of d₂, where the total number of daughter particles is n (see Fig. 3). Within the weathering mechanism framework, the largest daughter particle can be viewed as the remaining portion of the original parent particle, whereas the other particles are fragments that were separated from the parent particle during the weathering process.

Fig. 3.

Illustrative diagram of SSSPAM’s weathering model (Source: Welivitiya and Hancock 2024).

For simplicity, it is assumed that all particles are spherical. The principle of mass conservation indicates:

(2)d3=d13+(n−1)d23

If the larger daughter particle with a diameter of d₁ represents an α fraction of the parent particle, then:

(3)d1=α13d

Changing the α fraction value and the number of daughter particles, n, enables the simulation of different fracture behaviours, including symmetric fragmentation, asymmetric fragmentation, and granular disintegration (Wells et al. 2008):

(4)d2=( 1−αn−1 )13d

For example, when α = 0.5 and n = 2, it represents two daughter particles with symmetric fragmentation. When α = 0.99 and n = 11, it models a fracture mechanism similar to granular disintegration, where the largest daughter particle retains 99% of the original volume, and the remaining 10 smaller daughter particles together account for 1% of the original volume. In theory, the weathering module in SSSPAM can simulate any combination of α and n. However, the distribution of the daughter particles resulting from the weathering process is also influenced by the discretisation of the particle size distribution.

The range of particles resulting from the weathering of a single soil grading size class is computed by using Eqns 3 and 4. Let d_(i) represent the upper bound and d_(i−1) represent the lower bound of the ith soil grading class. The diameter range of larger particle can be determined as [α 1 3d(i−1), α 1 3d(i)] using Eqn 3. Similarly, Eqn 4 can be used to calculate the diameter range for the smaller particles, as [ ( 1−αn−1 ) 1 3d(i−1), ( 1−αn−1 ) 1 3d(i),]. The fragmentation process and the distribution of the resulting daughter particles are in Fig. 4.

Fig. 4.

Distribution of daughter particles resulting from the weathering of larger size classes. (Left) Transition from the parent to daughter grading classes, where M_i denotes the fraction of material undergoing weathering. (Right) Distribution of the larger daughter particles across different size classes (Source: Welivitiya and Hancock (2024) after Cohen et al. (2009)).

Within the context of modelling, it is useful to express the mass of each size class as mass fractions. The mass fractions remaining in the parent material grading size class and those transferred to smaller grading size classes are determined through an interpolation technique. Fig. 5 shows the process of distributing mass fractions to the smaller grading classes.

Fig. 5.

Interpolation of the contributing grading fractions of daughter products into smaller grading classes (Source: Welivitiya and Hancock 2024).

Depth dependent weathering functions

SSSPAM uses two primary pedogenesis functions: (1) the exponential decline, referred to as ‘Exponential’ (Humphreys and Wilkinson 2007); and (2) the humped exponential decline, referred to as ‘Humped Exponential’ (see Fig. 2) (Ahnert 1977; Minasny and McBratney 2006). Additionally, the SSSPAM modelling framework has been enhanced by the inclusion of a theoretical weathering function called the Dynamic Reversed Exponential function. In this function, the weathering rate peaks at the soil-bedrock interface and decreases exponentially in both upward towards the surface and downward towards the bedrock. In contrast to the Exponential and Humped functions, the Dynamic Reversed Exponential function allows the peak weathering rate to vary with changes in the depth of the soil-bedrock interface. The rationale for this function is that the transformation from bedrock to soil at the soil-bedrock interface may happen more rapidly than in other layers.

In the Exponential function, the weathering rate at a depth of h (m) below the soil surface, denoted as w_h, is expressed as:

where β is a constant that defines the maximum weathering rate, and δ₁ is the depth scaling factor.

For the Humped Exponential function, the weathering rate at a depth of h (m) below the soil surface, denoted as w_h, is expressed as:

where P₀ and P_a is the maximum and steady-state weathering rates, respectively, δ₂ and δ₃ are constants that define the function’s shape, and W_max is the maximum weathering rate at the hump’s peak, which is employed for normalising the function.

In the Dynamic Reversed Exponential function, the weathering rate of a layer h (m) beneath the soil surface, denoted as w_h, is expressed as:

where H (m) is the depth from the surface to the soil-bedrock interface, calculated based on the soil grading distribution at each iteration throughout the simulation.

Landscape and digital elevation model

A typical post-mining landform was used that was representative of many constructed landforms and provides an ideal test case for examining erosion, pedogenesis, and landform evolution (Fig. 1). The landform was 200 m long and 100 m wide (Fig. 6). At the top of the landform, a flat area of 50 m long gives way to a gentle slope of 20%. Halfway to the bottom of the hillslope, a back-sloping area (5% slope) of 10 m represents a break in slope, which is commonly constructed on post-mining landforms. High elevation areas bounded the landform on the right and left sides, acted as a catchment boundary for the modelling. A digital elevation model (DEM) with 1 m grid resolution (Fig. 6) was created to represent this landform. This size and grid dimension was sufficiently large to examine field scale processes.

Fig. 6.

Digital elevation model of the landform (1 m grid). This landform is used for all simulations to evaluate erosion, pedogenesis and landform evolution.

Post-construction, the site had been seeded with a mix of pasture, shrub, and tree species with the goal for the landscape to sustain cattle grazing. However, vegetation establishment was not successful, with the surface devoid of vegetation for approximately 20 years (Fig. 1). Hence the site is an ideal test site to examine pedogenesis without including the increased complexity of vegetation on soil production and function.

SSSPAM parameters

SSSPAM parameters were developed for mine material in the Bowen Basin, Queensland. The fluvial erosion rate in SSSPAM is calculated based on discharge per unit width and slope where e is the erodibility factor and α₁ and α₂ are exponents on discharge and slope, respectively. The sediment transport equation is described and parameter set developed (named 1SLV) using runoff and erosion data gathered using a laboratory flume, following the procedure described by Welivitiya and Hancock (2022) (Table 1).

Soil particle size distribution and rainfall data

Daily rainfall data obtained from the Bureau of Meteorology (BOM) weather station at Iffley, QLD (station no: 034100), the closest to the study site, was used here. The available daily rainfall record was repeated end-to-end, in order to generate a complete daily rainfall time-series of 3000 years. A 3000-year period was considered sufficient time for an equilibrium or steady state erosion output to be achieved.

Particle size distribution data of the material is needed as an input to the SSSPAM model. The particle size distribution of the mine material was determined by sieve analysis of the 1SLV material (Table 2, Welivitiya and Hancock 2022).

Influence of depth dependent weathering function

The three soil production functions all produce different temporal erosion patterns and rates (Fig. 7). Therefore, the choice of DDWF influences surface evolution (Temme and Vanwalleghem 2016; van der Meij et al. 2018), In the case of ‘no weathering’, initially (approximately within the first 600 years), the erosion rate reduces and then increases afterwards (Fig. 7). The erosion rate variation of the ‘no weathering’ simulation behaves somewhat differently compared to the Exponential and Humped Exponential pedogenesis function landscapes.

Fig. 7.

Erosion rate variation of the landform under different depth dependent weathering functions.

The rates of erosion for the Exponential and Humped Exponential DDWFs (Fig. 7) have different characteristics. Unlike the ‘no weathering’ scenario, the erosion rate of the Exponential DDWF does not have an initial region where the erosion rate is reduced. It starts from a non-zero value, continues to increase for some time, and stabilises near 2000 years. However, the erosion rate using the Humped Exponential DDWF shows a slight reduction in erosion rate in the initial stage, similar to the ‘no weathering’ simulation results, and then an increase, similar to the Exponential DDWF results. Interestingly after the initial stage, the erosion rate curve of Humped Exponential DDWF is almost parallel to the erosion rate curve of the Exponential DDWF.

Fig. 8 shows the final landforms after 3000 years of evolution without weathering, and the landform develops numerous gullies. Interestingly, the number of gullies observed above the mid-bench area is much higher than those below the mid-bench area. This observation suggests that the mid bench has acted as an obstacle to the gully propagation process.

Fig. 8.

Geomorphology of the landform after 3000 years of evolution without weathering and surface coarseness (represented by d₅₀).

Fig. 9 shows the final landforms after 3000 years of evolution with Exponential DDWF and Humped Exponential DDWF. At first glance, both landforms seem similar, yet subtle differences exist. The most notable difference is that the area below the mid-bench has more gullies at closer spacing using the Humped Exponential weathering function. In contrast, more severe erosion can be seen in the same area using the Exponential DDWF function. Both Exponential and Humped Exponential DDWF functions display higher erosion compared with the results of ‘no weathering’ simulations (Fig. 8).

Fig. 9.

Top: geomorphology of the landform after 3000 years of evolution with (a) Exponential DDWF, (b) Humped Exponential DDWF. Bottom: surface coarseness (represented by d₅₀) of the landform employing the (a) Exponential DDWF after 200 years, (b) Humped Exponential DDWF after 200 years, (c) Exponential DDWF after 3000 years, and (d) Humped Exponential DDWF after 3000 years.

Armouring is well understood to influence erosion (Sharmeen and Willgoose 2006, 2007; Yang et al. 2023). Hancock et al. (2017) found a reduction in both field-measured erosion and that predicted by the SIBERIA LEM at several mine-site erosion plots.

With ‘no weathering’, the initial reduction of the erosion rate can be associated with surface armouring. The armouring process predominantly removes fine material, retaining the coarser material. With time, the surface loses fine material and becomes enriched with coarse material. This is demonstrated in Fig. 9 where the surface coarseness of the resultant landform simulated without weathering is shown. The d₅₀ value of the initial soil grading (i.e. at the beginning of the simulation) is 12.5 mm. Since there is no weathering, deposition is the only process where an area of the landform can develop finer material. All other areas will have surface d₅₀ values greater than the initial value. This occurs as the armouring process removes fines and hence enriches the surface with coarse material. It is noteworthy to mention that the scale of the colour bar in Fig. 8 has a range of 0–30 mm, while the colour bar scale of Fig. 9 showing the surface coarseness of Exponential and Humped Exponential DDWS is 0–15 mm. For the Exponential and Humped Exponential DDWF in Fig. 9, the process of armouring coarsens the surface, while weathering breaks down soil particles, leading to low surface coarseness. Hence, the relative balance of armouring and weathering processes dictates the surface coarseness of any point in the landform.

Models such as CAESAR-Lisflood have the capability of preferentially eroding and depositing both coarse and fine sediment (Coulthard et al. 2013). Since the fine material available for erosion reduces, so does the erosion rate. Also, at the start of hillslope evolution, the landform does not have a developed drainage network. Hence, the drainage network evolves by eroding some landform regions and depositing the material in other areas. During this time, concurrently with armouring, most eroded material is deposited within the landform because the drainage system is not connected with the downstream landform outlets. Once the drainage system matures and connects, sediment eroded upstream can be routed to the hillslope outlet at the base of the slope, leading to increased erosion. This highlights the need for an integrated hillslope and/or catchment approach to predict both surface and subsurface evolution (Walker and Coventry 1976; van der Meij et al. 2018; Willgoose 2018).

The erosion rate difference between the Exponential and the Humped Exponential DDWF simulations can be explained by considering how the two DDWFs influence the weathering process. In the case of Exponential DDWF, the highest weathering rate occurs at the soil surface and, then exponentially reduces with soil depth. Consequently, the surface material weathers more vigorously, producing large amounts of fine particles. The breakdown of the coarse soil particles at the surface prevents armouring while keeping the flow transport limited (the surface retains enough fine particles, but the flow transport capacity is saturated, preventing the erosion of all the fine particles). As described earlier, as the drainage system matures, the flow’s transport capacity also increases, leading to the flow becoming weathering limited (the flow is able to entrain all the fine particles formed during the weathering process), stabilising the erosion rate in the long term. This demonstrates the role of hillslope connectivity in surface evolution (van der Meij et al. 2018; Willgoose 2018).

The surface layer coarseness of the early landform (around 200 years) is relatively high in the Humped Exponential DDWF (Fig. 9b) simulation compared with the Exponential DDWF simulation (Fig. 9a)). The highest weathering rate in the Humped Exponential DDWF is a finite depth below the surface (Humphreys and Wilkinson 2007; Stockmann et al. 2014). Therefore, in the simulation’s initial stage, the surface material’s weathering rate remains low, leading to a relatively strong armouring process. Similar to the erosion rate curve produced by the non-weathering scenario, the strong armouring of the surface at the initial stage reduces erosion. After some time, the surface material is weathered sufficiently, and the drainage system has matured, leading to the observed increase in the erosion rate in the Humped Exponential DDWF weathering rate curve. Following this stage, the erosion rate curve of the Humped Exponential DDWF landscape follows the same pattern as the Exponential DDWF landscape.

Influence of particle breakage geometry

Breakage pattern and geometry have a large influence on erosion and deposition processes, and landscape evolution (Wells et al. 2008; Willgoose 2018). Erosion commences with a non-zero erosion rate, then increases with time and stabilises (Fig. 10). The Granular Disintegration mechanism has the highest overall erosion, followed by Even Breakage, Symmetric Fragmentation, with the Spalling mechanism exhibiting the lowest total erosion.

Fig. 10.

Erosion rate variation of the landform under different weathering geometries.

Comparing Fig. 9 with Fig. 11, it is evident that the weathering mechanism influences landform geomorphology. Qualitatively, it can be observed that when a weathering mechanism produces relatively large daughter particles compared with the parent particle, the resultant landform develops fewer deep gullies and leads to lower erosion rates. Conversely, when a weathering mechanism produces smaller daughter particles, the landform develops a large number of shallow gullies and higher erosion rates (Yang et al. 2023).

Fig. 11.

Top: the landform’s geomorphology after 3000 years of evolution with (a) Spalling, (b) Symmetric Fragmentation, (c) Even Breakage, and (d) Granular Disintegration. Bottom: surface coarseness of the landform (represented by d₅₀) after 3000 years of evolution with (a) Spalling, (b) Symmetric Fragmentation, (c) Even Breakage, and (d) Granular Disintegration.

In all landscapes with weathering, the process of armouring is weakened by the breakdown of the coarse material, and the erosion rate continues to increase without any initial reduction throughout all the weathering simulations (Wells et al. 2008; Willgoose 2018; Yang et al. 2023). The relative size of the daughter particles (at least the largest) compared to the parent particle can explain the different magnitudes of erosion rates observed for each weathering mechanism simulation. From all the weathering functions, spalling-like weathering (where the larger daughter particle retains 90% of the parent particle) has the lowest erosion rate. Although fine particles are produced in the Spalling mechanism, most of the parent particle volume is retained in the large daughter particle. Although the weathering process interferes with the armouring process, the retention of coarse particles lessens the severity of this interference, leading to some armouring and a low erosion rate.

Granular Disintegration weathering, where the entire parent particle disintegrates into 100 daughter particles with 1% of the parent particle volume, exhibits the highest erosion rate. In this case, no coarse material is retained through the process of weathering, whereas the armouring process is severely hindered, leading to the highest erosion rate.

The erosion rate variation employing Symmetric Fragmentation and Even Breakage mechanisms can also be understood by taking the size of the daughter particles into account. In Symmetric Fragmentation, the parent particle breaks into two identical daughter particles, each keeping 50% of the parent volume. Therefore, the resultant daughter particles are relatively coarse, leading to strong armouring and relatively low erosion rates. The Even Breakage mechanism has features from both Symmetric Fragmentation and Granular Disintegration mechanisms. It produces a relatively large number of identical daughter particles, similar to the process of Granular Disintegration. However, the daughter particles’ sizes are larger than those of Symmetric Fragmentation. Therefore, the armouring process is weaker compared to the Symmetric Fragmentation simulation and stronger than the Granular Disintegration simulation. Consequently, Even Breakage exhibits a higher erosion rate than the Symmetric Fragmentation while showing a lower erosion rate than Granular Disintegration.

As discussed earlier, the erosion rate variation of the different weathering mechanisms depends on the competing processes of armouring and weathering (Sharmeen and Willgoose 2006, 2007). Armouring tends to coarsen the surface by preferential removal of fine particles, while weathering produces fine particles through material breakdown. Fig. 11 (bottom) displays the surface coarseness of the final landform after 3000 years under different weathering mechanisms. Spalling produces the landform with the highest surface coarseness, followed by Symmetric Fragmentation and Even Breakage mechanisms. Granular Disintegration produces the landform with the finest soil surface. By comparing the erosion rate (Fig. 10) and the final surface d₅₀, it is clear that a strong correlation exists between the landform’s surface coarseness and rate of erosion. The higher the surface coarseness, the lower the erosion rate of the landform.

It can be observed that Symmetric Fragmentation and Spalling weathering mechanisms produce relatively coarse surface material brought forward by strong armouring. Armouring mainly occurs in the gully bottoms and the gully banks. On the sidewall, the d₅₀ increases toward the gully’s bottom and reaches a maximum at the bottom. This armouring effect is much more substantial for Spalling as this process produces coarser daughter particles than the Symmetric Fragmentation process. Due to armouring, the erosion rate of Symmetric Fragmentation and Spalling weathering mechanism landscapes are lower compared to the other two simulations.

Employing Even Breakage and Granular Disintegration, the armouring effect is weak, leading to high erosion. For these two landscapes, armouring mostly occurs within the gully bottoms. Interestingly, the d₅₀ of gully bottoms of Granular Disintegration landscape appears higher than that of Even Breakage. This observation is somewhat counterintuitive as Granular Disintegration produces daughter products smaller than the Even Breakage mechanism.

However, with Even Breakage, some daughter products remain on the surface of the gully bottoms as armour since they are larger than the daughter products produced by Granular Disintegration. Since the daughter products of the Granular Disintegration mechanism are smaller, all the daughter products are readily eroded from the surface of gully bottoms. This process exposes the unweathered parent particles, which leads to a higher apparent d₅₀ in the gully bottoms.

The simulation results suggest that a fragmentation mechanism producing a mixture of particles of different sizes (or at least some large particles in the case of Symmetric Fragmentation) will lead to armouring and reduced erosion over time (Parker and Klingeman 1982; Sharmeen and Willgoose 2006, 2007). Also, if the breakage mechanism is known, then the likelihood of gullying (as opposed to erosion processes such as sheetwash) can be better predicted (Kirkby 1971).

Influence of the largest daughter particle in the weathering mechanism

The findings here demonstrate that the size of the weathered particle has a direct influence on erosion (Sharmeen and Willgoose 2006, 2007; van der Meij et al. 2018; Willgoose 2018). Fig. 12 demonstrates that the erosion rate increases as the α parameter of the Spalling weathering mechanism decreases. Interestingly, the difference between each erosion rate curve reduces as the α parameter is reduced. The two reasons for this are: (1) as the size of the largest daughter particle decreases, the armouring process is weakened, leading to increased erosion; and (2) although the size of the largest daughter particle is reduced when the α parameter is reduced, the sizes of other daughter particles increase. This size increase positively enhances the influence of armouring. Hence, although the erosion rate increases with a decreasing α parameter, the difference between the Spalling parameters visually reduces (Fig. 12).

Fig. 12.

Temporal variation in erosion rate using different Spalling weathering mechanism parameters.

The change in particle size has long-term consequences on landform evolution and the surface coarseness of the resultant landform. Fig. 13 displays the difference in surface morphology for the landforms employing different α values in the Spalling weathering mechanism. Comparing the morphology of the two landforms, it is evident that the landform developed through a lower α parameter has a higher number of gullies than the simulation employing a high α parameter. Also, the surface material is coarser in the high α parameter simulation compared to the low α parameter simulation. As described earlier, the different rates of armouring occurring in the landform due to different size daughter products produced by different weathering mechanisms are the reason behind this geomorphological difference. Hancock et al. (2017), when using the SIBERIA LEM and different armouring model inputs, found that geological timescale landscapes varied.

Fig. 13.

Geomorphology comparison between different parameters for the Spalling weathering mechanism, (a) α = 0.9 and (b) α = 0.4 (top). Surface coarseness comparison between different parameters for the Spalling weathering mechanism, (a) α = 0.9 and (b) α = 0.4 (bottom).

Modelling outcomes and future work

Using a LEM with a pedogenesis capability highlights the importance of quantitatively understanding both surface and subsurface landform evolution (Evans 2000; Sheridan et al. 2000; van der Meij et al. 2018). The findings demonstrate that different weathering and pedogenesis functions produce different erosion rates and that at centennial to millennial time scales different landforms can result. The results here at the hillslope scale supports the previous findings by Humphreys and Wilkinson (2007). The model can be easily employed to provide possible landscape outcome not just at the hillslope scale, but also at the catchment scale (Willgoose 2018).

A strength of this work is that the surface erosion parameters have been determined for the material, with a reasonable level of confidence in their reliability over a 10–20-year period, as demonstrated by Welivitiya and Hancock (2022). However, it remains uncertain which weathering and pedogenesis model is the most accurate. Field and laboratory data from both the surface and subsurface, combined with modelling efforts, are necessary to constrain model predictions (Legros and Pedro 1985; Wells et al. 2008; Sanchidrián et al. 2014). Ongoing fieldwork aims to determine weathering and armouring pathways on post-mining landforms.

The landscape employed here is a representative post-mining landform. The materials and parameters derived are sourced from a Queensland coal mine (Welivitiya and Hancock 2022). However, other landforms with different materials, hillslopes, slopes and climate are likely to perform differently (Simmons and Fityus 2016; Yang et al. 2023). The SSSPAM model provides the ability to assess potential landform influences and optimise designs. Other processes such as chemical weathering (and dissolution) (Willgoose 2018) can be examined in the SSSPAM modelling framework. Alternatively, if the weathering process is dominated by a chemical process, a model such as Lorica may be appropriate (Temme and Vanwalleghem 2016). However, assessment of the material here found that chemical weathering or dissolution was not a major process of material breakdown (Hazelton and Murphy 2016; Welivitiya and Hancock 2022).

The study site was ideal as it became a ‘badland’ with no vegetation. The material is highly erodible. Hence, the lack of vegetation and its influence over pedogenesis in the short-term can be neglected. However, the role of vegetation, and its influence on pedogenesis, cannot be ignored (Barzegar et al. 2002; Peterman et al. 2014; Ritschel and Totsche 2019). There are no LEMs at present that have a coupled erosion, deposition and vegetation growth model that actively evolves a landform. Some models such as CAESAR-Lisflood (Coulthard et al. 2013) have a model that evolves vegetation to climax cover but it does not respond to changes in rainfall once fully evolved, nor does it have a subsurface-vegetation component. The SIBERIA LEM integrates vegetation by assuming a set of representative erosion and hydrology parameters determined for the surface. It does not dynamically respond to vegetation growth or decay. Nor does it have a model for below ground biomass. However, the modelling here was for an approximately 20-year-old landform with an absence of vegetation (Fig. 1). A complete vegetation and decadal to centennial time scales is unlikely without human intervention. Hence, the modelling here assuming no vegetation provides somewhat series of plausible landscape outcomes.

An issue not highlighted here is the differences in erosion rates. Employing ‘No Weathering’ parameters produces the lowest erosion rate of <0.5 mm per year. This corresponds to an erosion rate of 6.5 t ha⁻¹ per year, assuming a bulk density of 1.3 tonnes m⁻³ for the soil. An erosion rate of <5 t ha⁻¹ per year approximates that of agricultural erosion rates for many Australian landscapes (see Bui et al. (2011) for a review of tolerable soil erosion rates in Australia). However, for all other parameter sets, the rate of erosion is much higher. Therefore, using the ‘no weathering’ parameters may underestimate erosion and highlights the need to better constrain pedogenesis and weathering processes.

In particular for coal mines, waste rock has been shown to rapidly weather with topsoil and vegetation and or competent rock needed to reduce erosion (Welivitiya and Hancock 2022). Field and model data from other sites, have also found that erosion rates are reduced by the armouring as well as the inclusion of surface roughness by ripping (Saynor et al. 2019). Other models such as SIBERIA (Willgoose 2018) and CAESAR Lisflood (Coulthard et al. 2013) can model the effect of armouring.

As described above, for simplicity, the modelling has assumed no vegetation. Vegetation is an essential component of successful, long-term mine reconstruction. However, the establishment of a consistent long-term vegetation cover (in this case, grazing pasture) is yet to be achieved for this site. It should be recognised that an absence of vegetation is not an acceptable post-mining outcome. It also should be noted that for many mine sites vegetation will disturbed by both grazing and fire (McKenna et al. 2024). Further, the pedogenesis functions assumed to be applicable here (Fig. 2) have been developed from natural (non-disturbed) landscapes, assuming vegetation is a component in pedogenesis. Hence, correctly or otherwise, vegetation is incorporated into the model framework for the subsurface used here.

Post-mining landforms will exist for geological time once they are constructed. The pedogenic trajectory for post-mining landforms is unknown (Emmerton 2019; Cox et al. 2021). If the pedogenic potential of mine overburden (or waste rock) as it is being removed is known, then there is potential to optimise material placement (Slukovskaya et al. 2019) such that the best soil development material can be placed at the surface with the least amenable material placed at depth. The model and process presented here demonstrate that we possess the technology to forecast landscape evolution. With knowledge and planning there is the potential to reconstruct a landform that has at least equal and potentially greater ecological function than that of the pre-mine landscape.