Forecast of AMD Quantity by a Series Tank Model in Three Stages: Case Studies in Two Closed Japanese Mines

Abstract

There are about 100 sites of acid mine drainage (AMD) from abandoned/closed mines in Japan. For their sustainable treatment, future prediction of AMD quantity is crucial. In this study, AMD quantity was predicted for two closed mines in Japan based on a series tank model in three stages. The tank model parameters were determined from the relationship between the observed AMD quantity and the inflow of rainfall and snowmelt by using the Kalman filter and particle swarm optimization methods. The Automated Meteorological Data Acquisition System (AMeDAS) data of rainfall were corrected for elevation and by the statistical daily fluctuation model. The snowmelt was estimated from the AMeDAS data of rainfall, temperature, and sunshine duration by using mass and heat balance of snow. Fitting with one year of daily data was sufficient to obtain the AMD quantity model. Future AMD quantity was predicted by the constructed model using the forecast data of rainfall and temperature proposed by the Max Planck Institute–Earth System Model (MPI–ESM), based on the Intergovernmental Panel on Climate Change (IPCC) representative concentration pathway (RCP) 2.6 and RCP8.5 scenarios. The results showed that global warming causes an increase in the quantity and fluctuation of AMD, especially for large reservoirs and residence time of AMD. There is a concern that for mines with large AMD quantities, AMD treatment will be unstable due to future global warming.

1. Introduction

Japan has more than 5000 abandoned/closed mines, and about 100 of their sites produce acid mine drainage (AMD) due to the presence of sulfide mineralization [1]. The general treatment for AMD is neutralization and sedimentation by addition of a neutralizer, such as lime, calcium carbonate, and sodium hydroxide [2], and solid/liquid separation [3] of the produced sludge from the neutralized effluents. In these treatments, all toxic elements are concentrated into the sludge by precipitation [4,5,6] and adsorption [7,8,9,10,11,12], and the sludge is controlled in a tailing pond at a mine site or final disposal site. For these last several decades, AMD has been treated properly in Japan and has not caused severe pollution. However, since our results of the statistical calculation (details are shown below) suggested that some mines have required AMD treatment for over 150 years [13,14] and other groups suggested that more than 1000 years of treatment will be necessary [15] in the current situation, more sustainable treatment to reduce both AMD generation [16,17] and treatment cost [18] is needed. To reduce the treatment cost of the addition of chemicals and of sludge generation, for example, a passive treatment that utilizes the natural environment of mines, such as topography, plants, and microorganisms, has attracted attention as a sustainable AMD treatment based on new concepts [19,20,21]. Several researchers are trying to successfully reuse this sludge as an industrial material [22,23].

Unlike industrial wastewater, AMD quantity and quality differ significantly in mines due to regional, geological, mineralogical, and biological factors. Therefore, it is necessary to customize an appropriate treatment for each mine. To select an optimal treatment method from the various treatment technologies, including those based on both active and passive concepts, an accurate understanding of the current potential for AMD generation [24,25] and the future forecast of AMD quantity and quality are essential.

The objective of this study was to determine a forecast for AMD quantity. To accomplish this, we constructed a model that reproduces the current AMD quantity using previous monitoring data of AMD, and then extrapolated it to the future. To reproduce AMD quantity, there are two ways: one is a hydraulic simulation [26,27,28,29,30,31] and the other is a tank model. A hydraulic simulation provides detailed information on the origin and distribution of AMD and can be a powerful tool for discussing AMD generation countermeasures, but it requires, in addition to meteorological data, detailed geological, mineralogical, and hydraulic data, which are generally difficult to obtain, especially for abandoned and closed mines. On the other hand, a tank model is a blackbox to determine the relationship between inflow and outflow, but just inflow data of rainfall and snowmelt and outflow of AMD quantity are sufficient for the model [32,33]. In this study, the tank model was selected for AMD quantity modeling, and rainfall and snowmelt were used as inflow. The rainfall data were corrected for elevation and adjusted using the statistical daily fluctuation model to suit each AMD site. Snowmelt was also estimated from rainfall by considering mass and heat balance by using temperature and sunshine duration data. We did not select a hydraulic simulation but chose a statistical model because our target mines are closed and it was difficult to obtain detailed monitoring, geological, and hydraulic data for this study.

