Onshore Spill Simulator

1. Purpose

The plugin simulates what happens to oil or fuel released on land during an onshore drilling campaign: where it flows, how far and how fast it travels, how much evaporates, how much soaks into the ground, whether it reaches the drainage network, and which communities, water points and facilities it reaches and when. Inputs are supplied as QGIS layers and parameters. Applications at an onshore site require a suitable DEM and independently checked release, surface, weather and receptor inputs.

Version 2 adds the features of temperature-dependent overland oil models such as OilFlow2D (heat transfer between the oil, the sun, the air and the ground; temperature-dependent density, viscosity and yield stress; roughness zones; barriers; recovery points and free-outflow cells; several sources), a hydrogen sulphide (H2S) dispersion model for sour-gas leakages, and a movable tanker that can be placed anywhere on the map and run with one click.

The results are intended for emergency response planning: sizing spill-response equipment and response times, placing containment, siting water-supply boreholes and camp facilities, briefing communities, and preparing the oil spill contingency section of the ERP and the ESIA. They are a screening simulation, not a forecast of a particular event.

2. How the model is built

The model was developed in nine steps. Each step corresponds to a module of the plugin (engine/ folder), so the method can be checked and changed step by step.

3. Process equations

3.1 Overland flow

Oil moves from a cell to its lower neighbours according to the slope of the oil free surface η = z + h (ground plus oil), so oil fills hollows and spills over their lowest side. The outflow is shared among the downhill neighbours in proportion to (Δη/d)^1.1 (multiple flow direction, Freeman, 1991). Only the depth above the surface retention depth hr can move. Its speed is the lesser of the viscous thin-film (gravity-current) speed (Huppert, 1982) and the Manning speed:

where ρ and μ are the oil density and viscosity, S the free-surface slope and n the Manning roughness (Engman, 1986). On sloping natural ground the oil runs in narrow threads, so only the rivulet width (default 6 m of a 30 m cell) holds retained oil and soaks it up; on flat ground and on pads the whole cell is wetted.

3.2 Evaporation

Fingas (2004) generic equations (t in minutes, T in °C, %D = mass % distilled at 180 °C):

The equations are applied cell by cell with the pseudo-time method, so parcels of different age that mix are weighted by mass. Evaporation stops at the product's evaporable limit. Viscosity rises with the evaporated fraction F_e: μ = μ0 exp(C Fe) (Mackay et al., 1980), C = 10 for crude oils and 1 for distillates. In heat-transfer mode the evaporation coefficient is evaluated with the local oil temperature, so warm oil evaporates faster (Section 3.7).

3.3 Infiltration and groundwater

Green-Ampt infiltration (Green and Ampt, 1911) with parameters for oil: the hydraulic conductivity is scaled by (ρo/ρw)(μw/μo) and reduced by the water already in the pores; the capillary suction is scaled by the ratio of surface tensions (about 0.4, Leverett scaling). The depth reached by the oil front is the infiltrated volume divided by the oil the soil can retain, Sor φ (default 0.10 x porosity, i.e. about 40 L/m³ for silty soils, consistent with CONCAWE 3/79). Where this depth reaches the depth to groundwater the cell is flagged.

3.4 Well spray ("oil rain")

Droplets (log-normal, median 250 µm) are released at an effective height H (default 40 m: derrick height plus jet rise), carried at the wind speed at that height (power law, exponent 0.14) and settle at their terminal velocity (Atlas et al., 1973 for drops above 0.1 mm, Stokes law below). Each size lands at x = u(H) H / vt downwind, spread by the wind meander and by the range of release heights. The deposit is recomputed every hour as the wind changes.

3.5 Rain wash-off

Storm runoff is calculated with the SCS curve number (USDA-NRCS, 1986). While runoff lasts, half of the retention depth is flushed and the oil is carried at half of the kinematic sheet-flow speed of the runoff, which grows with the upstream catchment area. Stream discharge also rises with the storm.

3.6 Streams

In a flowing stream the oil floats and moves down the D8 path at the surface speed of the water, u = 1.1 × 0.7 Q0.1 (hydraulic geometry, Leopold and Maddock, 1953), where Q is the stream discharge. Oil is lost to the banks at an exponential rate (default 15% per km) and stays there as stranded oil, where it continues to weather. Oil that reaches the edge of the model area is counted as having left it.

3.7 Heat transfer and temperature-dependent rheology (OilFlow2D-type)

