Machine Learning to Predict Area Fugitive Emission Fluxes of GHGs from Open-Pit Mines

Abstract

Greenhouse gas (GHG) emissions from open-pit mines pose a global climate challenge, which necessitates appropriate quantification to support effective mitigation measures. This study considers the area-fugitive methane advective flux (as a proxy for emission flux) released from a tailings pond and two open-pit mines, denominated “old” and “new”, within a facility in northern Canada. To estimate the emission fluxes of methane from these sources, this research employed near-surface observations and modeling using the weather research and forecasting (WRF) passive tracer dispersion method. Various machine learning (ML) methods were trained and tested on these data for the operational forecasting of emissions. Predicted emission fluxes and meteorological variables from the WRF model were used as training and input datasets for ML algorithms. A series of 10 ML algorithms were evaluated. The four models that generated the most accurate forecasts were selected. These ML models are the multi-layer perception (MLP) artificial neural network, the gradient boosting (GBR), XGBOOST (XGB), and support vector machines (SVM). Overall, the simulations predicted the emission fluxes with R${}^2$ (-) values higher than 0.8 (-). Considering the bias (Tonnes h${}^{-1}$), the ML predicted the emission fluxes within 6.3%, 3.3%, and 0.3% of WRF predictions for the old mine, new mine, and the pond, respectively.

1. Introduction

In developed and developing countries, energy is a critical component supporting the development of the economy and society [1]. Demand for energy is forecasted to expand globally by 33% from the 2010s to the 2030s, and the demand will lead to an increase in the production of more greenhouse gas (GHG) emissions [2]. In the United States, more than 25% of the methane emitted into the atmosphere is a result of the oil and gas industry. Abating such emissions from the oil and gas sector is a critical component of mitigating climate change [3,4,5]. The warming trend can be slowed down by reducing the amount of methane emissions [6,7]. Methane has a global warming potential (GWP) 28 times greater than carbon dioxide for a duration of 100 years. Approximately one quarter of total global warming is related to the methane emissions from natural and anthropogenic sources. The short-term warming trends are affected by the methane lifetime, which is approximately 10 years [8].

The oil sands in Canada are proven to be the third largest reservoir of oil in the world [9]. The exploitation of such a reserve demands rigorous environmental management, given the energy intensity required to extract the crude oil from the oil sands [10,11]. In 2015, Alberta oil sands accounted for the largest share of GHG emissions from the oil and gas sector, 37% and 10%, respectively, of GHG emissions in Canada [9]. The unconventional oil resources industry, such as oil sands, is more energy-intensive than the conventional oil, with higher GHG emissions associated with resource extraction, which shows the need for more attention to the emissions from these industries [12]. Oil sands operations emit GHGs due to various activities such as drilling, land development, extraction, chemical and physical processing, venting, flaring, and the release from large area sources (also known as fugitive emissions) [2].

For the development of adequate mitigation measures, a reliable quantification of GHG emission fluxes from open-pit mines is needed [13]. Large estimation differences of methane emission fluxes from such facilities are the result of ignoring the meteorological effects that complicate pollutant transport in complex topographies [14]. While orographic complexities confound atmospheric boundary layer (ABL) flows compared to ABL flows over flat and homogeneous terrains [15,16,17,18,19,20,21,22,23], only a handful of investigations are dedicated to ABL flows over open-pit mines [24,25], primarily due to operational, risk, access, economic, and other concerns [25,26]. Open-pit mines represent an unusual topographic feature. The meteorological patterns within and near these mines can be quite different from flat terrain [20,21,22,27,28], and understanding these patterns is important given the environmental impact of mining activities. The release of GHGs from such environments is associated with multiple steps of the resource extraction process, starting from land modification to the exportation of the mined goods [4,14]. Mine activities often create large amounts of fugitive dust, odorous compounds, and GHGs [29], and the flow circulations, shear layers, and plume meandering created by mine pits [30] impact the transport of materials to downwind environments. Approved emissions calculation methods for regulatory reporting fail in quantifying area fugitive emission fluxes from complex terrains. For instance, concerning open-pit mines, Liggio et al. (2019) [31] investigated that area-fugitive emissions of carbon dioxide probed in situ by an aircraft differed by 13–123% of values reported by emission inventory datasets. To address the common challenges of quantifying GHG emission fluxes from open-pit mines, measurement, modelling, and machine learning techniques of emission quantification are briefly reviewed in the subsequent sections.

1.1. Measurement Techniques to Quantify Emission Fluxes

The extrapolation of observations and point measurements is the basis for most of the open-pit mine GHG emission flux predictions by industrial reports. Flux chambers (FCs) are a simple but intrusive means of measuring gas emission fluxes. One of the difficulties with FC measurements is their spatially small measurement footprint, as well as their impracticality for measuring emission fluxes from vertical surfaces (e.g., mine faces). A spatially inhomogeneous and topographically complex source like an open-pit mine would require an extensive and difficult FC sampling survey [32,33]. Micro-meteorological measurement approaches, such as eddy-covariance or flux gradient, overcome some of the weaknesses of FCs (e.g., no interference with the emitting surface, larger measurement footprint, and capability for long-term monitoring), but they are fundamentally challenged by complex wind conditions, and the difficulty of interpreting their measurements for spatially inhomogeneous emission sources [34,35,36]. Another measurement paradigm is based on the concept of mass balance, which attempts to enclose the emitting facility in a hypothetical box in the atmosphere, while the emission flux is measured on the boundary of this box using aircraft or drone sampling [31,37,38,39]. While overcoming the limitations of the previous techniques, this approach has its own shortcomings. For instance, aircraft sampling misses the first few tens of meters near the ground, the measurement of which is crucial for the accurate estimation of the emission fluxes [37]. In addition, mass balance measurements are conducted for a facility only a few times a year, while a more continuous measurement technique is ideal to capture diurnal and seasonal variations of emission fluxes [22]. Because of the large footprint of open-pit mining facilities, precise measurements of their emission fluxes are not possible. Therefore, past research has attempted to offer modeling tools for the estimation of emission fluxes for operational use.