For the AMD quality model, we previously reported the geochemical calculation with first-order elution kinetics of sulfide minerals [13,14]. In the model, sulfides that should be the source of AMD were selected from the quality data, and their first-order elution rate and initial AMD generation potential were estimated by fitting to the time change of their elution amount obtained from the AMD quantity and quality data. The AMD quality could be estimated by the coupling of the kinetics for sulfide elution and oxidation, and the geochemical code for the chemical equilibrium calculation of precipitation and adsorption [34,35,36,37]. This means that accurate estimation of AMD quantity is crucial for the AMD quality model.

In this study, the AMD quantity model was constructed using two case studies of underground mines: a closed sulfur mine (Mine A) and a closed black-ore copper, lead, and zinc mine (Mine B). Mine A has a large quantity of AMD, averaging 18 m³ min⁻¹ with a small fluctuation, which is opposite to that from Mine B (1.5 m³ min⁻¹). From these case studies, the parameters of the model were estimated and the future AMD quantity for the next few decades was predicted using forecast data for rainfall and temperature based on the MPI–SEM (Max Planck Institute–Earth System Model) [37]. For this, we selected two kinds of global warming scenarios proposed by IPCC (Intergovernmental Panel on Climate Change): a low-stabilization scenario of RCP (representative concentration pathway) 2.6 and a high-level greenhouse gas emission scenario of RCP 8.5 [37]. We further discuss the effects of global warming on the forecast of AMD quantity stemming from the closed sulfide mines that were examined.

2. Materials and Methods

2.1. AMD Quantity Model

2.1.1. Tank Model

The AMD quantity model was constructed by using a series tank model in three stages, as shown in Figure 1. The first and second stages correspond to the nonpolluted water reservoir on the surface and inside the mine, respectively. The third stage corresponds to the polluted water reservoir in the ore deposit that causes AMD. Inflow, r (mm), is the summation of rainfall, r_w, and snowmelt, r_s. In each tank, a part of the inflow is distributed to the outflow (mm h⁻¹), q_oi (i = 1, 2, 3), and seepage flow to the next tank, q_si (i = 1, 2), according to the water reservoir height (mm), x_i (i = 1, 2, 3) and outflow height (mm), b_i (i = 1, 2). The water balance in each tank is as follows:

dx_i/dt = q_si–1 − q_si − q_oi, (i = 1, 2, 3)

(1)

where t is time, q_s0 = r = r_w + r_s, and q_s3 = 0. The outflow is calculated from:

q_oi = a_oi (x_i − b_i), (i = 1, 2, 3)

(2)

where a_oi (i = 1, 2, 3) is the outflow coefficient and b₃ = 0. The seepage flow is also calculated from:

q_si = a_si x_i, (i = 1, 2)

(3)

where a_si (i = 1, 2) is the seepage coefficient.

In this study, q_o3 corresponded to AMD quantity. The inflow, r, was set using the following procedure. The Kalman filter and particle swarm optimization methods were used for fitting q_o3 to the observed data of AMD quantity to estimate the x_i, b_i, a_si, and a_oi parameters [38].

2.1.2. Correction of Rainfall Data and Judgment of Snowfall

Rainfall data near each mine were derived from AMeDAS (Automated Meteorological Data Acquisition System) provided by the Japan Meteorological Agency [39]. To suit each mine situation, the daily data of rainfall and temperature obtained from AMeDAS were corrected for elevation and adjusted by using the statistical daily fluctuation model, as shown in Figure 2. Also, snowfall was estimated according to the corrected temperature, and if the temperature was under 2 °C, it was judged to be snowfall and not rainfall.

Since the capture rate of rainfall particles by the rain gauge decreases as the wind speed increases, known as the Jebons effect, Kondo et al. proposed the following statistical correction for daily rainfall data obtained from AMeDAS [40]:

r_w1 = (1.25 + 0.15 cos(w(d − 20)))r_w0

(4)

w = 2π/365

(5)

where r_w0 (mm) is the raw rainfall data, r_w1 (mm) is the corrected rainfall data by the model, d (days) refers to the days from 1 January. Furthermore, the amount of rainfall near mines depends greatly on elevation due to rapid updraft and adiabatic expansion, and is corrected as follows [40]:

r_w2 = (1 + c(h − h₀))r_w1

(6)

where r_w2 (mm) is the corrected rainfall data by elevation, h (km) is the elevation of mine site, h₀ (km) is the elevation of the AMeDAS observation point, and c is the coefficient (0.001 km⁻¹ for 5 °C or less and 0.00064 km⁻¹ for 5 °C or more). AMeDAS data of temperature, T, were also corrected according to following equation:

T₁ = −6(h − h₀) + T

(7)

where T₁ (°C) is the corrected temperature.

As shown in Figure 2, if the temperature at the mine site was above 2 °C, it was assumed that there was no snowfall and the inflow of rainfall was set as r_w2. On the other hand, if the temperature was below 2 °C, the rainfall data, r_w2, were judged to be equivalent to snowfall, r_f = r_w2, and r_w2 = 0.

2.1.3. Estimation of Snowmelt and Snow Cover

Daily snowmelt, r_s, and snow cover, r_c, were also estimated by following the mass and heat balance of snow using the AMeDAS data of rainfall, temperature, and sunshine duration, as shown in Figure 2. If the temperature was below 0 °C, no snowmelt was assumed and r_s = 0; otherwise, snowmelt was calculated according to the following procedure.

Snowmelt, r_s (mm), was calculated from the ratio of fusion heat, Q (J), and latent heat, L (334 J kg⁻¹), the catchment area, S (mm²), and the density of water, ρ (9.97 × 10⁻⁷ kg mm⁻³):

r_sSρ = Q/L

(8)

The fusion heat was calculated from the heat balance of snow exposure as follows:

Q = Q₁ + Q₂ + Q₃ + Q₄ + Q₅

(9)

where Q₁ is the short-wavelength radiation, Q₂ is the long-wavelength radiation, Q₃ is the sensible heat transfer, Q₄ is the latent heat transfer, and Q₅ is the transfer heat due to rainfall. Here, heat changes in the snow layer and heat transfer from the ground were assumed to be negligible [41].

The short-wavelength radiation was calculated from albedo, r, which is the ratio of reflected sunshine radiation to sunshine radiation on the earth’s surface and the daily average of solar irradiance, I (W m⁻²):

Q₁ = (1 − r)I

(10)

The average of solar irradiance was the function of the ratio of sunshine duration, N, and astronomical sunshine duration, N₀ [42]:

I/I₀ = 0.179 + 0.550N/N₀, for 0 ≤ N ≤ N₀

(11)

I/I₀ = 0.114, for N = 0

(12)

where I₀ is the solar irradiance at the top of atmosphere. The values of solar irradiance at the top of the atmosphere, I₀, astronomical sunshine duration, N₀, and albedo, r, are available from references [39,43], and sunshine duration data, N, are available from AMeDAS.

The long-wavelength radiation was the difference between the radiation from the atmosphere, Q_a, and the radiation from the snow surface, Q_s:

Q₂ = Q_a − Q_s

(13)

Q_a = σ(T₁ + 273.15)⁴(0.605 + 0.048e^0.5)

(14)

Q_s = 0.9σ(T_s + 273.15)⁴

(15)

e = 6.1078 × 10^{7.5T/(T+273.3)}

(16)

where σ is the Stefan–Boltzmann constant (5.67 × 10⁻⁸ W m⁻² K⁻⁴) and e is amount of saturated water vapor [44]. The temperature of the snow surface, T_s, was calculated from [45]:

T_s = 1.13T₁ − 1.67

(17)

