Diurnal and Seasonal Variation of Area-Fugitive Methane Advective Flux from an Open-Pit Mining Facility in Northern Canada Using WRF

Greenhouse Gas (GHG) emissions pose a global climate challenge and the mining sector is a large contributor. Diurnal and seasonal variations of area-fugitive methane advective flux, released from an open-pit mine and a tailings pond, from a facility in northern Canada, were simulated in spring 2018 and winter 2019, using the Weather Research and Forecasting (WRF) model. The methane mixing ratio boundary conditions for the WRF model were obtained from the in-situ field measurements, using Los Gatos Research Ultra-Portable Greenhouse Gas Analyzers (LGRs), placed in various locations surrounding the mine pit and a tailings pond. The simulated advective flux was influenced by local and synoptic weather conditions in spring and winter, respectively. Overall, the average total advective flux in the spring was greater than that in the winter by 36% and 75%, for the mine and pond, respectively. Diurnal variations of flux were notable in the spring, characterized by low flux during thermally stable (nighttime) and high flux during thermally unstable (daytime) conditions. The model predictions of the methane mixing ratio were in reasonable agreement with limited aircraft observations (R2 = 0.68). The findings shed new light in understanding the area-fugitive advective flux from complex terrains and call for more rigorous observations in support of the findings.