1.2. Modeling Techniques to Quantify Emission Fluxes

Access to high computational power has propelled methods for the calculation of area-fugitive emission fluxes by atmospheric models under two categories of modeling approaches. The first category groups models into diagnostic versus prognostic paradigms. In the diagnostic paradigm, wind field and pollutant dispersion modeling is achieved using empirical formulations, while in the prognostic paradigm, fundamental transport equations of atmospheric flow and dispersion are solved. Diagnostic models are computationally efficient, while they lack accuracy. On the other hand, prognostic models are computationally expensive, while they provide better accuracy. Using either of these paradigms, the second category groups models under either forward or backward dispersion modes. In the forward dispersion approach, the input pollutant sources are resolved both in space and time, while the model calculates the emission flux downwind. However, in the backward dispersion approach, the pollutant concentration is measured downwind, while by relying on a dispersion model, the emission source strength upwind is inferred. The backward dispersion modelling paradigm is also known as inverse dispersion modelling (IDM) [40,41].

The CALifornia PUFF Model (CALPUFF) is an example of a diagnostic model [42]. CALPUFF has been utilized to analyze emission fluxes from industrial and urban sites, as well as the associated air quality and health implications [43,44,45,46,47,48]. When applied to complex topography, CALPUFF has been shown to be less accurate when modeling atmospheric dispersion under the thermally stable condition [49]. In addition, CALPUFF is shown to overpredict the turbulent diffusion, and as a result, underpredict the concentration of pollutants when compared to observations or more accurate models [50]. CALPUFF is also shown to suffer from lack of adequate meteorological forcing for its wind field and dispersion predictions. For instance, when near-surface and upper-air grid points of the model are not sufficiently forced by appropriate meteorological conditions, the wind field predictions lose accuracy [51]. One remedy in operational CALPUFF modeling is to force it with meteorological fields obtained from a prognostic mesoscale weather model [52]. As far as emission flux estimation is concerned, CALPUFF is usually used in an IDM configuration, where a measurement of concentration of a pollutant downwind is carried out, while the emission flux at the source is estimated by relying on CALPUFF’s dispersion calculations. This paradigm in using CALPUFF can be seriously challenged, given its deficiencies when applied to non-conventional topographies, such as a mine pit.

The Weather Research and Forecasting (WRF) model is an example of a prognostic model. Such mesoscale models are employed to provide weather, air quality, and dispersion analyses. These models require physical parameterizations to account for the exchanges of momentum, heat, and atmospheric species (e.g., water) between the earth surface and the atmosphere [53]. Present-day mesoscale models allow for the selection of a rich spectrum of parameterizations, and choosing the effective parameterizations among the alternatives is a very complex task, informed by the desired outcome of a modeling investigation [54]. The accuracy of the actual terrain representation, which is under question, is the primary limitation of WRF at refined grid spacings [55]. In order to make accurate projections in weather, land-use and topographic data are needed, which are counted as significant inputs used for weather forecasting models. In addition, high-resolution terrain data are required for the elucidation of complex meteorological phenomena of different areas [56]. The complex terrains impact weather through changing different atmospheric variables, e.g., wind velocity vector components, radiative properties, and surface sensible and latent heat fluxes [57]. It has been demonstrated that an increase in input topographic and land-use resolutions, combined with refining grid spacing, plays a more significant role than simply using more refined grid spacings in the model toward simulations with more accuracy [58]. The weather patterns that are influenced by terrain features at the regional level play a crucial role in the numerical model by changing the distribution of rainfall, wind field, and atmospheric state variables [59,60]. The land-use information plays an integral part in weather alteration across the scales by influencing the exchanges of momentum, heat, and atmospheric species between the earth surface and the atmosphere [61]. Atmospheric state variables can vary, caused by the influences of changing aerodynamic roughness length scale when higher terrain and land-use resolution information are used [62]. The thermodynamic variables show improved agreement with measurements when detailed terrain information is utilized [20,63].

In regards to emission flux calculations, the WRF model is furnished with a built-in passive tracer dispersion method [22,64,65], which is accompanied by other chemistry modules that can be utilized to compute the temporal and spatial evolution of the mixing ratio of tracers, and consequently, emission fluxes in the forward dispersion paradigm of modeling [66,67,68,69]. The functionality of WRF for dispersion simulations in complex terrain and the land use environment of a mine pit has been investigated comprehensively [22]. The boundary condition for gases at the lowest model layer in mesoscale models is implemented either as fixed-value (fixed-mixing ratio) [22,67] or constant gradient (fixed-flux informed from inventory datasets) [68,70,71,72] (e.g., in accordance to calculations or emission records) representing spatial and temporal distributions. As proposed by Karion et al. (2019) [71], adopting the fixed-flux approach on the basis of emission inventories has been criticized, due to the model’s inability to conserve the emission flux near the surface and model’s exit boundaries, along with the spatio-temporal inaccuracies regarding the inventory datasets [68]. Based on the given observations, the use of the fixed-value (i.e., fixed-mixing ratio) boundary condition may be preferred. Despite the fact that fixed-value approaches do not provide an estimate for the discharged pollutants at the source, they aid the estimation of the emission flux at locations far away from the source, which is affected by weather modulations in the wind field and other state variables. This is named the advective flux, a proxy for the emission flux in the present study.