when T₁ ≤ 1.47.

The sensible heat was calculated from [46]:

Q₃ = K(1 − 0.0065h/(T₁ + 273.15+0.0065h))^5.257, for T₁ ≥ 0,

(18)

Q₃ = 0, for T₁ ≤ 0

(19)

where K is the transfer coefficient of the sensible heat and the latent heat; 3.5 was proposed for the area near the mines that were modeled in this case study [46]. The latent heat was calculated from:

Q₄ = 1.53K(e − 6.11), for T₁ ≥ 7

(20)

Q₄ = 0, for T₁ ≤ 7

(21)

The transfer heat due to rainfall was calculated from:

Q₅ = ρc_w(273.15 + T₁)r_wS

(22)

where c_w is the specific heat of the water (4.186 J kg⁻¹ K⁻¹).

The snowmelt, r_s, calculated from Equations (8) to (22), should be less than the snow cover, r_c. Snow cover was calculated from following summation of daily mass balance:

r_c = ∑r_f − ∑r_s

(23)

if r_s > r_c, then snowmelt should be r_s = r_c.

2.2. Forecast Data of Temperature, Rainfall, and Sunshine Duration

In the above-mentioned AMD quantity and quality models, daily data of rainfall, average temperature, and sunshine duration obtained from AMeDAS were used for model construction. Therefore, future forecasts of daily data of rainfall, average temperature, and sunshine duration were also necessary for the forecast of AMD quantity and quality in the future. For the AMD quality model, the geochemical calculation with first-order elution kinetics of sulfide minerals was used as mentioned above; the first-order elution rate and initial AMD generation potential were estimated by fitting to the time change of their elution amount obtained from the AMD quantity and quality data. The daily output data of MPI–ESM (Max Planck Institute–Earth System Model), which is an earth system model proposed by Max Planck Institute, were used in this study. Two kinds of IPCC RCPs for the greenhouse gas (GHG) concentration scenario were selected: RCP2.6 and RCP8.5. The former is the scenario with the lowest GHG emission to keep future temperature rise below 2 °C, and the latter is the scenario with the highest GHG emissions.

From the system, daily forecast data of rainfall, average temperature, and maximum and minimum temperatures were available. The daily forecast of the average of solar irradiance, I, was estimated from [47]:

I = 0.76I₀(1− exp(−AΔT^2.2))

(24)

A = 0.036exp(−0.154ΔT_ave)

(25)

where ΔT is the difference between the maximum temperature and the minimum temperature, and ΔT_ave is the monthly average of ΔT.

2.3. Case Studies in Two Closed Mines

In this study, two closed underground mines in the northern part of Japan were selected as a case study; Mine A has a large quantity of AMD, with small fluctuation of quantity and quality, and Mine B has a small quantity of AMD, with large fluctuation of quantity and quality. Snow is observed in winter at both of the closed mines. The locations of the mines are shown in Figure 3, and the AMD characteristics are shown in

Table 1. The AMD characteristics of Mine A and B.

Items	Mine A		Mine B		Effluent
	April 1982	March 2017	April 1972	March 2017	Standard
pH	2.0	2.3	3.2	4.7	5.8–8.6
Fe	547 mg L−1	178 mg L−1	N.D.	N.D.	10 mg L−1
As	3.3 mg L−1	0.93 mg L−1	N.D.	N.D.	0.01 mg L−1
Cd	N.D.	N.D.	0.360 mg L−1	0.024 mg L−1	0.03 mg L−1
Zn	N.D.	N.D.	90.4 mg L−1	14.9 mg L−1	2.0 mg L−1
Pb	N.D.	N.D.	0.700 mg L−1	0.117 mg L−1	0.1 mg L−1
Q	17.6 m3/min	14.8 m3 min−1	2.50 m3 min−1	0.77 m3 min−1

N.D.: Not detected; Q: Quantity.

.

