Phosphorus Migration in an Unconfined Aquifer Using Modflow and Mt3dms

Abstract The rapid rate of urbanization and increasing demand for water in agriculture and industry are the reasons for considering groundwater as a main source of water. This can be a prologue to contamination of groundwater. Phosphorus as a type of nutrient that is derived from fertilizers has adverse effect on surface and subsurface water. The aim of this study was to monitor groundwater quality and the fate of contamination via three dimensional finite-different groundwater flow simulation (i.e. Visual MODFLOW version 4.2.). The study area was the campus of University Putra Malaysia. The monitoring indicated that the concentration of phosphorus is higher than those imposed by the standard of the Malaysian Department of Environment (DOE). Results of contamination transport modelling revealed the different rates of phosphorus transport in layers at the end of simulation period.


Introduction
Overexploitation of groundwater, which results in contamination, is threatening the future of groundwater resources. Pollution of groundwater resource with harmful and poisonous substances endangers groundwater resource and human health. Recent studies have indicated that, human activities can alter the natural composition of groundwater through the disposal or distribution of chemicals and microbial matter over land surface and into the soil, or through injection of wastes directly into groundwater. According to Belanger and Mikutel (1985); Shaw et al. (1990); and Simmons and Lyons (1994) groundwater has been documented as an important source of transportation of nutrients (like phosphorus and nitrate) to different points and places. The transported phosphorus can lead to adverse effect on sustainability of groundwater as a source of potable water.
Contamination by phosphorus and phosphate is a common problem of groundwater around the world (Stollenwerk 1996). Knapp et al. (2002) reported that the range of dissolved inorganic phosphorus was found to be usually between 0.01 to 0.1 mg/l and mostly less than 0.2 mg/l. However human activities can increase the concentration of phosphorus in water. Long-term fate of phosphorus in aquifer was reported by few studies around the world (Daroub et al. 2000;Mercik et al. 2000).
In Malaysia, groundwater is recognized as one of the important sources of water. However widespread use of fertilizers and disposal of sewage threatens the quality of this resource. Contamination of groundwater has a negative impact on the ecosystem and human health (Ismail 2001). It is noted that accumulation and transportation of some nutrients such as phosphorus is recently increased. Excessive use of phosphorus for farming activities has adverse effect on quality of groundwater, which leads to various problems. Major problems caused by phosphorus include: i) eutrophication, which could kill fish and aquatic plants living in ponds, lakes, etc. ii) undesirable odour and high turbidity due to elimination of aquatic plants in water body (Vollenweider 1968;Khan, Ansari 2005).
In 2004, the Interim National Water Quality Standard (INWQSM) for Malaysia defined the maximum value of phosphorus in groundwater for Class IIA/IIB amounting to 0.1 mg/l (DOE 2004).
This paper describes a simulation of groundwater movement and contamination migration in an unconfined aquifer within an area of approximately 0.16 sq. km located around the Faculty of Engineering, University Putra Malaysia, Serdang, Selangor, Malaysia. The study area was used to be a paddy field where fertilizers were widely used. The regional groundwater flow within the aquifer is south-east to north-west. The regional groundwater flow is suspected to carry the pollutants into the study area and also downstream of the study area.
Groundwater flow and contaminant transport models are supplementary tools in the assessment and decision making process of groundwater pollution (Jakimaviciute-Maseliene et al. 2006). Groundwater field investigation projects -which focus on the flow properties and water quality -are time and money consuming tasks. Therefore, simulation by a numerical modelling software package, such as Visual MODFLOW, is a proper choice to determine the groundwater movement and contaminant transport rate.
The main objective of this study was to model the migration path and the concentration of phosphorus in the sludge, pond, and groundwater, as well as interaction between them in the study area by Visual MODFLOW (version 4.2.) for 10 and 50 years from now. It was based on the available data from the monitoring and the assessment of groundwater movement direction and the soil characteristics of the study area.
The 160 000 m 2 study area is located in a flat region with elevations ranging between 38 m and 39.5 m above the mean sea level. Grass, low bush and trees cover the land. A pond, which is located approximately in the centre of Fig. 1, and a sludge area, approximately 200 m to north-east of the pond, are major sources hydrologic elements in the area.
The climate of the area is the typical hot-wet equatorial that is characterized by high daily temperatures 33.01ْ and 23.18 ْ°C, humidity of 83 to 100 percent, and high rainfall. (During the last 10 years, the average maximum and minimum annual rainfall in the study area amounted to 3667.6 mm and 2137.9 mm per year).
The top part of the ground in the research area is covered by Kajang formation which include unconsolidated deposits generally underlain by limestone bedrock (Gobbett, Hutchiston 1973;Huat et al. 2004). Kuala Lumpur limestone belongs to the Upper Silurian and named the Central Selangor Limestone by Gobbett and Hutchiston (1973). Fig. 2 displays the fan diagram of the study area. The top layer of the aquifer is classified as clay loam with varying thickness (between 3 and 14.5 m) and the lower layer consists of limestone as bedrock. The total thickness of aquifer system under assessment amounted to 20 meters from ground surface including layers.

Field Measurement
The dominant form of phosphorus (orthophosphorus) was measured within the study area using the ascorbic acid method by Hach (1988). Orthophosphorus was measured at 7 locations -6 wells and 1 pond -to find concentration of this contamination in different parts of the region. The essential parameters used in the software for the simulation includes hydraulic conductivity, Porosity, dispersivity, thickness and type of soil, evapotranspiration, and recharge. The data were gathered between January 2007 and March 2008 (Saghravani 2009).