Motivated by the limitations of the diagnostic and prognostic models, machine learning (ML) may be employed to estimate emissions at a mining facility, given limited observed and simulated data. As opposed to regression techniques, which are linear models, the goal of neural networks is to elucidate non-linear patterns in information by constructing layers of nodes (neurons) in the model. The more popular and accurate machine learning methods, like support vector machines (SVM) [73], which engineer hyperplanes to separate instances, can be used as the main method for emission flux calculation from open-pits. SVM is optimized to create separating hyperplanes in the data among classes to increase the margin among classes. Another ML method is the multi-layer perceptron (MLP) artificial neural network, which has mostly been applied in hydrology [74], although the model is not preferred for present-day machine learning tasks. The gradient boosting (GBR) technique is an approach for classification, regression, and similar problems, offering a prediction among an ensemble of weak prediction models. These weak models are composed of decision trees [75]. XGBOOST (XGB) is a branch of the GBR model, which utilizes the Newton–Raphson method in a function space not similar to the gradient boosting that utilizes the gradient descent method. In the XGB model, a second order Taylor’s expansion is used in the loss function to construct the connection to the Newton–Raphson method. The long short-term memory (LSTM) technique [76] is used in deep learning problems. Not similar to standard feed-forward neural networks, LSTM includes feedback connections. Implementations of procedures such as LSTM allow network training to take place without having long-term parameters “explode” or “vanish” as a result of multiple learning updates [77,78]. ML models based on SVM, deep learning, LSTM, and more have been used in various facets of energy engineering predictions, such as power plant heat transfer rate, power plant emission reduction, fluidized adsorption bed processes, and generator power curves [79,80,81]. However, despite the importance of quantifying emissions from open-pits, the authors of this article have not been informed of any studies attempting to estimate methane emission fluxes from open-pit mining facilities during different diurnal times and seasons using ML models.

1.3. Research Gaps and Objectives

The above review uncovers deficiencies and strengths in all emission flux quantification methods applicable to open-pit mines. Techniques purely based on measurement can be accurate, but suffer from limited spatio-temporal coverage. Diagnostic models are computationally fast but struggle to provide accurate estimates of the emission flux. On the other hand, prognostic models are complicated to set up, and computationally expensive for operational use, but they offer more accurate estimates of the emission flux. Finally, machine learning techniques appear to require extensive training datasets and be specific for each emission source, without the possibility of extrapolating the use of one technique developed for one emission source to another emission source under different conditions.

This study is motivated by the deficiencies and strengths of each method to propose a surrogate method with operational capability to estimate the emission flux of GHGs from an open-pit mine fairly accurately and computationally fast. In the proposed method, meteorological variables and emission flux from a weather research and forecasting (WRF) mesoscale model were used as an input and training dataset for a few machine learning techniques to predict the emission flux. Two sets of models are used. The first set is composed of the GBR, XGBOOST, SVM, and MLP algorithms. The ML model learns in general. The second set includes LSTM, which is only considered as a research model, since in the operational sense, it is difficult to measure the emission flux at any given time (e.g., each hour) on a continuous basis to predict the emission flux in the next time interval (e.g., next hour).

The study is organized as follows. Section 2 describes the methodology. Section 2.1 introduces the open-pit mine site. The details of the WRF model are presented in Section 2.2, and the details of the machine learning models and the statistical analysis are presented in Section 2.3. In Section 3, the results of the WRF and machine learning model simulations for two mines (old and new) and a tailings pond are presented. Section 4 includes the main conclusions and recommendations.

2. Methodology

2.1. Open-Pit Mining Site

Figure 1a shows the site location of the present study, which is an open-pit mine in Canada. The mining site was composed of two open-pits, an old mine and a new mine, and a tailings pond. The old mine, with a depth of roughly 100-m and a span of roughly 2000-m, was mainly the place where the mining excavations were conducted. Another mine, the new mine, which was newly exploited, was smaller than the old mine.

