6 Verification of the Reactive Package within RUMT3D and Plausibility Tests
6.1 Acid Mine Drainage (AMD)
6.1.3 Setup of the Benchmark Problem and Model Comparison
To achieve accurate and numerically stable reactive simulation results using RUMT3D, this subsection examines appropriate spatial and temporal discretisations required for the present benchmark problem. Simulation results of RUMT3D are thereby compared with the PHREEQC-2 model. At the end of this subsection, once appropriate discretisations are found, minor differences in results between these two models and a non-physical effect occurring with both models are discussed.
Spatial scale. Initial tests for numerical stability, such as simulating a tracer breakthrough (chloride) revealed that a discretisation of 10 m per cell causes significant numerical dispersion when using RUMT3D/MT3DMS’s numerical solver TVD (third-order total-variation-diminishing scheme) for the simulation of advective transport. To ensure that TVD would not only handle the discretisation chosen but also different concentrations of the recharge, tailing and aquifer water, the following initial tracer concentrations were assumed: 0.01, 0.005, 0.001 mol L-1 in the recharge, tailing and aquifer solution, respectively. Fig. 10 illustrates the simulated tracer concentrations after one year for both, the RUMT3D and the PHREEQC-2 model in comparison for different spatial discretisations of 10, 2 and 1 m per cell, respectively. As the different solutions have different initial concentrations, two moving fronts can be seen in Fig. 10: the first front starts at the beginning of the tailing with the concentration of the recharge solution (0.01 mol L-1) and stops at the end of tailing (at 20 m) with the concentration of the tailing solution (0.005 mol L-1); the second front begins at the end of the tailing/beginning of the aquifer with the concentration of the tailing solution (0.005 mol L-1) and stops 20 m further downgradient in the aquifer with the concentration of the aquifer solution (0.001 mol L-1). Fig. 10 demonstrates that for conservative transport of chloride, the PHREEQC-2 model shows less numerical dispersion in comparison with RUMT3D for the simulations that use the coarse discretisation of 10 m per cell. However, significant numerical
H
+dispersion (flatter curve) becomes also apparent for PHREEQC-2 in the reactive case, especially for the early parts of the simulation when the most significant geochemical changes occur (Fig. 11). Fig. 11 shows the simulated calcium (a-d Figs.) and pH (e-h Figs.) profiles after 0, 1, 2 and 5 years for different discretisations of 10 and 1 m per cell, respectively. Though, as for the non-reactive case, the magnitude of numerical dispersion of PHREEQC-2 for the reactive case compared to the one of RUMT3D is not as large. When recharge solution reaches the right hand boundary of the aquifer after approximately 10 years, numerical dispersion becomes negligible.
For a discretisation of 1 m per cell, the simulation results of RUMT3D and PHREEQC-2 are in good agreement. Fig. 10 indicates that a discretisation of 2 m per cell would also be sufficient to generate satisfactory results. Nonetheless, since subsequently a conduit system (which might trigger numerical instability) is introduced into the original McNab problem, a discretisation of 1 m per cell was selected for the following investigations.
Fig. 10: Comparison of simulated tracer profiles of the McNab problem with the RUMT3D and PHREEQC-2 models using different spatial discretisations of 10, 2 and 1 m per cell after 1 year. The recharge, tailing and aquifer solution had initial concentrations of 0.01, 0.005, 0.001 mol L-1, respectively.
0 0.002 0.004 0.006 0.008 0.01
0 10 20 30 40 50
Distance downgradient (m)
Tracer concentration (mol L-1 ) 10 m - PHREEQC-2
2 m - PHREEQC-2 1 m - PHREEQC-2 10 m - RUMT3D 2 m - RUMT3D 1 m - RUMT3D
1 year
a) Calcium profile at the beginning of the simulation.
b) Simulated calcium profile after 1 year.
0 0.00005 0.0001 0.00015
0 20 40 60 80 100
Distance downgradient (m) Calcium (mol L-1 )
10 m - PHREEQC-2 1 m - PHREEQC-2 10 m - RUMT3D 1 m - RUMT3D 10 m
1 m
0 year
0 0.0005 0.001 0.0015 0.002 0.0025
0 20 40 60 80 100
Distance downgradient (m) Calcium (mol L-1 )
10 m - PHREEQC-2 1 m - PHREEQC-2 10 m - RUMT3D 1 m - RUMT3D
10 m - RUMT3D
0 0.0002 0.0004
10 20 30
10 m - PHREEQC-2
1 m 1 year
c) Simulated calcium profile after 2 years.
d) Simulated calcium profile after 5 years.
0 0.005 0.01 0.015
0 20 40 60 80 100
Distance downgradient (m) Calcium (mol L-1 )
10 m - PHREEQC-2 1 m - PHREEQC-2 10 m - RUMT3D 1 m - RUMT3D
10 m - RUMT3D
10 m - PHREEQC-2 1 m
2 years
0 0.005 0.01 0.015 0.02
0 20 40 60 80 100
Distance downgradient (m) Calcium (mol L-1 )
10 m - PHREEQC-2 1 m - PHREEQC-2 10 m - RUMT3D 1 m - RUMT3D
10 m - RUMT3D 10 m -
PHREEQC-2 1 m
5 years
e) pH profile at the beginning of the simulation.
f) Simulated pH profile after 1 year.
7 8 9 10
0 20 40 60 80 100
Distance downgradient (m)
pH
10 m - PHREEQC-2 1 m - PHREEQC-2 10 m - RUMT3D 1 m - RUMT3D 10 m
1 m
0 year
1.5 4.5 7.5 10.5
0 20 40 60 80 100
Distance downgradient (m)
pH
10 m - PHREEQC-2 1 m - PHREEQC-2 10 m - RUMT3D 1 m - RUMT3D 10 m - RUMT3D
10 m - PHREEQC-2 1 m
1 year
g) Simulated pH profile after 2 years.
h) Simulated pH profile after 5 years.
Fig. 11: Comparison of simulated calcium (a, c, e, g) and pH (b, d, f, h) profiles after 0, 1, 2 and 5 years, respectively, for the McNab problem using different discretisations of 10, 2 and 1 m per cell. Results are compared to the PHREEQC-2 model.
1.5 4.5 7.5 10.5
0 20 40 60 80 100
Distance downgradient (m)
pH
10 m - PHREEQC-2 1 m - PHREEQC-2 10 m - RUMT3D 1 m - RUMT3D 10 m -
RUMT3D 10 m -
PHREEQC-2
1 m
2 years
1.5 4.5 7.5 10.5
0 20 40 60 80 100
Distance downgradient (m)
pH
10 m - PHREEQC-2 1 m - PHREEQC-2 10 m - RUMT3D 1 m - RUMT3D
10 m - RUMT3D 10 m - PHREEQC-2
1 m
5 years
Reaction time step frequency. In RUMT3D, the number of transport but not the number of reaction time steps (see Chapters 2.5.1 and 3.3) is automatically controlled such that the appropriate stability criteria for physical transport are satisfied (Zheng & Wang 1999; Prommer et al. 2003). Therefore, the influence of the selected reaction time step length in RUMT3D was further investigated for the present simulation problem. For cases like the McNab problem, where all reactions are assumed to be in equilibrium, any error that occurs as a result of the (user-)selected reaction time-step length is based on temporal operator splitting errors (usual OS error, Walter et al. 1994; Steefel & MacQuarrie 1996), whereas in the case of kinetic reactions the solution accuracy might also depend, e.g., on the integration method used for the kinetic reactions (Parkhurst & Appelo 1999). It should be noted that within an investigated range of reaction step frequencies, an increase in the number of reaction time steps does not necessarily improve the accuracy of results once an inappropriate spatial discretisation is selected (Fig. 12). Fig. 12 displays the simulated calcium concentration profiles for RUMT3D after 1 year using a discretisation of 10 m per cell and reaction time step sizes of 100 and 1000 per year. Results obtained with the PHREEQC-2 model are shown in comparison.
When comparing simulated concentration profiles (after 1, 2, 5, 10, 15 and 20 years, respectively), a good agreement of results can be observed for both, the PHREEQC-2 and the RUMT3D solutions for all investigated cases with the following selection of reaction time steps per year: 100, 150 and 200 if a spatial discretisation of 1 m per cell is selected. Fig. 13a-d show the comparison of results for pH, sulphate, Fe(II) and goethite profiles after 20 years using either 10, 100 or 200 reaction time steps per year. Results from the PHREEQC-2 model are also shown in comparison. Based upon the good agreements of the RUMT3D model using 100, 150 and 200 reaction time steps per year with the PHREEQC-2 model, for all simulations in this section, 100 reaction time steps per year were used.
Fig. 12: Comparison of simulated calcium profiles with the PHREEQC and RUMT3D models after 1 year, for the McNab problem using discretisation of 10 m per cell. Two different reaction time step (ts) sizes were used for the RUMT3D model per year (100 and 1000).
a) Simulated pH profile.
0 0.0005 0.001 0.0015 0.002 0.0025
0 20 40 60 80 100
Distance downgradient (m) Calcium (mol L-1 )
10 m - PHREEQC-2 10 m - RUMT3D - 100ts 10 m - RUMT3D - 1000ts 100ts,
1000ts
PHREEQC-2
1 year
1 3 5 7
0 20 40 60 80 100
Distance downgradient (m)
pH
PHREEQC-2 RUMT3D - 10ts RUMT3D - 100ts RUMT3D - 200ts
1 3 5 7
17 20 23 26
10ts 20 years
b) Simulated sulphate profile.
c) Simulated Fe(II) profile.
0.014 0.018 0.022 0.026 0.03
0 20 40 60 80 100
Distance downgradient (m) Sulphate (mol L-1 )
PHREEQC-2 RUMT3D - 10ts RUMT3D - 100ts RUMT3D - 200ts
0.014 0.022 0.03
17 20 23 26
100ts, 200ts 10ts PHREEQC-2