Petroleum leakage has produced long-lived groundwater contamination in the Rey industrial area south of Tehran, where oil from refineries, storage facilities, pipelines, and other infrastructure has accumulated in and above the shallow aquifer. Oil present as a separate liquid phase on the groundwater table is called light non-aqueous phase liquid (LNAPL). A central remediation question is whether these free-phase oil lenses will spread substantially with groundwater flow or remain localized near their sources.
Agah et al. (2025) developed a field-informed two-dimensional finite element model to examine this behavior over two years. FEATool Multiphysics solved the time-dependent equations for oil-layer thickness and lateral movement. The simulations indicated limited migration: thick oil zones around the Tehran Oil Refining Company (TORC) and the seven oil companies in the Rey industrial area (REY7) contracted with time but remained localized. The authors recommended concentrating remediation on these source zones. The model tracked oil-layer redistribution; it did not simulate dissolution or biodegradation.
Illustrative composite based on the study-area map (Figure 8a), mesh (Figure 8b), and 7.3-day oil-layer result (Figure 11a) from Agah et al. (2025), published under CC BY 4.0. The lower image overlays an adapted plume and mesh on the orange study area. This redrawn image is not a georeferenced or quantitative reproduction. The aquifer scene above is illustrative.
Field data and oil-lens model
The modeled area covered 6000 by 6000 m in the Rey industrial zone. Field investigations had identified high-density contamination around TORC and REY7, while lower-concentration contamination was more widely distributed. Groundwater levels measured in monitoring wells and core borings indicated a general flow direction toward the southwest. These observations supplied the initial water-table distribution and the starting oil-layer thickness used by the numerical model.
The study represents the contaminant as an immiscible oil lens floating on the groundwater table. The formulation begins from oil-phase mass conservation and Darcy’s law, then reduces the problem to a depth-averaged equation for oil-lens thickness. In this form, plume motion combines advection with groundwater flow and gravitational spreading caused by the density difference between oil and water. The researchers assumed a homogeneous and isotropic aquifer, constant porosity, incompressible fluids, a sharp transition between saturation states, and horizontal flow consistent with the Dupuit-Forchheimer approximation.
Several additional assumptions narrow the model to lateral LNAPL migration. Air flow in the unsaturated zone is neglected, and the three-phase air-oil-water problem is simplified to two phases. Oil is treated as conservative, so dissolution, volatilization, soil reaction, and biodegradation are excluded from the governing formulation.
Finite element implementation
FEATool Multiphysics, a MATLAB FEA simulation toolbox, served as the numerical solver for the governing equations. The two-dimensional domain was discretized with finite elements, with greater mesh density in regions where oil saturation or capillary-pressure gradients were high. Time integration used an implicit stepping scheme. The model used site-specific hydraulic and fluid properties, including a soil porosity of 0.42, oil viscosity of 0.6 Pa·s, and oil and water densities of 800 and 1000 kg/m³, respectively.
Initial contaminant conditions were derived from measured maximum oil-layer thickness around TORC and REY7. The lateral boundaries were assigned zero oil-layer thickness. The paper also describes the upper and lower surfaces as impermeable, but Figure 8 shows a horizontal two-dimensional mesh and does not show how those vertical constraints were implemented numerically. Groundwater flow followed the measured site gradient. The calculation tracked changes in oil-lens thickness and lateral position rather than vertical migration through the full unsaturated and saturated subsurface.
Plume evolution over two years
The FEATool simulations were evaluated after 7.3 days, one year, and two years. After 7.3 days, the thickest parts of the plume remained centered on TORC and REY7. Oil-layer thickness exceeded 6.5 m in the thickest zones, while the 0.5 m oil-thickness contour remained close to its initial shape. The study attributes this low mobility to the combination of restrictive soil permeability and high hydrocarbon viscosity.
After one year, the area with thickness above 6.5 m had contracted substantially, but the enclosing 0.5 m oil-thickness contour changed much less. The modeled plume showed only slight displacement at its southern edge despite the groundwater-flow direction. By two years, thick oil zones became still more localized and the 0.5 m contour contracted further, yet the overall plume continued to show little lateral or southward migration.
Simulated oil-layer thickness after (a) 7.3 days and (b) one year. The color scale gives thickness in metres. Figure 11 from Agah et al. (2025), published under CC BY 4.0.
Limited lateral migration does not imply rapid cleanup. The simulations suggest that thick oil zones can persist near the original sources while groundwater flow redistributes the free-phase oil only slowly. The 0.5 m contour describes an oil-thickness threshold, not the extent of all groundwater contamination or the distribution of dissolved hydrocarbons.
Field comparison, interpretation, and limitations
The authors report agreement with LNAPL thickness measurements at wells 6 and 53, randomly selected from an oil collection plan containing 148 wells. Groundwater-level measurements also characterized the site hydraulics used by the model. The paper describes matching calculated and measured oil-layer thickness, then using the finite element model to estimate oil-layer behavior across the site. It does not establish an independent predictive test against withheld observations. No single error statistic is reported for the comparison, so the evidence is best read as a site-specific field comparison rather than a complete characterization of model uncertainty.
The model is explicitly two-dimensional and assumes negligible vertical variation. The authors note that buoyancy-driven flow, capillary retention, and heterogeneous subsurface layering could redistribute hydrocarbons vertically, and recommend three-dimensional modeling in future work. There is also a distinction between the mathematical formulation and the paper’s discussion of attenuation. The governing model treats oil as conservative and excludes dissolution and biodegradation, while the results section discusses contraction of the thick plume in terms of possible natural processes including dissolution and minor biodegradation. The simulated spatial evolution therefore supports the conclusion of slow, localized plume behavior, but the relative contribution of individual degradation mechanisms is not resolved by the FEATool formulation itself.
For the broader Rey area, the authors divide the problem into localized high-density zones and more widespread low-density contamination. They recommend concentrating active measures on TORC and REY7 and propose an additional oil collection well south of REY7 where another high-density zone may be present. Their wider remediation plan calls for free-phase oil removal with collection or pump-and-treatment systems before applying land-use-specific cleanup methods. Natural attenuation is presented as a possible option for lower-density areas after migration from the high-density sources has been controlled.
Related FEATool resources include the Flow in Porous Media tutorial, which demonstrates porous-domain flow modeling, the 1D time-dependent convection and diffusion example, and the Custom Equation physics mode for user-defined partial differential equation (PDE) formulations. These examples illustrate individual techniques relevant to the groundwater model, but do not reproduce the Rey hydrocarbon-plume simulation.
References
- A. Agah, F. Doulati Ardejani, M. Shehab, C. Butscher, R. Taherdangkoo. Modeling Hydrocarbon Plume Dynamics in Shallow Groundwater of the Rey Industrial Area, Iran: Implications for Remediation Planning, Water, 17(8), 1180, 2025, doi: 10.3390/w17081180.