The present study takes a forward dispersion modeling paradigm employing WRF 4.0 [82]. The WRF 4.0 source code is developed in Fortran (https://www2.mmm.ucar.edu/wrf/users/ (accessed 17 January 2020)). A tracer dispersion method was used and the field-measured methane (tracer) mixing ratio was forced at the model’s lowest level. This represented the spatio-temporal variations of methane (tracer) mixing ratio over the two mines and the pond. The advective methane flux (a proxy for emission flux in this study) was simulated using WRF and presented as it experienced changes diurnally and seasonally. At the model’s inner domain boundary, the dynamics of the surface-atmosphere interactions were not accounted for (e.g., biological processes within the pond, the intensity of excavations, or methane build-up inside the pit cavity at night); as an alternative, the meteorologically modulated emission flux downstream of the emission source was computed. Due to the presence of uncertainties regarding this method, it was not possible to accurately quantify the absolute flux of methane emissions from the entire facility. However, the quantification of diurnal and seasonal change in the emission flux could be performed with enough accuracy to motivate future research [22].

The supporting observations and modeling occurred in four campaigns. The first campaign lasted for twenty days in May 2018 (Summer’18 or S18). The second campaign lasted for twenty days in February–March 2019 (Winter’19 or W19). The third campaign lasted for 30 days in July–August 2019 (Summer’19 or S19). The fourth campaign lasted for 30 days in October–November 2019 (Fall’19 or F19). Overall, 100 (days) (2400 (hours)) of observations and modeling were conducted, which are summarized in

Table 1. Supporting observations and modeling in four emissions monitoring campaigns; the date durations for each of the emissions sources are stated.

Campaign	Year	Tailings Pond	Old Mine	New Mine
S18	2018	1 May–11 May	18 May–28 May	-
W19	2019	14 February–24 February	16 March–26 March	-
S19	2019	9 July–19 July	29 July–8 August	12 August–22 August
F19	2019	25 October–4 November	10 November–20 November	8 October–18 October

. The training and testing data for the ML models are provided at an hourly resolution, so the entire dataset includes 2400 records for each input or output variable. There are 960, 960, and 480 data records for each of the tailings pond, old mine, and new mine, respectively.

2.2. Weather Research and Forecasting Model

2.2.1. Model Configurations

Five nested domains were set based on one-way nesting (i.e., no feedback interaction between two adjacent domains), with the outermost domain (D1) including a large area of north-western America, and the innermost domain (D5) mainly encompassing the mining site, including the tailings pond and the two mines. Detailed information about horizontal and vertical grid spacings, time step, grid spacing and time step ratios between adjacent grids, and the time advance scheme are provided by Nambiar et al. (2020) [22]. Most notably, the inner most domain had a horizontal grid spacing of 0.51 (km), and the model included 12 pressure levels from 0.025 to 2 (km) above the surface. The map illustrating the domains is presented in Figure 1a. The WRF settings and physical parameterizations are summarized in

Table 2. WRF physical parameterizations.

Physics	Parameterization
Planetary Boundary Layer (PBL)	Mellor-Yamada-Janić [83]
Microphysics	Thompson Parameterization [84] (D1, D2, and D3 only)
Longwave Radiation	Rapid Radiative Transfer Model  [85]
Shortwave Radiation	Rapid Radiative Transfer Model [85]
Cumulus Parameterization	Tiedtke Parameterization [86] (D1, D2, and D3 only)
Surface Layer (SL)	Monin-Obukhov Eta Similarity Parameterization [83]
Land Surface (LS)	Noah Land Surface Parameterization

.

The standard Global 30 Arc-Second (GTOPO 30s) database was employed to define the latest topography of the site in the three largest domains with a horizontal resolution of 0.9 (km). For the two smallest domains, the Shuttle Radar Topography Mission (SRTM) 1s information was employed with a horizontal resolution of 0.030 (km) [62,87]. For the smallest domain, the extracted topographical resolution (SRTM 1s) was further overwritten with a LiDAR dataset with a horizontal and vertical resolution of 0.001 (km). Figure 1b presents the topography of the site updated with SRTM 1s and LiDAR datasets.

Nahian et al. (2020) [20] investigated that, by means of utilizing the latest modified land-use configurations and topography, the computations of atmospheric state variables at the same facility would agree better with the measurements. The representation of land-use according to the Moderate Resolution Imaging Spectroradiometer (MODIS) 30s dataset (updated 21-class IGBP-MODIS land use dataset) by Friedl et al. (2010) [88] was utilized with categorizations for the pond (21), the two mines (16), mine processing facilities and housing (13), and cleared forest (10). The lake models in WRF as parameterized by Subin et al. (2012) [89] and Gu et al. (2015) [90] could simulate water–atmospheric exchanges of momentum, heat, and humidity. A depth of 0.050 (km) was implemented for the water bodies at the site, and so the lake model included 25 layers of water, soil, and snow combined. Figure 1c presents the settings of the land use implemented for WRF calculations. The initial and boundary conditions in WRF were implemented, resorting to the Global Data Assimilation System (GDAS) dataset (from the National Centers for Environmental Prediction (NCEP)), with grid spacing of 0.25${}^{\circ }$ and a time stamp updated every 6 (hours).

2.2.2. Methane Transport and Flux Calculation

The campaigns provided field observations of methane mixing ratio near the surface at various locations of the mining facility. These data were acquired at heights from 2 to 10 (m) above the surface in four spots uniformly distributed in each of the areas of the tailings pond and the two mines, as depicted in Figure 1d. The instrument used was the Los Gatos Research Ultra-Portable Greenhouse Gas Analyzers (LGRs). Allen et al. (2019) [91] report that this instrument operates based on the principle of off-axis integrated cavity output spectroscopy and two near-infrared tunable diode lasers, which promptly scan a single strong (and isolated) absorption line of the target gas (here methane). The overall release areas, tiled by the rectangular patches, were equalized for the mines and the pond. This area, approximately 30 (km${}^2$), was the same for all simulations and release sources. Methane measurements were acquired every 15 (min). The measurements were averaged every 4 (hours) and forced in the model as a boundary condition at the model’s lowest layer [64]; this was necessitated by the WRF model, which required the recompilation of the program every 4 (hours). In summary, the program was customized to allow tracer release from a desired spatial area.

As proposed by Bhimireddy et al. (2018) [65], the passive tracer dispersion method was employed in WRF without accounting for any chemical reactions, and the tracer moved vertically and horizontally by the computed mean wind field and turbulent transport by means of the surface layer (SL) and planetary boundary layer (PBL) parameterizations. The Monin–Obukhov SL parameterization accounts for the vertical transport of the tracer by relying on similarity theory [83]. The Mellor–Yamada–Janić PBL parameterization accounts for turbulence kinetic energy ($TKE$) employing a 2.5 closure technique [92]. Both horizontal and vertical turbulent motions within PBL are modelled as diffusion using eddy diffusion functions that depend on $TKE$ and other model variables [93]. Twelve (hours) were counted as the spin-up time [22]. To discover the diurnal and seasonal variations in the emission flux, simulations were performed for 4-hour time intervals.

The tracer (methane) emission flux F was computed through a box placed around the emission sources by considering Gauss’ theorem. F can be given by [37]

$$F=F_a+F_t+F_s+F_m+F_c,$$

(1)

where $F_a$ is the advective flux and $F_t$ is the turbulent flux through the box surfaces, $F_s$ is the surface flux because of deposition or biogenic release of methane due to surface activity, $F_m$ is the flux due to density variations in the box, and $F_c$ is the flux because of chemical reactions encompassing methane as it interacts with other species in air. Large-eddy simulations (LES) predict a horizontal turbulent flux of atmospheric species upwind, which accounts for 10–20% of the advective flux [38,70]. However, this flux could not be computed in the present simulations due to the computational cost associated with overall simulations, extending a time span of 100 (days). $F_s$ was disregarded considering the fact that the site land was primarily barren or developed for industrial operations and housing [72]. In accordance with Oertel et al. (2016) [94], the boreal forest is not recognized as a source of methane either (on the contrary, this forest is suggested to absorb methane). Thus, the basic assumption is that the measured methane mixing ratio is only generated from mine pits and ponds, and is not affected by other processes. So, overall, only the advective flux $F_a$ and flux due to density variations in the box $F_m$ are considered to be the major contributors to the total flux, i.e., $F\approx F_a+F_m$. We further disregard $F_m$; after all, any such flux contribution would average to zero over many days. Therefore, $F_a$ was considered as the single major operational contributor to the total flux F. Considering these assumptions and former evidence for the same mining facility, it was suggested that the advective flux solely accounted for more than 95% of the overall emission flux [37]

$$F\approx F_a=\sum _{north}\overline{V}\overline{S}\Delta A+\sum _{east}\overline{U}\overline{S}\Delta A-\sum _{south}\overline{V}\overline{S}\Delta A-\sum _{west}\overline{U}\overline{S}\Delta A+\sum _{top}\overline{W}\overline{S}\Delta A,$$

(2)

where $\overline{U}$, $\overline{V}$, and $\overline{W}$ are the wind velocity vector components (m s${}^{-1}$) simulated at each grid point on the box enclosing the emission source, $\Delta A$ is area element (m${}^2$), and $\overline{S}$ is the concentration of the passive tracer (here methane) ($\mathsf{\mu }$ g m${}^{-3}$) for a specified grid point. It is worth mentioning that when the box height stretches to altitudes well beyond the PBL height, then the last term in the above equation can be neglected (equal to zero). By its formulation WRF simulates the tracer transport in units of (ppm); however to calculate the emission flux in (Tonne h${}^{-1}$), the tracer unit must first be converted to ($\mathsf{\mu }$ g m${}^{-3}$) for methane. Using thermodynamic considerations, the units for the tracer can be converted from (ppm) to ($\mathsf{\mu }$ g m${}^{-3}$). This is possible since all the required thermodynamic variables (e.g., temperature, potential temperature, and pressure) are available in WRF and can be extracted. Nambiar et al. (2020) [22] provide the detailed calculations. Overall, this procedure allows the operational calculation of the advective flux as a proxy for emission flux.

2.2.3. Uncertainty in WRF Estimate of the Methane Emission Flux

The ML model predictions are based on the WRF model estimates, so the success of the ML model predictions depends on the success of the WRF model in predicting the emissions flux of methane. If the accuracy of the WRF predictions of the methane emission flux can be quantified, then ML predictions will have practical relevance in the field.

For quantifying the WRF uncertainty of methane emission flux from the open-pit mine, the relevant uncertainties are bias in wind speed and bias in methane concentrations (or alternatively, mixing ratio). The bias (rather than root mean square error) is relevant because any emission quantification method attempts to report the cumulative emission flux accurately over an entire year, as opposed to a snapshot prediction.

In an earlier study, we quantified the WRF uncertainty in wind speed using airborne measurements in the actual mine field [20]. In that study, the WRF uncertainty and average wind speed from field observations were $\Delta U=0.0819$ and $\overline{U}=2.614$ (m s${}^{-1}$). In another study, we quantified the WRF uncertainty in the mixing ratio of methane using aircraft measurements in the same mine field [22]. In that study, the WRF uncertainty and average mixing ratio from field observations were $\Delta S=0.0543$ and $\overline{S}=0.417$ (ppm). Using these uncertainties and average values, we can estimate the uncertainty in WRF’s advective flux (proxy for emission flux) using

$$\frac{\Delta F}{F}\approx \sqrt{{\left(\frac{\Delta U}{\overline{U}}\right)}^2+{\left(\frac{\Delta S}{\overline{S}}\right)}^2}=\sqrt{{\left(\frac{0.0819}{2.614}\right)}^2+{\left(\frac{0.0543}{0.417}\right)}^2}=0.13,$$

(3)

which suggests that the uncertainty in WRF’s prediction of the methane emission flux is about 13%. Given the lack of alternative methods for prolonged measurements of the emission flux, this uncertainty justifies the use of the WRF dataset as training data for the ML model.

2.3. Machine Learning Model and Statistical Analysis

Under practical considerations, up to four measurement periods (fall, winter, spring, and summer) with no more than 10 days per season are required to quantify emission fluxes every year. This measurement period is sufficient because of the high cost of accurate observation methods, such as deploying instrumented aircraft or another technologically sophisticated system. The GBR, XGBOOST (XGB), SVM, and MLP are candidate ML models that can be used operationally, while predicting the emission flux in the next hour with LSTM, based on measuring the emission flux at the current hour, will remain as a research tool. All the available data for the four campaigns include 10 days of observations and WRF modelling for a given emission source (i.e., old mine, new mine, or pond).

Model tuning was achieved and optimal hyper-parameters were set by minimizing the error statistics with respect to the reference dataset (emission flux predicted by WRF). The settings can be provided as: GBR: number of estimators = 1000, minimum samples per leaf = 4, minimum weight fraction per leaf = 0.01; XGBOOST: gamma = $1\times {10}^{-4}$, max tree depth = 20, number of estimators = 100, learning rate = $5\times {10}^{-2}$, minimum child weight = 2, sub-samples = 0.5; SVM: gamma = $1\times {10}^{-2}$, epsilon = $1\times {10}^{-3}$, $C=1\times {10}^4$, kernel = rbf (radius basis function), cache = 1000 Mbytes; MLP: activation function = ReLu, loss function = MSE, learning rate = $1\times {10}^{-3}$, solver = ADAM, epsilon (for the ADA solver) = $1\times {10}^{-8}$, number of hidden layers = 1 with 100 nodes, alpha (L2 penalty to control over-fitting) = $1\times {10}^{-4}$, batch size = 100. The ML model source code is developed in Python 3.8 (https://www.python.org/ (accessed 10 August 2021)), and it utilized python’s basic data analysis and visualization packages [95,96]. Specific ML algorithms were imported from the Scikit-Learn (sklearn) version 1.0 [97,98] and xgboost version 1.6 [99] libraries.

Figure 2 shows the methodology diagram for the present analysis. Overall, up to 10 features are chosen from WRF as input to the ML models. These are meteorological variables: diurnal time (h), season, wind speed (m s${}^{-1}$) at 10 (m) above surface, temperature (K) at 2 (m) above surface, precipitation (mm h${}^{-1}$), PBL height (m), surface pressure (Pa), relative humidity (%), surface sensible heat flux (W m${}^{-2}$), and surface latent heat flux (W m${}^{-2}$). These variables are averaged over the footprint of each emission source. The WRF-generated emission flux is chosen as the training data for the ML models. The ML models generate the emission flux as output. All the data are shuffled and 80% of data are separated, as the training set from the remaining 20% of data, as the test set. Ten ML algorithms are tested to check for the accuracy. Since the models are fairly simple and run reasonably fast, it is needed to choose the models that best fit the purpose, ignoring the “complexity" of the models. The inclusion of the deep neural network (DNN), in the form of the MLP algorithm, indicated that the benefits of DNNs cannot be fully exploited here, as the training dataset is limited, and therefore, traditional ML algorithms perform better, with a reduced chance for over-fitting. In the present work the analysis is performed by adding one input variable (called feature in the ML method) at a time to check the accuracy of the prediction for a greater number of input variables.

Quantitative comparisons between the model prediction ($M_i$) and a reference dataset ($O_i$) (here methane emission flux from WRF) are performed by determining the bias and root mean square error (RMSE) (Tonnes h${}^{-1}$) defined by

$$Bias=\frac{{\sum }_{i=1}^n(M_i-O_i)}{n},$$

(4)

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{(M_i-O_i)}^2}{n}},$$

