Comparative Study of CALPUFF and CFD Modeling of Toxic Gas Dispersion in Mountainous Environments

Abstract

Verifying the pattern of toxic gas dispersion simulations under mountainous conditions is vital for emergency response and rescue. In this study, a comparative analysis is conducted between CALPUFF (California Puff Model) and CFD (Computational Fluid Dynamics) gas dispersion modeling focusing on the range of Semi-Lethal Concentration (LC₅₀) and Immediate Danger to Life and Health Concentration (IDLH). To identify general dispersion patterns, a hypothetical pipeline breakout accident in a mountainous area is simulated and thirteen groups of simulation conditions are set up for the experiments, including calm wind (velocity less than 0.5 m/s) and winds from the east (E), south (S), west (W), and north (N) at velocities of 1, 2, and 3 m/s with a 1 arc-second degree SRTM data as terrain data. Comparative experiments show the diffusion patterns of the two models are essentially consistent, and the overall dispersion range deviation between two methods is within 266 m. The evaluation of CALPUFF’s adaptability for microscale mountainous environments indicates its potential use for high-sulfur gas fields and gas dispersion simulations in emergency scenarios.

1. Introduction

China’s high hydrogen sulfide gas reservoirs are mainly located in the Sichuan Basin, producing raw gas with up to 14% hydrogen sulfide content. Risk assessment, leakage monitoring, safety control, and emergency response are vital components for ensuring Quality, Health, Safety, and Environment (QHSE) in high-sulfur gas field production. A toxic gas leakage accident will have a significant impact on the surrounding environment and may also cause widespread social impacts [1,2,3]. Therefore, it is an essential task to carry out gas dispersion simulation and establish emergency response and evacuation plans in case of toxic gas leakage accidents.

CALPUFF is a well-established environmental quality regulation guideline supported by the U.S. Environmental Protection Agency (EPA) and is particularly suitable for long-distance pollutant transport (over 50 km) [4,5]. In recent years, the CALPUFF atmospheric dispersion model has been successfully utilized for small-scale gas plant and pipeline scenarios in emergency situations [6]. This model offers faster calculation speeds and enhances calculation accuracy and the capability to handle complex terrain and diverse meteorological conditions [7,8]. However, there are still few studies validating the CALPUFF simulation results of toxic gas pipeline leakage under complex topography conditions. Therefore, we aim to compare the simulation results of CALPUFF and CFD models under complex terrain and meteorological conditions to verify the pattern of toxic gas dispersion simulations and evaluate the CALPUFF model for various pipeline leakage emergency scenarios.

Validating atmospheric dispersion simulation is a significant challenge. As for model evaluation, we strive to identify the model that most accurately estimates concentrations in the area of interest. Field experiments are considered the most reliable, as they closely replicate real-world conditions. Dedicated gas dispersion field experiments have been conducted, resulting in datasets that can be used for this purpose [9]. For instance, dispersion models have often compared observations from field sampling sites to model simulations [10,11,12]. However, the limited number of field experiments cannot provide a comprehensive and reliable evaluation of the models. There are natural fluctuations in the concentration fields due to random turbulence in the atmosphere [13].

A comparative study between wind tunnel experiments and computer simulations is a more feasible method. Atmospheric dispersion models have been validated using wind tunnel experimental data [14,15,16,17]. However, atmospheric dispersion is influenced by many complex factors such as meteorological conditions, topography, land use types, dry or wet deposition, leak source parameters, etc. Different experimental methods and models lead to different conclusions of environmental assessment, safety risks, and emergency evacuation.

In many studies, two or even more air dispersion simulation models are chosen for comparison and evaluation. For example, studies on the applicability of the AERMOD model and the CALPUFF model mostly focus on specific cases, such as urban areas [18,19], mountainous regions [20,21], industrial parks [22] and so on. Although it is hard to analyze model results from the perspective of data accuracy, comparisons are necessary to analyze the advantages, disadvantages, and adaptability of different models. The model performance can be evaluated from different aspects, such as the spatial distribution of gas dispersion, the diffusion time of specific locations, concentration peaks, the total concentration during the entire dispersion process, etc. If these different methods produce similar results, it indicates that these methods can be used to simulate the actual patterns of gas dispersion.

This study examines a hypothetical gas pipeline breakout incident that involves the release and dispersion of high-sulfur gas characterized by complex topography. A qualitative and quantitative comparison is conducted between the CALPUFF model and the CFD model to assess the adaptability and reliability of the CALPUFF model within the context of emergency response management systems.

2. Materials and Methods

2.1. Study Area

The study area, shown as the yellow line shaded area in Figure 1, is situated in the mountainous and hilly terrain, characterized by numerous gullies and significant changes in relative elevation. It covers approximately 110 ${\mathrm{km}}^2$ and spans 4 townships, Xiaba, Nanba, Tahe and Gaoqiao. There is a combined resident population of approximately 90,000. The well sites, pipelines, and valve chambers are close to high-density populated areas, posing significant safety and environmental risks and challenges. In Figure 1, the emergency management zones within 1.5 ${\mathrm{km}}^2$ (yellow line) and 2.5 ${\mathrm{km}}^2$ (purple line) around the pipeline are shown.

