Modelling of density-dependent flow systems: Sensitivity to spatial and temporal discretizations and numerical schemes using SEAWAT 2000

Abstract

Contaminant transport with significant density contrasts is of increasing interest in many subsurface-hydrology problems like seawater intrusion, saltwater up-coning in coastal aquifers, and dense contaminant migration. Several numerical codes have been developed to simulate these systems. Results from these codes are found to be sensitive to grid and time step sizes, and to the numerical schemes used to approximate the solution. Although the 3-D Finite Difference Method code, SEAWAT has been used over a wide range of problems, the sensitivity of the computed results to spatial and temporal discretization levels for different numerical schemes supported by the code has not been studied in detail. This study aims to: (1) provide an approximate guide for SEAWAT users to ensure proper selection of grid and time step sizes for a particular numerical scheme in order to minimize the numerical error (in 2-D) and (2) To investigate the ability of the SEAWAT code to simulate more complex 3-D systems. The Elder-Voss benchmark problem was selected for the study. Different levels of grid and time discretizations were found to produce significantly different results. A grid size with dx (0.38%) and dz each 0.6% of the total length and depth of the domain respectively is found to be fine enough to produce results with acceptable accuracy for most of the numerical schemes with Courant number (C) of 0.1. Some numerical schemes produce accurate results at coarser meshes compared with other schemes that produce similar accuracy at only an extremely fine mesh. For the frequently used C value of 1.0, the computed salinity-spread patterns varied significantly for different solution schemes compared with the case when C is 0.1 or 0.5. To ensure a high level of accuracy in the modeling results using SEAWAT, C should be ≤0.1 when the Peclet number is ≤ about 1. SEAWAT was able to capture the main physical features of the Elder-Voss convection pattern in 3-D at coarse mesh sizes compared with that produced by a Finite Element Method code “FEFLOW” at a very fine mesh. The 3-D modelling requires large numerical efforts compared with that of the 2-D case. Runtimes for the 2-D problem range from less than. hr for a coarser mesh to about 2.8 days in the case of a very fine mesh, whereas it takes around 6 days in the 3-D case with a relatively coarse resolution level. Similarity between 2-D and 3-D convection patterns was found to exist in terms of fingering, and up/down welling behaviours. Other factors like the level of accuracy required, computational expense and storage memory requirements should be considered along with the spatial and temporal resolution levels.