At Southern Colorado, the Lower Arkansas River (LARV) experienced a high amount of irrigation. Heavily irrigation increase water table led to salinization of soils, waterlogging of crops, and non-valuable consumptive use of water (Scanlon et al, 2007; Pimentel et al, 1997; Gates et al, 2012). Water quality one of the concerns that threaten the LARV at the Down Stream Study Region (DSR). The groundwater interacts with the geology which has oxidation-reduction (redox) processes cause trace element to be expelled (Shultz et al, 2017). This process between the groundwater and the geology release selenium (Se) and it occurs in the bedrock as seleno-pyrite (FeSe2) and leach to the surface water system. This oxidation process accelerates the presence of dissolved oxygen (DO) and Nitrate (NO3) in groundwater. Se, NO3, Sulfate (SO4), are highly affecting the hydrological system and its increase the level of contamination in both stream and aquifer. The aim of this paper is to develop a computational model that relate and link the surface water with the groundwater, and calibrate it manually and automated with the observed data until reach the best fits between the model and the observed data so that it can be useful at the Lower Arkansas River Valley.
* Corresponding author. Tel.: +1 970 581 1099
E- mail address: [email protected] (I.A. Qurban).
1. Study Region:
The Lower Arkansas River Valley (LARV) is located in southern Colorado between Pueblo at Colorado and the border of Colorado – Kansas State line. For this paper, studying and modeling was conduct for the DSR. The DSR located in downstream of John Martin Reservoir and around 30 km east of the eastern edge of the upstream reservoir as shown in the figure 1 (Tummalapenta, 2015). DSR bounded city of Lamar from the east to the eastward of the border of Colorado – Kansas State line. The area of the DSR around 55,200 ha (552 km2), and about 60 % of the area are irrigated, which is approximately 33,000 ha (330 km2), from alluvial pumping wells or from diverted water canals from the river in DSR. Around 92 groundwater observation wells for selenium and 86 observation wells for nitrate.
Figure 1 DSR area, cities, main canals, pumping wells, and main streams (Tummalapenta, 2015)
Selenium (Se) covered numerous of inorganic and biochemical processes on earth widely. Naturally, selenium and sulfur similarly conduct and follows much of that element’s geochemistry (Severson et al, 1990, USGS). SeO4 can be released from shale formation that consists seleno-pyrite through reduction of O2 or NO3 (Bailey et al, 2013a), the following equations illustrate SeO4:
2FeSe2 + 7O2 + 2H2O 2Fe2+ + 4SeO4-2 + 4H+ 1
5FeSe2 + 14NO3- + 4H+ 5Fe2+ + 10SeO4-2 + 7N2 + 2H2O 2
The residuum Se in booth bedrock or outcropped is oxide by O2 that leaching groundwater and soil water, as well as, the NO3-laden water that from applied fertilization over the cultivated irrigation area (Bailey et al, 2013a).
Figure 2 Root zones process and chemical reactions transformation (Bailey et al,2014)
Figure 2 from (Bailey et al, 2014) illustrate the soil organic matter that composed of humus (slow-decomposing), litter (fast-decomposing), and manure. In this paper, NH4-N, NO3-N, SeO4-Se, SeO3-Se will be replaced by NH4, NO3, SeO4, and SeO3 for simplicity.
2.1.Groundwater Flow Model:
Aquifer parameter, type of the aquifer, and source or sink must be known in order to solve for groundwater. MODFLOW became mostly used groundwater flow. MODFLOW is a model that solves for three-dimensional finite-difference groundwater flow equation which is provided by the United States Geological Survey (USGS). The DSR steady region has been studied as an unconfined aquifer, heterogeneous, antiseptic, and transient flow which applied a partial 3D flow equation:
Where, x, y, and z are coordinates axis (L); Kx, Ky, Kz are Hydraulics conductivity (L/T); h is hydraulic head (L); Q is flow rate that represent source/sink (L3/T). If Q have negative value that mean the flow is out from the system. For Q have positive value that mean the flow is entering the system; Ss is Specific storage of the aquifer (L-1); and t is time (T). MODFLOW is linked with the groundwater flow equation to solve for head and discharge of each cell (i,j,k) in the 3D model domain. In this project, MODFLOW-NWT has been used that linked with unsaturated zone flow (UZF) and surface water routing (SFR). UZF package takes into account the infiltration and evaporation (ET). SFR take into account the interaction between the stream and groundwater to solve the parameter of the aquifer with differencing the hydraulic head between stream and aquifer. These packages are important in the study region because of its link with the reactive transport model RT3D.
2.2.MODFLOW-NWT Model for DSR:
The study area was divided into 22,134 cells, with a total of 18,600 active cell and 3,534 inactive cells (outside the domain). The grid cell was identified by 102 rows and 217 columns, and the dimension of the grid cell is 250 m x 250 m. In addition, MODFLOW simulates 261 stress periods, a stress period is a time step where the cells stress by source/sink represented by weeks. The observation data that available in this model from a period of 2002 to 2007. The model simulated for 40 years known as spin up instead of 5 years in order to reach the steady state fluctuation in species concentration. The number of stress periods has been increased from 261 weeks to 2085 weeks in order to match with the 40 years spin up. Thus, from MODFLOW two main files have been obtained from the output files which are heads of each grid cell in each layer (hff) and stream discharge (fort.11). These two files are the new input file fort UZF-RT3D model.
2.3. Reactive Transport Model (UZF-RT3D):
The Se and N reactive transport model (Bailey et al, 2013a) has been constructed by combining Se and N reaction modules into the base reactive transport model UZF-RT3D (Bailey et al, 2013b), that connected with MODFLOW-NWT at (UZF1) package (Niswonger et al,2011). Similarly, with MODFLOW the initial input files have been changed and the model spins up for 40 years which is equal to 2085 weeks. So, there are five main input files has been modified which are: 1. Basic transport package 2. Source-sink mixing (ssm). 3. Irrigation (irg). 4. Agriculture (agr). 5. Surface water transport (swt).
Increased the number of stress period (NPER) from 252 to 2085 weeks. The number of days has been changed as well from 1764 to 14595 days to match the stress period. Finally, loop the number of stress period length and the number of flow time steps 7 times to match the simulation spin up which is 40 years.
The point of sources/sinks for each stress period has been extended from 252 to 261 weeks by match and repeat the same data in that file of weeks 201, 202, 203, 204, 205, 206, 207, 209 at the year 2006. This step has been done in order to match the observed data which is from Jan 2003 to December 2007. Finally, the source/sinks have been looped 7 times until reach the stress period number 2085 weeks.
The initial list of irrigation canals has been changed such as the NPER and the sampling data from sampling event from 29 to 240. This change has been done after discussing with the committee Dr. Bailey and Professor. Gates, so, based on some assumption the number of sampling data has been changed from 29 sample dates to 240 sample dates. This data has been taken from the observation filed at the from the period 2003 to 2007 (i.e. day 115, 149, etc.) so, in order to match the data with the 40 years spin up the assumption dates has been agreed with the committee to run the model. Finally, the infiltrated portions of each grid cell have been looped 7 times widely in order to achieve the 40 years spin up.
The number of years of simulation has been increased from 5 years to 40 years so that the number increased until 2042. The Cell crop information has been increased first from 252 stress periods to 261 stress periods similarly with section 2.3.2, after that the data has been looped for 40 years to each week number 2085. The temperature data for the year before the beginning of simulation (2002) has been looped 7 times as well as the section 2.3.3.
The number of observation times has been changed from 252 to 2085 weeks. The time of general concentration output has been changed, it is in hours, so the total number of hours has been increased to 350280 hours by adding 168 hours which is equal one week. The total timing information was 42336 hours and it has been changed as well to be 350280 hours. The number of stress period at the upstream boundary condition has changed similarly with section 2.3.2. and then complete the missing data at that file to complete the total hours of the simulation. The daily temperature and daylight hours, and the hourly solar radiation have been increased similarly with section 2.3.2. Finally, the algae concentrations to be included with groundwater mass discharge to streams has been increased as the section 2.3.2, then looped until the day 14595.
After all, modifications that have been done on the initial input files three main files has obtained from MOFLOW which are (frt.11, hff., and sfr.) in order to run the UZF-RT3D for 40 years simulation. So that, the model has run and the output data will discuss the result.
2.4. UZF-RT3D Calibration and testing:
is important to ensure that the model and the observation data are matched and fit adequately to each other, (Bailey et al, 2014) calibration and adjustment both the stream and groundwater by selecting the groundwater reactive transport parameter. Two main processes of calibration the model which are manual and automated calibration using the Parameter Estimation (PEST) software package (Doherty, 2016). In this paper, the manual calibration has been started. The observation and file data of Se and NO3 at the LARV are available from the period of 2003 – 2007 so that this period is designed for the testing period.
2.4.1. Stream Calibration:
Based at the discussion with the committee the LARV is divided into streams and tributary as a grid cells, the Arkansas river stream is divided into five streams which are, stream1, stream4, stream5, stream7, and stream9 for calibration purpose and four main tributaries affecting the area which are, Clay Creek, Big Sandy Creek, Buffalo Creek, and Wild Horse Creek as shown in the figure (3). After that, each stream and tributary will be calibrated with the observed data separately.
Figure 3 Number of the streams, and the main tributaries.
2.4.2. Groundwater Calibration:
Groundwater calibration is similar to the one that described by Bailey et al, (2014) which is divide the area into segments that include groundwater concentration in the observation wells. In this paper, the area was divided into 15 segments based on the geology of the area, the main canals Amity, Buffalo, Fort Lyon, Lamar, the sub-basins, stream segments, and observation wells location (figure 4). So, by computing the average concentration values of Se and NO3 of each observation well and applied in the map and assigned them in each segment the calibration and the comparison between the model and the observed data can be noticed.
Figure 4 (a): Se monitoring wells location, main stream, and tributaries. (b): NO3 monitoring wells location, main stream, and tributaries. (c): Divide the study area into 15 segments, location of the average concentration of each monitoring well, and main stream with the tributaries.
The output files from the RT3D model has been used for the calibration. The concentration of each species in the stream and the concentration each species in the groundwater are the main output files from RT3D. Both stream and tributaries have been compared with the observed data. Thus, the observation data that is available in this project is from 2003 – 2007, and in order to calibrated with the model this paper include the last five years of the spin up simulation which is from 2037 – 2042 which is the data that the model is wormed up and reach the steady state condition, so it can be compared. The results of stream 5 and Big Sandy creek as tributaries for both Se and NO3 are shown in the figure (5, 6). Initially, the model results did not match with the observation data and this is because there is no chemical reaction at this area due to the riparian zones, which indicate that there is no source of energy for the bacteria to carry autotrophic reactions, so that by modifying the RT3D code and by adjusting the concentration values of each species the results might get better as shown in figures (5b, 5c, 6b, and 6c).
Figure 5 (a): Total Se concentration between the model and the observation data at Big Sandy Creek (the model is not react and has no data). (b): Total Se concentration between the model and the observation data at Big Sandy Creek (when the concentration has been changed in RT3D model). (c): Total Se concentration between the model and the observation data at Stream 5 (the model is not react and has no data). (d): Total Se concentration between the model and the observation data at stream 5 (when the concentration has been changed in RT3D model).
Figure 6 (a): Total NO3 concentration between the model and the observation data at Big Sandy Creek (the model is not react and has no data). (b): Total NO3concentration between the model and the observation data at Big Sandy Creek (when the concentration has been changed in RT3D model). (c): Total NO3concentration between the model and the observation data at Stream 5 (the model is not react and has no data). (d): Total NO3concentration between the model and the observation data at stream 5 (when the concentration has been changed in RT3D model).
On the other hand, the groundwater calibration seems it’s got a better result than the surface water. The model is working great and the results of the all 15 segments fit with each other for both Se and NO3. The results are shown in the figure (7).
Figure 7 (a) Total Se concentration in groundwater compared with the observed data. (b) Total NO3 concentration in groundwater compared with the observed data.
The results obtained in section 3 are not the final results since the manual calibration between the model and the observed data did not fit well. It can be seen from the figure (7) there are kinks at segment 1 for both Se and NO3 with a concentration of 3760 µg/L 685 mg/L respectively, and this is because the screen of the observation well contact directly with the bedrock that the Se and NO3 are reacted and release high concentration. Figures (5 and 6) illustrates that the model is not working well so the code has been changed to solve the chemical reaction at the riparian zones which is affecting the stream, if there is no reaction at the riparian zone that mean the stream will not have any concentration of Se or NO3. In addition, the reaction file at the RT3D which maintain the concentration rate of each species has been changed by using the trial and error method until reach the best match between the model and the observation data. Once manual calibration will be done, the next step would be the automated calibration using parameter estimation model (PEST).
In conclusion, the effect of Se and NO3 in the Lower Arkansas River Valley is noticeable. The excessive irrigation in this area through the years put the LARV under the threat of water quality for the farmers and the crops. So that the RT3D model is has been constructed and devolved for some species. Also, the model has been started to calibrate manually to get best matches between the observation data and the mode. In the end, this project is an ongoing research and the result that discussed in section 5 need to be improved by changing the reaction rates and then run PEST. For future work, after the automated calibration (PEST) the Best Management Practice BMP should have applied in the area in order to decrees the contamination of Se and NO3