The terrain generally ascends from southwest to northeast, marking the transition from the basin hills and low mountains to the middle mountain. The rising and stripping of low mountains and hills dominate the area. The terrain is undulating with altitudes ranging between 340 m and 1000 m.

The region experiences a subtropical humid monsoon climate, which is characterized by an average temperature of 16.8 °C, annual rainfall of 1239.4 mm, and an average wind speed of 1.5 m/s. In spring, summer and autumn, the prevailing wind is from the northeast. In winter, the prevailing wind is from the east. The dominant wind direction is northeast with a relative humidity of 77%. The area has a long frost-free period, four distinct seasons, and limited sunshine.

As for land use, forest land accounts for the largest portion of the total land area at 37.34%, followed by arable land at 25.55%, and unused land at 22.05%.

This hypothetical leakage accident locates near Qili Village (31.367° N, 108.213° E) in Kaizhou District, Chongqing, as shown in Figure 2. The blue circle represents a radius of approximately 500 m from the leakage point. The terrain in this area is very rugged with many mountains and forests and little flat land. It mainly consists of low mountains and hills with deep valleys and narrow ravines. Roads run along those ravines, and there are dense houses around them. Qili Village is surrounded by various residential areas, including kindergartens, village committees, clinics, primary schools, junior high schools, and nursing homes. These sensitive points are located less than 1000 m away from the potential leak source, which is shown as a red circle. In the event of a leakage accident, the consequences could be catastrophic.

2.2. Assumed Scenario

Assuming that the pipeline leakage accident occurred on the hillside between the two valve chambers Q and S, in an extreme scenario, the gathering pipeline undergoes a total break. Due to the pressure drop fast in the pipeline, the Emergency Shutdown Device (ESD) system monitored the leakage signal of the pipeline within 26 s, the valves closed automatically, and the separation of the pipeline and other pipelines was completed in 20 s. In view of the fact that a certain amount of time is remaining, the isolation of the upstream and downstream valve chambers of the pipeline was completed within 50 s after the total breakage accident occurs. The total leakage volume is the leakage volume before the emergency cut-off valve of the pipeline section adds the residual volume of the pipeline section between the cut-off valves after 50 s of leakage. The emission source of the pipeline can be detected by the sensors. The emission strength, volume, location and time can be obtained precisely.

According to the above settings, the phase model developed by Det Norske Veritas (DNV) is used, the instantaneous emission is used to calculate the leakage rate from the pipeline, and the curve of the leakage emission rate is obtained, as shown in Figure 3. In this figure, the risk assessment chooses the most extreme scenario of a full pipeline breakout, and the leakage source emission rate between the two valve chambers Q and S is calculated.

The overall assumed leakage parameters are as follows:

The pipeline length is 5 km, the pipe diameter is 500 mm; the pressure is 7.3 MPa; The volume percentage of $C{H}_4$ is 82.71%, and that of $H_{2}S$ is 9.73%; The total estimated leaking volume of $H_{2}S$ is 60,566 kg, and the total leakage time is 185 s.

2.3. Modeling Principles

2.3.1. CALPUFF Model

The CALPUFF model is an environmental quality regulation model that has long been supported by the U.S. Environmental Protection Agency (EPA) [4]. As an important air quality simulation tool, it has played a significant role in environmental protection and atmospheric science research. The CALPUFF model is an improved Lagrangian Gaussian model. It is assumed that the concentration contribution of puff at a receptor point follows Gaussian Equations (1) and (2).

$$C={\displaystyle \frac{Q}{2\pi {\sigma }_x{\sigma }_y}}gexp[-d_a^2/\left(2{\sigma }_x^2\right)]exp[-d_c^2/\left(2{\sigma }_y^2\right)]$$

(1)

$$g={\displaystyle \frac{2}{{\left(2\pi \right)}^{1/2}{\sigma }_z}}\sum _{n=-\infty }^{\infty }exp\left[-{\left(H_e+2nh\right)}^2/\left(2{\sigma }_z^2\right)\right]$$

(2)

where C is the ground level pollutant concentration at the receptor site (g/m³); Q is the mass of the pollutant slug (g); ${\sigma }_x$, ${\sigma }_y$ and ${\sigma }_z$ are the standard deviation of the Gaussian distribution of the pollutants in the X, Y, and Z directions, respectively (m); $d_a$ is the distance from the center of the smoke mass to the receptor point in the X direction (m); $d_c$ is the distance from the center of the smoke mass to the receptor point in the Y direction (m); g is the vertical term of the Gaussian equation (m); $H_e$ is the effective height of the center of the smoke mass above the ground (m); and h is the height of the mixed layer (m). The vertical term g explains the multiple reflections between the mixed layer and the ground. If ${\sigma }_z>1.6$ h, g will be simplified as $1/h$. Smoke masses in the general convective boundary layer satisfy this condition several hours after release.

Equation (1) shows that the CALPUFF model uses a Gaussian approach to treat the interior of the puff. The CALPUFF model describes the diffusion of pollutants as a normal distribution along the y-axis and z-axis directions. The diffusion process is divided into puff clusters over multiple time intervals. For a specific reception point, the total concentration is calculated by summing the averages of all the sampling intervals of the surrounding smoke plumes over the basic time step. The model’s outcomes are comparable to those of a Gaussian model, but instead of a simple straight line, the center line follows the trajectory curve of the wind, which changes in real time. Therefore, the model combines the characteristics of the Lagrangian and Gaussian models, making it an improved version of the Gaussian model that uses the Lagrangian method. It is also known as the Lagrangian–Gaussian puff model.