Groundwater Flow Modelling
The numerical model of the groundwater was created using Visual MODFLOW software package (version 4.2.) based on the available data to simulate the groundwater flow properties such as direction, loadings, dynamics, fate, and transport of phosphorus in the groundwater. The following equation describes the mathematical representation of three dimensional groundwater flow equation (Harbaugh 2005): where K xx , K yy , K zz are values of hydraulic conductivity along the x, y, and z directions that are assumed to be parallel to the major axes of hydraulic conductivity; h is piezometric head; W is volumetric flux per unit volume that is represented for pumping, recharge, or other sources and sinks; S s is specific storage coefficient; x, y, z are coordinate direction, and t is time.
The rate and transport of contaminations in a 3-D groundwater flow systems are described as follows (Zheng 2006): Where, ω is subsurface porosity; C k is concentration of parameter in dissolved form (M L -3 ); k is Parameter; t is time (T); X i,j is distance (L); D ij is diffusion and dispersion coefficient (L 2 T -1 ); v i is seepage (L T -1 ); q s is Volumetric flow rate per unit volume of aquifer (T -1 ); C k s is concentration of source for k (M L -3 ) and n R ∑ is the chemical reaction (M L -3 T -1 ). The MT3DMS transport model follows the same spatial discretization convention as MODFLOW.

Model Boundaries
Some boundary conditions that were assigned include, i) Constant Head Boundary (CHB), which was placed in the north of the study area and assumed to be constant throughout the simulation, ii) Horizontal Flow Boundary (HFB) or Wall Boundary that is used to illustrate any obstacle in path of groundwater, iii) recharge boundary and evapotranspiration boundary, which were used to create the model.

Model Calibrations
The study area was modelled using 45 columns and 30 rows in the software package and each column and row of the model represents 13-m wide strip of the aquifer within the study area. The model was initially calibrated by adjusting the input data. The result of normalized Root Mean Square for concentration of contaminant was 28.6%, root mean squared was 0.14 mg/l, residual mean was 0.059 mg/l, and the standard error of the estimate was 0.057 mg/l.

Model Limitations
Despite the high accuracy and capability of Visual MODFLOW to simulate the groundwater flow and the prediction of contamination fate, some limitations lead to decrease in the precision of the simulation. One of the important factors that directly affected the result was the availability of accurate input data because precision of contamination prediction would be reduced whenever the data is determined by average and/or estimation. Evapotranspiration and recharge were assumed constant whereas in the actual fact, these amounts are subjected to changes during months and years.

Simulation of Groundwater Flow and Phosphorus Migration
Research area was assumed as heterogeneous and isotropic unconfined aquifer. Pollution transport simulation was applied to the Visual MODFLOW. Simulations were performed for two periods of 3650 and 18250 days and implemented in steady state type due to constant hydraulic head. Long 50-year period was selected due to i) low hydraulic conductivity of Layer No. 1, and because ii) phosphorus is adsorbed by and spread in soil very slowly, and in short periods of time, however, once phosphorus adsorption capacity reaches the maximum amount the phosphorus starts to flow in groundwater (Knapp et al. 2002). These two steady state scenarios were applied to assess migration of contamination in groundwater. The result of the simulations is as follows:

Case 1: The Ten-Year Simulation Period
Groundwater flow direction was evaluated through similarity between interpolation method in Fig. 4

and the direction of arrows in simulation results in Figs 3 and 4.
Contamination migrations through the layers are different. It should be mentioned that migration of phosphorus within 10 years is not substantial and phosphorus cannot move long distances from the sludge as illustrated in Fig. 3. It seems that a longer period of time is needed for the migration of contamination to further places. In other parts of the study area, like the pond, contamination has not changed. However, in the second layer, the pollution migrates faster when compared to Layer No. 1. In case of the present study, the source of contamination is located in the first layer. The contaminant has to flow through medium in Layer No. 1 before it can reach Layer No. 2. It should be noted that X and Y axis in figures are expressed in metres.  After the simulation period has reached 10 years, it was observed that the concentration of the phosphorus within the sludge in Layer No. 1 was 0.6 mg/l, whereas the concentration of the phosphorus in Layer No. 2 below the sludge was still 0.25 mg/l. Within Layer No. 2, the contamination is spreading faster than in Layer No. 1 because of the higher hydraulic conductivity of this layer. After the simulation period had reached 10 years, the minimum phosphorus concentration within Layer No. 1 was 0.2 mg/l. This value of phosphorus concentration is at the upper limit of Class IIA/IIB of INWQSM (2004). This means that the groundwater within Layer No. 1 is approaching to the upper limit of the permitted phosphorus concentration in water used for drinking. It should be noticed that vertical migration is responsible for the occurrence of phosphorus in the study area and also its migration to further places during the simulation period. Fig. 7 clearly shows movement of contamination in two layers and vertical migration.

Conclusions
This research was focused on the future of groundwater pollution by phosphorus in the area around the Faculty of Engineering, University Putra Malaysia. The study consists of numerical simulation of the groundwater flow and phosphorus dispersion throughout the study area using Visual MODFLOW at two time interval, namely, ten and fifty years. Determination of contamination by phosphorus in different parts of groundwater and demonstration of pollution transport in aquifer were the key points of this research.
Simulation of phosphorus dissipation within the first layer showed that it cannot migrate widely from the sludge due to soil characteristics of Layer No. 1. In other words, the spread of contaminants from the sludge in Layer No. 1 is restricted. On the other hand, contamination migration by groundwater in Layer No. 2 was different and the process of spreading of phosphorus via this layer is much faster.
To study this aspect further, sources of contamination, which have effects on groundwater quality, sludge, and pond in wide range should be examined to protect the source of water in the research area and other reservoirs situated downstream.