(5)

where n is the number of data points accounted for in the error statistic calculation. The coefficient of determination R-squared (R${}^2$ (-)) is the statistical measure of the correlation between a model’s performance and a reference dataset defined by

$$R^2=\frac{{[n{\sum }_{i=1}^n\left(M_i{O}_i\right)-\left({\sum }_{i=1}^n{M}_i\right)\left({\sum }_{i=1}^n{O}_i\right)]}^2}{[n{\sum }_{i=1}^n{\left(O_i\right)}^2-{\left({\sum }_{i=1}^n{O}_i\right)}^2][n{\sum }_{i=1}^n{\left(M_i\right)}^2-{\left({\sum }_{i=1}^n{M}_i\right)}^2]}.$$

(6)

The Pearson’s correlation coefficient ($PCC$ (-)) between $M_i$ and $O_i$ is often used to assess the ability of a model to simulate temporal variation and is given by

$$PCC=\frac{{\sum }_{i=1}^n(M_i-\overline{M})(O_i-\overline{O})}{\sqrt{{\sum }_{i=1}^n{(M_i-\overline{M})}^2}\sqrt{{\sum }_{i=1}^n{(O_i-\overline{O})}^2}},$$

(7)

where overbars signify averages and n is the sample size.

3. Results and Discussion