The CALPUFF modeling system consists of four parts: (1) geographic and meteorological pre-processors; (2) the CALMET meteorological module; (3) the CALPUFF dispersion module; and (4) the CALPOST post-processor. The workflow is shown as Figure 4.

Unlike AERMOD or other Gaussinian models, the CALMET/CALPUFF modeling system has the ability to model calm hours by simulating stagnant puffs. In meteorology, calm wind refers to wind speeds less than or equal to 0.5 m/s. In this situation, the wind has very little impact on environmental temperature and air quality, and smoke and pollutants may linger in the air for a longer time. Stagnant puffs are not dispersed via advection (since the wind speed is zero), but they may still undergo turbulence-related dispersion. Furthermore, even if the measured wind speed is zero, CALPUFF accounts for other possible flow components, e.g., puff transport caused by divergence or slope flow. Therefore, the model will calculate concentrations during calm periods.

Generally, complex terrain affects the 3D wind field in the CALMET meteorological module. The terrain effect will be manifest in the boundary conditions for the flow field. GEO.DAT in CALMET stores various geographic data required for model calculations, including land use types, topographic data, and some surface parameters (such as surface roughness, albedo, Bowen ratio, soil heat flux parameters, and vegetation leaf area index) as well as human activity heat flux. The topographic data are recorded in a gridded format. CALMET performs diagnostic adjustments to the wind field based on the topographic data to account for terrain effects on a fine grid. Specifically, the method proposed by [23] is used to assess the dynamic terrain effects. The regional scale wind field is utilized to calculate the vertical wind speed affected by the topography and meets the exponential stability decay function. The dynamic terrain effects on the horizontal wind components are assessed by applying a least-squares minimization method to the initial guessed wind field. The least-squares minimization method is iteratively applied until the three-dimensional deviation is below a threshold.

What’s more, the complex terrain module in CALPUFF is based on the approach used in the Complex Terrain Dispersion Model (CTDPLUS). This method is to simulate the puff which is too small to influence the large-scale wind field. Smaller-scale terrain features encountered by a puff in this flow can be simulated by a separate subroutine. For some special geomorphological conditions, such as isolated mountains, the model can be set up with secondary terrain to simulate obstacles to the airflow.

2.3.2. CFD Principles

Computational Fluid Dynamics (CFD) is a commonly used method for studying the atmospheric dispersion pattern following hydrogen sulfide leakage from a natural gas pipeline. The primary objective is to elucidate the temporal changes in the density, pressure, velocity, and temperature of hydrogen sulfide gas. The following assumptions are made: after the leakage of a natural gas pipeline, hydrogen sulfide gas is in a state of high pressure and high velocity flow. However, only the exchange and dispersion with the surface air are considered, which is consistent with the standard model of turbulence distribution. Furthermore, the gas flow in the three-dimensional space is considered to be a constant temperature transport.

(1) Mass conservation equation

$${\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial w}{\partial z}}=0$$

(3)

where u, v, and w are the components of the velocity vector in the x, y, and z directions.

(2) Turbulent pulsation kinetic energy equation

$${\displaystyle \frac{\partial \left(\rho k\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\rho k{v}_i\right)}{\partial x_i}}={\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\beta }_k}}\right){\displaystyle \frac{\partial k}{\partial x_j}}\right]+G_k-\rho \epsilon -Y_m+S_k$$

(4)

Turbulent kinetic energy dissipation equation

$${\displaystyle \frac{\partial \left(\rho \epsilon \right)}{\partial t}}+{\displaystyle \frac{\partial \left(\rho \epsilon v_i\right)}{\partial x_i}}={\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\beta }_{\epsilon }}}\right){\displaystyle \frac{\partial \epsilon }{\partial x_j}}\right]+{\displaystyle \frac{\epsilon C_1{G}_k}{k}}-{\displaystyle \frac{C_2\rho {\epsilon }^2}{k}}+S_t$$

(5)

$${\mu }_t=\rho C_0{\displaystyle \frac{k^2}{\epsilon }}$$

(6)

where $\rho$ is the fluid density, kg/m³; k is the turbulent kinetic energy, J/kg; $\epsilon$ is the dispersion rate, m²/s³; $S_k$ is the source term of the turbulent kinetic energy equation; $S_t$ is the source term of the dissipation rate equation; $v_i$ is the hydrogen sulfide air flow velocity in the direction of grid i, m/s; $x_i$ and $x_j$ correspond to the i, j directions of the grid, respectively; $G_k$ is the turbulent kinetic energy generated by the velocity gradient; ${\mu }_k$ is the dynamic viscosity of the fluid; ${\beta }_k$, ${\beta }_{\epsilon }$, is the corresponding value for turbulent kinetic energy and turbulent dissipation rate, Planck’s number; Finally, $C_0$, $C_1$, and $C_2$ are empirical constants.

(3) Conservation equation for components