With Heat-transfer mode ticked, every cell carries an oil temperature T as well as a thickness. Heat moves with the oil (the heat content h T is transported with the volume), and each time step the temperature of the layer is advanced with the energy balance used by temperature-dependent overland oil models (Garcia et al., 2023; Echeverribar et al., 2023):

• S - global solar radiation on a horizontal surface: clear-sky value (Haurwitz, 1945) from the sun's elevation for the site latitude, date and hour, reduced for cloud and atmospheric haze (Kasten and Czeplak, 1980). Albedo of oil a = 0.05.
• Long-wave exchange - oil emissivity ε = 0.95; sky emissivity ε_a from air temperature and humidity (Brutsaert, 1975), raised for cloud.
• Convection - hc = 5.7 + 3.8 U W/m²K with the hourly wind U (Watmuff et al., 1977); strong wind cools the oil and shortens its runout.
• Ground conduction - h_g = 10 W/m²K to a ground surface at the mean air temperature plus 2 °C.
• Air temperature follows a daily cycle between the monthly mean minimum and maximum (minimum near 03:00, maximum near 15:00); relative humidity and cloud follow the month.
• Source temperature - the oil leaves the wellhead or tank at the given temperature (e.g. 60 °C for produced fluid); spray droplets reach air temperature before they land.

The balance is linearised about the current temperature and integrated exactly over the step, so it is stable for any thickness. The oil properties then follow the temperature:

The oil moves as a Bingham film: where the driving stress τb = ρ g hm S is below the yield stress it does not move; above it the laminar speed is reduced by 1 - 1.5 r + 0.5 r3, r = τy/τb. A waxy crude therefore spreads while hot, gels as it cools at night, and moves again when the afternoon sun warms it - the behaviour shown by OilFlow2D for high pour point oils. In flowing streams gelled oil is carried as floating lumps. With heat transfer off, the properties are fixed at the mean air temperature (constant-property mode, as in the basic OilFlow2D module).

3.8 Control features

• Roughness zones - polygons whose Manning n replaces the surface-class value (a field manning_n or one value for the layer): crops, dense grass, gravel, paving.
• Barriers - lines or polygons (earth berms, bund walls, sandbag lines, road embankments) raised by the given height on the flow surface. Oil ponds behind them and overtops when it is deeper than the barrier.
• Recovery points and free-outflow cells - points, lines or polygons where oil is removed (sump, interceptor trench, skimming point, culvert outlet); the volume is reported as recovered in the mass balance.
• Several sources - additional point releases with their own volume, start time, duration and temperature (configuration key extra_sources).

3.9 H2S dispersion from a sour-gas leakage

The gas phase of a well spillage is modelled separately with a Gaussian puff model (engine/h2s.py).

• Source - gas rate G (MMscf/d at 60 °F) and H2S content x (ppm by volume) give the H2S mass rate q = G 106 0.0283 / 86400 x P/(R Ts) x 34.08 g/s. Two release modes: a free vertical jet from the rig floor (plume rise from momentum and buoyancy, Briggs, 1975, 1984) or an impinged leak, where the jet is stopped by the substructure or BOP stack and the gas spreads from a low height with no plume rise (the conservative case, used for the reference runs). A mixture heavier than air is flagged, because a passive Gaussian model is then not conservative near the source.
• Transport - puffs are released every 90 s and carried by the hourly wind at plume height (power-law profile), turning smoothly as the wind veers; each puff is spread along the path to the next puff so the plume is continuous.
• Spread - Pasquill-Gifford stability class from wind speed, sun elevation and cloud (Turner, 1964); Briggs (1973) open-country dispersion curves; reflection at the ground and at the top of the mixed layer (150 m on stable nights to 1,500 m on unstable afternoons); the plume height is reduced by half the rise of the ground above the wellhead.
• Concentration - ground-level 10-minute means, converted to ppm at 25 °C; 60-minute means are running hourly averages.
• Thresholds - EPA Acute Exposure Guideline Levels for H2S (60 min): AEGL-1 0.51 ppm (odour, discomfort), AEGL-2 27 ppm (irreversible or escape-impairing effects), AEGL-3 50 ppm (life-threatening); NIOSH IDLH 100 ppm; 500 ppm (rapid loss of consciousness). AEGL zones use the highest 60-minute mean; IDLH and 500 ppm zones use the highest 10-minute mean. The 10-minute AEGLs (0.75, 41 and 76 ppm) are higher, so the 60-minute values are the stricter test.
• Outputs - H2S time series and worst 10-min / 60-min rasters, time at which AEGL-2 is first exceeded, threshold zones with their furthest reach, the planned H2S zone for comparison, receptor exposure (worst class, first time above AEGL-1/2/3 and IDLH), HTML player and GIF/MP4.