Time series of wind speed at 10 (m) above surface, potential temperature at 2 (m) above surface, surface heat flux, and emission flux of methane during 10 days of F19 campaign for the old/new mines and the pond are presented in Figure 3. All data shown are extracted from WRF simulations. In the time period from 48 to 96 (hours), higher wind speed and surface heat flux are predicted over the pond compared to the mines (Figure 3a,c). The higher wind speed and surface heat flux caused a higher emission flux from the pond (Figure 3d). Overall, the pond shows greater variability in the wind speed, surface heat flux, and emission flux noticed by the peaks from Figure 3a,c,d. The reason for such variations could be related to meteorological effects over water bodies. For the mines, the emission flux drops during the night as the wind speed decreases over night in the mine cavity. The change in the diurnal variation of temperature for the mines are higher than the pond as the heat capacity of water bodies are higher than soil, which could damp the amplitude of temperature oscillations in the atmosphere above the pond. Figure 3d shows that as the new mine is more active, the emission flux from it is higher. The new mine emission flux exhibits a notable diurnal pattern, possibly due to the regular schedule of the mining activity. Meanwhile, such diurnal pattern is not as notable for the pond, which emits methane given the natural variation in the meteorological conditions. Likewise, such diurnal variability in the emission flux is weaker for the old mine.