$${\displaystyle \frac{\partial c_{\eta }}{\partial t}}+▽·J_{\eta }+v_i·▽c_{\eta }=0$$

(7)

$$J_{\eta }=-D_i·▽c_{\eta }$$

(8)

where $c_{\eta }$ is the molar concentration of component $\eta$ (mol/m³); $J_{\eta }$ is the diffusive molar concentration of component $\eta$ (mol/(m³·s)); and $D_{\eta }$ is the diffusion coefficient of the gas $\eta$ (m²/s).

The most commonly used software for CFD calculations is COMSOL (COMSOL), Fluent (ANSYS), FLACS (Gexcon), and other similar software. In this study, COMSOL Multiphysics 6.2 is selected for the simulation and solution processes. The platform product for creating physics-based models and simulation applications. The main steps are as follows, as shown in Figure 5:

(1) Construct the physical model. The topographic map is used to create a Digital Elevation Model (DEM) with a range of 64 ${\mathrm{km}}^2$. These data then are imported into Rhino, a three-dimensional modeling software, to generate a three-dimensional surface physical model and a leakage pipeline model. Finally, the model is imported into COMSOL software for Boolean operations to create a 6 km × 6 km × 2.5 km cubic model, which represents the surface above the earth’s surface. This will result in a vertical space of approximately 300–2200 m above the ground surface, which is established as the simulation range for the experiment.

(2) Set the physical field. The standard turbulence Navier–Stokes model is selected and coupled with the component models. The mixed materials of hydrogen sulfide and air condition are set, as are the properties of the inlet, outlet and wall of the fluid region.

(3) Setting Grid Division. The tetrahedral grid is employed as the fundamental division unit with local encryption of the grid at the mouth of the leakage pipe.

(4) Model Simulating. The simulation is conducted using a pressure, concentration solver multi-field coupling for transient solution.

(5) Export and post-process simulation results. The calculation results can be displayed in a variety of formats, including contour, gradient, vector, particle flow traces, three-dimensional slicing, transparent and semi-transparent, and other graphical methods through the post-processing module.

The principle of CFD model involves establishing continuity equations, turbulence models, and component conservation equations under various conditions, combined with initial and boundary conditions, applying the finite element theory to solve the Navier-Stokes equations, and thus achieving the prediction of the spatial and temporal patterns of the hydrogen sulfide-air field during leakage. In this paper we employs the Reynolds-Averaged Navier-Stokes (RANS) model. The RANS model introduces a roughness wall function, subsequently modifying the traditional k-$\epsilon$ two-equation model, to fully consider wind field variations in complex terrain and landforms. Generally, for urban areas with a few trees or buildings, the surface roughness is set to 0.03 m; for snowy mountainous areas, it is set to 0.001 m; and for densely forested mountains, it ranges from 15 to 25 m.

2.4. Qualitative and Quantitative Comparative Analysis

2.4.1. Ranges of the LC50 and IDLH

Hydrogen sulfide $H_{2}S$ is a potent neurotoxicant that can cause irritation to skin and mucous membranes. It is also an explosive gas with an explosion limit range of 4.0% to 46% (volume ratio). The natural gas in our research contains 8.28–14.25% hydrogen sulfide. $H_{2}S$ has a significant impact on the human body. Once it reaches lethal concentration, the human body quickly loses consciousness, and without prompt rescue, respiration will cease, leading to death.

In the field of environmental impact assessment, LC₅₀ is defined as the concentration of a substance that causes unconsciousness after short-term exposure (30 min), cessation of breathing without prompt treatment, dizziness, loss of consciousness and balance, and the need for immediate artificial respiration and/or CPR techniques. IDLH (Immediately Dangerous to Life or Health) concentration is defined as follows: if a worker does not wear a mask or lacks experience in escaping, and the gas concentration in the work environment reaches a dangerous level, a 30-min exposure will cause permanent damage to the human body or reduce the level of human health.

For this research, we assume an emergency response time of 30 min. The LC₅₀ threshold of $H_{2}S$is 986 ${\mathrm{mg}/\mathrm{m}}^3$, which is consistent with toxicity data. However, for consistency with the predicted evaluation results of the Environment Impact Assessment Report, we retain the LC₅₀ threshold of $H_2$S (618 ${\mathrm{mg}/\mathrm{m}}^3$) and the concentration value of serious injury, IDLH (432 ${\mathrm{mg}/\mathrm{m}}^3$), with an exposure time of 30 min.

To evaluate the consistency of concentrations obtained from the CALPUFF model and CFD, we compare the contours of surface layer concentration results output from both models. The ranges of gas dispersion for different hazardous levels are compared, specifically the maximum ranges of the LC₅₀ and the IDLH.

2.4.2. Statistical Performance Measures:FB, NMSE, R, MAE, PBIAS

To evaluate the gas concentration obtained from CALPUFF model simulation and CFD experiments, we use the quantitative indexes proposed by Hanna and Chang [13]. Additionally, a resampling method for predicted data is proposed, and the predicted concentration values at discrete points are quantitatively analyzed. We select Fractional Bias (FB) and Normalized Mean Square Error (NMSE) as indicators.