Not modelled: ignition of the release (an ignited leak burns H2S to SO2), dense-gas slumping, building and rig-structure wakes, and calm winds below 1 m/s.

4. Default parameters

4.1 Products (oils.py)

Temperature dependence (all products): Andrade constant B 1,500-6,500 K, density change -0.65 kg/m³ per °C, specific heat 2,000 J/kgK, yield stress 1-5 Pa at the pour point rising tenfold every 4-5 °C of cooling.

These are generic product-class defaults, not measurements of a particular reservoir fluid. Replace the values with suitable product assay data when available.

4.2 Surface classes (model.py)

Soil suction values follow Rawls et al. (1983) for sandy loam to clay loam; retention depths are indicative holding capacities of 2-10 L/m² for compacted to vegetated ground.

4.3 Monthly climate defaults (climate.py)

Synthetic example only. No station observations or regional climate profile are bundled. All twelve months use invented, constant values: mean temperature 20 C, maximum 25 C, humidity 50%, wind u = 3 m/s and v = 0 m/s, cloud fraction 0.5, no rain or stream discharge, soil water saturation 0.5, and groundwater depth 5 m. These are demonstration inputs, not an environmental assessment. Supply reviewed weather, hydrology, soil and groundwater values for the study area using the wind CSV, Processing overrides and explicit storm inputs. Numerical results will differ from releases with other defaults.

5. Running the plugin

1. Install: Plugins > Manage and Install Plugins > Install from ZIP, choose onshore_spill_sim.zip. The tool appears under Plugins > Onshore Spill Simulator, on the toolbar and in the Processing Toolbox.
2. Open a project with a DEM in a projected CRS (metres). Arrange and style the layers the way the animation should look; only the layers that are switched on are drawn.
3. Run Simulate onshore spill. Set the DEM, click the spill point on the map, choose the release type, the product and the release size, the month, the start time and any storms (e.g. 24,3,50; 54,2,25).
4. Add the ground layers (pads, site plan, roads, dugouts, woodland), the mapped rivers and the receptor layers.
5. Choose an output folder and run. The log shows each stage. Runtime depends on the grid size, simulated duration, output settings and available computing resources.
6. To redraw the animation after restyling the map or changing the extent, run Render animation from a finished run and select the *_run.json file.

Advanced settings (expand Advanced parameters): spray share, release height and droplet size; model area; cell size; DEM smoothing and depression tolerance; stream threshold; snapping the spill point onto a stream (use it for road crossings); temperature, stream flow, soil wetness and depth to groundwater overrides; curve number; rivulet width; bank stranding; animation size and speed. Heat-transfer settings (source temperature, cloud, humidity) and the control layers (roughness zones, barriers, recovery points) are in the main form.

5.0 Simulation panel (set up, run and play inside QGIS)

Click Simulation panel on the toolbar (green play button) or in Plugins > Onshore Spill Simulator. The panel docks on the right of the map and works with the study project; it reads 00_Config/model_inputs.json, so a run made in the panel uses the same DEM, ground layers, rivers and receptors as the reference scenarios.