In Mine A, native sulfur and pyrite were mined during its operation. The size of the ore deposit was about 1500 m East–West, about 1500 m North–South, and 25–150 m of thickness, and ore reserves were about 230 million tons per year. The mine produced about 1 million tons of ore and one-third of Japan’s sulfur demand, but it closed in 1971 due to the market influence of sulfur recovered from oil refining. The quantity of AMD is about 18 m³ min⁻¹ on average annually, which is one of the largest AMD values in Japan [48]. The AMeDAS point is located 11 km east and 825 m below the mine.

In Mine B, copper, lead, and zinc were mined during its operation. The ore deposit was a black ore type, which has changed from lower to yellow ore, black ore, and quartz band. The mine produced a maximum of about 25,000 tons per year, but closed in 1985 due to ore depletion. The annual average of the AMD quantity was 1.72 m³ min⁻¹ in 2017 and increased from 5 to 7 m³ min⁻¹ during the snowmelt season [49]. The AMeDAS point is located 18 km northwest and 465 m below the mine.

3. Results and Discussion

3.1. AMD Quantity Model Construction

The relation between the input data of AMeDAS rainfall (upper side) and the observed data and the calculated value of AMD quantity (underside) are shown in Figure 4. In this calculation, the daily observed data of AMD quantity in the prior one year were used for fitting, and later one-year data were used for the model validation. The fitting period was also changed from half a year to two years and the correlation coefficients between the observed and calculated values were compared, as shown in

Table 2. Relationship between the fitting period and the correlation coefficients.

Mine	Fitting Period	Correlation Coefficient in the Fitting Period	Correlation Coefficient in the Validation Period 7
A	Half-year 1	0.89	0.65
	One year 2	0.81	0.69
	Two years 3	0.80	0.67
B	Half-year 4	0.86	0.77
	One year 5	0.86	0.83
	Two years 6	0.87	0.85

¹ 2015/10/1–2016/3/31; ² 2015/4/1–2016/3/31; ³ 2014/4/1–2016/3/31; ⁴ 2016/10/1–2017/4/1; ⁵ 2016/4/1–2017/3/31; ⁶ 2015/4/1–2017/3/1; ⁷ 2016/4/1–2017/3/31 for Mine A; 2017/4/1–2018/3/31 for Mine B.

and Supplementary Figure S1. Of course, the longer the fitting period, the higher the correlation coefficient in the validation period, but a fitting period of one year seemed to be generally sufficient for the reproduction of AMD quantity in the next one year. As shown in Supplementary Figure S2, when the correction for elevation and the statistical daily fluctuation and the snowmelt estimation were not conducted, the reproducibility of AMD quantity became worse, especially for Mine A. This is because Mine A is located at a higher elevation and the effects of elevation correction and snowfall are larger than for Mine B.

The parameters obtained for the tank model are shown in

Table 3. Parameters of the tank model obtained for Mines A and B.

Mine	Tank Stage	Outflow Coefficient ao (day−1)	Seepage Coefficient as (day−1)	Outflow Height b (mm)
A	First stage	0.768	0.568	0.00484
	Second stage	0.936	0.352	0.00140
	Third stage	0.00282
B	First stage	0.738	0.863	0.0358
	Second stage	0.292	0.607	0.0205
	Third stage	0.0500

. Mine A has a smaller outflow height and a larger AMD reservoir than Mine B. Additionally, Mine A has the smaller value of outflow coefficient in the third stage, which directly affects AMD generation, compared to Mine B. This trend means that Mine A has the bigger reservoir and the longer residence time of AMD, which resulted in the smaller fluctuation of AMD, compared to Mine B. As we mentioned in the previous section, ore production was 230 million tons per year in Mine A and 25,000 tons per year in Mine B. This difference in scale should directly affect the difference in reservoir and residence time of AMD.

3.2. Forecast of AMD Quantity

The forecast of AMD quantity (underside) is shown in Figure 5, with the forecast of rainfall (upperside) proposed by MPI–ESM. According to MPI–ESM in RCP2.6, the forecast for temperature rise around the mines is about 2 °C by 2050, but remaining at about 1.0–2.5 °C after 2050. On the other hand, in RCP8.5, the forecast for temperature continues to rise and reaches +5.9 °C in 2100.