$$FB={\displaystyle \frac{2({\overline{X}}_o-{\overline{X}}_p)}{{\overline{X}}_o+{\overline{X}}_p}}$$

(9)

$$NMSE={\displaystyle \frac{\overline{{(X_o-X_p)}^2}}{{\overline{X}}_o{\overline{X}}_p}}$$

(10)

where P and O represent the CALPUFF model simulation and CFD experiment, respectively. They denote the mean values of the respective sampling points ${\overline{X}}_o$ and ${\overline{X}}_p$. The FB and NMSE indicators reflect the overall systematic bias of the data. In the ideal case, when the model and the experimental results are in perfect accordance with each other, the values of FB and NMSE are equal to 0.

Additionally, a correlation coefficient R is selected to reflect the correlation between the two groups of data. When $\left|R\right|$ is closer to 1, it indicates a stronger linear correlation between the two sets of data. Conversely, when $\left|R\right|$ is closer to 0, it indicates a weaker linear correlation.

$$R={\displaystyle \frac{\overline{(X_o-\overline{X_o})(X_p-\overline{X_p}})}{{\sigma }_o{\sigma }_p}}$$

(11)

where ${\sigma }_o$ and ${\sigma }_p$ are the standard deviations of the sample points for the CALPUFF and CFD simulations, respectively.

$$MAE=\overline{\left|\overline{X_o}-\overline{X_p}\right|}$$

(12)

$$PBIAS={\displaystyle \frac{\overline{\left|\overline{X_o}-\overline{X_p}\right|}}{\overline{\left|\overline{X_p}\right|}}}\times 100$$

(13)

MAE and PBIAS also are selected as indicators for quantitative analysis. The Mean Absolute Error (MAE) represents the average absolute error between the values predicted by CALPUFF and CFD. It measures the mean of the absolute differences between the two sets of values. Percent Bias (PBIAS) reflects the cumulative deviation between CALPUFF and CFD models. When the models align well, PBIAS provides a more accurate assessment of the overall model performance. An optimal PBIAS value of 0 indicates that the model is in good overall agreement.

3. Results

3.1. CALPUFF Simulation

Table 1. CALPUFF simulation parameters.

Parameter Type	Parameter Name	Value
Source attributes	UTM coordinates (km)	(234.877 3473.580)
	Pollutant	$C{H}_4+H_{2}S$
	Pollutant volume concentration	9.73% $H_{2}S$
	Chimney height (m)	0
	Outlet initial velocity (m/s)	1
	Emission rate (g/s)	0.49 s: 300,000 50–184 s: 190,000 after 185 s: 0
	Diffusion type	Point source
Meteorological parameters	Air velocity (m/s)	0, 1, 2, 3
	Wind direction (°)	0, 90, 180, 270
	Temperature (°C)	7
Simulation time setting	Wind field time step (s)	3600
	Concentration field time step (s)	60
	Total duration (h)	2
Simulation range setting	X-direction length (km)	10
	Y-direction length (km)	10
	Grid size (m)	50

shows the detailed parameters for the dispersion simulation of CALPUFF, including pollution sources, meteorological conditions, and simulation settings. Ground observation and MM5 (Fifth-Generation Penn State/NCAR Mesoscale Model) are selected for meteorological data, and the horizontal grid spacing, vertical stratification, and number of grid points of CALPUFF were the same as those of CALMET. The meteorological module CALMET can generate three-dimensional gridded meteorological data through sophisticated treatment and the assimilation of surface/upper air/precipitation observations, including prognostic wind field data from the mesoscale models such as MM5. In our study, surface meteorological data such as wind speed, wind direction, temperature, cloud cover and relative humidity are directly recorded as ‘SURF.DAT’. Upper meteorological data such as wind speed, wind direction, temperature, pressure, and elevation are set directly as ‘UP.DAT’. Pseudo-upper station data can be extracted for a set of user-defined MM5 grid points of interest. CALMM5 is a processor that extracts and interprets data in the output file from MM5 and creates a file of meteorological data named ‘MM5.DAT’. The MM5 mesoscale numerical model employs a double-nested (DOMAIN1, DOMAIN2) configuration with horizontal grid spacings of 56 km and 18 km, respectively. Vertically, the research MM5 data are divided into 23 layers, centered at a latitude and longitude of 104° E and 30° N, with grid point counts of 60 × 60 and 72 × 82, respectively. The integration step is set to 90 s. MM5 background data are derived from the T231 forecast products issued daily for the corresponding period and every six hours. For objective analysis, data are from the National Meteorological Center’s twice-daily upper air and eight-times-daily surface original bulletin.

The CALPUFF receptor points had a grid spacing of 100 m, a time step of 1 min, and a simulation dispersion time of 60 min. The simulation area was 10 km × 10 km × 3 km with a grid spacing of 50 m. The terrain file used the SRTM1 dataset with a resolution of 30 m, and the land use file used the LUCC dataset with a resolution of 1 km. Thirteen groups of simulation conditions were set up for the experiments, including calm winds from the E, S, W, and N at speeds of 1, 2, and 3 m/s. Choosing four wind directions instead of eight can significantly reduce the computational load, allowing for more efficient numerical calculations and simulations.

3.2. CFD Numerical Simulation