Introduction
Atmospheric methane (CH 4 ) is a Greenhouse Gas (GHG) with about 28-35 times more global warming potential (over a time period of 100 years) than carbon dioxide (CO 2 ) [1 -4]. Methane contributes to the generation of further GHGs such as tropospheric ozone and stratospheric water vapour [4][5][6]. Methane emissions from natural and anthropogenic sources account for approximately one-quarter of today's total warming [7]. In comparison to the longer atmospheric lifetime of CO 2 , which implicates long-term warming trends, the lifetime of methane is approximately 10 years, which influences short-term warming trends. Mitigating the amount of methane emissions can possibly cause an immediate slowdown of the warming trends [8,9]. Therefore, it is very important to better understand, measure, and quantify the magnitude and variations of methane emissions from different sources.
The fossil fuel industry is one of the major contributors of anthropogenic methane emissions, globally (∼100−180 Tg yr −1 ) [10,11]. For instance, the United States Environmental Protection Agency mining production is very inaccurate for area-fugitive emissions, since such emission fluxes do not scale with production [29]. Satellite observations have numerous drawbacks: (1) they cannot observe GHGs at night or in high latitudes with very oblique solar zenith angles; (2) the spatial and temporal resolutions are low such that they cannot map a facility continuously with high spatial detail [35,36]. The rudimentary IDM methods suffer in complex topographies and land-use environments where horizontal homogeneity in meteorological fields cannot be assumed [37,38].
Recent research advancements and availability of computational power have enabled alternative estimation of area-fugitive emission fluxes of GHGs using forward (or backward for advanced IDM) dispersion modelling employing mesoscale meteorological models that ingest near-surface-level observations of GHGs (or downwind observations of GHGs for advanced IDM) and predict the advective fluxes of GHGs at the model exit boundaries (or near-surface-level emission fluxes of GHGs at the source for advanced IDM). The advective flux can be different from the emission flux due to various mechanisms in storing, reacting, and depositing the GHGs in the atmospheric domain, to be discussed in the methodology section. However, if the advective flux is monitored or predicted reliably over long time durations of multiple days, it can be used as an alternative measure of the amount of GHGs that enter the atmosphere. This motivates a paradigm shift in GHG emission monitoring and regulations. This approach accounts for complex meteorological transport phenomena, particularly over complex topographies, and enable operational advective flux estimation as it varies diurnally and seasonally. For instance the Weather Research and Forecasting (WRF) model has a built-in passive tracer dispersion model [39,40] as well as other chemistry plugin options such as WRF-Chem and WRF-GHG that can be used to simulate the mixing ratio of GHGs and subsequently advective fluxes forward in time [41][42][43][44]. In such models, the boundary condition for near-surface-level gases is either implemented as a fixed-mixing ratio [42] or fixed-flux [43,[45][46][47][48][49] (e.g., based on measurements or emission inventories) accounting for spatio-temporal distributions. The fixed-flux approach is not used in some cases where spatio-temporal inaccuracies exist in the inventory datasets for surface emission fluxes [43]. Also, over complex terrains with topography and land use variations, the near-surface-level emission fluxes are very heterogeneous and difficult to measure due to flow variations near the surface, while mixing ratios of GHGs are less spatially variable [50,51]. This may justify use of the fixed-mixing ratio boundary condition based on observations. Although such methods do not estimate the emission flux at the source (e.g., an open-pit mine or a tailings pond), they enable estimation of the advective flux downwind of the source as influenced by meteorological effects.

Objectives
While most methane emission estimates in the literature are only concerned with midday conditions, this study aims to quantify diurnal and seasonal variation of methane advective flux from an open-pit mining facility. This objective is inspired by and based on the premise that varying earth surface (temperature, moisture, land cover, land use, etc.) [52] and meteorological (thermal stability, wind speed, etc.) [17] conditions can influence methane advective fluxes that are currently not reported by previous research. As shown in Figure 1a, the open-pit mining facility is located in northern Canada near the Wood Buffalo National park. The facility includes a tailings pond, which is an area of refused mining waste where the waterborne refuse material is pumped. The mining excavations are primarily conducted over the mine area, which is approximately 100 m deep, with a width-to-depth aspect ratio of greater than 20. The supporting observations and modelling occurred for twenty days in May 2018 (late spring) and twenty days in February-March 2019 (winter and early spring). This study utilized a forward dispersion modelling approach using WRF 4.0 with a passive tracer dispersion option supplied by near-surface-level boundary conditions with field-measured methane mixing ratios accounting for spatio-temporal variations over the mine and the pond. The methane advective flux was then calculated at the model's inner domain boundary and reported as it varied diurnally and seasonally and as a function of meteorological conditions. This methodology to determine the advective flux did not reflect the dynamics of the surface-atmosphere boundary (e.g., biogenic activity inside a tailings pond, rate of exposure of the freshly-excavated mining surface to the air, or storage of methane in a thermally stable surface layer inside the mine at nighttime); instead the focus was on the meteorologically-modulated advective flux downstream of an emission source. Given the uncertainties of this method, accurate quantification of real-time emission flux of methane was not possible; however, quantification of the emission flux over long time durations is possible by monitoring or predicting the advective flux over multiple days (assuming no reactions, surface deposition, or biogenic generation of GHGs). This is true since the net stored GHG in the model inner domain will be zero over extended time periods. In addition, the relative change in the advective flux as it varied diurnally and seasonally could be provided with sufficient accuracy to inspire further investigation.

WRF Model Set-Up
This study used the Weather Research and Forecasting (WRF) 4.0 model with the Advanced Research WRF (ARW) dynamical core [53] and with a passive tracer dispersion option. The model was configured with five nested domains using one-way nesting, with the outermost domain (d01) covering most of the western and northern Canada and parts of the northern United States, while the innermost domain (d05) mainly covered the mining facility, comprised of the tailings pond and the mine. As in Nahian et al. [17], a ratio of 1/3 was used for both time step and horizontal grid spacing in successive nested grids. The time step for d01 and d05 were ∼60 and ∼1 s, respectively. The third order Runge-Kutta split-explicit scheme was used for time integration. In the horizontal directions, the grid spacings for d01 and d05 were 41,000 and 506.17 m, respectively. In the vertical direction the model had 90 eta levels from the surface up to an altitude of 20 km, with the lowest level located at 25 m above ground and the first twelve levels below 2 km above ground (approximate height of the daytime Planetary Boundary Layer (PBL)). The map illustrating the position of the domains is shown in Figure 1a. The WRF configuration and physics scheme details are summarized in Table 1. To construct the most recent topography of the mining facility, in d01, d02, and d03 the standard Global 30 Arc-Second (GTOPO 30s) dataset in WRF was used, which provided a horizontal resolution of 900 m. For d04 and d05, the Shuttle Radar Topography Mission (SRTM) 1s dataset with a horizontal resolution of 30 m was used [59,60]. For d05, the background resolution (SRTM 1s) was overwritten with an up-to-date Light Detection And Ranging (LiDAR) dataset with a horizontal and vertical resolution of 1 m to reflect the presence of the tailings pond and the mine. Figure 1b shows the topographic map updated with SRTM 1s and LiDAR datasets.
For a realistic simulation of the atmospheric dynamics and dispersion, the finest details of the land-use configuration should be accounted for in the WRF model [17]. Nahian et al. [17] found that WRF simulations of meteorological variables, at the same mining facility, were closer to the observations when the latest mine topography and modified land-use configurations were used. As a result, in that study the model biases reduced from 1.10 to 0.08 m s −1 , from 1.04 to 0.50 m s −1 , from 0.98 to 0.32 K, and from 45.7 to 17.3 W m −2 , for near-surface-level wind speed, boundary-layer wind speed, near-surface-level potential temperature, and turbulent sensible heat flux, respectively. Land-use configuration of the Moderate Resolution Imaging Spectroradiometer (MODIS) 30s data product (Modified 21-category IGBP-MODIS land use) by Friedl et al. [61] was used with classifications for the tailings pond (lake = 21), mine (barren = 16), plant and lodge (urban and built-up = 13), and other (grassland = 10). The Subin et al. [62] and Gu et al. [63] lake model in WRF was enabled to properly simulate the thermodynamics and aerodynamics effects of water bodies on the atmosphere. For the tailings pond and other water bodies, a depth of 50 m was considered, and the lake model accounted for 25 model layers including up to 5 layers for snow, 10 layers for water, and 10 layers for soil. The land-use configuration used for WRF simulations is shown in Figure 1c.
The meteorological initial and lateral boundary conditions were obtained from the National Centers for Environmental Prediction (NCEP) Global Data Assimilation System (GDAS) dataset with a spatial resolution of 0.25 • and a temporal resolution of 6 h [64].

Methane Transport in WRF
The near-surface-level methane mixing ratios were obtained from field observation campaigns at the mining facility. Mixing ratios of methane were measured at an elevation of ∼2-10 m above ground in various locations surrounding the mine and the tailings pond using Los Gatos Research Ultra-Portable Greenhouse Gas Analyzers (LGRs). The LGRs were placed at four locations surrounding the mine and four locations surrounding the tailings pond. As shown in Figure 1d, in WRF, four rectangular patches were defined, each centred around the measurement point, corresponding to both the mine and the tailings pond. The total methane source areas covered by the patches were made equal for both the mine and the tailings pond and were kept the same for all simulations. The combined surface area of four source locations was 29.72 km 2 . The frequency of measurement was 15 minutes. The mixing ratio of the measurement was averaged every four hours and updated in the WRF model as the boundary condition; the four-hour averaging was required by the limitations of running WRF, which mandated a recompilation of the source code every four hours with updated near-surface-level boundary conditions for methane. The mixing ratios were applied at the model's lowest level at 25 m above ground [39]. The WRF source code was modified to enable the methane release from a specified location of interest. This allowed analyzing relative advective fluxes between the mine and the tailings pond and reporting relative advective fluxes as they varied diurnally and seasonally.
In order to explore the methane mixing ratio variation, and the advective flux, in different meteorological conditions, WRF model simulations were conducted during May 2018 (late spring) and February-March 2019 (winter and early spring) corresponding to methane observations for both the mine and the tailings pond at the mining facility. Mine simulations were conducted for ten days on 18-27 May 2018 (late spring) and on 16-25 March 2019 (late winter and early spring), respectively. Tailings pond simulations were conducted for ten days on 01-10 May 2018 (late spring) and on 14-23 February 2019 (winter), respectively. The operations of the mining facility remained the same during the study period, except for a short duration of shift change for the operation staff ocurring daily at 0600 and 1800 local standard time (LST = UTC − 7).
In WRF the passive tracer mixing options were enabled. The passive tracer had no chemical properties and was treated by the model as a scalar variable. The passive tracers were released at every time step and they were transported horizontally and vertically by the model-simulated wind field and turbulent parameterization using the Surface Layer (SL) and PBL schemes [40]. The main transport mechanisms influencing the spatio-temporal distribution of a passive tracer mixing ratio in the PBL were horizontal/vertical advection, horizontal/vertical turbulent transport, and SL fluxes [65]. The Monin-Obukhov SL scheme provided vertical flux of the tracer using similarity relations [54,66]. The MYJ PBL scheme parameterized the Turbulence Kinetic Energy (TKE) using a 2.5 closure scheme [67]. For horizontal turbulent transport, a second-order diffusion equation was solved employing an eddy diffusion coefficient as a function of TKE and other model parameters. For vertical turbulent transport, the tracer vertical flux was computed by employing an eddy diffusion coefficient, which was also formulated using TKE and other model parameters [68].
WRF was spun-up for 12 h before releasing the passive tracer from the specified surface grid points. The tracer output after the first 3 h of release was used as the initial and lateral boundary conditions, except for the near-surface-level boundary condition, where the observed mixing ratio was used. In order to find the variation in advective flux under different diurnal and seasonal times, simulations were conducted for four-hour time intervals of 0000-0400, 0400-0800, 0800-1200, 1200-1600, 1600-2000 and 2000-2400 LST, when measured near-surface-level mixing ratios were used as constant boundary conditions, while meteorological fields were updated according to the GDAS dataset at a much higher time resolution of ∼1-60 s. The advective methane flux at the model inner domain exit boundary was calculated and is reported with a one-hour time resolution.

Flux Calculation
Using the divergence theorem, the methane emission flux F from a control volume box surrounding an open-pit mining facility can be given by [26] where F adv and F turb are the advective and turbulent fluxes through the box walls and top lid, F sur f is the surface flux due to deposition or biogenic generation of methane, F mass is the flux associated with increasing mass within the control volume due to density changes, and F chem is the flux due to chemical reactions involving methane. Accounting for turbulent flux in WRF requires Large-Eddy Simulation (LES) of methane transport with a high sampling frequency of wind and tracer solutions every 5 min or less [27,45]. However, in this case a PBL scheme was used at coarser resolutions to be able to simulate an overall long time-span of 40 days for the estimation of the flux with the sampling frequency of the solution every 1 h. Using LES and wind tunnel experiments, it was estimated that the horizontal turbulent flux of the tracer was upwind and 10-20% of the magnitude of the advective flux [27,45]. Over short ranges in the inner WRF domain, the flux due to deposition, biogenic generation, and chemical reactions of methane may be ignored given the fact the site land was mostly bare (except for the emission source areas) or developed for residence and industrial operations [49]. Further, the surrounding boreal forest is known not to be a source for methane flux (in fact the boreal forest is a sink for methane flux) [52], so the measured methane mixing ratio could be assumed to be generated from area-fugitive sources and not influenced by background methane levels. Therefore, the divergence theorem reduces to F = F adv + F mass . If this flux is monitored or predicted over a long time duration, F mass should average to zero since the built-up and release of the stored methane in the model inner domain cancel over many days, so a prolonged monitoring or prediction of F adv provides an alternative measure of the amount of GHGs entering the atmosphere. In fact, these assumptions could be justified using ground-based and aircraft measurements of the methane emission flux at the same facility, where it was estimated that the turbulent, surface, mass, and chemical fluxes of methane were at most up to ∼3-5% of the advective flux, combined [26]. With these assumptions and previous evidence, the advective flux alone contributed to ∼95-97% of the total emission flux. So for operational purposes, we could calculate the advective flux, to provide an estimate of the total flux, as where U, V, and W were the average components of the wind velocity vector in (m s −1 ) evaluated at each grid cell on a bounding box containing the emission source, and S was the average concentration of the methane tracer in (µg m −3 ) for a given grid cell on the bounding box. Note that if a tall box is considered that stretches to altitudes many multiples of the PBL height, then integration of the advective flux on the box top is not necessary, since it would be zero. The tracer units of (ppm) could be converted to (µg m −3 ) having the molecular weight of methane M CH 4 = 16.04 g mole −1 , total pressure for a grid cell as a function of height P(z) in (Pa), average temperature for a grid cell as a function of height T(z) in (K), and the gas constant R = 8.314 J mole −1 K −1 by where MR was the average tracer mixing ratio for the grid cell in (ppm). Having the average potential temperature for the grid cell as a function of height Θ(z) in (K) as a model output, we could calculate the temperature at any given height for the grid cell using where P 0 was the reference sea-level pressure in (Pa). This allowed calculation of the advective flux of methane in units of mass per time, for example, in (tonne hr −1 ). Even though there were inherent errors for neglecting some flux terms, this simplified way enabled the method to quantify diurnal and seasonal variations of flux for operational purposes.

Statistical Analysis
The advective fluxes computed from the methodology of the previous section associated with different emission sources (e.g., mine or pond), diurnal time, or seasons, can be compared to one another using a statistical test. The two-sample t-test was performed to find if the difference between two calculated average advective fluxes (e.g., F CH 4 ,1 and F CH 4 ,2 ) associated with two simulations was statistically significant with a confidence level of 95%. The t-test rejected (R for p < 0.05) or accepted (A for p > 0.05) the Null Hypothesis (NH) that there was no difference in the mean advective fluxes between two cases in question, i.e., µ 1 = µ 2 . The t statistic is given by where overbars represented averages, n was the sample size, and S 1 and S 2 were the sample standard deviations for the fluxes. A degree of freedom was also considered for the t-test according to Nahian et al. [17]. The Normalized Average Difference (N AD) was used to represent the shift between two sets of model simulations, given by [69] The Pearson's Correlation Coefficient (PCC) between variables F CH 4 ,1 and F CH 4 ,2 is given by where overbars represented averages and n was the sample size. A PCC value close to 1 indicates that both variables are highly correlated, while a PCC value close to zero shows that there is no correlation. A negative PCC value close to −1 indicates anti-correlation.

Diurnal and Seasonal Variation of Advective Flux
The near-surface-level mixing ratios of methane were measured at four locations surrounding the mine and four locations surrounding the tailings pond using Los Gatos Research Ultra-Portable Greenhouse Gas Analyzers (LGRs). The LGR measurements were made during field campaigns in late spring 2018 and winter and early spring 2019. These field measurements were, for the mine from 18 to 27 May 2018 (M18), the tailings pond from 1 to 10 May 2018 (P18), the mine from 16 to 25 March 2019 (M19), and the tailings pond from 14 to 23 February 2019 (P19). Figure 2 shows the time series of four-hour-averaged near-surface-level mixing ratio of methane at four observation locations during the M18, P18, M19 and P19 field campaigns; the four-hour averaging was required by the limitations of running WRF, which mandated a recompilation of the source code every four hours with updated near-surface-level boundary conditions for methane. Figure 3 shows the time series of hourly advective methane flux for all the simulations. These figures can be studied in conjunction with Figure 4 showing the PCC matrix calculated for the near-surface-level mixing ratio, normalized advective flux, and near-surface-level meteorological variables. The mine showed more diurnal variations in the mixing ratios as can be seen from the many observed elevated peaks for both the M18 and M19 measurements. For the M18 observations in late spring, the elevated mixing ratios were mostly in the nighttime or early morning hours. One reason for diurnal variation could be meteorological and thermal stability conditions. This is also in agreement with the negative correlation between the average Mixing Ratio (MR) of four locations with WRF-simulated wind speed at 10 m (S10)(PCC = −0.32) and upward surface Sensible Heat Flux (SHF) (PCC = −0.44). M19 observations of MR also showed negative correlation with S10 (PCC = −0.33). In comparison to the mine observations, the tailings pond mixing ratio showed less diurnal variation. P18 observations showed more variability in comparison to P19 observations.  The hourly-simulated methane advective flux F CH 4 from the inner domain (05) boundaries, surrounding the open-pit mining facility, was calculated. There were four groups of WRF simulations, each 10 days long, following the same time periods of M18, P18, M19 and P19 LGR observations of methane mixing ratio. Figure 3a shows the times series of hourly-normalized advective flux of methane (F CH 4 ) for 10 days associated with M18, P18, M19 and P19 simulations, respectively. The advective flux was normalized by the average hourly advective flux (F CH 4 ) calculated for the duration of 10 days in each case. Here the value of one on the vertical axis represents the average advective flux. This figure can be studied in conjunction with Figure 4 showing the PCC matrix calculated for the near-surface-level mixing ratio, normalized advective flux, and near-surface-level meteorological variables. It must be noted that reporting absolute advective fluxes with high accuracy using this method is only possible if the mixing ratio of near-surface-level methane is rigorously measured in numerous points. For this study, data from only four measurement locations were available at a time, therefore, the advective flux was reported in the normalized form. Nevertheless, this analysis was informative to provide the behaviour of diurnal and seasonal variations of advective flux with some confidence.  As can be seen in Figure 3a for the mine simulations M18 in the late spring, during the thermally stable nighttime and early morning periods, the advective flux could be as low as 20% of the average, while during the thermally unstable mid afternoon hours it could be as high as 400% of the average. These percentage values could change from day to day as the meteorological conditions over the mining facility changed. In M18 simulations, the advective flux was highly correlated with S10 (PCC = 0.82) and SHF (PCC = 0.74) (Figure 4a). This can also be seen from the time series of S10 and SHF ( Figure A1a in Appendix A) and colour plots ofF CH 4 at given S10 and SHF values ( Figure A4a in Appendix A). The lower advective flux at nighttime corresponded to low S10 and SHF, while the higher advective flux during daytime corresponded to high S10 and SHF (Figures A1a and A4a in Appendix A). Diurnal changes in the advective flux from the tailings pond simulations P18 could be seen in Figure 3b. In comparison to M18, there were attenuated advective flux peaks during the thermally unstable mid afternoon hours, while the same trend in advective flux was witnessed at the thermally stable nighttime and early morning hours. In P18 simulations, the advective flux was highly correlated with S10 (PCC = 0.81) but anti-correlated with SHF (PCC = −0.43) (Figure 4b and Figures A1b and A4b in Appendix A). The anti-correlation was expected given the difference in surface temperatures between the water bodies and the surrounding earth surfaces [17]. Per Figure 4a,b the correlation coefficients among other meteorological variables were indicative of diurnal variations associated with local weather patterns.  For mine simulations M19 in late winter and early spring, as shown in Figure 3c, for the first few days a higher advective flux was predicted during the nighttime and early morning hours, while a lower advective flux was predicted during the mid afternoon hours. This seasonal change in the advective flux could be caused by the meteorological conditions in winter that were associated with synoptic events. In M19 simulations, the advective flux was highly correlated with S10 (PCC = 0.62) but not SHF (Figure 4c and Figures A1c and A4c in Appendix A). The Sea Level Pressure (SLP) for M19 simulations ( Figure A3c in Appendix A) in winter was higher than SLP for M18 simulations in the late spring ( Figure A3a in Appendix A). Nevertheless, the last few days of M19 simulations showed a diurnal variation in the advective flux similar to the M18 simulations (Figure 3a) possibly as local weather patterns emerged in the early spring. The methane advective flux from the tailings pond simulations P19 in winter is shown in Figure 3d. The advective flux exhibited a totally different pattern from the P18 simulations in the late spring (Figure 3b). A clear diurnal pattern in variation of advective flux could not be seen in P19 simulations suggesting that synoptic weather conditions could be prominent over local weather conditions. In this case the advective flux was only highly correlated with S10 (PCC = 0.79), similar to the first few days of M19 simulations. The SLP for P19 simulations ( Figure A3d in Appendix A) in winter was higher than SLP for P18 simulations in the late spring ( Figure A3b in Appendix A). The findings were in agreement with previous studies that showed, in winter, there were high correlations between the synoptic scale pressure systems and pollutant mixing ratios [70,71]. The high pressure system in the lower elevation combined with low wind speed and temperature inversion could prevent the vertical dispersion of the GHGs in the air and therefore the GHGs were not transported significantly away from the source. In M19 and P19 simulations no notable correlations were found between the advective flux and any other meteorological variable. Overall, among all the meteorological variables considered, the wind speed was most highly correlated with the advective flux for all cases, in agreement with previous findings concerning short-range landfill methane emission fluxes [72].
To better understand the diurnal and seasonal variation of the advective flux statistically, the simulated hourly advective fluxes were grouped in four-hourly time intervals 0000-0400, 0400-0800, 0800-1200, 1200-1600, 1600-2000, and 2000-2400 local standard time (LST = UTC − 7) over the 10 days of M18, P18, M19 and P19 simulations. The Normalized Average Difference (N AD) was calculated to determine the difference in average advective fluxes (F CH 4 ) between four sets of model simulations, i.e., M18-P18, M18-M19, M19-P19, and P18-P19. The two-sample t-test was employed to investigate if the differences between the average total (or the four-hourly) advective fluxes (F CH 4 ) were statistically significant. Statistical p values were used to check the statistical significance for a confidence level of 95%. The t-test rejected (R for p < 0.05) or accepted (A for p ≥ 0.05) the Null Hypothesis (NH) that there was no difference in advective fluxes between two cases in question. Table 2 shows the values of N AD, and the two-sample t-test results.
As can be seen in Table 2, the average total (combining all the time intervals) advective flux (F CH 4 ) from the mine for the M18 simulation was slightly higher than that from the pond for the P18 simulation, but the difference was not statistically significant (N AD = 0.131, NH = A). The average total advective flux from the mine for the M18 simulation was significantly higher than that of M19 simulations (N AD = 0.356, NH = R). The average total advective flux from the mine for M19 simulation was significantly greater than that from the pond P19 simulations (N AD = 0.542, NH = R). The average total advective flux from the pond for the P18 simulation was significantly higher than that of P19 simulations (N AD = 0.747, NH = R). Overall, it was found that the average total advective flux in late spring was statistically significantly greater than those in winter and early spring by 36% and 75%, for the mine and pond, respectively. The four-hourly analysis revealed that there could be a statistically significant difference (NH = R) in the average advective flux from the mine versus the pond in the same season or from the same source in different seasons. Concerning M18 and P18 simulations, it can be seen that the mine advective flux could be less than the pond advective flux under thermally stable conditions by up to 62% (0000-0400 LST) while being greater under thermally unstable conditions by up to 31% (1200-1600 LST). Concerning M19 and P19 simulations, on the contrary it can be seen that the mine advective flux could be greater than the pond advective flux for most diurnal times by up to 89% (all time intervals except for 1200-1600 LST). This could be explained by dominance of synoptic weather conditions in the winter and early spring as opposed to local weather conditions in the late spring. Concerning M18 and M19 simulations, winter and early spring advective flux from the mine could be greater than that in late spring for some diurnal times by up to 124% (0000-0400, 0400-0800, 2000-2400 LST) but less at other times up to 126% (0800-1200, 1200-1600 LST). Concerning P18 and P19 simulations, the late spring advective flux from the pond could be greater than that in the winter and early spring for some diurnal times by up to 112% (1200-1600, 1600-2000 LST). Table 2. Normalized Average Difference (N AD) of WRF simulated advective flux (F CH 4 ); two-sample t-test to check the null hypothesis (NH) that there was no difference in average advective fluxes (F CH 4 ) between two cases in question; accept (A) meant there was no difference, while reject (R) meant there was a difference with a confidence of 95%; n is the sample size; the time intervals in local standard time (LST = UTC − 7).      Figure 5a shows the predicted total column mixing ratio of methane released from the mine during the afternoon at 1400 LST on 27 May 2018. Influenced by westerly winds (Figure 6a), during this period the methane plume appearance was primarily on the eastern part of the mining facility. During the mid afternoon hours in late spring the boundary layer was thermally unstable, characterized by higher wind speeds and temperatures ( Figures A1a and A2a in Appendix A), and the methane plume was well mixed within the convective boundary layer up to 1-2 km Above Grade Level (AGL), in agreement with previous aircraft observations [26,28,29], as can be seen in the eastern boundary of Figure 7a. Figure 6a shows that temperature and wind speed over the mine pit were higher than those over the tailings pond, which can be explained on the basis of higher heat capacity of the water bodies in comparison to the surrounding land [17,73]. Figure 5b shows the predicted total column mixing ratio of methane released from the mine during the nighttime at 0200 LST on 20 May 2018. Influenced by southerly winds (Figure 6b), during this period the methane plume was seen on the northern part of the mining facility. At this time the boundary layer was thermally stable, characterized by the lower vertical height of the plume under 150 m AGL (Figure 7b) and lower wind speed and temperatures (Figure 6b). At this time, the temperature over the pond was higher than that over the mine [17,73].  (Figures 7g,h). Under synoptic weather conditions in winter and early spring, it was noted that the plume rose up to no more than 400 m AGL. This was supported by the fact that the boundary layer exhibited conditions closer to a thermally stable state as evidenced by snow cover on the ground (not shown), lower temperatures ( Figure 6 and Figure A2 in Appendix A), and lower Sensible Heat Flux (SHF) (Figures A1 and A4 in Appendix A) compared to late spring.

Model Evaluation against Aircraft Observations
Limited aircraft observations of the methane mixing ratio in the PBL were available to compare with the WRF predictions. The Convair-580 aircraft belongs to the National Research Council (NRC) of Canada, and methane measurements were performed using a Picarro model G2401-m instrument. Aircraft measurement details (including those of methane and state parameters) were detailed elsewhere [28,29]. The aircraft flew in a box flight pattern around the facility on 31 May 2018 and measured the methane mixing ratio with a sampling frequency of 0.5 Hz. On this day the methane plume was measured on the west face of the box flight and the measurement was centred around 1100 LST (the box flight typically took 1 to 2 h to complete). The aircraft minimum altitude on this flight was about 550 m AGL and the maximum altitude was 1850 m AGL. For model validation, the WRF simulations were run using four-hourly-averaged methane mixing ratios over both the mine and the pond in May 2018 as boundary conditions, which were updated in the model every four hours. This was needed since in reality there would be simultaneous emissions from both locations although they were not measured simultaneously. The WRF solution was extracted at 1100 LST in grid cells closest to the aircraft latitude, longitude, and altitude positions. In other words, for every aircraft measurement at a frequency of 0.5 Hz, a solution was sampled from WRF. Figure 8 shows the comparison of aircraft versus WRF predictions of methane mixing ratio after removing the background mixing ratio from the methane measurements of the aircraft. The bias and root mean square error for this comparison were 0.0543 and 0.0530 ppm, respectively. The coefficient of determination R 2 for this comparison was 0.68. This evaluation provided confidence in the ability of the WRF model to simulate the plume transport adequately.

Sources of Uncertainty
Although the methodology in this study enabled simulation of diurnal and seasonal variations of area-fugitive methane advective flux from the open-pit mining facility with some confidence, various sources of uncertainty could be proposed for the current approach. The meteorological predictions of the WRF model for this mining facility were evaluated in detail against near-surface-level and PBL observations in a previous study [17]. The methane mixing ratio predictions of the WRF model was also compared against limited aircraft observations for the current facility in this study. Nevertheless, the authors did not possess rigorous methane mixing ratio observations within PBL over many diurnal times and seasons to evaluate the advective flux predictions in more detail. Most satellite products considered (e.g., TROPOMI, GOSAT, and GHGSat) did not provide total column methane mixing ratio measurements at nights or co-incident with near-surface-level observations using LGRs for the open-pit mining facility. Finding methane mixing ratio observations within PBL in the winter was even more difficult due to lack of satellite data products or aircraft measurements. In the current model it was assumed that each of the four measured methane mixing ratios near surface was constant within a box area, however, in reality the mixing ratio could vary at finer spatial scales, so deployment of more near-surface-level instruments measuring the methane mixing ratio could enhance the accuracy of the flux predictions. Current research effort for this mining facility aims to deploy many tens of such near-surface-level instruments for future predictions. The WRF model is sensitive to grid spacing changes [17,74]. It was beyond our scope to provide sensitivity results for changing the model grid spacing, although for the PBL scheme used, the horizontal grid spacing of 506.7 m in the inner domain and use of 90 vertical levels were appropriate and also frequently used by other relevant studies in the field [68].

Conclusions
To the authors' knowledge, this is the first modelling study to quantify the diurnal and seasonal variation of area-fugitive advective flux of methane from a complex open-pit mining facility in northern Canada using WRF. The methane mixing ratio boundary conditions for the WRF model were obtained from the in-situ field measurements in various locations surrounding the mine pit and a tailings pond for two seasons in late spring in 2018 and winter and early spring in 2019. Methane transport was simulated using passive tracer dispersion in WRF. To find the variation in the advective flux under different diurnal times, the simulations were conducted for blocks of four-hour time intervals, where methane mixing ratio boundary conditions were updated in the model according to the observations, while meteorological fields were updated according to the Global Data Assimilation System (GDAS) dataset. The dispersion and advective flux patterns from an active mining facility were confounded by the complex topography, variation in land use, and meteorological effects at micro and mesoscales.
It was found that there were significant (p < 0.05 or with 95% confidence level) diurnal and seasonal differences in the methane advective flux from the mine and the pond. In late spring, advective flux was influenced by local weather conditions, while in winter and early spring, the advective flux was influenced by synoptic weather conditions, in agreement with previous findings [70,71]. Among the meteorological variables considered, the advective flux was most highly correlated with wind speed (winter and spring), in agreement with previous findings [72], but the advective flux could also be moderately correlated with other variables such as sensible heat flux near surface (late spring). Overall, it was found that the average total advective flux in late spring was statistically significantly greater than those in winter and early spring by 36% and 75%, for the mine and pond, respectively. The diurnal variation of the advective flux was more notable in late spring, characterized by lower advective flux at nighttime under thermally stable conditions and higher advective flux during daytime under thermally unstable conditions. In late spring, methane plume visualizations indicated a plume rise up to the top of the PBL height of 1-2 km AGL under thermally unstable conditions, in agreement with previous findings [26,28,29], while the visualizations indicated a plume rise up to 150 m and 300 m AGL under thermally stable conditions for the mine and pond, respectively. In winter and early spring, methane plume visualizations indicated a lower plume rise up to 400 m AGL due to the boundary layer being closer to the thermally stable condition. The WRF model predictions of methane mixing ratio were evaluated against aircraft observations on midday on 31 May 2018. The bias and root mean square error for this comparison were 0.0543 and 0.0530 ppm, respectively. The coefficient of determination R 2 for this comparison was 0.68.

Implications
Although this approach did not estimate the emission flux at the source (e.g., an open-pit mine or a tailings pond), it enabled estimation of the advective flux downwind of the source as influenced by meteorological effects. This approach motivates a paradigm shift in GHG flux monitoring and prediction from area-fugitive sources. Given that measuring the emission flux near the surface over complex heterogeneous terrains is inherently difficult [50,51], alternatively the advective flux can be estimated downwind using near-surface-level observations of GHG mixing ratios at the source with the aid of mesoscale modelling. This can be justified due to the fact that the storage and release of GHGs in the mesoscale modelling domain (enclosing the source) average to zero over multiple days, so the advective flux (in the absence of GHG chemical reactions, surface deposition, or biogenic generation) ultimately provides a proxy or estimate for the amount of GHGs that enter the atmosphere. The methodology used can be applied to other situations where large scale land surface modifications require micro and mesoscale Numerical Weather Prediction (NWP) in dispersion modelling of atmospheric pollutants and quantification of area-fugitive advective fluxes of GHGs. This approach addresses a need for utilization of mesoscale models in GHG flux quantification. Most existing mesoscale models are not currently utilized to quantify flux. Rather, they take the emission flux as input (via emission inventory datasets) and predict the mixing ratio of pollutants throughout the atmospheric domain [43,[45][46][47][48][49]. While this is useful in air quality studies, it does not address the issue of quantifying flux. Instead, the proposed method in this study provides a mechanism for quantification of GHG fluxes, as they vary diurnally and seasonally, using mesoscale models, forced by near-surface-level observations of the GHG mixing ratio. This approach may also eliminate the need for resource-intensive top-down approaches in quantifying GHG fluxes using satellites, aircraft, and drones [23,25,26,29,31,32]. The proposed method only relies on near-surface-level measurement of GHG mixing ratio, while it relies on mesoscale modelling to derive the advective flux of GHGs.

Future Work
The method can be improved by ingesting boundary conditions for the methane mixing ratio at higher spatio-temporal resolutions. To achieve this, either high-resolution surface-level or satellite-based observations can be used. The method can also be improved by more rigorous evaluations against PBL observations of the methane mixing ratio as it varies diurnally and seasonally.