1. Set up and run. Choose the hazard (oil well spillage, H2S leakage, rig diesel tank, road tanker, continuous leak) and the location (wellhead, rig fuel storage, stream crossing, or Pick a point on the map). Set the date and time the release starts. The season box underneath updates with the date: season name, mean and afternoon temperature, humidity, monthly rain and rain days, the resultant wind and its steadiness, stream flow, soil wetness and groundwater depth for that month (synthetic example values unless overridden), with notes on what this means for oil and gas. Typical for the month fills in the storms expected over the run (mean number of rain days in the period, each at the mean rain-day depth, as a late-afternoon storm). Adjust the release figures and press Run simulation. The model runs in the background; the map stays usable.
2. Play. When the run ends the panel switches to Play. The result is added to the Animation player group at the top of the Layers panel and stepped through frame by frame with the slider, the play / pause buttons and the speed list. Any earlier run in 04_Outputs can be chosen from the list. For oil runs the oil temperature and the yield stress of gelled oil can be shown instead of the thickness. Subtitles describe each frame: the clock and date, whether the release is still running, the wind and rain or the atmospheric stability, the oil temperature and gelled volume, how far the oil or gas has reached and in which direction, the receptors reached and the nearest settlement. Make GIF / MP4 from this map draws the animation from the layers that are switched on, including an online basemap such as Google Hybrid; Open the HTML page opens the interactive page of the run.
3. Settlements. Lists the settlements around a well within a chosen radius with their distance, direction and district, and how each is affected in the run being played (oil arrival time, or the H2S hazard class and peak concentration). Rows within 1 km are shaded red and within 2 km orange. Add a settlement lets a compound or village be clicked on the map (with the Google Hybrid basemap on); it is saved in the settlements_digitised layer and used as a receptor in every later run.

5.1 Movable tanker

1. Click Place a movable tanker and simulate (orange truck on the toolbar), then click the map where the tanker loses its load - a road, a stream crossing, a laydown area.
2. A small form asks for the product, the volume lost, the time to empty, the month, the start time, any storm and whether to snap the point onto a stream. The point is added to the Movable tanker layer.
3. The model inputs (DEM, ground layers, rivers, receptors, output folder) are read from 00_Config/project_config.json beside the project file, so the Processing form opens already filled in - press Run, or tick Run straight away in the form. Results go to 04_Outputs/Tanker_Runs/<name>, are added to the project and the HTML player opens.
4. To run several positions later: move or add points in the Movable tanker layer with the editing tools (fields name, volume_m3, product, release_min, month, start, storms, snap_m), then run Simulate spills at movable tanker points; each point becomes a run.
5. Every HTML player also has a Sandbox tanker button: drag the orange marker anywhere within 6 km of the source and the page traces the downhill path, the distance to the first watercourse, an indicative travel time and the receptors within 120 m of the path. It is a quick screening aid; use the QGIS tool for volumes and timing.

5.2 H2S release

Run Simulate H2S leakage at a well: click the well, enter the gas rate and H2S content, choose the release mode, the month and start time, the receptor layers and the planned H2S zone radius to compare with.

6. Outputs

7. Using the results

• Response times: the arrival isochrones and the receptor table give the time available to deploy booms, berms and absorbents before oil reaches a stream, a dugout or a settlement.
• Containment: the footprint and the first entry into the drainage network show where interceptor trenches and earth berms should be pre-planned (pad perimeter drain, culverts on the access road).
• Water supply: water points reached by oil or by spray should not be used for drinking water for the camp; new boreholes should be sited upwind and up-gradient of the rig, outside the pad drainage.
• Season: compare dry- and wet-season runs - in the wet season storms spread the oil and flowing streams carry it kilometres in hours; in the dry season oil stays near the source and soaks deeper into dry soil.
• Mass balance: the volumes on the ground and stranded on banks are the volumes to be recovered and the basis for waste handling estimates.

8. Calibration and limitations

• The DEM controls the result. With a 30 m global DEM, ditches, culverts, bunds and the pad drainage are not resolved. A drone or RTK survey DTM of the site (1-5 m) should be used for the final assessment; reduce the depression tolerance to 0.2-0.5 m with such data.
• Release figures are placeholders until the drilling engineering basis gives the worst-case discharge, the duration to capping and the rig's tank inventories and positions.
• Surface retention, roughness, soil conductivity, stream flow and bank stranding are literature-based defaults; they should be checked against site soil tests (infiltration, texture), stream gauging and, where possible, field spill-trial observations.
• Preferential flow through cracks, root channels and fractured laterite is not modelled; actual penetration can be deeper than the computed front.
• Fire is not modelled; H2S is modelled with a passive Gaussian puff model (Section 3.9). Response actions are applied only where barriers and recovery points are given; the reference runs show the unmitigated case.
• The heat balance treats the oil layer as well mixed; for layers thicker than a few centimetres the surface cools faster than the core. Solar radiation and cloud are monthly typical values, not a forecast.
• Default gas rate and H2S content are placeholders and do not describe a particular well.
• Dissolution and biodegradation are neglected over the simulated periods (hours to days).

9. References