In this experiment, we choose the 6 km × 6 km × 2.5 km area around the pollution source for simulation. In the CALPUFF model, CALMET generates a 3D wind field diagnosis based on the principles of mass conservation and the continuity equation. In CFD simulations, the meteorological parameters, such as wind speed, wind direction, and temperature, are kept consistent with those used in the CALPUFF model. Although CALMET data are not directly utilized in the CFD model, the mass conservation equation and the species transport equation are configured with the same meteorological parameters to simulate the gas flow behavior and the diffusion process of sulfur hexafluoride pollutants. Mesoscale MM5 data have few impacts on the microscale emergency gas dispersion model. There are few studies combining CFD models with mesoscale meteorological models for gas dispersion simulations. Consequently, MM5 data are not involved in COMSOL directly.

The 3D surface physical model and leakage pipeline model were generated by the 3D modeling software Rhinoceros 7 (Robert McNeel and Associates). The Boolean operation was imported into the COMSOL software to cut the 6 km × 6 km × 2.5 km cube to obtain a cubic area for simulation. The area comprises six boundaries: the ground surface, the northwest–southeast air inlet, and the top. The source of the pollutant leakage is situated in the center of the ground, which is a hexahedron with a length of 10 m, width of 10 m and height of 1 m.

Table 2. CFD experimental boundary conditions and parameter settings.

Boundary Name	Boundary Type	Parameter Value
Source vent	Concentration and speed	Concentration was controlled using a segmented function with a $H_{2}S$ leak rate of 1 m/s. The leak concentration was $20.54\times {10}^{-3}$ mol/m3 from 0 to 50 s, and $5.58\times {10}^{-3}$ from 50 to 185 s. After 185 s, there was no leak.
Ground	Zero flux	No slip and no flow boundary, soil material
Outer wall of the source	Zero flux	No slip
West inlet	Speed and Pressure	For the velocity boundary, use 0.1, 1, 2, or 3 m/s. For the pressure boundary, use 0 gauge pressure.
North inlet	Speed and pressure	Same as above
South inlet	Speed and pressure	Same as above
East inlet	Speed and pressure	Same as above
Top	Pressure	Gauge pressure is 0

shows the detailed parameters for the dispersion simulation of COMSOL.

We use tetrahedral grids as the basic unit of division. The grid scale is encrypted with an ultra-detailed grid (cell size from 40.6 to 215 m) used from the ground to the 550 m height and a regular encrypted grid (cell size from 108 to 602 m) used from 550 m to 2500 m height. The model is divided into a total of 694,690 domain tetrahedral cells. The experiment has a total simulation time of 60 min with a 6 s time step.

3.3. Comparative Analysis of LC50 and IDLH Impact Ranges

We measured the maximum distance from the leakage source to the $L{C}_{50}$ and IDLH boundary under 13 different conditions. The influence ranges of $L{C}_{50}$ and IDLH measurements, shown in

Table 3. Comparison of maximum impact distance between Max LC50 and Max IDLH.

Scenario	CALPUFF		COMSOL		Error (CALPUFF-COMSOL)
	${\mathit{LC}}_{\mathbf{50}}$ (m)	IDLH (m)	${\mathit{LC}}_{\mathbf{50}}$ (m)	IDLH (m)	${\mathit{LC}}_{\mathbf{50}}$ (m)	IDIH (m)
Calm	766	847	682	914	84	−67
North 1 m/s	1042	1121	938	1196	104	−75
East 1 m/s	1744	2022	1861	2147	−117	−125
South 1 m/s	2542	2698	2473	2713	69	−15
West 1 m/s	1052	1391	1214	1365	162	35
North 2 m/s	1052	1283	1132	1369	−80	−86
East 2 m/s	1532	1756	1681	1981	−149	−225
South 2 m/s	2002	2604	2069	2338	−67	266
West 2 m/s	1781	1926	1801	2003	−20	−77
North 3 m/s	821	1031	855	1050	−34	−19
East 3 m/s	1323	1572	1167	1422	−156	−150
South 3 m/s	1919	2279	1945	2242	26	37
West 3 m/s	1610	1677	1582	1752	28	−75

, are determined using CALPUFF and CFD simulation results.

Figure 6 shows concentration predictions from the CALPUFF and CFD models of the 30-min slip-averaged maximum concentration of each grid under four wind conditions. This value indicates that the duration of this maximum concentration in each grid is 30 min. The reason for choosing this value is that the time for the semi-lethal concentration ${\mathrm{LC}}_{50}$ of $H_{2}S$ is 30–60 min, and the time for immediate danger to life health (IDLH) time is also 30 min. The black line shows the range of ${\mathrm{LC}}_{50}$ (618 ${\mathrm{mg}/\mathrm{m}}^3$), and the red line refers to the IDLH (432 ${\mathrm{mg}/\mathrm{m}}^3$).

From the topographic map, it can be observed that there are steeper mountains in the north and west of the leakage point than in the south and east, causing toxic gases to diffuse toward the northwest and southwest. It is evident that the steep mountain has a significant impact on the diffusion direction.