Figure 4, Figure 5 and Figure 6 show the $PCC$ matrix computed for measured tracer (methane) mixing ratio, normalized emission flux (computed by WRF), and other atmospheric state variables (computed by WRF) for the old mine, new mine, and pond, respectively. For this analysis all the data is used from the WRF simulations, i.e., the entire 2400 records, with 960, 960, and 480 records related to the tailings pond, old mine, and new mine, respectively. As can be seen for the old mine in Figure 4 for S18, W19, S19, and F19 cases, the emission flux was greatly correlated with wind speed at 10 [m] above surface ($S10$) in all seasons but highly correlated with Surface Heat Flux ($SHF$) (W m${}^{-2}$) only during the warm conditions (S18 and S19 campaigns). For example, in the S18 case, the $PCC$ s for $S10$ and $SHF$ are 0.82 (-) and 0.74 (-), respectively, for the old mine (Figure 4a). The same correlation for the new mine and pond can be seen for $S10$ and emission flux in warm conditions (S18 and S19 campaigns) (Figure 5a and Figure 6a,c). The emission flux of the pond in S18 and S19 campaigns was highly correlated with $S10$ ($PCC=0.81,0.62$ (-)) but anti-correlated with $SHF$ ($PCC=-0.43,-0.37$(-)) (Figure 6a,c). The anti-correlation was foreseen having the knowledge that the water temperature would be different than that of the nearby earth surfaces. This is itself due to larger thermal capacity of water compared to soil [20].

In S19 and F19 campaigns, the emission fluxes from the old and new mines were highly correlated with $S10$ for the old mine ($PCC=0.74,0.84$ (-)) and for the new mine ($PCC=0.89,0.85$ (-)). In S19 and F19 campaigns, the emission fluxes were highly correlated with $SHF$ for the old mine ($PCC=0.35,0$ (-)) and for the new mine ($PCC=0.4,0.68$ (-)). However, in S19 and F19 campaigns, the emission flux from the pond was highly correlated with $S10$ ($PCC=0.62,0.77$ (-)), but either anti-correlated or correlated with $SHF$ ($PCC=-0.37,0.73$ (-)), possibly due to transitioning from the warm to the cold conditions. The greatest correlation with the emission flux is with the wind speed. The higher the wind speed, the greater the emission flux. For the case of the mines, surface sensible heat flux is positively correlated with the emission flux. On the other hand, for the case of the pond, surface sensible heat flux is negatively correlated with the emission flux in the warm conditions ($PCC=-0.43,-0.37$ (-) for S18 and S19, respectively) but positively correlated with the emission flux in the cold conditions ($PCC=0.12,0.73$ (-) for W19 and F19, respectively). This difference can be explained given the meteorological processes. In the warm conditions (S18 and S19), the high emission flux is accompanied with high wind speed during the daytime, where the pond surface is cooling the atmosphere due to being at a lower temperature (negative heat flux). While in the cold conditions (W19 and F19), the high emission flux is accompanied with high wind speed during synoptic events, where the pond surface is warming the atmosphere due to being at a higher temperature (positive heat flux).

Figure 7, Figure 8 and Figure 9 show the scatter plots of emission fluxes predicted by the MLP, GBR, XGB, and SVM models versus the WRF predictions. This data is associated with the 20% testing data, for which the ML models were not trained. On the same figures some statistics for the comparison are calculated and reported, which include the Bias (Tonnes h${}^{-1}$), RMSE (Tonnes h${}^{-1}$), and R${}^2$ (-) for the old mine, new mine, and the pond.

Table 3. Bias (Tonnes h${}^{-1}$), RMSE (Tonnes h${}^{-1}$), and R${}^2$ (-) values for the four machine learning methods (MLP, GBR, XGB, SVM) of old and new mines and the pond for S18, W19, S19, and F19.

		Old Mine (S18, W19, S19, F19)		New Mine (S19, F19)		Pond (S18, W19, S19, F19)
Bias (Tonnes h${}^{-1}$)	MLP	−0.47	Average = −0.62	−0.51	Average = −0.37	0.25	Average = 0.03
	GBR	−0.46		−0.17		0.11
	XGB	−0.63		−0.13		−0.06
	SVM	−0.91		−0.65		−0.20
RMSE (Tonnes h${}^{-1}$)	MLP	3.93	Average = 3.96	2.87	Average = 3.20	3.75	Average = 3.73
	GBR	3.85		3.57		3.96
	XGB	4.04		3.36		3.81
	SVM	4.02		2.98		3.38