• Andrade, E.N. da C. (1930). The viscosity of liquids. Nature, 125, 309-310.
• Atlas, D., Srivastava, R.C. and Sekhon, R.S. (1973). Doppler radar characteristics of precipitation at vertical incidence. Reviews of Geophysics and Space Physics, 11(1), 1-35.
• Barnes, R., Lehman, C. and Mulla, D. (2014). Priority-flood: an optimal depression-filling and watershed-labeling algorithm for digital elevation models. Computers & Geosciences, 62, 117-127.
• Briggs, G.A. (1973). Diffusion estimation for small emissions. ATDL Contribution File No. 79, Atmospheric Turbulence and Diffusion Laboratory, Oak Ridge.
• Briggs, G.A. (1975). Plume rise predictions. In: Lectures on Air Pollution and Environmental Impact Analyses, American Meteorological Society, Boston, 59-111.
• Briggs, G.A. (1984). Plume rise and buoyancy effects. In: Atmospheric Science and Power Production, DOE/TIC-27601, 327-366.
• Brutsaert, W. (1975). On a derivable formula for long-wave radiation from clear skies. Water Resources Research, 11(5), 742-744.
• CONCAWE (1979). Protection of groundwater from oil pollution. Report 3/79, The Hague.
• Echeverribar, I., Fernandez-Pato, J., Martinez-Aranda, S., Garcia, R. and Garcia-Navarro, P. (2023). A model for computing thermally-driven shallow flows. IOP Conference Series: Earth and Environmental Science, 1136, 012037.
• Engman, E.T. (1986). Roughness coefficients for routing surface runoff. Journal of Irrigation and Drainage Engineering, 112(1), 39-53.
• Fingas, M.F. (2004). Modeling evaporation using models that are not boundary-layer regulated. Journal of Hazardous Materials, 107(1-2), 27-36.
• Freeman, T.G. (1991). Calculating catchment area with divergent flow based on a regular grid. Computers & Geosciences, 17(3), 413-422.
• Garcia, R., Echeverribar, I., Martinez-Aranda, S., Fernandez-Pato, J., Garcia-Navarro, P., Calero, N., Cordova, J.R. and Artigas, J. (2023). Overland oil spill modelling considering heat transfer for high pour point oils. Proceedings of the 45th AMOP Technical Seminar on Environmental Contamination and Response, Environment and Climate Change Canada, 39-50.
• Garcia, R., Lacasta, A. and Garcia-Navarro, P. (2015). Two-dimensional flexible-mesh model to simulate oil spills over complex terrain. Conference proceedings, 517-525.
• Green, W.H. and Ampt, G.A. (1911). Studies on soil physics, part 1: the flow of air and water through soils. Journal of Agricultural Science, 4(1), 1-24.
• Haurwitz, B. (1945). Insolation in relation to cloudiness and cloud density. Journal of Meteorology, 2, 154-166.
• Huppert, H.E. (1982). The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface. Journal of Fluid Mechanics, 121, 43-58.
• Kasten, F. and Czeplak, G. (1980). Solar and terrestrial radiation dependent on the amount and type of cloud. Solar Energy, 24, 177-189.
• Leopold, L.B. and Maddock, T. (1953). The hydraulic geometry of stream channels and some physiographic implications. U.S. Geological Survey Professional Paper 252.
• Leverett, M.C. (1941). Capillary behavior in porous solids. Transactions of the AIME, 142, 152-169.
• Mackay, D., Buist, I., Mascarenhas, R. and Paterson, S. (1980). Oil spill processes and models. Environment Canada, Report EE-8.
• National Research Council (2010). Acute Exposure Guideline Levels for Selected Airborne Chemicals, Volume 9 (hydrogen sulfide). National Academies Press, Washington DC; EPA AEGL values for hydrogen sulfide.
• NIOSH (1994). Immediately Dangerous to Life or Health (IDLH) values: hydrogen sulfide, 100 ppm.
• Rawls, W.J., Brakensiek, D.L. and Miller, N. (1983). Green-Ampt infiltration parameters from soils data. Journal of Hydraulic Engineering, 109(1), 62-70.
• Turner, D.B. (1964). A diffusion model for an urban area. Journal of Applied Meteorology, 3, 83-91.
• USDA-NRCS (1986). Urban hydrology for small watersheds. Technical Release 55.
• Watmuff, J.H., Charters, W.W.S. and Proctor, D. (1977). Solar and wind induced external coefficients for solar collectors. COMPLES, 2, 56.