Figure 6a shows the CALPUFF concentration range of ${\mathrm{LC}}_{50}$ and IDLH under calm wind. It can be seen that the diffusion range of IDLH under calm wind conditions is basically circular with a diameter of about 847 m. Two independent high concentrations of ${\mathrm{LC}}_{50}$ also accumulate in the valley southwest of the leakage source. Figure 6b shows the CFD diffusion range of $L{C}_{50}$ and IDLH under calm wind conditions by using COMSOL. It can be seen that the diffusion range under calm wind conditions is basically an east–west elongated, north–south shortened ellipse with the maximum IDLH range reaching 914 m, which is similar to the CALPUFF results. Under calm wind conditions, the maximum range error for the two models is 84 m.

Figure 6c shows the CALPUFF concentration range of ${\mathrm{LC}}_{50}$ and IDLH under wind speeds of 1 m/s in four directions. Figure 6d shows the CFD diffusion range under wind speeds of 1 m/s in four directions. It can be observed that the diffusion trends of the two models are similar: both of them are spreading along the valleys, and the farthest diffusion direction is S wind. Under the S wind direction, the toxic gas extends northwest along the valley with the IDLH range reaching a maximum of 2700 m. The maximum range error for the two models is 125 m.

Figure 6e shows the CALPUFF diffusion range of ${\mathrm{LC}}_{50}$ and IDLH under wind speeds of 2 m/s. Figure 6f shows the CFD diffusion range under wind speeds of 2 m/s. The trend of the two models is similar with the largest diffusion remaining S wind. Under S wind, the diffusion extends along the northwestern valley with the IDLH range reaching a maximum of 2600 m. In the N wind, due to the obstruction of the mountains, the diffusion range of IDLH is about 1300 m, which is half that of the S wind. The maximum range error for the two models is 266 m.

Figure 6g shows the diffusion range of ${\mathrm{LC}}_{50}$ and IDLH under wind speeds of 3 m/s. Figure 6h shows the CFD diffusion range under wind speeds of 3 m/s. Due to the relatively high wind speed, the range of concentration of toxic gas is quickly diluted, reducing the diffusion range about 300 m in four directions. The maximum IDLH range under south wind conditions is about 2300 m, extending to the northwest. In the E and W wind, the IDLH ranges also gradually decrease to 165 m and 1670 m, respectively. Due to the obstruction of the mountains in the north wind direction, the diffusion range is only 1050 m. The maximum range error for the two models is 156 m.

From Table 3, it can be seen that the results of the CALPUFF and CFD models exhibit high agreement. In the comparison of maximum diffusion range, under calm wind, N, E, S and W, both models show good consistency with the IDLH range difference about 260 m. As for environmental assessment and emergency management, the maximum range of IDLH and ${\mathrm{LC}}_{50}$ for a leak point using four-wind-direction or eight-wind-direction calculations is very meaningful.

The two maximum diffusion errors, 225 m and 266 m, of IDLH occur under S and E wind at 2 m/s. The error is due to the differences in the principles of the two models considering complex terrain. When calculating pollutant diffusion, CALPUFF modifies the three-dimensional wind field model by setting terrain roughness. In contrast, the CFD method generates a three-dimensional finite element grid and uses a terrain roughness function for calculations. During the calculation process, different methods cause a slight difference in the gas dispersion patterns.

3.4. Quantitative Analysis of CALPUFF and CFD Experiments

Table 4. Quantitative evaluation of CALPUFF and CFD experiments.

Simulation Scheme	FB	NMSE	R	MAE	PBIAS
Calm	1.74	3.54	0.44	0.68	1364.58
North 1 m/s	0.09	0.44	0.05	0.14	9.22
East 1 m/s	−0.44	0.08	0.22	0.38	−35.78
South 1 m/s	1.72	25.13	0.17	3.21	1248.93
West 1 m/s	1.89	9.67	0.57	2.03	3515.17
North 2 m/s	0.51	1.97	0.02	0.20	69.62
East 2 m/s	−0.29	0.12	0.26	0.25	−25.31
South 2 m/s	1.46	10.78	0.21	1.53	539.76
West 2 m/s	1.58	6.03	0.05	0.96	750.78
North 3 m/s	0.78	3.41	0.16	0.20	127.88
East 3 m/s	0.28	0.23	0.29	0.21	32.08
South 3 m/s	1.14	10.18	0.23	1.18	479.84
West 3 m/s	1.70	11.94	0.29	1.21	1167.98

shows the quantitative evaluation of CALPUFF and CFD experiments under different leakage scenarios. The gas concentrations obtained from the CFD experiments and CALPUFF simulations cannot be compared directly. To achieve quantitative data analysis and comparison, the first step involves converting concentration values, which are the units used by the CALPUFF model, to ensure consistency in data units. Next, we employ Kriging interpolation to resample the data. This process will cover an X-coordinate range of 3 km and a Y-coordinate range of 3 km with a grid spacing of 0.05 km × 0.05 km, resulting in the creation of two grid files to maintain the spatial consistency of the data. Finally, Equations (9)–(13) will be used to perform statistical analysis of the FB, NMSE, R, MAE and PBIAS values on the grid points of the resulting data.

The FB and NMSE indices reflect the overall magnitude deviation between the two sets of data with smaller values indicating a higher degree of fit. Ideally, when the model and experimental results perfectly match, the values of FB and NMSE are both equal to 0. As shown in Table 4, the overall magnitude deviation of FB ranges in the E direction, indicating a moderate fit. Other cases show no good fit.