R${}^2$ (-)	MLP	0.89	Average = 0.89	0.94	Average = 0.92	0.85	Average = 0.85
	GBR	0.89		0.90		0.83
	XGB	0.88		0.91		0.84
	SVM	0.89		0.94		0.88

shows the Bias (Tonnes h${}^{-1}$), RMSE (Tonnes h${}^{-1}$), and R${}^2$ (-) values for the four machine learning methods (MLP, GBR, XGB, SVM) applied to the old and new mines and the pond during the individual campaigns. The error statistics are associated with the 20% testing data, for which the ML models were not trained. The averaged Bias (Tonnes h${}^{-1}$) shows lower values for the pond and new mine, while the old mine has the highest Bias (Tonnes h${}^{-1}$) and RMSE (Tonnes h${}^{-1}$). This may suggest a more successful prediction of emission flux by the ML methods for the pond and the new mine compared to the old mine. The average R${}^2$ (-) of prediction by the deployed ML methods shows better forecast for the new mine although the two other predictions have R${}^2$ (-) higher than 0.8 (-). This may be due to more regular and scheduled mining activity in the new mine, compared to the old mine and the pond. Comparing the Bias (Tonnes h${}^{-1}$), RMSE (Tonnes h${}^{-1}$), and R${}^2$ (-) shows that the GBR and XGB models performed better predictions for all the different emission sources and seasons. The average WRF emission flux for the old mine, new mine, and the pond are 9.90, 11.15, and 9.19 (Tonnes h${}^{-1}$), respectively. Considering the Bias (Tonnes h${}^{-1}$), the ML-predicted emission fluxes are on average within 6.3%, 3.3%, and 0.3% of WRF prediction for the old mine, new mine, and the pond, respectively. The same numbers considering RMSE (Tonnes h${}^{-1}$) would be around 30% of WRF predictions.

4. Conclusions and Recommendations

Open-pit mines represent an unusual topographic feature. The meteorological patterns within and near these mines can be quite different from those of flat terrain, and understanding these patterns is important given the environmental impact of mining activities. Mine activities often create large amounts of fugitive dust, odorous compounds, and greenhouse gases (GHGs), the quantification of which is confounded by flow circulations, shear layers, and plume meandering created by mine pits. Open-pit mine emissions measurement techniques have limitation of spatio-temporal coverage. Although the diagnostic models are computationally fast, the models have difficulty to estimates the emission flux accurately. On the other hand, more accurate prognostic models are complicated to set up and computationally expensive for operational use. The machine learning techniques are fast, but require extensive training datasets to be accurate.

This work proposed a surrogate method with operational capability to estimate the emission flux of GHGs from a mine pit in Canada fairly accurately and computationally fast, identifying diurnal and seasonal variations of the flux. The study used Weather Research and Forecasting (WRF) and machine learning (ML) methods. The in situ field measurements of methane mixing ratio were conducted in various spots of two mine pits, denominated old and new, as well as a tailings pond in May 2018 (S18), February–March 2019 (W19), July–August 2019 (S19), and October–November 2019 (F19) to provide tracer boundary conditions for WRF. Methane transport in WRF was simulated using a passive tracer dispersion method. The meteorological variables and emission flux from the WRF model were used as training and test datasets for various ML models. Four ML techniques (multi-layer perceptron (MLP) artificial neural network, gradient boosting (GBR), XGBOOST (XGB), and support vector machines (SVM)) were used as the main modeling methods. The long short-term memory (LSTM) technique was considered as a research model since, in the operational sense, it is difficult to measure the emission flux at any given time interval on a continuous basis to predict the emissions in the next time interval.

As the new mine was more active than the old mine, the emission flux out of the new mine was higher than the old mine and showed significant diurnal variations. The emission flux was greatly correlated with wind speed at 10 (m) above the surface ($S10$), but it exhibited lower or anti-correlations with other meteorological variables. The ML methods were more successful at the prediction of emission flux for the pond and the new mine compared to the old mine. The average bias and RMSE (Tonnes h${}^{-1}$) were lower for the pond and new mine, while the old mine had the highest bias (Tonnes h${}^{-1}$) and RMSE (Tonnes h${}^{-1}$). Because of the more scheduled mining activity in the new mine, compared to the old mine and the pond, the average R${}^2$ (-) of emission flux prediction by the deployed ML methods was higher than the two other predictions; however, the old mine and pond still exhibit an R${}^2$ (-) higher than 0.8 (-). Considering the bias (Tonnes h${}^{-1}$), RMSE (Tonnes h${}^{-1}$), and R${}^2$ (-), the GBR and XGB algorithms performed better predictions for all the different emission sources and seasons. Combining all emission sources (mines and pond), seasons, and the four ML methods, the new WRF-ML modeling paradigm can predict the emission flux with an average Bias, RMSE, R${}^2$ of −0.32 (Tonnes h${}^{-1}$), 3.63 (Tonnes h${}^{-1}$), and 0.88 (-), which reveals a potential interest in use of ML methods for predicting emission fluxes from the complex terrains of open-pit mines.

Some future work is possible. More accurate observations of methane mixing ratio with higher spatial resolution near the surface and over prolonged periods can potentially improve the WRF-ML prediction results. One of the main tasks which remains is whether WRF can conserve the mass flux of methane near the surface and at the exit boundary of its inner domain. If such conservation can be established, the proposed method can offer an accurate proxy for the emission flux from open-pit mines.