In Mine A, the MPI–ESM shows that both rainfall and heavy rain frequency, which is the number of days per year with greater than 50 mm of rainfall, increase due to the temperature rise. In 2100, the forecast of rainfall is +21% for both RCP2.6 and RCP8.5, and forecast of heavy rain frequency increases 4 days/year for RCP2.6 and 8 days/year for RCP8.5, compared to the present. According to these trends, the AMD quantity calculated from the constructed model increases, as shown in Figure 5. The forecast for AMD quantity in 2100 is +27% for RCP2.6 and +31% for RCP8.5.

In Mine B, the MPI–ESM shows that the temperature rise of around 2 °C in the RCP2.6 scenario does not have much effect on rainfall and heavy rain frequency. In 2100, the forecast for rainfall decreases 1.5% and heavy rain frequency decreases 2 days/year, which results in a 0.55% increase for the forecast of AMD quantity, compared to the present. However, the temperature rise of 5.9 °C in the RCP8.5 scenario affects the forecast of rainfall and heavy rain frequency as much as for Mine A. In 2100, the rainfall forecast increases 22% and heavy rain frequency increases 5.5 days/year, which results in a 25% increase for the forecast of AMD quantity, compared to the present.

The forecast of the standard deviation of AMD quantity is shown in Figure 6. A comparison of the coefficient of variation of AMD quantity between the present and the future is shown in

Table 4. The coefficient of variation for AMD quantity at present and in the future.

Mine Name	Present	Future
Mine A	0.083	RCP2.6: 0.092 RCP8.5: 0.090
Mine B	0.56	RCP2.6: 0.49 RCP8.5: 0.54

. Here, the coefficient of variation was calculated for 10 years from 2010 to 2020 for the present and from 2100 to 2110 for the future. The temperature rise due to global warming caused larger fluctuations in the AMD quantity for Mine A. In the case of Mine B, since the AMD reservoir is small and the AMD residence time is short, even if the rainfall fluctuation increases in the future due to global warming, the AMD fluctuation will remain largely as it is now. On the other hand, in the case of Mine A, since the AMD reservoir is larger and the AMD residence time is longer, AMD fluctuation tends to increase gradually in the future, affected by increases in rainfall fluctuation due to global warming. This trend suggests that AMD treatment might be unstable because of global warming in the future, especially for mines with larger AMD quantities.

In general, AMD quality tends to deteriorate as the AMD quantity increases. This is because when the AMD quantity increases, AMD comes in contact with a new pollution source in the mine. Actually, at present, the fluctuation of AMD quantity is larger in Mine B than in Mine A, and the fluctuation in AMD quality tends to be larger in Mine B as well. This suggests that fluctuations in AMD quantity due to global warming will cause large fluctuations in AMD quality.

4. Conclusions

The AMD quantity model was constructed for two closed mines in Japan. The model was constructed with a series tank model, and fitted by using daily data for one year, which were enough to obtain adequate parameters. The results showed that Mine B has a smaller AMD reservoir and a shorter AMD residence time than Mine A, resulting in a large fluctuation of AMD quantity in Mine B. The forecast of AMD quantity was also estimated based on the forecast of rainfall and temperature proposed by the MPI–ESM with IPCC RCP2.6 and RCP8.5 scenarios. The forecast results showed that temperature rise due to global warming will cause an increase in rainfall, resulting in increased AMD quantity. The fluctuation of rainfall will also increase due to global warming, increasing the fluctuation of AMD quantity in Mine A. The effect of global warming in Mine A will be bigger than in Mine B due to its larger reservoir and longer residence time of AMD.

In this study, it is expected that the quantity and fluctuation of AMD might increase due to global warming. This suggests that fluctuations in AMD quality might also increase. Therefore, when selecting future treatment methods, careful consideration should be given to whether or not the AMD fluctuation can be sufficiently dealt with in the future, especially for passive treatment.