The R index indicates the correlation between the two sets of data with higher values representing greater correlation. Generally, when 0 < |R| < 0.3, there is a weak positive linear correlation between variables; when 0.3 < |R| < 0.7, there is a moderate positive linear correlation; and when 0.7 < |R| < 1, there is a strong positive linear correlation. As seen in Table 4, all cases show weaker correlations. This might be due to the fundamentally different principles of the two models, resulting in weaker data correlations.

Based on the NMSE, MAE and PBIAS values in the results of CALPUFF and CFD, the models for the south and west wind directions are quite large. Although from Table 3, the diffusion direction and maximum diffusion distance are similar, the diffusion patterns are different. Therefore, there is considerable variation under the south and west wind conditions. We need a further analysis and improvement in the model’s consistency.

From the quantitative data analysis above, it is clear that the CFD model and the CALPUFF model do not precisely match in terms of data. However, this situation is understandable. The literature review indicates that model comparison studies are primarily conducted between CALPUFF, AERMOD, and other Gaussian models [6,12], all of which are Gaussian models, leading to smaller data errors and better consistency in metrics such as FB, NMSE, and R. When comparing the CALPUFF model with wind tunnel tests, field measurements, or CFD models, the fundamental principles of these methods are entirely different, and comparisons are generally made using only a subset of sampling points [24]. In this study, we resample and compare models by using a X–Y spatial range, which has made the results difficult to interpret. Reference [25] also shows concentration predictions from the CALPUFF and CFD-LS models of gas emissions from an open pit, which generally disagree with the data in all simulations.

3.5. Discussions

The above results and analysis show that the CALPUFF model and CDF model have similar characteristics when simulating gas diffusion cases in mountainous areas:

(1) It can be seen that the results of the CALPUFF and CFD models exhibit high agreement in the diffusion direction and maximum diffusion distance.The statistical indicators are not good due to the different principles of wind field calculation. There are differences in the wind field simulation methods between diagnostic modeling (CALMET) and numerical modeling (CFD). The reasons for the poor statistical data are mainly attributed to the complex structure of the flow under terrain condition.

(2) Many researchers indicate that CALPUFF modeling with calm or very light wind conditions is often associated with high impacts and, for a given application, may result in the overall highest impacts. However, in our experiments, the CALPUFF model results are not significantly greater than the CFD model results in terms of diffusion distance and direction.

(3) Atmosphere numerical models can be categorized as regional or mesoscale (<200 km), urban scale (5–50 km), microscale (<2 km), building scale (<100 m), indoor or component scale (<10 m) and human scale (<l m). Our study focuses on the simulation of microscale toxic gas dispersion, with a spatial range of 100 m to 2000 m, and a time duration within 30–120 min. The resolution of the mesoscale meteorological data is approximately 2 to 200 km with a temporal resolution 0.5 days. Therefore, using mesoscale data has few impacts on the local dispersion model in this study due to the different scale.

(4) According to our research, the evacuation ranges can be designated with 1.5 km for the N and W wind direction, 2.0 km for the E direction and 2.5 km for the S direction. In environmental assessment and emergency management, calculating the worst-case scenario range of IDLH and ${\mathrm{LC}}_{50}$ for a leak point using four wind directions or eight wind directions is meaningful. Regardless of whether it is the CALPUFF or CFD model, the simulation results show that the gases diffused spread through the valley, where there are numerous human activity hubs such as schools, resident and shops. Therefore, a specific emergency response plan needs to be formulated for this particular scenario. We choose the maximum value from the two calculation results as the reference for the emergency evacuation range.

(5) Compared to the CALPUFF model, which has minute-level computational efficiency, the CFD simulation method has heavier computational demands to produce application results. The CALPUFF model is able to reflect the influence of complex topography on airflow and pollutant distribution. The CALPUFF modeling for complex terrain for emergency scenarios can be used as a dispersion model for high-sulfur gas fields and to perform gas dispersion simulation for emergency situations.

4. Conclusions

Many of China’s high hydrogen sulfide gas reservoirs are located near mountainous and densely populated areas. The research objective is to simulate a hypothetical pipeline breakout in mountainous environments using two different models and verify the pattern of toxic gas dispersion. The research results can assist QHSE professionals with understanding a general dispersion pattern under mountainous regions and developing more effective emergency and evacuation plans and emergency management systems.

In this study, a comparative study between CALPUFF and CFD methods of toxic high-sulfur gas dispersion simulation was carried out. The experimental results show that the ${\mathrm{LC}}_{50}$ and IDLH ranges in the simulation results of the emergency-oriented CALPUFF model and the CFD model are highly consistent. The overall dispersion range deviation between the CFD model calculation results and the CALPUFF calculation results is not more than 266 m. The comparative experiments conducted in this study suggest the utility of the CALPUFF model in the emergency response platform instead of the use of CFD models, which have incompatible computational demands with the rapid response required in emergency events. Based on this study, the Emergency Response/Management System (ERMS) for the high-sulfur gas field can be developed. Further analysis would be required to verify the computation of the wind field by the two models and to provide more consistent results for complex topography scenarios.