Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Hydrogeologic modeling of submarine groundwater discharge in the Gulf of Mexico near southeastern Louisiana Thompson, Craig E. 2005

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2005-0665.pdf [ 17.19MB ]
Metadata
JSON: 831-1.0052556.json
JSON-LD: 831-1.0052556-ld.json
RDF/XML (Pretty): 831-1.0052556-rdf.xml
RDF/JSON: 831-1.0052556-rdf.json
Turtle: 831-1.0052556-turtle.txt
N-Triples: 831-1.0052556-rdf-ntriples.txt
Original Record: 831-1.0052556-source.json
Full Text
831-1.0052556-fulltext.txt
Citation
831-1.0052556.ris

Full Text

HYDROGEOLOGIC MODELING OF SUBMARINE GROUNDWATER DISCHARGE IN THE GULF OF MEXICO NEAR SOUTHEASTERN LOUISIANA B Y CRAIG E. THOMPSON B.Sc, The University of British Columbia, 2002 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES Geological Sciences THE UNIVERSITY OF BRITISH COLUMBIA June 2005 © Craig E. Thompson, 2005 A B S T R A C T Submarine groundwater discharge (SGD) has been the subject of increasing amounts of research in the past decade. It is defined as the total mixture of seawater and fresh groundwater flowing out from the seabed, into the coastal water, through the underlying sediments. Since the early recognition of its ecological significance in the mid-1960's, numerous studies have been performed whose goal was its quantification. It has been studied using a variety of techniques. These include direct methods, geochemical tracers, and numerical models. A number of studies have incorporated multiple techniques with limited succes. These efforts have found that results obtained using geochemical tracers and direct methods are in general agreement; while those obtained through modelling methods vary to some degree. The goal o f this study was to quantify the rate o f offshore groundwater discharge in the coastal zone of southeastern Louisiana near the Mississippi delta. A previous study indicated that S G D may be important within the region, with the volume discharged reaching approximately 10% of the Mississippi River discharge. Numerical groundwater flow modelling and the geochemical tracer radium were used in this study to examine the interpretation of the previous results. Quantifying the volume of groundwater discharged in this region is of particular significance due to the presence of hypoxic bottom waters that develop each summer. Two cross-sectional models and a three-dimensional model were developed in order to understand the groundwater flow systems of the region. The model region extended inland from the upland recharge zone out to the edge of the continental shelf. The hydrogeologic units present in the study region are part of the Coastal Lowlands n aquifer system. The system is comprised of a thick package of unconsolidated sediments that extend kilometres into the subsurface. The models were simulated using the density-dependent numerical code F R A C 3 D V S . A number of two-dimensional simulations were performed using a variety of parameter and boundary conditions to determine i f there was any scenario under which a significant volume o f S G D may be expected. A single three-dimensional model was run. In all simulations, it was found that the groundwater flow regime could be described by three flow systems. These were an onshore system, a seawater recirculation system, and an offshore convection system. Results of all o f the simulations showed that little S G D is expected within the study region. N o groundwater discharge was predicted adjacent to the coastline. Insufficient topographic relief is present within the system to drive fresh groundwater discharge offshore. Instead, saline groundwater recharge was predicted due to the seawater recirculation system. This system resulted in saline groundwater being present in the deeper parts of the aquifer system over a hundred kilometres inland from the shore. Further offshore, a low rate of discharge was predicted due to the offshore convection system. The predicted volume of S G D was approximately 0 . 1 % o f Mississippi River discharge. A second research group from the University of East Carolina, Tulane University, and the University of M i a m i performed a companion study aimed at quantifying S G D using the geochemical tracer radon. Results from their work showed again that little offshore discharge was occurring within the study region. Computed discharge rates were i i i less than 1 cm/d during all sampling periods. These results confirmed the results of the numerical modelling. The results from these studies are significantly different to what was previously found. This may be due to the low groundwater concentrations of radium that were assumed when calculating the rate of offshore discharge. I f higher radium concentrations had been used, as found in deep saline formation waters, then results from both sets of research would be match more closely. Further research is necessary to determine offshore groundwater pathways in order to quantify the concentration of tracer present in any offshore discharge. iv T A B L E OF CONTENTS Abstract i i Table of Contents v List o f Tables ix List of Figures x Acknowledgements xv 1. I N T R O D U C T I O N 1 1.1. S U B M A R I N E G R O U N D W A T E R D I S C H A R G E : D E F I N I T I O N A N D H I S T O R I C A L P E R S P E C T I V E 1 1.2. R E S E A R C H D E F I N I T I O N A N D P U R P O S E 4 1.3. S T U D Y R E G I O N 5 1.4. A D D I T I O N A L R E A S O N S F O R Q U A N T I F Y I N G S U B M A R I N E G R O U N D W A T E R D I S C H A R G E : 7 1.5. P R E V I O U S S U B M A R I N E G R O U N D W A T E R D I S C H A R G E S T U D I E S 8 1.5.1. D I R E C T M E T H O D S 8 1.5.2. T R A C E R T E C H N I Q U E S 9 1.5.3 M O D E L I N G M E T H O D S 12 1.6. G R O U N D W A T E R C O N V E C T I O N I N T H E C O N T I N E N T A L S H E L F 15 2. H Y D R O L O G I C U N I T S O F T H E C O A S T A L L O W L A N D S A Q U I F E R S Y S T E M 22 2.1. I N T R O D U C T I O N 22 2.2. L O C A T I O N O F T H E G E O P R E S S U R E D Z O N E 23 2.3. D I V I S I O N OF A Q U I F E R S Y S T E M I N T O H Y D R O G E O L O G I C U N I T S 24 2.4. H Y D R O G E O L O G I C U N I T S 26 2.5. S U M M A R Y 29 3. N U M E R I C A L F O R M U L A T I O N OF F R A C 3 D V S : A T H R E E - D I M E N S I O N A L , D E N S I T Y D E P E N D E N T G R O U N D W A T E R F L O W M O D E L 34 3.1. I N T R O D U C T I O N 34 3.2. N U M E R I C A L F O R M U L A T I O N 34 3.2.1. F L O W 34 3.2.2. T R A N S P O R T 37 3.3. S O L U T I O N P R O C E D U R E 38 4. D E V E L O P M E N T O F T H E T W O - A N D T H R E E - D I M E N S I O N A L M O D E L S 39 4.1 M O D E L R E G I O N 39 4.2. G R I D G E N E R A T I O N 40 4.2.1. T H R E E - D I M E N S I O N A L M O D E L 40 4.2.2. T W O - D I M E N S I O N A L M O D E L S 43 4.2.2.A. C R O S S S E C T I O N 1 43 4.2.2.B. C R O S S S E C T I O N 2 44 4.3. D I V I S I O N O F T H E D O M A I N I N T O H Y D R O G E O L O G I C U N I T S 44 4.4. P A R A M E T E R V A L U E S 46 v 4.4.1. F L O W 47 4.4.2. T R A N S P O R T 49 4.5. F L U I D D E N S I T Y 50 4.6. T E M P E R A T U R E 52 4.7. I N I T I A L C O N D I T I O N S 53 4.8. R A Y L E I G H N U M B E R C A L C U L A T I O N S 55 5. T W O - D I M E N S I O N A L S I M U L A T I O N S OF S U B M A R I N E G R O U N D W A T E R D I S C H A R G E I N T H E G U L F OF M E X I C O N E A R S O U T H E A S T E R N L O U I S I A N A 77 5.1. I N T R O D U C T I O N 77 5.2. B O U N D A R Y C O N D I T I O N S 79 5.2.1. C R O S S - S E C T I O N 1 79 5.2.2. C R O S S - S E C T I O N 2 81 5.3. R E S U L T S F O R C R O S S - S E C T I O N 1 81 5.3.1. B A S E C A S E S C E N A R I O 81 5.3.1. A . M O D E L V E R I F I C A T I O N 89 5.3.2. H Y D R A U L I C C O N D U C T I V I T Y 92 5.3.2. A . S I M U L A T I O N 1: I N C R E A S E D H Y D R A U L I C C O N D U C T I V I T Y .... 93 5.3.2.B. S I M U L A T I O N 2: R E M O V A L OF T H E N E A R S H O R E M U D L A Y E R 95 5.3.2.C. S I M U L A T I O N 3: L O W P E R M E A B I L I T Y S U R F A C E U N I T 97 5.3.3. D I S P E R S I V I T Y 100 5.3.4. O N S H O R E B O U N D A R Y C O N D I T I O N 101 5.3.5. S H A L L O W F L O W S Y S T E M 104 5.3.6. S E A F L O O R B O U N D A R Y C O N D I T I O N 106 5.3.7. B R I N E S 107 5.4. R E S U L T S F O R C R O S S - S E C T I O N 2 110 5.4.1. B A S E C A S E S C E N A R I O 110 5.4.2. D I S P E R S I V I T Y 113 5.5. D I M E N S I O N L E S S P A R A M E T E R S 115 5.4.1. P E C L E T N U M B E R 115 5.4.1. C O U R A N T N U M B E R 116 5.6. S U M M A R Y 118 6. T H R E E - D I M E N S I O N A L S I M U L A T I O N O F S U B M A R I N E G R O U N D W A T E R D I S C H A R G E I N T H E G U L F OF M E X I C O N E A R S O U T H E A S T E R N L O U I S I A N A 149 6.1. I N T R O D U C T I O N 149 6.2. B O U N D A R Y C O N D I T I O N S 151 6.3. R E S U L T S OF T H E T H R E E - D I M E N S I O N A L S I M U L A T I O N 152 6.3.1. R E S U L T S N E A R C R O S S S E C T I O N 1 152 6.3.2. R E S U L T S N E A R C R O S S S E C T I O N 2 154 6.4. T H R E E - D I M E N S I O N A L D I S T R I B U T I O N O F R E C H A R G E A N D D I S C H A R G E 155 6.5. D I M E N S I O N L E S S P A R A M E T E R S 157 vi 6.5.1. P E C L E T N U M B E R 157 6.5.2. C O U R A N T N U M B E R 158 6.6. S U M M A R Y 159 7. C O M P A R I S O N O F H Y D R O G E O L O G I C M O D E L S U B M A R I N E G R O U N D W A T E R D I S C H A R G E E S T I M A T E S W I T H T R A C E R - B A S E D E S T I M A T E S 167 7.1. I N T R O D U C T I O N 167 7.2. G E O C H E M I C A L T R A C E R S T H A T H A V E B E E N U S E D T O M E A S U R E S U B M A R I N E G R O U N D W A T E R D I S C H A R G E I N T H E S T U D Y R E G I O N 168 7.2.1. R A D I U M (Ra) 168 7.2.2. R A D O N (Rn) 170 7.3. R E S U L T S OF T H E T R A C E R - B A S E D E S T I M A T E S 171 7.3.1. R A D I U M (Ra) 171 7.3.2. R A D O N (Rn) 173 7.4. C O M P A R I S O N B E T W E E N M E T H O D S 174 7.4.1. N U M E R I C A L M O D E L L I N G A N D R A D I U M (Ra) 174 7.3.2. N U M E R I C A L M O D E L L I N G A N D R A D O N (Rn) 175 8. S U M M A R Y A N D C O N C L U S I O N S 184 A P P E N D I X A l . D E P O S I T I O N A L H I S T O R Y O F T H E S T U D Y A R E A 187 A L L I N T R O D U C T I O N 187 A1.2 . S T R U C T U R A L F E A T U R E S 187 A1.3 . D E P O S I T I O N A L P A T T E R N S 189 A1.4. S T R A T I G R A P H I C U N I T S 190 A1.5 . O T H E R F E A T U R E S 193 A 1.6. S U M M A R Y 194 A 2 . B E N C H M A R K P R O B L E M S 199 A2 .1 . I N T R O D U C T I O N 199 A2.2 . B O X P R O B L E M 199 A2.3 . H E N R Y P R O B L E M 201 A2.4 . E L D E R P R O B L E M 202 A2.5 . H Y D R O C O I N P R O B L E M , L E V E L 1, C A S E 5 203 A2.6 . S A L T L A K E P R O B L E M 204 A2.7 . S A L T P O O L P R O B L E M 204 A2.8. R O T A T I O N A L F L O W OF T H R E E I M M I S C I B L E F L U I D S 205 A2.9 . S U M M A R Y 206 A 3 . M E A S U R E D S O L U T E C O N C E N T R A T I O N A N D D E N S I T Y D A T A F R O M T H E S H A L L O W G R O U N D W A T E R W E L L S 214 v i i REFERENCES 217 viii LIST OF T A B L E S Table 2.1. Previously defined hydrogeologic units along with those used for this research 31 Table 4.1. Horizontal hydraulic conductivity estimates reported by Prudic(1991). . 75 Table 4.2. F low and transport parameters used in the base case scenarios of the two-dimenional simulations 75 Table 4.3. Parameter values used for calculation of the solutal and thermal Rayleigh numbers 76 Table 5.1. We l l depth and average fluid density values of groundwater wells that were sampled 125 Table 5.2. Summary of the results obtained from the two-dimensional models 148 Table A L L Stratigraphic units present in Louisiana and Mississippi 196 Table A3 .1 . Measured solute concentration and fluid density values of shallow groundwater wells 215 ix LIST OF FIGURES Figure 1.1. Schematic diagram of processes associated with submarine groundwater discharge 18 Figure 1.2. Map of the southeastern United States showing the location of the study area 19 Figure 1.3. Extent of the hypoxic zone in the current study region for the summers of 1986, and 1996 20 Figure 1.4. Plot of streamlines for results obtained by Wilson (2005) 21 Figure 2.1. Extent of the Coastal Lowlands aquifer system, southeastern United States, and the area of the system that is bounded by the geopressured zone 30 Figure 2.2. Distribution of hydrogeologic layers for section A - A ' located on Figure 2.1 33 Figure 4.1. Comparison of the model region to the tracer study region 57 Figure 4.2. Plan view map of the three-dimensional grid 58 Figure 4.3. Elevation of the surface of the Coastal Lowlands aquifer system 59 Figure 4.4. Elevation of the base of the Coastal Lowlands aquifer system 60 Figure 4.5. Elevation of the land surface for (A) cross-section one and (B) cross-section 2 61 Figure 4.6. Location of cross-sectional models and seafloor sediment samples taken 62 Figure 4.7. Magnified views of the grid for the shoreline region of cross-section 1 63 Figure 4.8. Gr id employed for cross-section 1 65 Figure 4.9. Gr id employed for cross-section 2 66 Figure 4.10. Elevation of the surface of Permeable Zone B 67 Figure 4.11. Elevation of the surface of Permeable Zone C 68 Figure 4.12. Elevation of the surface of Permeable Zone D 69 Figure 4.13. Elevation of the surface of Confining Unit E 70 Figure 4.14. Elevation of the surface of Permeable Zone D 71 Figure 4.15. Outcrop Pattern of the hydrologic units within the modeled region 72 Figure 4.16. Distribution of hydrogeologic layers for cross-section 1 73 Figure 4.17. Distribution of hydrogeologic layers for cross-section 2 74 Figure 5.1. Schematic diagram of the boundary conditions employed for the two-dimensional models 119 Figure 5.2. F lu id density and streamlines for the base case scenario 120 Figure 5.3. Groundwater fluxes through the surface of the model for the base case scenario 121 Figure 5.4. Magnified view of the region surrounding the shoreline for the base case scenario 122 Figure 5.5. Measured fluid density near the location of cross-section 1 123 Figure 5.6. Locations of shallow groundwater wells that were sampled 124 Figure 5.7. Fluid density and streamlines for the case where the hydraulic conductivity of Permeable Zones A , B , and C was increased by a factor of 6 126 Figure 5.8. Groundwater fluxes through the surface o f the model for the case where the hydraulic conductivity of Permeable Zones A , B , and C was increased by a factor of 6 127 Figure 5.9. F lu id density and streamlines for the case where the mud layer was removed from the first 15 km of seafloor 128 Figure 5.10. Groundwater fluxes through the surface of the model for the case where the mud layer was removed from the first 15 km of seafloor 129 Figure 5.11. Fluid density and streamlines for the case where Permeable Zone A was assigned a hydraulic conductivity equal to the mud layer 130 Figure 5.12. Groundwater fluxes through the surface of the model for the case where Permeable Zone A was assigned a hydraulic conductivity equal to the mud layer 131 x i Figure 5.13. Fluid density and streamlines for the case where the dispersivity was increased an order of magnitude 132 Figure 5.14. Groundwater fluxes through the surface of the model for the case where the dispersivity was increased an order of" magnitude 133 Figure 5.15. Fluid density and streamlines for the case where a third type exit boundary was assigned to the onshore nodes 134 Figure 5.16. Groundwater fluxes through the surface of the model for the case where a third type exit boundary was assigned to the onshore nodes 135 Figure 5.17. Fluid density and streamlines for the case where only Permeable Zones A and B were active 136 Figure 5.18. Groundwater fluxes through the surface of the model for the case where only Permeable Zones A and B were active 137 Figure 5.19. Fluid density and streamlines for the case where a specified concentration boundary was applied to the seafloor 138 Figure 5.20. Groundwater fluxes through the surface of the model for the case where a specified concentration boundary was applied to the seafloor 139 Figure 5.21. Fluid density and streamlines for the case that included brines 140 Figure 5.22. Groundwater fluxes through the surface of the model for the case that included brines 141 Figure 5.23. Fluid density and streamlines for the base case scenario of cross-section 2 142 Figure 5.24. Groundwater fluxes through the surface of the model for the base case scenario of cross-section 2 143 Figure 5.25. Magnified view of the groundwater fluxes through the surface of the model for the base case scenario of cross-section 2 for (A) the northern shore of Lake Ponchartrain and (B) the shore of the Gulf of Mexico ..144 Figure 5.26. Fluid density and streamlines for cross-section 2 for the case where dispersivity was increased by an order of magnitude 145 Figure 5.27. Groundwater fluxes through the surface of the model for cross-section 2 for the case where the dispersivity was increased an order of magnitude 146 xii Figure 5.28. Magnified view of the groundwater fluxes through the surface of the model for cross-section 2 for the case where the dispersivity was increased an order of magnitdue for (A) the northern shore of Lake Ponchartrain and (B) the shore of the Gulf of Mexico 147 Figure 6.1. Schematic diagram of the boundary conditions employed for the three-dimensional simulation 160 Figure 6.2. Fluid density and streamlines for the three-dimensional model near cross-section 1 161 Figure 6.3. Groundwater fluxes through the surface of the three-dimensional model near cross-section 1 162 Figure 6.4. Fluid density and streamlines for the three-dimensional model near cross-section 2 163 Figure 6.5. Groundwater fluxes through the surface of the three-dimensional model near cross-section 2 164 Figure 6.6. Predicted recharge and discharge areas for the three-dimensional model 165 Figure 6.7. Predicted recharge and discharge rates for the three-dimensional model 166 Figure 7.1. Generalized diagram (or box model) depicting the sources and sinks of geochemical tracers within the study region 178 Figure 7.2. Study region for the radium tracer research performed by Krest a/(1999) 179 Figure 7.3. Average calculated rate of SGD using radon from the fall of 2003 180 Figure 7.4. Average calculated rate of SGD using radon from the spring of 2004 181 Figure 7.5. Average calculated rate of SGD using radon from all sampling periods 182 Figure 7.6. Comparison of predicted discharge rates for the three-dimensional model and 2 2 2radon 183 Figure A L L Structural features that controlled depositional patterns in the study area 195 Xlll Figure A1.2 . Surficial geologic map of Louisiana and Mississippi 197 Figure A1.3 . Location of salt basins present within the study area 198 Figure A2 .1 . Relative concentration profiles for the Henry Problem 207 Figure A2.2 . Relative concentration plots for the Elder Problems for (A) 10 years, and (B) 20 years 208 Figure A2 .3 . Vector plots for the Elder Problem at 20 years for different levels of discretization 209 Figure A2.4 . Relative concentration and vector plots for the H Y D R O C O L N Problem 210 Figure A2.5 . Salinity values obtained for the Salt Lake Problem 211 Figure A2.6. Salinity breakthrough curves obtained at the outlet for the Salt Pool Problem 212 Figure A2.7 . Evolution of the salinity interfaces for the Rotational F low Problem 213 xiv ACKNOWLEDGEMENTS I would like to thank all those who made this thesis work possible. First and foremost, I would like to thank my advisor, Dr. Leslie Smith. He provided thoughtful insight and support throughout this project, and it was a privilege to work with him. Thanks to the members of my committee, Dr. Roger Beckie and Dr. Uli Mayer. Through courses and responding to numerous questions they provided invaluable assistance. I would like to thank Clay McCoy and Dr. Reide Corbett from East Carolina University, and Dr. Brent McKee from Tulane University for providing the data from their work and allowing me to assist on a sampling cruise. Thanks also to Rob McClaren from the University of Waterloo for assistance with FRAC3DVS and adding the necessary coding. Finally, I would like to thank my colleagues at the University of British Columbia for their comments and support during my time at UBC. A special thanks to Randi Williams for all her help and support along the way. xv 1. INTRODUCTION 1.1. SUBMARINE GROUNDWATER DISCHARGE: DEFINITION AND HISTORICAL PERSPECTIVE Submarine groundwater discharge (SGD) has been the focus of increasing amounts of research in the past decade (Burnett et al, 2001). It has become recognized that groundwater, along with river water, may be an important nutrient source to the coastal marine environment if the former contains elevated nutrient concentrations. Submarine springs and seeps have been recognized for many centuries, but their role in the marine ecosystem was not studied due to the difficulty of quantification and the perception that it was unimportant (Burnett et al, 2003). Simmons and Netherton (1986) state that Kohout (1966) and Kohout and Kolipinski (1967) were among the first researchers to recognize the ecological significance of SGD. Although seepage rates may vary over an area, when averaged out SGD may be important both chemically (due to elevated nutrient concentrations) and volumetrically (due to sea water dilution) in coastal bays and estuaries (Johannes, 1980). A recent report also suggests that groundwater discharge may provide a pathway for bacterial contamination of the coastal environment (Paytan et al, 2004). The definition of SGD is ambiguous. Some groups have defined it to include only that component of fresh groundwater that discharges through the sea floor. This part of SGD is driven by hydraulic gradients that exist onshore. The supply of fresh groundwater essentially comes from recharge to coastal aquifers. Other groups have tended to include recirculated seawater in the definition of SGD, which includes discharge due to tidal pumping, wave set-up, and convection processes. For this research, SGD will be defined 1 as the total mixture of seawater and fresh groundwater flowing out from the seabed, into the coastal water, through the underlying sediments (Destouni and Prieto, 2003). Figure 1.1 provides a generalized diagram of the components of SGD listed above. Taniguchi et al (2002) compiled a list of published SGD data on a worldwide basis. They found that, in general, discharge decreased with increasing surface water depth, or by inference, distance away from shore. This result may be expected, as typical beach faces form topographical lows that induce groundwater discharge. Taniguchi et al (2002) state that one must be cautious about drawing too many conclusions from these results. Most studies have focussed on the coastal zone, and too few studies have been performed to conclusively state that SGD decreases away from the shoreline. Some researchers have found that SGD may be occurring up to distances of kilometres offshore (Burnett et al, 2003). For this to occur, an aquitard that extends offshore must be present to force the freshwater to discharge farther out (see Figure 1.1). Additional reasons for higher SGD rates being measured near the shoreline may be due to tidal pumping and wave-setup. Li et al (1999) presented a theoretical model that accounted for the various sources of groundwater discharge. From their calculations, L i et al (1999) estimate that up to 96% of SGD may be due to groundwater circulation and oscillating flow caused by wave-setup and tides. Discharge of fresh groundwater due to onshore hydraulic gradients may then be expected to make up a much lower proportion of the total SGD. As the potential significance of SGD has become recognized, studies have begun to focus on its quantification. Several methods have been employed, including the use of geochemical tracers (e.g. radon), direct methods (i.e. seepage meters), and 2 hydrogeological models. In most studies, only one method has been employed. Thus, the confidence that can be placed in the results obtained is directly related to the uncertainties associated with the measurement technique. Due to the often diffuse nature of S G D and the limitations of each technique, the potential for erroneous results is quite high. Recently, the need for the use of multiple measurement techniques at a given site has become apparent (Burnett, 1999; Burnett and Turner, 2001). A working group composed of researchers from a variety of disciplines was set up by the Land-Ocean Interactions in the Coastal Zone (LOICZ) and the Scientific Committee on Oceanic Research (SCOR) groups to address this topic. To date, they have performed several intercomparison experiments whose aim was not only to quantify S G D , but also to understand the discrepancies that occur between techniques. A n issue of the journal Biogeochemistry (Volume 66, Issue 1-2, 2003) was published whose whole content was devoted to this topic. In general, they observed that results obtained using seepage meters and geochemical tracers were in relatively good agreement, while those obtained through the use of hydrogeological modelling typically estimated S G D values less than the other two methods. Smith and Zawadzki (2003) state that transient processes that were not included in their model, such as tidal pumping, may be causing the difference in results. More recently, a special issue of Ground Water (Volume 42, Number 7, Special Issue 2004) was published that was entirely dedicated to S G D studies on the Atlantic Coastal Plain of the United States. Results from two of the articles found in this issue confirmed the observation noted above by Smith and Zawadzki (2003). Martin et al (2004) and Cable et al (2004) found that the discrepancy between previous land-based estimates of S G D and offshore measurements of S G D could be explained by mixing in 3 the shallow (< lm) sediment of the seabed. The mixing of seawater with shallow groundwater provided an additional source of SGD that was measured by offshore tracer techniques, but was not accounted for using a numerical model. The shallow mixing may be driven by wave-setup, tidal pumping, and convection. The relative importance of these driving forces was found to be site-specific. 1.2. RESEARCH DEFINITION AND PURPOSE The goal of this project is to examine the interpretation of SGD estimates obtained in the region surrounding the Mississippi River delta by Krest et al (1999). They derived estimates of SGD using radium ( Ra and Ra) isotopes as geochemical tracers. Numerous ocean water column measurements were taken in order to determine the relative Ra contributions to the water column from sources such as river discharge, sediment desorption, and the background concentration present in the ocean waters of the Gulf of Mexico. Using these values, the rate of groundwater discharge required to be entering into their study region to achieve a mass balance was calculated. They estimated that an advective groundwater velocity of about 1 cm/day may be discharging into the Gulf between the mouths of the Mississippi and Atchafalaya Rivers (see Figure 1.2). If this value was assumed to represent the average discharge rate of groundwater through the seabed throughout their study region, than SGD may be a significant source of nutrients to the estuarine environment of southern Louisiana, depending on nutrient concentrations within the groundwater. In order to build upon the work of Krest et al (1999) and to place more confidence in the results obtained, a multiple technique approach was applied to this project. A regional scale hydrogeological model was developed using existing data on the geology, 4 hydrogeology and salinity within the study region. The density-dependent groundwater flow model FRAC3DVS (Therien et al, 2003) was utilized in order to incorporate density effects. A series of simulations were performed in order to examine the effect that boundary conditions and hydrogeologic parameters had on the model-derived estimates of SGD. Estimates of SGD were computed for pre-development (i.e. no internal sources or sinks) conditions. Simulations that included the effects of pumping wells were not performed for reasons that will be discussed later. A second group, with members from the University of East Carolina, Tulane University and the University of Miami (Fla.), obtained an independent estimate of SGD using radon (Rn) and helium (He) as geochemical tracers. They used a series of box models to calculate a mass balance for the tracers, and determine estimates of SGD. Using the tracer-derived estimates, the hydrogeological model was then examined and updated in order to evaluate the reasons for any discrepancies that were found between the two methods. 1.3. STUDY REGION The study area consists of the coastal region west of the Mississippi River delta from 89.0 to 90.8° longitude. It extends seaward just pass the 100 m isobath, and can be seen in Figure 1.2. The onshore region consists of a broad, low-lying coastal plain underlain by thick sedimentary deposits. The Mississippi River flows through the central portion of the study region, and terminates in a typical bird's foot delta at the coastline. The Mississippi River has the third largest drainage basin area, the sixth largest water discharge and the 5 seventh largest suspended sediment load among world rivers (Milliman, 1991). The river has an average freshwater discharge of 380 km3/yr, and an average dissolved nitrogen load of 1.82 Tg N/yr. Low-lying hills surround the Mississippi River to its east, while the western flank exhibits a gently sloping alluvial plain. Topography reaches an elevation of about one hundred metres at the northern extent of the study area. Along with the Mississippi, numerous surface water bodies, including lakes and canals or bayous, dissect the land. The climate of southeastern Louisiana is warm and humid. It is characterized by long, hot, humid summers and relatively short, mild winters (Martin and Whiteman, 1999). The average annual temperature is about 22 °C. Average annual precipitation within the region is 1575 mm/year. Precipitation generally exceeds potential evaporation throughout most of the year, providing an abundant water source to the surface water bodies (Williamson and Grubb, 2001). It is thought that little precipitation is effective in recharging the underlying aquifers, as most of it is lost to surface run-off. In their model results, Williamson and Grubb (2001) estimated that net regional recharge throughout their study area was only 0.6% of the total precipitation, although it did vary significantly throughout their large study region. For example, in southwestern Mississippi and adjacent Louisiana, maximum recharge rates reached 152.4 mm/yr, or about 10% of the average annual precipitation rate. 6 1.4. ADDITIONAL REASONS FOR QUANTIFYING SUBMARINE GROUNDWATER DISCHARGE Quantification of SGD in the region is of significant importance due to the presence of hypoxic waters in the bottom waters of the Gulf of Mexico during the summer season. The areal extent of this zone has been as large as 21 700 km , as determined in 2001 by Turner et al (2005). They define hypoxic waters to be those which contain less than 2 mg/L of dissolved oxygen. It is currently thought that the cause of hypoxia in the region is due to elevated levels of nitrates that have been transported by the Mississippi and Atchafalaya Rivers to the Gulf of Mexico. The source of nitrate is believed to be nitrogen-rich fertilizers, which are applied to farmlands in the Mississippi watershed. The enriched nitrate levels lead to increased production of organic material. When it sinks, the organic material consumes oxygen faster than it can be replenished through vertical mixing in the stratified water column, causing depleted dissolved oxygen concentrations in the bottom waters (Turner et al, 2005). These authors make no mention of the role that groundwater may play in this system, however. If a significant volume of groundwater is being discharged to the Gulf, than characterizing its nitrogen concentration would be important in understanding the source of the problem. Figure 1.3 shows the distribution of the hypoxic zone in the current study area for the years 1986 and 1996. As can be seen from the figure, the zone of hypoxia grew significantly over that time period. Along with its ecological importance, the presence of hypoxic waters in the Gulf of Mexico has had a direct effect on economic activity in the region. Commercial fishing is a primary employer in the region, and has been deleteriously affected by the presence 7 of the hypoxic zone. Over the past three decades, the size, and therefore the price of shrimp in the Gulf of Mexico have been falling (Raloff, 2004). Decreased fish stocks have also led to a decline in sport fishing, which is another source of revenue in coastal Louisiana. An action plan has been adopted to minimize the effects of hypoxia. By 2015, the goal of the multi-agency collaborative plan is to reduce the size of the hypoxic zone to 5000 km 2 by reducing the riverine nitrogen load by 30% of current levels (Mississippi River/Gulf of Mexico Watershed Nutrient Task Force, 2001). Obviously, if SGD in the region is high and nitrate concentrations are elevated, than groundwater could play an important role in the success of the plan. 1.5. PREVIOUS SUBMARINE GROUNDWATER DISCHARGE STUDIES Methods of quantifying SGD fall into three categories: direct methods, tracer techniques, and modelling (Burnett et al, 2001). Direct measurements are limited to the use of seepage meters. Tracer techniques make use of chemical species present naturally in the environment or those that have been introduced anthropogenically. Modelling studies can take several forms, ranging from simple water balance calculations to more complex three-dimensional hydrogeologic models. 1.5.1. DIRECT METHODS The most basic method for measuring SGD is through the use of seepage meters. In their simplest form, they consist of an open chamber that is inserted into the seafloor. A bag is attached to the chamber to allow collection of discharging water. SGD estimates can then be derived based on the volume of water that was collected over a known time 8 period. Seepage meters provide point measurements of SGD, which can give an integrated value of seepage with distance from shoreline when placed in transect (Oberdorder, 2003). Discharge collected may range from fresh to saline. Based on the salinity of the water, it should be possible to calculate the fresh and recirculated seawater components of SGD (Oberdorder, 2003). Several problems have been found to occur when using seepage meters. Anomalously high initial inflows may occur when the meter is first installed due to mechanical properties of the collection bag (Shaw and Prepas, 1989; Shaw and Prepas, 1990). Cable et al (1997) found that these effects could be removed when the bag is pre-filled prior to deployment. Other problems include leakage around the base of the meter, and interference from waves (Burnett et al, 2001). Newer seepage meters that employ heat pulse methodology (see Taniguchi and Fukuo, 1993, for example) may avoid some of these problems. This type of meter also offers the advantage of providing automated seepage measurements with time. 1.5.2. TRACER TECHNIQUES Environmental tracers have been used in numerous studies for quantifying SGD (e.g. Krest et al, 1999). Chemical species that have been employed include radium (e.g. Moore, 1996), radon (e.g. Cable et al, 1996), helium (e.g. Price et al, 2003), and barium (e.g. Shaw et al, 1998). The basic premise behind using a geochemical tracer is to quantify all of its sources and sinks in the coastal water column, then to deduce the amount that is present due to SGD. Quantifying sources and sinks of tracer has proven to be a difficult task, as uncertainties related to the properties of each tracer are present. The 9 choice of tracer used is limited to those that are enriched in groundwater relative to seawater. Ideally, the tracer should also be conservative and easy to measure (Burnett et al, 2001). Tracer-based SGD estimates provide an integrated value of SGD over a larger, diffuse seepage zone, rather than the discrete point measurements provided by seepage meters. This is due to the fact that as the tracer is discharged into the water column it becomes mixed with the ocean water, making it impossible to discern the location of the discharge point. Radium (Ra) isotopes have been the principal species used as geochemical tracers to quantify SGD (e.g. Burnett et al, 1990; Webster et al, 1994; Rama and Moore, 1996; Charette et al, 2003). It is essentially immobile in fresh (non-saline) groundwater because it adsorbs to the surface of sediments. The ratio of Ra adsorbed on sediment to Ra dissolved in the surrounding water varies inversely with the ionic strength of the solution (Krest et al, 1999). Thus, it is released through desorption when brackish and saline waters flush through the sediments (Moore, 1996). Estimates derived through the use of Ra represent the mixture of fresh and salt water that is discharged. Radon (Rn) has been used as a tracer for groundwater discharge studies in both marine and freshwater systems (Cable et al, 1996; Corbett et al, 1997; Corbett et al, 1999; Purkhl and Eisenhauer, 2004). 2 2 2 R n (half life of 3.8 days) is produced by the decay of 2 2 6 R a that is present in the aquifer matrix. As groundwater flows through the sediments it becomes enriched in the dissolved gas. This enrichment, relative to seawater, allows for its use as a tracer. Rn may be discharged by fresh groundwater or through re-circulated seawater, or some combination of the two. Therefore, estimates of SGD obtained through 10 Rn provide an overall value, which includes both fresh and recirculated seawater (Oberdorfer, 2003). The potential use of helium (He) as a tracer of groundwater flow was first recognized by Carter et al (1959). They injected 4He enriched water for experiments at both the laboratory and field scale and found that it was a viable tracer. Although they recognized the presence of naturally occurring 4He, they considered it a background level that had to be subtracted from their measurements. 4He is produced through the decay of uranium and thorium chain elements. In groundwater, its concentration is dependent on the residence time of the groundwater; while in seawater, it is in equilibrium with the atmosphere. Recently, Price et al (2003) suggested the use of naturally occurring 4He as a tracer of groundwater in the coastal environment. They found elevated concentrations of the tracer in brackish groundwater samples taken in Everglades National Park, Florida. Elevated concentrations were also found in the coastal water column by Top et al (2001), suggesting that SGD is occurring into Florida Bay. Several studies have examined the flux of barium (Ba) into coastal environments (e.g. Coffey et al, 1997; Moore 1997). These studies focussed on riverine input of both dissolved and particulate components. Ba adsorbs reversibly onto particulate matter in fresh water rivers. Similar to Ra, it desorbs upon encountering saline waters. Shaw et al (1998) proposed the use of this species as a geochemical tracer of SGD, where the adsorbed source of Ba was the aquifer solids. Shaw et al (1998) found enriched concentrations of Ba in the South Atlantic Bight off the coast of South Carolina, which could not be explained by riverine input and coastal upwelling alone. Based on the 11 concentrations found in groundwater wells, they concluded that SGD must be the source of Ba enrichment. 1.5.3 MODELING METHODS The simplest of the modelling methods are water budget calculations. When applied to SGD studies, the model assumes that all water entering the ground discharges into the marine environment. Water entering the ground is calculated as the difference between precipitation and the sum of surface runoff and evapotranspiration. This form of calculation is typically done on a regional scale, and only attempts to quantify the fresh water component of SGD. Oberdorfer et al (1990) estimated SGD into Tomales Bay, California using this approach. From their calculations, they found that groundwater could be an important component of the nutrient budget of the bay. Hydrograph separation has been used to estimate SGD at the continental and global scales by Zekster et al (1973) and Zekster and Dzhamalov (1981). This technique determines the component of flow contributed to streams and rivers from groundwater, termed baseflow. It is then assumed that groundwater discharge along the shoreline is similar to adjacent stream or river baseflow, and SGD can be estimated from this relation. Only freshwater flows from the shallowest aquifers can be estimated using this method. The use of numerical models in coastal areas has generally been restricted to saltwater intrusion studies. But several studies have used numerical models to estimate SGD in recent years. These range from two-dimensional representations of three-dimensional systems (e.g. Smith and Zawadski, 2003) to fully three-dimensional flow models (e.g. Langevin, 2001, 2003). Numerical models provide estimates of only the 12 freshwater component of SGD if a uniform density fluid is assumed. If density effects are taken into consideration by solving the density-dependent flow equation, than an estimate of total SGD may be obtained. It should be noted, however, that no models to date have attempted to account for the portion of SGD due to wave set-up and tidal pumping. This is probably due to difficulties associated with incorporating these small scale features into larger scale hydrogeological models. Destouni and Prieto (2003) performed a series of hypothetical two-dimensional simulations using SUTRA (Voss, 1984), in which water management scenarios were varied. Three different coastal aquifers were simulated. They found that a linear relation existed between the model predicted amount of SGD and the net land-determined groundwater drainage (i.e. the difference between all groundwater inputs and outputs) at steady state. The relation was present for all three systems, suggesting that it is independent of site specific parameters. Where SGD predictions were high, they found that it was dominated by fresh groundwater. In cases where it was low, discharge was predicted to contain mostly recirculated seawater. They also ran simulations where outputs equalled inputs, and found that a nonzero value of SGD was still predicted. This was due to recharge downstream of pumping wells and fresh groundwater flow paths evading well capture, as well as the recirculation of seawater. Smith and Zawadski (2003) simulated two-dimensional flow into the Gulf of Mexico at the Florida State University Marine Lab using FEFLOW (Diersch 1996). Their model incorporated density effects, accounting for fresh and saline components of SGD. Due to the lack of agreement between their model-derived estimates with those obtained from seepage meters and tracers, they performed a series of sensitivity simulations. 13 Hydraulic conductivity values of the surficial aquifer and the underlying confining unit were varied to a maximum reasonable limit, providing a range of SGD estimates. They found that their steady state models were unable to match the SGD values predicted using seepage meters and geochemical tracers. The results for all of their simulations provided values that were less than the other two methods. Smith and Zawadzki (2003) suggest that the reason for the discrepancy between measurement techniques may be due to transient processes, such as tidal pumping, that were not incorporated into the model. Kaleris et al (2002) modelled SGD into Ekcenforde Bay in the western Baltic Sea using three different scale models, in which only the smallest scale model incorporated density effects. Two scenarios were used for their regional scale model. In the first scenario, groundwater was allowed to discharge only through seeps termed pockmarks. A detailed survey of the sea bottom in the region identified these features. The rest of the seafloor was specified to be impermeable. For the second scenario, groundwater was allowed to discharge everywhere along the sea bottom. Here, the hydraulic conductivity of the sediments varied spatially. Results of the simulations showed that SGD values were similar for both cases. In the first scenario, discharge was forced to occur only through the pockmarks. Groundwater discharge in the second scenario only occurred through areas where the hydraulic conductivity was elevated relative to the rest of the seafloor, providing a similar result to the first scenario. Their smaller scale models were used to show how flow rates varied across the pockmarks. Where density effects were included, they found that discharge was nonuniform and concentrated at the edges of the pockmarks. 14 SGD into Biscayne Bay, Florida was simulated using SEA W A T (Guo and Langevin, 2002) by Langevin (2001, 2003). The large-scale density-dependent three-dimensional model covered an area of approximately 6300 square kilometres. It was calibrated using extensive field data, collected over a period often years. Langevin found that during dry seasons, SGD may have exceeded surface water flow into Biscayne Bay. But throughout the entire simulation period, average discharge was only about ten percent of surface water input. SGD estimates were found to be quite sensitive to the boundary condition used to represent the seafloor. This feature will be discussed later in the thesis. 1.6. GROUNDWATER CONVECTION IN THE CONTINENTAL SHELF As mentioned earlier, SGD is not restricted to near shore regions, and is not driven exclusively by fresh onshore groundwater flow. Groundwater discharge across the seabed may also be driven by recirculation and convection of seawater due to density differences. Convection-driven flow processes range from pore-scale exchange (Burnett et al, 2003) to kilometre-scale convection cells that extend thousands of meters into the subsurface (Wilson, 2003). However, at the time this thesis was prepared, few articles were found within the literature that addressed SGD through the entire continental shelf, nor large scale convection within continental margins. A summary of some of the findings of previous authors is provided below. Burnett et al (2003) discuss the potential for shallow advective pore water exchange in the continental shelf. They state that where surface water depths are less than 100 m, currents due to waves reach the seabed, shearing away loose sediments at the 15 surface. This process can eventually lead to removal of the fine material along the seafloor, thereby increasing its permeability and forming wave ripples. The ripples create micro-scale pressure gradients that lead to local recharge through the upstream side of individual ripples. The water follows a curved path through the ripple before discharging through the downstream side. Burnett et al (2003) state that vertical groundwater velocities may reach up to cm/hr in the uppermost centimetre of rippled sediment. The modelling results of Destouni and Prieto (2003) show that recirculation and convection of seawater may occur on a larger scale than described above. Results of their numerical simulations showed that a series of convection cells developed beneath the seafloor under all groundwater management scenarios that were examined. The cells were quite narrow, and extended the entire aquifer thickness of about 60 m. They formed as a result of density differences along the seabed. SGD from the convection cells was predicted to increase as the volume of water flowing from the onshore portion of their model decreased. Sharp et al (2001) showed that salinity-driven groundwater convection cells may be operating within the continental shelf sediments of the Texas part of the Gulf of Mexico. In a series of scenarios, they found that measured salinity inversions could be formed by the convection of groundwater through thin, alternating sand and shale beds. They assumed that thermal effects would be negligible, as solute-induced density variations were far greater than thermal variations. The salinity distribution within the study region of this project is similar to the Texas region (Williamson and Grubb, 2001). Therefore, the same groundwater flow patterns may be expected within the sediments of the southeastern Louisiana region. 16 Numerical simulations that were performed by Wilson (2003) showed that SGD may be expected hundreds of kilometres offshore due to large-scale convection cells. She ran a series of cross-sectional simulations that represented two-dimensional transects across the Atlantic continental shelf of eastern North America. In the simulations, only the offshore system was modelled. Groundwater salinity was assumed constant and equal to seawater. Convection cells formed as a result of density differences induced by thermal gradients, and were centered on slope breaks in the topography of the seafloor. The temperature varied from approximately 0 °C at the seafloor to almost 200 °C at depths greater than 6 km. Groundwater flow velocities reached cm/yr, although they were low enough that the conductive thermal regime was not disturbed. In 2005, Wilson furthered the study of SGD through the entire continental margin. She ran a series of hypothetical simulations that incorporated fluid density variations due to both solute and temperature gradients. She found that fresh SGD may be expected up to 10 km offshore where sufficient topographic gradients are present (see Figure 1.4 A). However, in simulations where the elevation of the water table did not exceed 1 m up to a distance of 10 km inland from shore, she found that no fresh SGD was predicted. In this case, the sole source of SGD was due to deep recirculation of seawater further offshore driven by convection of saline groundwater through the continental shelf (Figure 1.4 B). Therefore, simulation of the entire continental margin, or at least a portion, may be required to evaluate the total volume of SGD released into a given coastal region. 17 Figure 1.1. Schematic diagram of processes associated with submarine groundwater discharge. Arrows indicate fluid movement. Modif ied from Burnett et al (2003). Edge of the Continental Shelf Figure 1.2. Map of the southeastern United States showing the location of the study area. 1 9 Figure 1.3. Extent of the hypoxic zone in the current study region for the summers of 1986, and 1996. Modified from Rabalais et al, 2002. 20 ^Topographically-driven discharge to the ocean X(km) Figure 1.4. Plot of streamlines for results obtained by Wilson (2005). A. Case with greater onshore topographic relief. B. Case with lesser topographic relief. Streamline interval is 5000 kg/m yr onshore, and 500 kg/m yr for offshore. Modified from Wilson (2005). 2. HYDROLOGIC UNITS OF THE C O A S T A L LOWLANDS AQUIFER SYSTEM 2.1. INTRODUCTION Hydrogeologic units underlying the study region that were incorporated into the regional scale groundwater flow model belong to the Coastal Lowlands aquifer system. The Coastal Lowlands aquifer system (Figure 2.1) is composed of sediments that were deposited from the Oligocene to present. It has been defined by the U.S. Geological Survey to contain all sediments overlying the Vicksburg and Jackson stratigraphic units (Weiss, 1992). The description provided here is a summary of the information found in that report. The Vicksburg and Jackson units were grouped together to form a regional confining unit, generally several hundred meters thick, composed predominantly of dense marine clays that were deposited during the Eocene and Oligocene (Hosman, 1996). This unit separates the Coastal Lowlands aquifer system from the underlying Texas Coastal Uplands aquifer system and the Mississippi Embayment aquifer system to the north. The Coastal Lowlands and Mississippi Embayment aquifer systems are, however, connected through a narrow band that is 50 to 80 kilometres wide in northern Louisiana. The aquifer system extends west to the Rio Grande River. To the east, the sediments undergo a facies change to less permeable limestone in western Florida where the system was truncated. The southern extent of the Coastal Lowlands system is the edge of the Continental Shelf, which is roughly equivalent to the 100 metre bathymetric contour. Here, the sediments grade to less permeable clays. The Coastal Lowlands aquifer system is hydraulically connected to the Mississippi Embayment aquifer system at its northern extent. In this location, Permeable Zone A exists as a thin blanket of 22 permeable sediments with thickness generally less than 100 m. The Vicksburg-Jackson confining unit forms the base of the aquifer system. 2.2. LOCATION OF THE GEOPRESSURED ZONE In the southern portion of the study region, a zone of geopressure is found at depth. Geopressured sediments have been defined by Dickinson (1953) to have pressures exceeding the hydrostatic pressure of a column of water containing 80 000 mg/1 of dissolved solids. This would be equivalent to a pressure of approximately 10.5 KPa per meter (Jones, 1969). The geopressured zone has likely formed as a result of sediment compaction combined with the entrapment of fluid. Sediment deposition occurred throughout the Tertiary and Quaternary periods, often at rapid rates. As sediments accumulated at the surface, the underlying sediments began to compact under the increasing weight. Fluid expulsion was resisted by the presence of compacted clays. Trapped fluid lead to undercompaction of the sediments. Growth faults further reduced fluid movement by isolating individual sand beds. The geopressured zone has formed as result. Although it is likely that some fluid is being expelled across the geopressure boundary, the volume must be small to allow the pressure to persist on a geologic time scale. Where the zone of geopressure rises above the Vicksburg-Jackson unit, the base of the Coastal Lowlands aquifer system is defined where it occurs. The portion of the aquifer system where the geopressured zone forms the base can be seen in Figure 2.1. 23 2.3. DIVISION OF AQUIFER SYSTEM INTO HYDROGEOLOGIC UNITS Sediments comprising the Coastal Lowlands aquifer system were deposited during the Cenozoic. They were laid down as complexly interbedded sands, silts and clays, with lesser amounts of gravel, lignite and limestone (Hosman, 1996). The sediments thicken toward the Gulf of Mexico, where they reach thousands of metres in depth. A brief discussion of the depositional history of the study area can be found in Appendix 1. Although small, individual confining beds are present throughout the Coastal Lowlands aquifer system, they are not regionally extensive in nature. This leads to difficulties in separating the basin into a traditional aquifer and confining unit system. Numerous localized definitions of hydrogeologic units in the study region have been proposed. In southeastern Louisiana, individual beds have been mapped in detail. But, none of the beds have been extended regionally across the entire study area. An example comes from Meyer and Turcan (1955), who described ten different aquifers beneath Baton Rouge. The aquifers were named according to the depth at which they were encountered in the industrial district of the city. They recognized several clay layers between the aquifers, but did not label these as formal confining layers. Their naming scheme was used with the knowledge that it could not be extrapolated to other areas, because the aquifers' depth does not remain constant spatially. Similar naming conventions have been used throughout the study area, some of which can be found in Table 2.1. For the reasons discussed above, the Coastal Lowlands aquifer system defined by the US Geological Survey was adopted for this work. The aquifer system contains five permeable zones and two regional confining units. The average thickness of the system is 24 about 1800 metres, with a maximum thickness of greater than 5 000 metres occurring in an area offshore from southern Louisiana. Layer divisions were first defined at pumping centers (e.g. Baton Rouge), where the greatest amount of data was available. Layers were then extended using information from approximately 550 geophysical logs (Wilson and Hosman, 1988) that were located across the aquifer system (about 100 from offshore locations). The relation of the permeable zones to previously defined hydrologic units can also be found in Table 2.1. The first criterion used to define unit divisions was the presence of areally extensive, low permeability sediments (Weiss and Williamson, 1985). These were designated as regional confining units. Where none were present, large hydraulic conductivity contrasts between permeable sediments were used to assign layers. In locations where permeability contrasts were not present, unit divisions were located to minimize vertical hydraulic head changes within layers. The vertical gradients are due to the presence of localized clay beds, which restrict fluid movement in the vertical direction. If none of the above criteria could be used, the hydrologic units were extended by maintaining their thickness as a constant proportion of the total aquifer system thickness. The criteria listed above were not strictly adhered to at geophysical log locations. At these locations, minor adjustments to layer thickness were made if deemed appropriate. Also, divisions were adjusted in regions where the geopressured zone rose above the Vicksburg-Jackson confining unit. Boundaries between layers were set so that hydrologic continuity of each layer was maintained. Towards the Gulf of Mexico, each layer was gradually pinched out by allowing it to dip below the geopressured zone, until 25 only the top two units remained. In the updip area of the aquifer system, layers were pinched out successively by extending the units to the surface to form outcrop regions. 2.4. HYDROGEOLOGIC UNITS The upper unit of the Coastal Lowlands aquifer system is permeable zone A. It is roughly equivalent to the Chicot aquifer in southwest Louisiana, the "200", "400" and "800" foot sands of Baton Rouge, and the "200", "400", and "700" foot sands of New Orleans. It is composed of sediments deposited during the upper Pleistocene to Holocene. Permeable zone A extends north to the boundary of the aquifer system, where it is connected with the Mississippi River Valley alluvial aquifer of the Mississippi embayment system. The aquifer thins towards the continental shelf, as the dip of the seafloor is greater than the dip of the base of the aquifer. Permeable zone A has an average thickness of 200 metres and a maximum thickness of more than 370 metres in the Terrebonne embayment in offshore Louisiana. The location of the embayment can be found in Figure A 1.1 in Appendix 1. The boundary between permeable zone A and the underlying permeable zone B was defined based on changes in vertical hydraulic gradient and lithologic variations. It is equivalent to the upper part of the Evangeline aquifer and as much as about half of the Chicot aquifer in southwest Louisiana, and the "1200", "1500", and "1700" foot sands in Baton Rouge. It is composed of sediments of upper Pliocene to lower Pleistocene in age. The average thickness of permeable zone B is 580 metres. Maximum thickness of 1720 metres of this unit also occurs in part of the Terrebonne embayment of offshore Louisiana. The unit gradually thickens towards the Gulf of Mexico, where it is truncated by the edge of the continental shelf. 26 Upper Miocene to lower Pliocene sediments make up permeable zone C. It was differentiated from permeable zone B by large changes in vertical hydraulic gradient and local lithology of the sediments. The most substantial gradient changes were found in the pumping centers of Baton Rouge, Louisiana, and Houston, Texas. The permeable zone is equivalent to at least the lower half of the Evangeline aquifer (up to 80% in some locations), and the "2000" and "2400" foot sands of Baton Rouge. Within the study area, it thickens towards the Gulf of Mexico, reaching its maximum thickness of 1920 metres in the Terrebonne embayment offshore from Louisiana. The downdip limit of permeable zone C is located offshore where the top of the geopressured zone rises above the unit's surface. Permeable zones C and D are separated by a regional confining unit in parts of Texas and offshore of western Louisiana. The confining unit, called confining unit D, is not present within the study region of this thesis and thus was not used to place a division between the two permeable zones. Instead, sharp changes in vertical hydraulic gradient were used to place the division. The largest changes in gradient were found at Baton Rouge and Houston. Permeable zone D is approximately equivalent to the Jasper aquifer in southwest Louisiana and includes the "2800" foot sands of Baton Rouge. It is composed of sediments deposited during the mid Miocene series. The average thickness of the unit is about 550 metres, with its maximum thickness occurring in the Terrebonne embayment. The unit is truncated downdip by the geopressured zone in offshore Louisiana. The zone E confining unit lies directly beneath permeable zone D in a portion of Louisiana and Mississippi. It does not extend to the land surface. It is composed of 27 predominantly low permeability marine clays, and is the only regional confining unit present within the study area. It reaches its maximum thickness of about 1200 metres near the Louisiana - Texas state line, and generally thins towards the east. The average thickness of the zone E confining unit is about 300 metres. Beneath the zone E confining unit lies permeable zone E. It is composed of sediments of lower Miocene to upper Oligocene in age. Updip of the confining unit, permeable zone E lies directly beneath permeable zone D. In these areas, the layers were divided based on changes in vertical hydraulic gradient (at Baton Rouge and Houston) and differing Ethology. The average thickness of the unit is 425 metres, and it reaches a thickness of about 1200 metres near Baton Rouge. Included in the permeable zone is the Catahoula sandstone, the uppermost stratigraphic unit of the Vicksburg group. Figure 2.2 displays a cross section through section A - A ' located on Figure 2.1. All hydrogeologic units present within the study region are present in the cross section. From this figure, it is seen that the units comprising the Coastal Lowlands aquifer system form a thick groundwater flow system. However, relative to its horizontal extent, the aquifer system appears quite thin. 28 2.5. SUMMARY The Coastal Lowlands aquifer system forms a horizontally extensive flow system, covering an area of almost 415 000 km 2. It underlies the coastal areas in parts of Texas, Louisiana, Mississippi, Alabama, and Florida, as well as nearby offshore areas. The aquifer system is composed of a thick, interbedded sequence of sediments deposited during the Cenozoic. The maximum thickness reaches more than 5000 meters offshore from southern Louisiana. The aquifer system was divided into hydrogeologic units by the U.S. Geological Survey (Weiss, 1992) using well and geophysical log data from numerous onshore and offshore locations. The description was used to provide a hydrogeological framework for the numerical flow model that was prepared for this thesis. 29 o Permeable Zone A Permeable Zone B Permeable Zone C Permeable Zone D Permeable Zone E Vicksburg - Jackson Confining Unit Line where the top of the geopressured zone intersects the top of the Vicksburg-Jackson confining unit Alabama 2 0 0 K m Figure 2.1. Extent of the Coastal Lowlands aquifer system, southeastern United States, and the area of the system that is bounded by the geopressured zone. Modif ied from Weiss (1992) and Will iamson and Grubb (2001). them stem Series Jones et al(1956), Turcan et al (1966) Meyer a n d Turcan (1955) Rollo (1960) Weiss (1992) Era >-C O Series Southwestern Louisiana Baton Rouge area, Louisiana Baton Rouge area, Louisiana Gulf Coast Quaternary ene Holocene Prairie Formation Montgomery Formation Bentley Formation Williana Formation Chicot aquifer "400 foot" sand "600 foot" sand "800 foot" sand "1000 foot" sand "1200 foot" sand "1500 foot" sand "1 700 foot" sand "400 foot" sand "600 foot" sand Permeable Zone A Quaternary Pleistoo Chicot aquifer Permeable Zone B >cene Foley Formation aquifer (No sands assigned to "800 foot" sand "1000 foot" sand "1200 foot" sand u El lie aquiclude Evangeline Pliocene) "1500 foot" sand "1700 foot" sand Permeable Cenozoic Fleming Formation lie aquiclude Evangeline "2000 foot" sand "2400 foot" sand "2000 foot" sand "2400 foot" sand Zone C ocene of Fisk (1940) Jasper aquifer Burkevi "2800 foot" sand "2800 foot" sand Zone D Confining Unit Tertiary Jasper aquifer Burkevi Lena Member of Fisk (1940) d aquiclude Zone E Confining Unit Unname Permeable Oligocene And Eocene Catahoula Sandstone Vicksburg-Jackson Confining Unit Table 2.1. Previously defined hydrogeologic units along with those used for this research. 31 E CD .C O CD % C O CO C O Nyman and Fayard (1978) Tangipahoa and St. Tammany Parishes Louisiana Rollo (1966) New Orleans area, Louisiana Gandl (1982) Southern Mississippi area Weiss (1992) Gulf Coast b c CD O a o "o N o c CD U g t .0) Shallow Aquifer Upper Ponchatoula Aquifer Lower Ponchatoula Aquifer Big Branch Aquifer Kentwood Aquifer Albita Aquifer Covington Aquifer Slidell Aquifer Shallow Aquifer 200 foot' sand "500 foot' sand "700 foot' sand "1200 foot-sand Tchefuncta Aquifer Hammond Aquifer Amite Aquifer Ramsay Aquifer Franklinton Aquifer CD c 0 O o o Alluvium Loess Terrace deposits Citronelle Formation Miocene Aquifer System piigocene Aquifer System Permeable Zone A Permeable Zone B Permeable Zone C Zone D [Confining Unit Zone E Confining! Unit Permeable Zone E Vicksburg-Jackson Confining Unit Table 2.1 continued. Previously defined hydrogeologic units along with those used for this research. 32 X ( k m ) Figure 2.2. Distribution of hydrogeologic layers for section A - A ' located on Figure 2 . 1 . Point S L denotes the location of the shoreline. 3. NUMERICAL FORMULATION OF FRAC3DVS: A THREE-DIMENSIONAL, DENSITY DEPENDENT GROUNDWATER FLOW MODEL 3.1. INTRODUCTION The numerical simulator FRAC3DVS (Therien et al, 2003) was chosen for this project. It is a three-dimensional finite element model with numerous capabilities. For this project, the density-dependent flow component of the model was employed. Prior to its use here, this aspect of FRAC3DVS had received only limited testing. Therefore, a series of benchmark problems were simulated in order to verify its accuracy. A discussion of the problems simulated and the results obtained can be found in Appendix 2 at the end of this thesis. 3.2. NUMERICAL FORMULATION The numerical formulation used for the density-dependent flow component is based on that described by Frind (1982), and will only be summarized here. At the time this thesis was prepared, the FRAC3DVS user's manual had not been updated to include any documentation. Therefore, the following description has been taken from Frind (1982). 3.2.1. FLOW The input variable in FRAC3DVS describing flow in a non-uniform density fluid is equivalent freshwater head, hf. It is defined by equation (3.1): h r = - ? - + z (3.1) Peg 34 2 3 where p is fluid pressure (M/LT ), p 0 is the reference (freshwater) fluid density (M/L ), g is acceleration due to gravity (L/T 2) and z is the elevation (L) relative to a chosen datum. The Darcy equation for density-dependent flow is: dp_ KdXj (3.2) where qi is the flux in the i-direction (L/T), ky- is the permeability tensor (L 2), p. is the dynamic viscosity (M/L/T), p is the density (M/L 3) and r\j = 0 or 1 for the horizontal and vertical directions, respectively. Substituting equation (3.1) into (3.2) yields: dhf KdXJ (3.3) where pr is the relative density, defined as: Po (3.4) and KM is the freshwater hydraulic conductivity, defined as: PoS (3.5) assuming an isothermal system (i.e. p. = constant). Note that equation (3.5) shows that under this formulation, hydraulic conductivity is independent of fluid density. The relationship between fluid density and solute concentration employed in the model is linear. It is expressed as: p = Po(l + jC) (3.6) where y is a constant derived from the maximum fluid density, p m a x , given by: 35 r Pmax _ j Po (3.7) The concentration term, C (dimensionless), in equation (3.6) is written in relative terms, and can be defined by: (3.8) c max where c is the nodal solute concentration (M/L 3), and c m a x is the maximum solute concentration (M/L ) present in the system. Substituting equation (3.7) into equation (3.3) allows for density to be removed, yielding: dhf KdxJ + yCr?j (3.9) The governing equation for flow can than be written as: _d_ 3X; dhf K8XJ + yCnj 8hf ~dt (3.10) S s in equation (3.10) is the specific storage (1/L), which is defined as: Ss=Pog{a + nj3) (3.11) where a is the compressibility of the porous medium (LT 2 /M), n is the porosity (dimensionless), and p is the compressibility of the fluid (LT 2 /M). Note that fluid density has been assumed to be solely a function of solute concentration (assuming isothermal conditions, and p « a ) . Under the formulation of equation (3.11), S s is assumed to be independent of fluid density. 36 3.2.2. TRANSPORT The solute continuity equation, written in terms of relative concentration, is given by: _d_ 8X; D,. dC u dx V OXJ J dx, dt (3.12) where Vj is the fluid velocity (L/T), and D y is the dispersion tensor (L 2/T). The velocity dependent terms of Dy are defined as: at\v\dij+{al-at)V-p-\ + (DmTi^ (3.13) 1 v where ai is the longitudinal dispersivity (L), a t is the transverse dispersivity (L), D m is the molecular diffusion coefficient (L 2/T), T y is the tortuosity tensor (unitless; assumed equal to 1), |v| is the magnitude of the groundwater velocity vector (L/T), and Sy is the Kronecker delta (dimensionless), which has the properties: 8„ = 0 for / * j ij J 5,, = 1 for i - j (3.14) Equation (3.12) can be simplified by assuming Cdvjdxi « v. dC/dxi. The final form of the equation governing solute transport is then given by dC 8C dx, 11 dx J J ' dx: dt (3.15) 37 3.3. SOLUTION PROCEDURE Equations (3.10) and (3.15) are solved using a Picard iteration scheme. This method linearizes the system of equations, and solves them via successive substitution until the convergence criterion is satisfied. Central spatial weighting and fully implicit temporal weighting schemes were employed for all simulations. A full description of the solution procedure can be found in Diersch and Kolditz (2002). If convergence is not obtained, the current time-step is halved, and the iteration scheme begins again with the concentration values from the previous time. The transient term on the right hand side of equations (3.10) and (3.15) is solved using a standard finite difference formulation. The solution is marched forward in time using an adaptive time-step procedure. Starting from the specified initial step size, time-steps are either increased or decreased based on the maximum percent change in chosen parameters. Parameters that may be used include head and concentration. It is also possible to set a maximum timestep value prior to commencing a simulation. Other time-stepping procedures are available but were not used and thus will not be discussed here. 38 4. DEVELOPMENT OF THE TWO- AND THREE-DIMENSIONAL MODELS 4.1 MODEL REGION The offshore study region for the tracer portion of the project is located from 89° to 90.8° from the Louisiana coastline out to approximately the 100 m isobath. The model domain was extended outside of this region in order to include onshore flow processes that will directly influence the predicted amounts of SGD, and to prevent the model's boundary conditions from determining the generated solution within the tracer study region. Figure 4.1 depicts both the tracer study region and the modeled region. The western and eastern boundaries of the model were placed at approximately 92.1° and 88.6°, respectively. These locations were chosen to follow approximate topographic divides. Groundwater flow paths likely follow the topography of the land surface at shallow depths. However, the large vertical extent of the Coastal Lowlands aquifer system makes it unlikely that groundwater flow paths do not cross these boundaries at depth. The northern extent of the model was placed at about 31.4°. At this location the elevation of the land surface reaches almost 100 m, providing a driving force for groundwater flow to the coast. The southern boundary of the model was placed at the continental shelf. Here, the sediments of Permeable Zones A and B undergo a facies change to less permeable marine clays (Williamson and Grubb, 2001). Permeable Zones C, D, and E are truncated by the zone of geopressure prior to the edge of the shelf. Little hydrologic information was available to extend the model further. It was expected that the majority of groundwater discharged through the seafloor would occur closer to shore, and the presence of this 39 boundary should have little effect on the model-predicted values of SGD. Overall, the modelled region extends 320 km in the west to east direction, and 400 km in the north to south direction along the western boundary (the southern boundary is variable owing to the geometry of the continental shelf). 4.2. GRID GENERATION 4.2.1. THREE-DIMENSIONAL MODEL The three-dimensional grid was generated by discretizing the model region into variably sized, rectangular elements. The grid was oriented such that the origin was placed at the northwest corner of the model region, with the X dimension increasing towards the shoreline. Figure 4.2 gives a plan view of the grid over the study region. As can be seen from the figure, the grid extends past the edge of the continental shelf in the southern part of the study area. Elements lying outside of this boundary were deactivated within FRAC3DVS. Nodal spacing in the X-direction was set to a uniform value of 5 kilometers from the origin to X = 80 kilometers. From X = 80 kilometers to X = 280 kilometers, a uniform spacing of 2 kilometers was employed. This area was given an increased level of discretization due to the presence of the shoreline and the freshwater/saltwater interface. Early simulations showed that in this region nodal spacing needed to be refined in order to achieve a realistic solution due to sharp salinity gradients. From X = 280 kilometers to the continental shelf boundary, nodal spacing was returned to 5 kilometers. Node spacing in the Y-direction was uniform at 5 kilometers. The total number of nodes per layer was then 9165. 40 The flow system was discretized in the vertical direction using 50 variably spaced nodal layers. The distance between nodes was kept small near the surface of the model domain, because the highest groundwater velocities should be found there due to the presence of more permeable units. Nodal spacing was gradationally increased with depth, as it was expected that less vigorous groundwater flow would occur towards the base of the system. The model domain was divided into five layers. Each layer was divided using a number of sublayers. Nodal spacing within each layer was constant for a given lateral position. The uppermost layer (layer 1) had a uniform thickness of 5 m. It was subdivided using four layers of nodes with 1 m separation between each. The bottom of layer 2 was set at five percent of the total thickness of the domain (i.e. surface elevation minus 5% of the total thickness). It was subdivided using nine layers of nodes. Maximum and minimum nodal spacing was 27 and 1 m, respectively. The bottom of layer 3 was placed at twenty percent of the total domain thickness. It was subdivided using nine nodal layers, yielding maximum and minimum node spacing of 83 and 5 m, respectively. The bottom of layer 4 was set at 50% of the total domain thickness. Nine layers of nodes were used to subdivide it, resulting in maximum and minimum nodal spacing of 167 and 9 m. Layer 5 comprised the lower fifty percent of the flow domain. It was subdivided into 14 layers using 13 sheets of nodes. The maximum and minimum vertical node spacing for this layer was 199 and 11 m. In total, the three-dimensional grid was made up of fifty sheets of nodes. The resulting total number of nodes was 458 250, comprising 439 040 rectangular elements. 41 Nodes at both the surface and the base of the grid were assigned spatially varying elevations. The elevation of the land and seafloor surface was taken from Williamson et al (1990), while the elevation of the base of the aquifer system was taken from Weiss (1992). Surface files were created by overlaying the model grid onto maps of the surface and base. Control points were determined and input into the grid-building program GridBuilder (McLaren, 2003). Nodal elevations were calculated through a kriging procedure available in the program. The generated layers were then read into FRAC3DVS in order to define the surface and base of the aquifer system. Figures 4.3 and 4.4 give contour maps of the generated surface files which display the elevation of the land and seafloor surface and the elevation of the base of the Coastal Lowlands aquifer system, respectively. The zero meter contour in Figure 4.3 denotes the coastline of the model. Little topographic relief is present in this area. A broad, flat coastal plain occupies the majority of the model. The surficial topography does not reach one meter in elevation up to a distance of 10 - 15 km landward from shore. A plot of the elevation of the land surface with distance for cross-sections one and two (described below) is provided in Figure 4.5 to show how little topographic relief is present within the majority of the model region. At the northern edge of the model region, topographic relief reaches just over one hundred meters. Nodes with elevation less than zero meters were assigned to the seafloor. 42 4.2.2. TWO-DIMENSIONAL MODELS A series of quasi two-dimensional models (i.e. unit thickness in the transverse horizontal direction) were run due to the extended length of time required to simulate the three-dimensional model. Run times for the three-dimensional model were on the order of months. Two cross sections were chosen, whose location can be seen in Figure 4.6. The first section, denoted as A - A' in Figure 4.6, was chosen to approximately follow a flow line as determined from initial three-dimensional modeling. B - B' in Figure 4.6 depicts the second section, which was chosen to show the effect that Lake Ponchartrain would have on SGD into the study region. The lake is actually a shallow saline inlet that is directly connected to the Gulf of Mexico. 4.2.2.A. CROSS SECTION 1 Simulating the flow system in two dimensions allowed for increased nodal discretization. Horizontal discretization was again set to be variable. Nodal spacing was set at 1 km for the first 235 km from the northern boundary to a distance of 3 km from the shoreline. Here, the grid was refined for the next 3 km, with nodal spacing gradually decreasing to 500, 100, 50, 20 and 10 m. Nodal spacing was increased at the same interval for the first 3 km offshore in the same manner. Figure 4.7 shows the gradational change in horizontal node spacing near the coastline at a magnified view. Discretization was returned to 1 km for the remainder of the domain. A unit thickness was used for the third dimension. The result was 798 nodes per layer (i.e. 399 nodes times 2). Vertical discretization was kept the same as for the three-dimensional model. Using this level of discretization, the resulting total number of nodes was 39 900. A diagram of the grid employed for this series of simulations can be found in Figure 4.8. Darker areas on the 43 figure represent regions where the level of discretization was too fine to be seen at the scale used. 4.2.2.B. CROSS SECTION 2 The pattern of horizontal discretization for this model was similar to the first cross section. Nodal spacing was gradually increased approaching each shoreline, and than decreased moving landward. The same level of refinement used in the first cross section was again used here (see Figure 4.7), resulting in 880 nodes per layer. Vertical discretization was again set at 50 variably spaced nodal layers. The total number of nodes was 44 000. A diagram of the grid used can be found in Figure 4.9. 4.3. DIVISION OF THE DOMAIN INTO HYDROGEOLOGIC UNITS The model domain was divided into hydrogeologic layers through a similar procedure to that used for dividing the surface and base of the flow domain. Information on the depth of occurrence and extent of hydrogeologic layers can be found in Weiss (1992). The model grid was overlain on maps of the surface elevation of each layer, and control points were determined. The control points were then input into GridBuilder, and surface files were created through the kriging procedure. Where a hydrogeologic layer was not present, the elevation assigned to the nodes was set below the base of the system. In the outcrop region of each layer, the nodal elevation was set above the top of the system. The generated surface files were read into FRAC3DVS to divide the aquifer system into hydrogeologic layers. All elements falling between two surface files were assigned to a single unit. Where individual hydrogeologic layers were found to pinch out 44 artificially due to the level of vertical discretization, the surface files were modified to preserve hydrologic continuity. Due to the generally thick nature of most units, it is not thought that this process should have a significant effect on the overall flow system. Figures 4.10 - 4.14 provide contour maps of the surface elevation and extent of each hydrogeologic unit. Figure 4.15 provides a map depicting the outcrop pattern of the hydrogeologic layers. Note that Permeable Zone E and Confining Unit E do not outcrop within the modeled region. Core (up to 3 m depth) and seafloor sediment samples were taken at fifteen offshore locations. The location of the samples taken can be seen in Figure 4.6. In all samples, the sediments were found to contain only fine-grained muds. For this reason, an extra layer was added to the model to simulate a mud unit on the sea bottom. The layer was placed along the sea bottom, and given a uniform thickness of 5 m. The location and extent of the hydrogeologic layers for the two-dimensional models was determined from the surface files created for the three-dimensional model. Diagrams depicting the distribution of hydrogeologic layer for cross sections 1 and 2 can be found in Figures 4.16 and 4.17, respectively. Numerous growth faults have developed in the thick sedimentary sequence that makes up the Coastal Lowlands aquifer system due to the settling of deposited sediments. FRAC3DVS uses a fixed coordinate system; therefore any compaction-driven flow components were neglected. Most of these faults parallel the shoreline, with some of them believed to be currently active. The effect the faults have on groundwater flow in the region is not well understood (Williamson and Grubb, 2001). A report by Whiteman (1979) describes one fault as a leaky flow barrier. Williamson and Grubb (2001) note that 45 in most cases fault throws are not great enough to offset individual sand beds. Therefore, their effect on groundwater flow should be minimal. In a practical sense, it would have been difficult to include the faults in the model. Although many of the faults have been mapped at the surface, the depth the faults extend is often unknown. Assuming depths would have added increased uncertainty into the model results. For these reasons, they were not included in the model. 4.4. PARAMETER VALUES No field measurements of hydrogeologic parameters were made as a part of this study. Values used in the model were obtained from modeling studies carried out by Martin and Whiteman (1999) and Williamson and Grubb (2001). Martin and Whiteman (1999) modeled the freshwater portion of the eastern Coastal Lowlands aquifer system in Louisiana, Mississippi, Alabama, and Florida. The model of Williamson and Grubb (2001) included the entire Coastal Lowlands aquifer system, as well as the underlying Mississippi Embayment and Texas Coastal Uplands aquifer systems. Parameter values reported by both groups of authors were determined from their calibrated models. Both of these models were prepared as part of the US Geological Survey's Regional Aquifer System Analysis (RASA) project, whose goal was to characterize pre- and post-development groundwater flow in regional aquifer systems. The regions of both studies contained the model region that was simulated here. Only those values that were used in the base case scenario of the two-dimensional models will be described here. Parameters that were varied in subsequent simulations will be described in the relevant sections of the results. 46 4.4.1. FLOW Field measurements from more than 1500 aquifer tests and specific capacity data from nearly 6000 wells were used to obtain estimates for horizontal hydraulic conductivity in the model simulated by Williamson and Grubb (2001). Of these, 3077 were taken from within the Coastal Lowlands aquifer system (Martin and Whiteman, 1999). Prudic (1991) summarized the estimates obtained from the analyses of the test data. He found that the conductivity values could not be contoured due to large variations that were found to occur over short distances. Therefore, the estimates used for the initial model runs of the U.S. Geological Survey's Regional Aquifer System Analysis were derived from statistical analyses performed by Prudic (1991). He grouped the data by hydrogeologic layer, geographic area and a combination of the two to calculate average horizontal conductivity values. The geometric mean values he reported can be found in Table 4.1. The values that were used for the final models of both U.S Geological Survey groups were obtained through model calibration. Williamson and Grubb (2001) provide contour plots of hydraulic conductivity values by layer determined through model calibration. The plots, in general, show that conductivity increases from western Texas to eastern Louisiana, and decreases with depth, but there is significant deviation from this trend. Hydraulic conductivity values used for this study were derived from the contour plots and values given by Martin and Whiteman (1999) and can be found in Table 4.2. A single value was assigned to each hydrogeologic unit. The hydraulic conductivity values in Table 4.2 are similar to those reported by Prudic (1991) for Permeable Zones A and B. The conductivity values for Permeable Zones C, D, and E are somewhat lower. However, the majority of well tests 47 that were performed in these layers were taken from the more permeable outcrop regions (Wilson and Hosman, 1988). In calibrating their flow model, Williamson and Grubb (2001) found that the permeability of these units was much lower in southern Louisiana. Hydraulic conductivity values of the seafloor mud unit and Confining Unit E were set to values that were deemed representative of low permeability sedimentary materials. The values that were used for vertical hydraulic conductivity in the U.S. Geological Survey models were not provided by either group of authors. Therefore, a uniform anisotropy value of 10:1 was assigned to each layer. It is thought that this value should be appropriate in the sedimentary setting found in the model region. No mention was made by either group of authors regarding the porosity values that were used for their simulations. Therefore, a uniform value of 0.3 was assigned to all nodes. This value falls within the range given for sedimentary materials by Domenico and Schwartz (1997). The specific storage used in the U.S. Geological survey models was derived through multiple simulations by Kuiper (1994). Using a parameter estimation method, he found that his models were quite sensitive to the value of storage used. But, in all cases, his models indicated that a very small value for the storage coefficient yielded the best results. A specific storage of 1.55 x 10"7/m was deemed to give the best result. This value was also employed for the simulations described here. A summary of the flow parameters used for the base case scenario may be found in Table 4.2. 48 4.4.2. TRANSPORT In the models prepared by the U.S. Geological Survey, all nodes were assigned a spatially variable, but temporally constant value for solute density. They assumed that over the time of their simulations, the spatial distribution of solute concentration would not change (Williamson and Grubb, 2001). This allowed them to neglect an equation for solute transport. A uniform value of longitudinal dispersivity of 5 m was assigned to all cells within the model domain. This value falls within the range given by Gelhar et al (1992) for data classified to have high reliability at the scale of interest. Higher values for this parameter were tried in early simulations, but it was found that an unrealistic amount of dispersion occurred. Lower values, although potentially more appropriate in a physical sense, were not used for the base case scenario due to the level of discretization and associated longer simulation time that would have been required to obtain a convergent solution. Transverse horizontal dispersivity and vertical dispersivity were set at 0.5 m for all cells. Both of these values fall just outside the upper range for data classified to have high reliability by Gelhar et al (1992). But, due to the reason mentioned above, lower values for these parameters were not used. It should be noted, however, that the ratios between the longitudinal dispersivity and transverse horizontal and vertical dispersivities are within the range given by Gelhar et al (1992) for data of intermediate reliability. For the three-dimensional simulation, the dispersivity values were increased to accommodate the coarser grid spacing that was required. For this simulation, the dispersivity values listed above were increased by an order of magnitude. 49 The diffusion coefficient for the solute was set at 8.64 x 10"5 m/d for all cells, similar to the value employed by Smith and Zawadzki (2003) for their models in a coastal setting. A summary of the transport parameters used for the base case scenario can also be found in Table 4.2. 4.5. FLUID DENSITY Groundwater salinity data collected by research group members from E C U and UBC within the study region was only available for a limited number of shallow wells. The deepest well had a depth of 130 m, as compared to the maximum depth of the flow system of about 6 km. Therefore, data files compiled for the U.S. Geological Survey's RASA study (Williamson and Grubb, 2001) along with maps of the distribution of dissolved solids found in Pettijohn (1996) were used to determine fluid density within the model region. Salinity values ranged from freshwater (density of 1000 kg/m ) in the onshore regions of the system to seawater (density of 1025 kg/m3) beneath the seafloor of the Gulf of Mexico. Groundwater composed of brines has lead to fluid densities of up to 1150 kg/m3 being present at the offshore base of the system. The elevated density at the base of the system is due to the presence of numerous salt domes. The majority of the salt domes are thought to occur within a few hundred of meters of the seabed (Beckman and Williamson, 1990). Fluid density greater than 1065 kg/m3 was found only to occur at isolated locations. These locales are presumably associated with salt domes. Due to the fact that this feature was not accounted for within the model, it would not be possible to maintain the pockets of elevated density as the model was marched through time. Therefore, the maximum density was set to 1065 kg/m3 in simulations that included brines. 50 At this point, it is necessary to discuss several assumptions related to FRAC3DVS. The first assumption is the linear equation of state relating fluid density to solute concentration (see equation 3.6). Several authors have employed a nonlinear equation of state to account for volume changes due to the addition of solute, especially where fluid densities greater than seawater were present (see Diersch and Kolditz, 2002, for example). However, Williamson and Grubb (2001) calculated fluid densities in the Gulf region using a linear multiple regression equation that incorporated dissolved solids concentrations, temperature and hydrostatic pressure. Using the regression equation, a good linear fit was found to exist between dissolved solids concentration and the calculated fluid density. The relation proved to be valid for waters of all densities encountered, and was present for all layers. Therefore, the linear equation of state in FRAC3DVS was not modified for this work. The second assumption is the use of traditional Fickian dispersion (see equation 3.12). Traditional Fickian dispersion assumes that the dispersive mass flux of a solute is linearly proportional to the concentration gradient. However, several authors have suggested that an extended, nonlinear version of the dispersion coefficient may be more appropriate, particularly where large density contrasts are present (see, Watson et al, 2002, for example). These authors argue that gravity effects become non-negligible in these cases, and act to reduce the amount of dispersion that occurs. There have been results both supporting and refuting the extended version. Schotting et al (1999) found that the extended version was required to simulate freshwater being displaced by brine (maximum density of 1064 kg/m ) in a column experiment. The extended version was only required when density contrasts between the resident and displacing fluids were 51 high. Experiments involving high absolute concentration, but low density gradients could be modeled using traditional Fickian diffusion. Johanssen et al (2002) found that traditional Fickian dispersion was sufficient to simulate a saltwater upcoming process with fluid densities up to 1070 kg/m3. Due to the lack of definitive evidence that a nonlinear form of dispersion produces more appropriate results, the traditional Fickian form was employed. 4.6. TEMPERATURE Groundwater temperatures within the model domain vary from about 20 °C near the land surface to more than 150 °C at depths greater than 3000 m (Williamson et al, 1990) where groundwater of elevated density is present. At the time this study was being performed, no equation for heat transport was present in FRAC3DVS. Therefore, temperature dependence was not included. Ranganathan and Hanor (1988) performed several simulations in sediments surrounding a generic salt dome. They found that flow velocities induced by salinity gradients were up to 20 times greater than those generated due to thermal gradients. They did not, however, perform any simulations in which both salinity and thermal effects were considered. Evans et al (1991) ran a series of similar models of groundwater flow near a generic salt dome, in which both salinity and thermal effects were taken into account. They found that salinity gradients were the dominant flow driving mechanism as well, except in cases where the background salinity was high. Where the background density was set to 1100 kg/m3 or higher, thermally induced flow became apparent. It is believed that neglecting heat transport should not greatly affect the results as temperature effects near the land surface should be minimal. Groundwater in the region 52 has a relatively uniform temperature near the surface, and the majority of the salt domes present in the region are thought to occur hundreds of meters into the subsurface (Beckman and Williamson, 1990). 4.7. INITIAL CONDITIONS Initial salinities were assigned to each node using the data files provided in the report by Williamson and Grubb (2001). This allowed the model to march forward in time starting from an initial salinity distribution that was close to what has been measured in the field. Values were assigned in the same manner as the subdivision of the system into hydrogeologic layers. Salinity values and the depth of their occurrence from the U.S. Geological Survey's data files were input into GridBuilder. The values were smoothed using the kriging procedure to create another set of surface files. Files containing initial salinities corresponding to fluid densities of 1000, 1015, 1025, 1035, 1045, 1055, and 1065 kg/m3 were created. Initial nodal salinities were then assigned by choosing all nodes falling between a set of surfaces. In simulations where the maximum density was set to seawater, nodes that previously had densities greater than 1025 kg/m3 were assigned an initial density of seawater. The resulting salinity distribution was a series of sloping layers of constant density, separated by abrupt interfaces. The interfaces were found to blur due to dispersion within the first few hundred days of simulation time. Initial salinity values for the two-dimensional models were derived from the surface files created for the three-dimensional model. The surface files were modified where nodes were present without values due to the increased level of discretization. Existing values were linearly interpolated to make the surfaces continuous. 53 Typically, initial freshwater head values would also be assigned to each node. This was done in initial model runs, where each node was assigned an equivalent freshwater head corresponding to its initial fluid density. But, it was found that the model experienced difficulties in achieving a convergent solution using the initial head distribution, causing unreasonably long simulation times. It is believed that the problem was due to the large density contrasts that were present initially. It was found that the model was able to march forward in time where initial heads were set to zero. Changing the initial freshwater head distribution did not affect the results that were obtained. Therefore, all subsequent simulations used this initial condition. A two-dimensional simulation was performed using cross-section one to verify that the imposed set of initial conditions did not affect the model outcome. If the same results could be obtained using a different starting point, then evidence would be provided that the initial conditions did not determine the final result. Frind (1982) performed a similar numerical experiment to verify his model of Henry's problem (see Appendix 2 for a description of the problem). For this simulation, all nodes were specified to have an initial concentration equal to seawater. All boundary conditions were held constant in order to isolate any variations due to the different starting point. Results from the simulation proved to be identical to those obtained previously, as determined through comparison of the density distribution and fluxes through the model surface. The sole difference was an increased amount of simulation time required to flush the salt out of the system and reach equilibrium. The result implies that the initial conditions that were imposed did not dictate the final result of the model. 54 4.8. RAYLEIGH NUMBER CALCULATIONS The Rayleigh number (Raw), defined as the ratio between buoyancy and viscous forces, is a dimensionless number that is commonly used to predict whether or not free convection may occur. It was derived under the assumption of homogeneous parameters and simple flow domain geometries (i.e. rectangular), and thus is not strictly applicable here. However, RaN may be used as a rough calculation to determine whether or not convective flow patterns may develop within the Coastal Lowlands aquifer system. Diersch and Kolditz (1998, 2002) define the Rayleigh number as: RaN = RaT +Ras (4.1) where Raj and Ras are the thermal and solutal Rayleigh numbers, respectively. Rax is given by: BATKd RaT = - (4.2) A where P is the thermal fluid expansion coefficient (1/°C), AT is the temperature difference (°C) measured across the thickness d (L), K is hydraulic conductivity (L/T), and A is the thermal diffusivity (L 2/T). Ras can be calculated using: R aACKl where AC is the concentration difference (M/L 3) measured across the thickness d, n is porosity (unitless), and D m is the molecular diffusion coefficient (L /T). The parameter a in equation (4.3) is the fluid density difference ratio normalized by the maximum concentration. It is given by: cc = (p™-p°"> /cmax (4.4) 55 where p m a x and p 0 are the maximum and reference fluid densities (M/L 3), and C m a x is the maximum solute concentration (M/L 3). Studies have shown that a convective flow regime may be expected when Raw exceeds a critical value of 4K 2, or about 40 (Diersch and Kolditz, 2002). Previous numerical studies have investigated the stability of flow systems which develop at high Raw (Diersch and Kolditz, 2002). These studies have found that, in general, a stable convergent solution exists between Raw of 40 and 240-300. Between 240 and 300, a second critical value is reached, above which unstable (transient) convection develops. Numerical simulations may converge to a single solution that is often dependent on the imposed initial conditions. Rax and Ras were calculated assuming an average system thickness of 3000 m, and an average hydraulic conductivity of 2 m/d. It was assumed that groundwater at the surface had a temperature of 20 °C and a density of 1000 kg/m3 at the surface, while at the base the groundwater had a temperature of 150 °C with a fluid density of 1065 kg/m3. All parameters used for calculation of Raj and Ras may be found in Table 4.3. Rax and Ras were found to have values of approximately 6.0 x 103 and 1.5 x 107, respectively. Clearly a convective flow regime will be expected within the Coastal Lowlands aquifer system, as both components of RaN are greater than the second critical value of RaN- Comparing the Rax with Ras, it appears as though fluid density variations due to solute concentration will be significantly more important than those due to thermal effects. This observation was previously reported by Sharp et al (2001) for sediments in the Gulf Coast region of Texas. Concerns regarding the stability of the simulations that were performed will be discussed in the next chapter. 56 Figure 4.1. Comparison of the model region to the tracer study region. 5 7 Louisiana Edge of the Continental Shelf t N O 1 0 0 2 0 0 km Figure 4.2. Plan view map of the three-dimensional grid. 58 0 100 200 300 400 X(km) Figure 4.3. Elevation of the surface of the Coastal Lowlands aquifer system. Contours are in meters. 0 100 200 300 400 X(km) Figure 4.4. Elevation of the base of the Coastal Lowlands aquifer system. Contours are in meters. Figure 4.5. Elevation of the land surface for (A) cross-section one and (B) cross-section 2. Figure 4.6. Location of cross-sectional models and seafloor sediment samples taken. 62 -5 -28.5 Z (ID) -52.1 -75.6 -99.2 dx = 1 km dx = 500 m dx = 100 m dx = 100 m dx = 500 m dx = 1 km \ i I \ i \ \. \ \ J A 233 238 243 X(km) Figure 4.7. Magnified views of the grid for the shoreline region of cross-section 1. The dark region in (A) can be seen at the greater magnification used in (B). X(km) -5 Z ( m ) -28.5 dx = \ = IOC i \ ) m r d) < = 5 J 3 r n ix = c  2 ix Or I = 1C m j d nr x = i = 2 O r nJ dx n = 50 I m dx = J = IOC i \ m r 1 BI 237 238 239 X(krm) Figure 4.7 continued. Magnified views of the grid for the shoreline region of cross-section 1. The dark region in (A) can be seen at the greater magnification used in (B). F i g u r e 4.9. Grid employed for cross section 2. Darker areas represent regions where the level of discretization was too fine to be viewed as separate grid lines at this scale. Figure 4.10. Elevation of the surface of Permeable Zone B. Contours are in meters. Updip limit 0 100 200 300 400 X(km) Figure 4.11. Elevation of the surface of Permeable Zone C. Contours are in meters. Updip limit 0 100 200 300 400 X(km) Figure 4.12. Elevation of the surface of Permeable Zone D. Contours are in meters. 0 100 200 300 400 X(km) Figure 4.13. Elevation of the surface of Confining Unit E. Contours are in meters. Y(km) X(km) Figure 4.14. Elevation of the surface of Permeable Zone E. Contours are in meters. Permeable Zone D 0 100 200 300 400 X(km) Figure 4.15. Outcrop pattern of the hydrologic units within the modeled region. Figure 4.16. Distribution of hydrogeologic layers for cross-section 1. Point S L denotes the location of the shoreline. X ( k m ) Figure 4.17. Distribution of hydrogeologic layers for cross section 2. Note that the depth of Lake Ponchartrain is not to scale. Point SL denotes the location of the shoreline. Hydraulic Conductivity Hydrogeologic Unit Geometric Mean of all Measurements (m/d) Geometric Mean of Measurements Within the eastern Coastal Lowlands Aquifer System (m/d) Permeable Zone A 47.5 40.2 Permeable Zone B 15.5 25.0 Permeable Zone C 15.2 25.9 Permeable Zone D 20.1 23.8 Permeable Zone E 18.0 17.1 Table 4.1. Horizontal hydraulic conductivity estimates reported by Prudic (1991). Hydrogeologic Unit Horizontal Hydraulic Conductivity (m/d) Vertical Hydraulic Conductivity (m/d) Porosity (-) Specific Storage (1/m) Longitudinal Dispersivity (m) Transverse Dispersivity (m) Diffusion Coefficient (m2/d) Seafloor Mud Layer 0.043 0.0043 0.3 1.55 x 10"' 5 0.5 8.64x10"b Permeable Zone A 43 4.3 0.3 1.55 x 10"8 5 0.5 8.64 x10"6 Permeable Zone B 15.2 1.52 0.3 1.55 x 10'9 5 0.5 8.64 x 10"7 Permeable Zone C 6.1 0.61 0.3 1.55 x 10"10 5 0.5 8.64 x 10"8 Permeable Zone D 1.52 0.152 0.3 1.55 x10 - 1 1 5 0.5 8.64 x10"9 Confining Unit E 0.043 0.0043 0.3 1.55 x 10"12 5 0.5 8.64 x10 - 1 0 Permeable Zone E 0.305 0.0305 0.3 1.55 x10 - 1 3 5 0.5 8.64 X10"11 Table 4.2. Flow and transport parameters used in the base case scenarios of the two-dimensional simulations. Parameter Value K (m/s) 2 x 10~5 d(m) 3000 A C (kg/m 3) 100 C m a x ( k g / m 3 ) 100 Pma X(kg/m 3) 1065 Po(kg/m 3) 1000 a 0.065 n 0.3 D m (m 2/s) 1 x 10" 9 AT CC) 130 P(1/°C) 4 x 10" 4 A (m 2/s) 6 x 10~7 Table 4.3. Parameter values used for calculation of the solutal and thermal Rayleigh numbers. 76 5. TWO-DIMENSIONAL SIMULATIONS OF SUBMARINE GROUNDWATER DISCHARGE IN THE GULF OF MEXICO NEAR SOUTHEASTERN LOUISIANA 5.1. INTRODUCTION Several two-dimensional simulations were performed to develop estimates of SGD within the study region. Two cross-sectional models were used, whose location can be found in Figure 4.5. Cross-section one (denoted A - A ' in Figure 4.5) was chosen to approximately follow a flow line as determined from initial three-dimensional modeling results. Cross-section two (denoted B - B' in Figure 4.5) was chosen to examine the effect that Lake Ponchartrain would have on SGD into the study region. The two-dimensional models were developed due to the extended length of time that was required to simulate the system in three dimensions. Simulation times for the three-dimensional models were on the order of months, whereas the two-dimensional models could be run to completion within four or five days. The two-dimensional models were run without any internal sources or sinks due to pumping stresses. Therefore, results obtained from this set of simulations represent predevelopment flow conditions. A series of two-dimensional simulations were performed to investigate the sensitivity of the system to the chosen set of parameters, as well as the imposed boundary conditions. The sensitivity analysis was performed using the two-dimensional models due to significantly less simulation times that were required. Extensive simulation times required for the three-dimensional model precluded its use for this purpose. The first simulation that will be described for each cross-section will be referred to as the base case scenario. These simulations used the parameter values that were described in section 4.4 and are listed in Table 4.2. The chosen set of parameters represent average values that 77 were found in the literature from previous studies or were chosen to accommodate grid constraints. The boundary conditions that were assigned to the base case simulations will be described in the next section. Changes to the parameters and boundary conditions that were made in the subsequent runs will be described in the relevant section. Each two-dimensional simulation was run until a dynamic equilibrium was attained. Simulations were deemed to have reached equilibrium once salinity fronts ceased to move and changes in the groundwater flux through the surface layer of nodes were less than one percent relative between output times. The total amount of simulation time required to reach equilibrium varied from 2 to 2.5 million days. In order to simplify the simulations presented here, the maximum salinity of groundwater was set equal to seawater. This simplification was found to reduce simulation times significantly. Due to the fact that salt domes were not incorporated into the model, it would not have been possible to simulate complex flow patterns that may be present near them. A simulation was run that included groundwater with density up to 1065 kg/m3. Because the salt domes and the associated brines were not included, flow patterns occurring at depth, particularly in the offshore region of the model, should be interpreted with care. In a field setting, flow near the salt domes will be influenced to some extent by the brines. It is not believed, however, that this simplification should significantly affect flow patterns near the surface of the aquifer system. Groundwater of elevated density is generally restricted to the deeper portions of the Coastal Lowlands aquifer system (i.e. greater than 1 km depth). Previous simulations performed by Kuiper (1994) found that removing the portion of the aquifer system where brines were present did not significantly affect groundwater flow in Permeable Zones A and B. 78 Thermally induced groundwater flow and compaction-driven flow were also not included in these simulations. Density variations due to solute gradients are much greater than thermally induced variations within the study region and are therefore probably more significant (Sharp et al, 2001). Relative to topographically-driven flow, compaction-driven flow rates are expected to be negligible. Person et al (1996) note that previous modeling results in sedimentary basins have found maximum groundwater flow rates of 1 cm/yr due to sediment compaction. In comparison to topographically-driven groundwater flow rates which can reach up to 10 m/yr, flow induced by sediment compaction appears quite small. 5.2. BOUNDARY CONDITIONS 5.2.1. CROSS-SECTION 1 The boundary conditions employed for the base case scenario are shown in Figure 5.1 A. Modifications made to the boundaries in subsequent simulations will be noted in the relevant section. For the onshore portion of the model domain, the nodes corresponding to the water table were assigned fixed values of freshwater head, corresponding to the elevations of the nodes. A constant concentration equal to freshwater was assigned to these nodes. Nodes representing the seabed were also assigned fixed freshwater heads adjusted to account for the salinity of seawater. Changes in head due to tidal fluctuations were not incorporated into the two-dimensional models. Instead, a constant sea level at an elevation of zero was assumed. Due to this simplification, the component of SGD due to 79 tidal pumping and wave-setup could not be simulated. At the regional scale of this analysis, it would not have been possible to model these small scale processes accurately. A third type exit boundary was specified as the boundary condition for solute transport across the seabed. Where groundwater was determined to be exiting the domain, nodal salinities were calculated by the model. Where the flux was into the domain, the inflowing water was assigned a concentration equal to seawater. Use of this boundary condition allows the model to compute the location of the freshwater/saltwater transition, rather than specifying it a priori using a specified concentration boundary. The northern limit of the model was assigned a constant head equal to the elevation of the topmost node, and a constant concentration equal to freshwater. The southern boundary, located at the edge of the continental shelf, was assigned a constant freshwater head calculated for the specified concentration of seawater and an assumed water level equal to mean sea level. For both of these boundaries, hydrostatic conditions were assumed. The base of the aquifer system was assumed to be impermeable to both flow and solute transport. Any flow across this boundary should be small due to the presence of the thick Vicksburg-Jackson regional confining unit or the geopressured zone. Kuiper (1994) investigated the effect that flow across the geopressured zone would have on flow in the entire Coastal Lowlands, Mississippi Embayment, and Texas Coastal Uplands aquifer systems. He found that upward flow from the geopressured zone into the overlying aquifers could be as high as 33 m/d without significantly affecting the hydraulic head distribution in his model, as determined primarily from shallow hydrogeologic units. Discharge from the geopressured zone would be expected to be much less than this value 80 in order to maintain the elevated fluid pressures that have been found. Based on the findings of Kuiper and a lack of field data available to estimate discharge rates, no simulations were performed in which flow was allowed across the base of the system. 5.2.2. CROSS-SECTION 2 A diagram of the boundaries employed for the base case scenario of cross-section 2 can be found in Figure 5.1 B. All boundary conditions were the same as for the base case scenario of cross-section 1, except in the area of Lake Ponchartrain. The seabed of the lake was assigned a third type exit boundary. Water predicted by the model to flow into the domain across the seabed was assigned a concentration equal to seawater. Nodes representing the sea bed were assigned fixed equivalent freshwater heads calculated for the salinity of seawater. 5.3. RESULTS FOR CROSS-SECTION 1 5.3.1. BASE CASE SCENARIO The parameter values and boundary conditions assigned for the simulation that will be referred to as the base case scenario are described in previous sections. Results of the simulation are plotted in Figures 5.2 and 5.3. Figure 5.2 illustrates the groundwater flow patterns and fluid density distribution of the aquifer system. Selected streamlines were plotted in order to show the dominant flow patterns within the aquifer system. The density of streamlines is not indicative of the relative magnitude of the groundwater fluxes within the system. 81 Three main flow systems can be identified in Figure 5.2. An onshore flow system extends from the northern boundary to approximately X = 130 km (at the water table). In the center of the domain, a seawater recirculation system is present. It spans approximately 130 km at the land surface, and extends 5 km offshore. Recirculating seawater is present beneath the onshore system throughout much of the coastal plain. The third system, located within the sediments on the continental shelf, is the result of density-driven convective circulation. Groundwater flow in the onshore system is controlled by throughflow across the northern boundary and recharge across the land surface. Groundwater flows across the northern boundary due to the lateral hydraulic gradient induced by the relief of the land. A volume of 0.039 m3/d per meter width crosses the boundary at an average flux rate of 0.0065 cm/d. Little recharge through the model surface is predicted near the boundary because the land surface is relatively flat and groundwater flow is essentially horizontal. Inflow across the specified concentration boundary has led to fresh groundwater occupying the entire vertical extent of the aquifer system near the boundary. The magnitude of groundwater fluxes through the upper model surface can be found in Figure 5.3. Flux values are plotted as a three point moving average to reduce the numerical oscillations in the calculated discharge rates. Some amount of oscillation is expected due to the rather coarse level of horizontal discretization that was employed. In a true physical setting, the transition from a recharge zone to a discharge zone would occur gradually. Here, the transition is forced to occur over distances up to 1 km. The result is that the plot appears to have jagged spikes as adjacent nodes are predicted to have groundwater flowing in opposite directions. 82 Freshwater recharge enters the onshore system through local topographic highs. Recharge rates reach a maximum of 0.06 cm/d, with an average of 0.014 cm/d, in the upland region. The average recharge rate corresponds to 3.3% of the average precipitation rate. The average precipitation rate for southeastern Louisiana is 1575 mm/yr, or 0.43 cm/d. Relative to this value, recharge rates to the aquifer system appear quite low. Williamson and Grubb (2001) state that the amount of precipitation recharging the Coastal Lowlands aquifer system is low in most regions, and not directly dependent upon the rate of precipitation. In most areas, aquifers are fully saturated, causing the majority of precipitation to be lost to the atmosphere through evapotranspiration or flow over the land surface as runoff. In their model results, Williamson and Grubb (2001) estimated that net regional recharge throughout their study area was only 0.6% of the total precipitation under predevelopment conditions, although recharge rates did vary significantly. In the area surrounding the upland region of cross-section 1, estimated recharge rates reached a maximum of 0.042 cm/d (Williamson and Grubb, 2001). In this context, recharge rates in this region appear reasonable, as they are almost identical to the previous modeling results. All water entering the domain in the onshore system is discharged at the northern extent of the coastal plain. Here, fresh groundwater flows up over more dense seawater. Discharge rates reach a maximum of 0.06 cm/d. This value is equal to the maximum recharge rate that occurs to the north. Therefore, discharge rates at this location are directly dependent on recharge rates in the upland region. The discharge zone is horizontally extensive, and begins where the topographically higher uplands meet the broad, flat coastal plain. No groundwater originating in the onshore system is observed to 83 flow past the discharge zone. It is expected that the water discharged in the zone provides recharge to the numerous surface water bodies that are present. Relative to the surface water inventory, the volume of water released in the discharge zone is not substantial. According to Williamson and Grubb (2001), groundwater-surface water exchange is less than one percent of the total surface water inventory, making this small quantity difficult to measure using conventional methods. Groundwater discharged through the coastal plain may also be expected to seep into the extensive marsh regions that are present. The marshes span the majority of the nearshore region adjacent to the coastline. In many areas, marshland extends several hundred metres inland from shore. The remainder of the coastal plain is occupied by a second discharge zone. The zone starts at the southern extent of the onshore system, located at approximately X = 130 km, and extends to X = 237.5 km (500 m from the shoreline). Groundwater discharged in the second zone enters the aquifer system near the coastline. The recharged seawater follows a deep flow path before discharging through the coastal plain. The assignment of the specified freshwater concentration boundary at the model surface has allowed a narrow lense of low salinity groundwater to be present directly below the land surface. Discharge values in this region are an order of magnitude less as compared to the discharge zone from the onshore flow system located to the north. The inset plot in Figure 5.3 provides a magnified view of the flux plot at the coastline (located at X = 238 km). As seen from the plot, no SGD is predicted to occur near the shoreline. A downward flow component is seen in the final 500 m of land prior to the coast. 84 A magnified view of streamlines in the vicinity of the coastline can be found in Figure 5.4. This figure shows how flow patterns differ from traditional shoreline models (see Figure 1.1, for example). In traditional models, a narrow transition zone from freshwater to seawater is expected at the mean tide line. At this location, fresh groundwater, driven by onshore hydraulic gradients, may be expected to discharge across the seabed. SGD of intermediate to seawater salinity would also occur at the seaward extent of the transition zone due to the convection of seawater. Flow patterns predicted by this simulation differ significantly from the traditional representation. From 500 m inland of the shore seaward past X = 241 km, recharge is driven by the extensive seawater recirculation system beneath the coastal plain. This system also acts as the driving force for groundwater flow to the surface of the coastal plain up to 100 km landward of shore. Low hydraulic gradients induced by the flat coastal plain are insufficient to drive groundwater discharge through the seafloor sediments adjacent to the coastline. A very gradual increase in topography is found from the shoreline to the northern edge of the coastal plain (see Figure 4.5). The elevation change of the plain is less than a meter over a distance of 13 km from shore. Shallow groundwater flow within 500 m of the shoreline is essentially horizontal, with low rates of recharge crossing the land surface up to the shore. At the shoreline, seawater enters the system and flows vertically downward through the mud layer, before becoming entrained in the convective flow pattern of the recirculation system. The results obtained here are similar to what was found by Wilson (2005) for a similar coastal setting. In simulations that were performed where the elevation of the water table did not rise above 1 m up 10 km inland from shore, no fresh SGD was 85 predicted. In her simulations, all freshwater flowing from inland regions of higher elevation was predicted to discharge prior to 10 km from the coast. A deep seawater recirculation system that was 35 km wide was adjacent to the discharge zone. The shoreline recharge zone extends 20 km offshore. Seawater entering the first 5 km of seafloor follows a deep flow path before being discharged at the northern extent of the recirculation system. The deepest flow line that is plotted in Figure 5.2 mimics the topography of the base of the Coastal Lowlands aquifer system, except where it encounters Confining Unit E. At this point, streamlines divert around the low permeability unit before returning to the base and eventually discharging through the coastal plain. Recharge at the shoreline has resulted in the formation of a finger of groundwater with salinity slightly less than seawater. The finger extends about 1 km vertically, and spans a horizontal distance of approximately 10 km. It has formed as a result of recharge water entering the domain through the constant freshwater concentration boundary between X = 237 km and X = 238 km. In a field setting, verifying whether or not it is present would be difficult because the salinity of groundwater within the finger is 90% of seawater. Typical field instruments would be able to discern between this salinity and that of seawater. However, obtaining a data set that was spatially dense enough to delineate the finger could prove difficult. A future simulation will show that the finger is not a numerical artifact created by the constant concentration boundary condition. The remaining 15 km of the near shore recharge zone contributes water to the convection system on the continental shelf. Two convection cells are present in the system, with the first extending from 5 km offshore to approximately X = 280 km. It 86 extends vertically to the base of the domain. The convection cell produces SGD at an average rate of 0.0095 cm/d. Groundwater is discharged over a horizontal extent of 22 km, yielding 2.21 m /d per meter width. The second convection cell spans the remainder of the continental shelf. It produces SGD at an average rate of 0.025 cm/d, discharging a total volume of 9.12 m /d per meter width over a distance of 36 km. Discharge due to this convection cell is greater because a wider zone of recharge is contributing water to the flow system. To place the volume of groundwater discharged by the convection cells in the base simulation into context, the total volume of SGD may be computed across the entire study area as a percentage of Mississippi River discharge. If the rate of groundwater discharge is assumed constant across the study area that is approximately 90 km wide, then SGD is predicted to be only 0.1% of river discharge (assuming an average annual river discharge of 380 knrVyr). Therefore, relative to the magnitude of river discharge, SGD is predicted to be relatively insignificant. As noted in the previous chapter, the Rayleigh number that was calculated (RaN = 1.5 x 107) suggests that unstable (transient) convection cells may develop within the aquifer system. According to Diersch and Kolditz (2002), unstable convection will occur where RaN is greater than 240 - 300, although numerical simulations may converge to a single solution that is generally dependent on the imposed initial conditions. In order to check that the observed flow patterns did not change with time, the base case simulation was run for an additional 2.5 million days. Results of the simulation showed that the convective flow patterns did not change. Another simulation was run where the hydraulic conductivity of all hydrogeologic layers was decreased by four 87 orders of magnitude. In this simulation, the Ran would be decreased to 1.5 x 103. The convective flow patterns were not significantly altered in this simulation either, although the magnitude of the fluid flux was reduced correspondingly. Also, the simulation that was run using different initial conditions (see section 4.8) generated identical flow patterns as compared to the base case. These results suggest that a stable convective regime that is not dependent on the numerical formulation of the model has developed within the aquifer system. Several authors have noted that the Ran may not be a good indicator of instability (e.g. Schincariol et al, 1994; Simmons et al, 1999). Also, the assumptions that were built into RaN limit its applicability to systems with complex geometries such as the one described here. Wilson (2005) also notes that there is no Rayleigh number threshold for convection where dipping isopycnals (constant density surfaces) exist, which is the case here. The convective flow patterns must still be interpreted with caution. The cross-sectional model results are a two-dimensional representation of a three-dimensional process. The three-dimensional model will help to show how convection patterns may be altered in three-dimensional space. Some of the recharged water exits the domain through the southern continental shelf boundary. Williamson and Grubb (2001) note that near the edge of the continental shelf, the sediments in Permeable Zones A and B undergo a facies change to less permeable marine clays. Therefore, assigning the elements located next to the boundary low permeability values, or making the boundary entirely impermeable may provide a more realistic representation of the boundary condition. A simulation was performed where this boundary was impermeable. Results of this simulation showed that water 88 formerly crossing the boundary was forced to discharge through the seafloor, increasing the SGD due to the convection cell. No other differences from the base case were present; therefore graphical results of the simulation are not presented. Results of the base case simulation show that little SGD is expected along the extent of cross-section one. Groundwater driven by freshwater recharge and lateral hydraulic gradients in the upland region is predicted to discharge in a region over 100 km north of the coastline. From 500 m landward of the coastline out to X = 258 km, a large recharge zone is predicted by the model. No SGD is present near the shore. Water entering the subsurface in this zone becomes entrained in either the seawater recirculation system, or the offshore convection system. The recirculation system drives groundwater discharge in the coastal plain, in a zone that extends to within 500 m of the coastline. The goal of the sensitivity study described below was to determine if there are any situations under which SGD at the shoreline may be expected, and how certain boundary conditions and parameter values affect the simulation results. 5.3.1. A . M O D E L V E R I F I C A T I O N Data available to test the model results are salinity values provided in the report by Williamson and Grubb (2001) and salinity values obtained from shallow groundwater samples collected as part of the tracer study. A diagram of the salinity distribution (converted to fluid density) reported by Williamson and Grubb (2001) near the location of cross-section 1 is provided in Figure 5.5. The contour lines on Figure 5.5 were manually placed by interpolating between discrete measurement points. The measurement points spanned the vertical extent of the aquifer system. Williamson and Grubb (2001) compiled an extensive data set as part of their study. However, only sparse 89 data is available for the shoreline and offshore portion of their study region. They focused primarily on describing onshore flow patterns for the purpose of groundwater management issues. As seen from Figure 5.5, the model results match the measured salinity profile quite well. The locations of contour lines with values less than seawater appear similar to those obtained from the base case simulation. The thickness of fresh-water bearing sediments is somewhat greater in Figure 5.5, suggesting that the model has underestimated the thickness of this zone. This may be due to the coarse vertical node spacing that was used. At 5 m below the surface of the land, nodal spacing increased to around 30 m. The actual value was dependent on the horizontal position of the node. Employing decreased nodal spacing may have allowed the model to simulate the depth of freshwater-bearing sediments more accurately. Also, it is not possible to discern whether or not the finger of intermediate salinity is present at the shoreline. Figure 5.5 also shows that nearly two thirds of the Coastal Lowlands aquifer system is occupied by groundwater with density greater than seawater. The density distribution within this region could not be contoured because salinity values vary greatly over short horizontal and vertical distances. Instead, a number of point values have been plotted to demonstrate the complexity of the salinity distribution, which reflects the presence of salt domes. Because the maximum fluid density was set to that of seawater, flow processes that may be occurring at depth cannot be fully understood using the simulations described here. Therefore, the offshore convection system should only be interpreted in a semi-quantitative sense, as the salt domes may exert an influence that is not included in the 90 models. Also, thermal effects and compaction-driven flow may become more significant in the deeper regions of the aquifer system. However, it is believed that discharge through the upper boundary of the model domain should not be significantly affected by this simplification. Previous simulations performed by Kuiper (1994) tested the effects that the brine part of the system had on flow patterns within the Gulf Coast regional aquifer system. He found that the system was not very sensitive to two types of model truncation.The first truncation was at the shallowest depth where water density was 1005 kg/m3.The second was at a fixed depth of 915 m, which effectively removed all of the brines from the flow system. A no flow boundary was placed at the base of the truncations. In each case, results from the overall model did not change significantly. Differences occurred near the bottom no flow boundary, where, as expected, flow patterns were substantially different. Changes in near-surface groundwater flow were minimal, suggesting that the brines had little effect in this region. Figure 5.6 shows the location of wells that were sampled within the study region. All of the wells were completed within Permeable Zone A. The depth of each well, as well as average fluid density values that were measured are provided in Table 5.1. Solute concentration and fluid density measurements from each sampling period can be found in Table A3.1 in Appendix 3. Wells M W 4 through 9, as well as M W 1, 13, and 14 lie directly on or are near the location of cross-section 1. Salinity data from these wells shows that groundwater with density near that of freshwater is present up to a depth of at least 130 m beneath the coastal plain. The well data reinforce the fact that the model has underestimated the depth of freshwater-bearing sediments. However, the large-scale structure of the base case scenario matches the measured salinity values rather well. 91 Although no measurements of hydraulic head were taken at the monitoring wells, it should be noted that the majority were found to be artesian (i.e. flowing). The piezometric surface at these locations was found to be up to 1 m above ground surface, confirming the fact that the coastal plain is serving as a discharge zone for groundwater. Whether or not all of the wells that were sampled are artesian could not be confirmed due to pump completion. These wells may be expected to be flowing as well. A list of the wells that were found to be artesian is provided in Table 5.1. 5.3.2. HYDRAULIC CONDUCTIVITY Data files provided with the report by Williamson and Grubb (2001) indicate that the hydraulic conductivity of the hydrogeologic layers varies significantly within the study region. Values were found to vary over two orders of magnitude for all layers. Three simulations were performed to evaluate the role that the chosen permeability values had on the model-predicted SGD. Differences in hydraulic conductivity values that were used in these simulations were as follows: • Simulation 1: The hydraulic conductivity of Permeable Zones A, B, and C was increased by a factor of six. • Simulation 2: The seafloor mud unit was removed from the first 15 km of the seabed. • Simulation 3: Permeable Zone A was assigned a low permeability value. Results of these simulations will be discussed in the following sections. 92 5.3.2. A . SIMULATION 1: INCREASED H Y D R A U L I C C O N D U C T I V I T Y It was hypothesized that increasing the permeability of the top three hydrogeological layers may result in near shore SGD. The presence of more conductive surficial units would increase groundwater recharge rates for the given set of boundary conditions which would lead to a corresponding increase in horizontal groundwater flow velocities near the model surface. This may allow the onshore flow system to displace the recirculation system further seaward. Discharge due to the recirculation system would then be pushed further seaward, potentially resulting in SGD being predicted through the seafloor in the shoreline region. Simulation results can be found in Figures 5.7 and 5.8. Viewing the streamline plot in Figure 5.7, it is apparent that the onshore system did not displace the recirculation system further seaward. Recharge to the uplands has increased relative to the base case scenario, but the onshore system still occupies the same horizontal extent. Recharge fluxes reach 0.35 cm/d (see Figure 5.8), which is less than the average precipitation rate of the region. The average recharge rate is 0.08 cm/d, or 19% of the average precipitation rate. This value is higher than the recharge rate that was predicted in the U.S. Geological Survey's regional scale model (Williamson and Grubb, 2001). Flux rates in the discharge zone of the onshore system have increased as well. The discharge zone occupies the same horizontal extent as in the base case scenario. Increased permeability of the surficial units has altered flow patterns in the near shore part of the domain as well. Recharge through near shore sediments has increased. The recharge zone for this simulation extends 1.5 km north of the shoreline, with maximum flux rates reaching almost 0.78 cm/d at the shoreline. Increased recharge 93 through onshore sediments has allowed the finger of intermediate salinity groundwater to extend to a depth of 2 km. The depth of the finger does not match well with previously measured salinity values (see Figure 5.5). The depth of the finger suggests that the permeability of the system near the shoreline is probably too high in this simulation. Flow patterns on the continental shelf are also different in this simulation. Discharge from the convection cell near the shoreline occurs through a region that is only 8 km wide. As a result, the amount of SGD that is produced by the cell has been reduced to 0.21 m /d per meter width, with an average flux rate of 0.0025 cm/d. Convection occurring further offshore adjacent to the continental shelf results in more SGD in this simulation. A total volume of 31 m /d is discharged over 43 km in this region, with average flux rates of 0.046 cm/d. The volume of water discharged by this convection cell has increased because a wider offshore recharge zone is contributing water to the cell. The hydraulic conductivity values assigned to Permeable Zones A, B, and C in this simulation are probably maximum values that may be expected. Maximum recharge rates in the upland region are found to approach the average amount of precipitation in southeastern Louisiana, but are greater than what has been predicted in previous modeling work by Williamson and Grubb (2001). The depth of the finger of intermediate salinity is too great relative to measured salinity values (see Figure 5.5). Therefore, the permeability of the surficial units in the shoreline region is probably too high. However, simulation results do show that increasing the permeability of the sediments in the surficial layers does not significantly alter the dominant flow patterns. No SGD is predicted near the shoreline. Groundwater discharge is expected further offshore due to 94 convection processes in the shelf sediments. The value and location of the discharge is dependent on the permeability of the sediments. Increasing the hydraulic conductivity of the three upper hydrogeologic units further was found to increase maximum recharge rates to values greater than the average precipitation rate. Another simulation was run where the hydraulic conductivity of Permeable Zones A, B, and C was increased by a full order of magnitude. Maximum recharge rates through the upland region in this simulation were found to be 0.6 cm/d, a value greater than the average precipitation rate. Therefore, the order of magnitude increase in hydraulic conductivity probably does not provide an accurate representation of the flow system. The results were similar to those for the case where hydraulic conductivity was increased by a factor of 6, with no SGD predicted to occur at the shoreline. 5.3.2. B. SIMULATION 2: R E M O V A L O F T H E N E A R S H O R E M U D L A Y E R This simulation investigated the effect that the low permeability mud layer has on the model-derived SGD estimates near the shoreline. It has been suggested that higher permeability sand beds may extend from the onshore region of the study area out some distance into the mouths of the coastal bays (Meade Allison, personal communication, September 2003). The presence of a more conductive surficial unit into the Gulf may result in SGD entering the bay near the shoreline. Cross-section one passes through Terrebonne Bay. Along with Barataria Bay, it may be expected that sand beds are present some distance offshore. To simulate this situation, the low permeability mud layer at the surface was removed from the first 15 km 95 of the seafloor. Elements that were formerly part of the mud layer were assigned to Permeable Zone A. This allowed the seafloor sand beds to extend out to the mouth of Terrebonne Bay. Sediment samples were not collected within the bay, therefore the presence of the offshore sand beds could not be confirmed. Seafloor sediment samples that were taken just seaward of Terrebonne Bay (see Figure 4.5) were found to contain only fine-grained material, therefore the sand beds were not extended any further into the Gulf. Plots of the streamlines and fluxes through the upper boundary of the model domain are provided in Figures 5.9 and 5.10 respectively. The model results show that SGD is not predicted near the shoreline. Removal of the mud layer has resulted in increased flux rates in the near shore recharge zone, although it only extends 10 m landward of the mean tide line. Recharge rates reach a maximum of 1.4 cm/d a distance of 20 m seaward from shore. The influence of the increased recharge can be seen in both the seawater recirculation system and the near shore convection cell. Discharge through the coastal plain near the shoreline has increased. The majority of the increased recharge that enters the recirculation system follows a shallow flow path and discharges near the shoreline. Flux rates reach 0.4 cm/d a distance of 20 m landward from shore. Discharge due to the near shore convection cell has also increased relative to the base case scenario. A total volume of 4.31 m3/d of saline groundwater is discharged over a horizontal distance of 28 km (from X = 252 km to X = 279 km). A peak in the discharge rate is present on the northern side of the convection cell. Here, the recharged water follows a shallow flow path through the seafloor sediments and discharges at a 96 maximum rate of 0.03 cm/d. The average discharge rate of the convection cell is 0.015 cm/d. The volume of water discharged by the offshore convection cell remained unchanged. Simulation results show that no SGD is predicted near the shoreline despite the absence of the mud layer on the first 15 km of seafloor. The more permeable seafloor has allowed a greater volume of saltwater recharge to enter the system. Increased recharge has resulted in more discharge through the coastal plain near the shoreline. 5.3.2. C. SIMULATION 3: LOW PERMEABILITY SURFACE UNIT Another scenario that could be invoked to force groundwater discharge through the near shore sediments would be the presence of a low permeability unit at the surface that extends offshore. In this situation, groundwater flowing towards the coast through Permeable Zone B would be expected to seep through the overlying confining unit at decreased flux rates. In order for a similar quantity of water to be discharged (relative to the base case scenario), the discharge zone would be forced to extend further seaward, possibly reaching offshore. Field sampling of groundwater wells in the study region (see Figure 5.6) suggested that small confining units must be present near the land surface. Several of the wells were found to have natural artesian flow (see Table 5.1), which could only occur if a low permeability layer was present to restrict vertical groundwater movement. To examine this scenario, Permeable Zone A was assigned a hydraulic conductivity equal to that of the seafloor mud layer (i.e. K x = 0.043 m/d; K z = 0.0043 m/d). This scenario represents a physically unrealistic situation because the permeable 97 zone is known to be a conductive unit throughout the state of Louisiana. The majority of pumping wells in use in the State are completed within Permeable Zone A (Mesko et al, 1990). Any confining unit present, therefore, is not regionally extensive. No information was found that could be used to delineate the confining unit, therefore assigning all of Permeable Zone A to have low permeability was the simplest method to simulate this situation. Plots of the results are provided in Figures 5.11 and 5.12. From Figure 5.12, it is seen that flow patterns through the surface of the model are significantly altered. Recharge to the onshore flow system has been reduced from the updip limit of Permeable Zone A to the discharge zone at the northern extent of the coastal plain. Recharge flux rates are less than 10% of the average precipitation rate in the outcrop region of Permeable Zone B. Flux rates in the discharge zone of the onshore flow system are also reduced due to the presence of the low permeability surface layer. The freshwater discharge zone is calculated to extend from X = 77 km to X = 110 km. Discharge is more uniformly distributed within the zone, with an average flux of 0.006 cm/d. In Figure 5.11, it is seen that the reduced amount of recharge in the upland region has allowed the recirculation system to extend 20 km further inland. From X = 110 km to X = 230 km, a low rate of discharge is calculated for the coastal plain. Flux rates are less relative to the base case scenario because inflow rates in the near shore recharge zone have been reduced. Recharge to the seawater recirculation system extends 8 km inland from the coastline. SGD is again not predicted near the shoreline. 98 Convection patterns in the offshore region of the model also differ from the base case. The convection cell located adjacent to the shoreline is located deeper in the aquifer system within the more conductive sediments of Permeable Zone B (see Figure 5.12). Slow seepage of saline groundwater through the overlying confining unit does not occur. Instead, seawater percolates through the seabed and becomes entrained in the convection cell. As a result, no seabed discharge is produced by the cell. The convection cell located further offshore produces a diminished amount of SGD due to the thicker surficial confining unit. 0.73 m /d per meter width is discharged across the seafloor over a horizontal distance of 39 km, with an average flux rate of 0.0019 cm/d. Specifying Permeable Zone A to have a low hydraulic conductivity has not caused SGD to occur beyond the shoreline. Reduced recharge in the upland region of the model allowed the recirculation flow system to move further inland. Recharge at the shoreline has been reduced as well, although the near shore recharge zone has expanded further inland. The result is that less discharge is predicted in the two discharge zones on the coastal plain. 99 5.3.3. DISPERSIVITY The dispersivity values used in the model simulations were chosen based on data published by Gelhar et al (1992) and to accommodate grid constraints. Therefore, a simulation was performed to investigate the effect this parameter had on the model-derived SGD estimates. In this simulation, both the longitudinal and transverse dispersivities were increased by an order of magnitude to 50 m and 5 m, respectively. This maintained the ratio of ai: a t at 10:1. The values of dispersivity used in this simulation were also employed for the three-dimensional model to accommodate the larger node spacing necessary in that analysis. Therefore, results of this simulation may be used to compare the three-dimensional model and to elucidate any effects the higher dispersivity values may have on the results. Streamlines and fluid density values obtained from the simulation are plotted in Figure 5.13. Comparing the density contours with the base case scenario, it is seen that the width of the mixing zone between the freshwater and seawater endmembers has increased due to a greater amount of dispersion. The depth and extent of freshwater-bearing sediments appears essentially unchanged. The deeper mixing zone has increased the depth at which groundwater of seawater concentration is encountered. Greater dispersion has also reduced the depth of the finger of intermediate salinity and made it more rounded. At the shoreline, a wider mixing zone between fresh and saline groundwater is present due to increased dispersion. This has allowed onshore discharge driven by the seawater recirculation system to extend to within 300 m of the shoreline (see Figure 100 5.14). The recharge zone for the recirculation system extends 16 km offshore. Recharge rates at the shoreline have increased somewhat, as the system now has to move the same volume of water (relative to the base case) through a constricted zone. Groundwater flow in the offshore convection system has not changed significantly. The volume of SGD from the near shore convection cell has increased slightly to 2.54 m3/d per meter width due to increased recharge rates at the shoreline. The dispersivity values do not affect flow in the offshore region because groundwater has a uniform salinity. Discharge rates from the convection cell located further offshore are identical to the base case scenario. 5.3.4. O N S H O R E B O U N D A R Y CONDITION The boundary condition applied to the onshore portion of the model surface was changed for this simulation in order to evaluate how the specified freshwater concentration boundary was affecting shallow flow patterns. In previous simulations, a thin layer of freshwater was present beneath the coastal plain, despite the fact that saline groundwater is flowing towards the surface. The layer has persisted due to the imposed freshwater concentration boundary. For this simulation, all onshore nodes were assigned the same third-type exit boundary as the seafloor. For nodes that were predicted by the model to receive recharge, the concentration applied to the inflowing water was equal to freshwater. Where there was upward flow of recirculated seawater to the water table, it is expected that the depth of the layer of freshwater will be decreased, if not entirely removed. 101 From the plot of streamlines and the density distribution provided in Figure 5.15, it is seen that the layer of fresh groundwater is no longer present beneath the coastal plain. The upward flux of saline waters has increased the density of the groundwater that is present. Density values in the nodes located at the surface of the plain vary from 1005 kg/m at the northern extent of the seawater recirculation system to 1017.5 kg/m at the shoreline. In comparison to salinity values obtained from near surface groundwater samples (see Figure 5.6 for well locations) the model-predicted density values beneath the coastal plain are too high. In all wells that were sampled and in all sampling periods, the maximum groundwater density that was measured was 1006.7 kg/m . This value was recorded at M W 10, which is located 25 km inland from the coast. Other wells located near the cross-section had measured density values that were less than 1005 kg/m . Use of the first type concentration boundary allowed the model to calculate a more accurate density distribution beneath the coastal plain. This may be due to the relatively large transverse dispersivity value that was used, or the coarser vertical discretization that was employed below a depth of 5 m. Either of these factors could have resulted in the model calculating fluid densities beneath the plain that are too high. Differences in the fluxes through the model surface relative to the base case scenario can be found near the shoreline (see Figure 5.16). The concentration gradient at this location has been reduced, which has decreased the vertical gradient in equivalent freshwater head. The result is a diminished rate of recharge at the shoreline due to the seawater recirculation system. The maximum recharge rate is only 0.009 cm/d in this scenario. The zone of recharge has extended further inland, with downward flow 102 predicted up to 2 km inland from the coast. This has lead to a reduced volume of water being discharged through the coastal plain. Recharge through the onshore third type concentration nodes has allowed the finger of intermediate salinity to persist, suggesting that it did not form as a result of the imposed first type specified concentration boundary condition. The portion of the recharge zone contributing water to the near shore convection cell has moved further inland as well. This has increased the volume of water entering the convection system. The result is a slight increase in the volume of SGD to 2.66 m3/d per meter width produced by the near shore convection cell. Discharge rates due to the convection cell located further offshore have been unaffected. Results of this simulation show that a third type exit boundary condition cannot be used for the onshore nodes to accurately model the density distribution beneath the coastal plain with the given set of parameters. Upward flow of saline waters driven by the seawater reciruclation system has resulted in model-predicted density values being too high in this region. Field measurements from shallow groundwater wells have shown that fresh groundwater is present beneath the coastal plain. Changes in flux rates relative to the base case scenario are present at the shoreline, where reduced concentration gradients have decreased the rate at which water enters the domain. 103 5.3.5. SHALLOW FLOW SYSTEM Nearly two thirds of the Coastal Lowlands aquifer system in the study region contains groundwater with density greater than seawater. Brines are present at depth, occupying the majority of pore space in Permeable Zones C, D, and E, and all of Confining Unit E. Flow processes in the brine region cannot be modeled accurately without the inclusion of salt domes, which will add solute to the groundwater through dissolution and act as a barrier to the movement of fresher groundwater. Due to the fact that salt domes were not incorporated into the model, another simulation was run that included only Permeable Zones A and B. These two units primarily contain groundwater with density less than or equal to seawater, although the base of Permeable Zone B does contain fluid of elevated density. For this simulation, Permeable Zones C, D, and E and Confining Unit E were specified to be inactive. Inflow at the northern boundary will be removed from the flow system, creating an impermeable boundary along the bottom of Permeable Zone B throughout the domain. In the base case, the sum of inflow across the northern boundary and recharge through the outcrop area of Permeable Zone C represented only 0.64% of the total volume of water entering the onshore system, therefore it is not expected that removing this source of water will significantly affect the onshore system. All other boundary conditions remained the same as in the base case. Results of the simulation are plotted in Figures 5.17 (note the change in scale) and 5.18. The plot of fluxes through the model surface shows that both the recharge and discharge rates due to the onshore system have decreased relative to the base case, although the patterns are the same. This has occurred because a smaller volume of porous 104 material is available to conduct groundwater. The seawater recirculation system has not moved further inland because recharge rates at the shoreline have decreased as well. A greater proportion of the recharged seawater becomes entrained in the recirculation system, allowing discharges rates through the coastal plain to remain similar to the base case scenario. Flow in the system is forced to follow a shallower path, causing discharge to continue to within 200 m from shore. Removal of the deeper hydrogeologic units has the greatest impact on flux rates in the offshore convection system. Recharge rates through the entire seafloor portion of the model have decreased because deep convection patterns are now truncated at the base of Permeable Zone B. Two smaller convection cells evolve adjacent to the shoreline. They result in 1.86 m3/d per meter width of SGD, or a 16% reduction. SGD from the offshore convection cell has been reduced by 37% to 5.77 m3/d per meter width. Removing the deeper hydrogeological layers of the Coastal Lowlands aquifer system has had little effect on discharge rates through the onshore portion of the study region. More of the water recharged at the shoreline enters the seawater recirculation system, maintaining the same discharge rate through much of the coastal plain. Discharge rates through the offshore portion of the model have been affected as well. Offshore convection cells cannot extend as deep, and produce lesser amounts of SGD to the Gulf of Mexico. Future modeling efforts should focus on characterizing flow processes occurring in the deeper parts of the aquifer system, as they may indeed exert influence on discharge patterns through the seafloor. 105 5.3.6. SEAFLOOR BOUNDARY CONDITION Langevin (2001, 2003) discussed the sensitivity of his model-derived SGD results to the boundary condition applied to the seafloor. He found that the calculated volume of fresh SGD entering Biscayne Bay, Florida, was approximately 55% less when using a first type constant concentration boundary condition for the seafloor in comparison to a third type boundary because the model-calculated salinity beneath the bay was less. The reason for the difference in results is due to the way salt enters the system. With a first type boundary, salt may enter the system through both advection and dispersion.' Utilizing a third type boundary allows salt to enter only through advection. The latter scenario offers a more physically correct representation, as the concentration gradient across the boundary would be zero except in the region of the transition zone between fresh and saltwater. Also, the third type boundary does not force the freshwater/saltwater interface to be located at the shoreline. It allows the model to determine the location. A simulation was performed to determine the influence the seafloor boundary condition had on the groundwater system modeled here. All seafloor nodes were assigned a first type concentration boundary. Simulation results are provided in Figures 5.19 and 5.20. The plots show that the first type boundary condition has caused only small changes in the model results. Differences from the base case scenario may be found at the shoreline. Here, the thickness of the transition zone from freshwater to saltwater has decreased as greater quantities of solute have entered the system in this simulation. This has steepened the concentration gradient below the last onshore node, resulting in greater recharge rates at the shoreline. Discharge through the coastal plain is still predicted up to a distance of 106 500 m from shore. Increased recharge through the specified freshwater concentration nodes has broadened the extent of the finger of intermediate salinity. Most of the extra recharge water becomes entrained in the convection cell just offshore. SGD from the convection cell has increased to 2.72 m /d per meter width. The volume of water discharged by the convection cell located next to the continental shelf remained unchanged. Results of this simulation are consistent with the findings of Langevin (2001, 2003). The first type boundary has allowed more salt to enter the aquifer system, and decreased the thickness of freshwater-bearing sediments near the shore. No SGD is predicted through the near shore sediments, so a direct numerical comparison cannot be made. 5.3.7. BRINES Numerous salt domes are present within the study region. Dissolution of salt from the domes has resulted in nearly two thirds of the Coastal Lowlands aquifer system containing brines (see Figure 5.5). Brines are defined as water containing a dissolved solids concentration greater than 35 000 mg/1, which is equivalent to the density of seawater. The goal of the next simulation was to examine the influence that high density fluids may have on the flow system. It was not expected that the model would be able to reproduce the density profile found in Figure 5.5. Density inversions could not be simulated without explicitly incorporating salt domes into the flow system. However, not enough information on the location of the domes was available to input them into the 107 model. The simulation was run only to determine if the presence of brines would result in major changes to the flow system. For this simulation, the maximum fluid density was set to 1065 kg/m . This value was chosen as a representative average of fluid density in the more saline portion of the system. As noted earlier, density values greater than 1065 kg/m3 were found only at isolated locations, presumably due to the presence of salt domes. The initial salinity distribution was modified to include higher density fluids (see section 4.7). A source of salt was not added to the system; therefore, the distribution of the brines at steady state was determined by the imposed initial condition, and its interaction with the flow system. Different initial conditions required that the southern continental shelf boundary be modified. For this boundary, specified concentrations were assigned which varied from 1025 kg/m3 to 1035 kg/m3. Freshwater heads were then recalculated to maintain hydrostatic conditions given the assigned concentrations. Preliminary simulations which incorporated brines showed that numerical oscillations were present near the base of the system where node spacing was greatest. To eliminate the oscillations, dispersivity values were increased an order of magnitude from the base case scenario. Results from the simulation are plotted in Figures 5.21 and 5.22. Looking at the salinity distribution within the aquifer system, it is seen that the thickness of freshwater-bearing sediments has not changed. Brines are only present in the lower part of the domain in the offshore region of the flow system. As a result, the recharge/discharge patterns in the onshore system have not been altered. 108 Differences from the base case scenario are present in the recirculation system and the offshore convection system. Recharge rates at the shoreline have increased due to more vigorous convection that is occurring as a result of the brines. The layers of fluid with elevated density at depth have introduced a concentration gradient at depth that was not present in previous simulations. Most of the increased recharge enters the recirculation system, resulting in more water being discharged through the coastal plain near the shoreline. Discharge extends to within 300 m of the shoreline, similar to the scenario described above where dispersivity values were increased an order of magnitude. SGD from the offshore convection cells has decreased in this scenario. The volumes of water discharged are 1.00 m3/d per meter width and 4.27 m 3/d per meter width for the near shore and offshore convection cells, respectively. Discharge across the seabed has decreased because groundwater flow paths tend to follow the density contours and discharge through the continental shelf. Upon comparing the model predicted density distribution with Figure 5.5, a relatively poor match is found in the offshore region of the model. Vertically downward flow driven by the seawater recirculation system and offshore convection cell has pushed the brines deeper in the aquifer system than has been measured. This result does not necessarily indicate that the downward flow is not occurring in the field, though. Dissolution of salt from any salt domes that may be present could maintain the elevated density values that have been found. Future modeling efforts may be able to determine whether or not this is occurring. In any case, the results of this simulation are not sufficiently different to warrant the inclusion of brines in future simulations. The 109 dominant flow patterns present within the aquifer system are the same as for scenarios where the maximum fluid density was set to seawater. 5.4. RESULTS FOR CROSS-SECTION 2 5.4.1. BASE CASE SCENARIO Cross-section two is located in the region where groundwater flow paths to the Gulf of Mexico are potentially intercepted by discharge to Lake Ponchartrain (located from X = 140 km to X = 162 km). Flow patterns occurring in this section are similar to cross-section one (see Figure 5.23). An onshore flow system extends from the northern boundary to approximately X = 125 km. Flow in this region is controlled by the topography of the land surface. 0.031 m3/d per meter width of groundwater enters the domain across the northern boundary with an average flux of 0.0086 cm/d, driven by lateral hydraulic gradients. This value is slightly less than the volume that was found to enter the system through the northern boundary in cross-section 1 because the depth of the aquifer system is less. Aquifer recharge enters the system through local topographic highs. The average rate of recharge is 0.025 cm/d, with maximum values just over 0.07 cm/d. The average value is 5.8% of the average precipitation rate in the region. These values are similar to what was found in the upland region of cross-section one. Figure 5.24 depicts the model calculated fluxes at the water table. Inflow into the Coastal Lowlands aquifer system across the specified concentration boundaries located to the north and at the surface has resulted in freshwater occupying the entire vertical extent at the boundary. 110 A large discharge zone extends from X = 100 km to approximately X = 125 km. The location of the discharge zone corresponds to a major break in slope, similar to cross-section one. In this case, the slope break occurs near the northern shore of Lake Ponchartrain, rather than at the coastal plain (see Figure 4.5). The elevation at the surface continues to decline towards the lake, with discharge rates to approaching 0.1 cm/d. Fresh groundwater at the base of the Coastal Lowlands aquifer system flows up over the denser saline wedge and discharges at the southern extent of the system. Similar to cross-section one, all groundwater entering the domain in the onshore system is predicted to discharge prior to the shoreline of the lake. Flow in the central portion of the model region is controlled by a horizontally and vertically extensive seawater recirculation system. Recharge waters enter the sediments at the shoreline to the Gulf of Mexico and flow north to the coastal plain. The deepest flow lines reach the base of the aquifer system, traveling under Lake Ponchartrain before discharging near its northern shore. A shallow recirculation system is also present beneath the lake itself. Water recharging near the southern shore of Lake Ponchartrain flows to the north and discharges near its northern shore. The maximum flux rate on the northern shore occurs 20 m inland from the shoreline with a value of 0.49 cm/d (see Figure 5.25 A). SGD is predicted in Lake Ponchartrain due to the shallow recirculation of seawater through the lake bed. A volume of 0.25 m3/d is predicted by the model to seep through the floor of the lake over a horizontal extent of 9.4 km. Discharge rates vary over an order of magnitude from 1.3 x 10"4 cm/d to 3.8 x 10"3 cm/d, with the maximum flux rate occurring at the center of the lake. Lake Ponchartrain lies outside the study region of 111 the tracer portion of the project. Therefore the rate at which SGD may be occurring cannot be confirmed. Also, the rate at which groundwater is observed to be seeping through the lake floor is probably too small to be confirmed by tracer techniques or through the use of seepage meters. Groundwater discharge driven by the large seawater recirculation system occurs through that part of the coastal plain located between Lake Ponchartrain and the Gulf of Mexico. It is calculated to continue up until 20 m from the Gulf of Mexico shoreline (see Figure 5.24 B). A thin layer of fresh groundwater is present beneath the plain as a result of the imposed freshwater concentration boundary condition. In comparison to field data that was collected in groundwater wells near the location of the cross-section, predicted density values are again too high beneath the coastal plain. Data from M W 2 and M W 3 show that relatively freshwater may be present up to a depth of 125 m below ground surface. As in cross-section 1, no SGD is predicted near the shoreline. A zone of recharge extends from the shoreline out to 15 km offshore, supplying water to both the seawater recirculation system and an offshore convection system. The maximum recharge rate at the shoreline appears quite high, reaching a value of 0.98 cm/d. In the offshore portion of this cross-sectional model, the offshore convection system is composed of three cells. Each cell extends to the base of the aquifer system, and spans an approximate distance of 30 km. Flux rates from the first two convection cells offshore 0.02 cm/d producing 6.32 m3/d of SGD per meter width. Offshore recharge rates are significantly less than at the shoreline, therefore SGD from the convection cell further offshore is significantly less at 0.44 m3/d per meter width. If the volume of SGD predicted in this simulation is assumed constant across the study area, than groundwater 112 discharge into the study region may be expected to make up to be only 0.06% of Mississippi River discharge. This percentage is even smaller than what was predicted for cross-section 1. The results obtained from the base case scenario of cross-section two are similar to cross-section one. No SGD is predicted near the shoreline of the Gulf of Mexico. Topographically-driven groundwater flow in the northern portion of the model discharges at a major break in slope near the northern shore of Lake Ponchartrain. From this point to the Gulf of Mexico shoreline, flow is controlled by an extensive seawater recirculation system. Recharged seawater enters the system at the shoreline, discharging through the coastal plain and near the northern shore of the lake. A shallow recirculation system present beneath Lake Ponchartrain results in slow seepage of groundwater through the lake bed. Flow in the offshore part of the Coastal Lowlands aquifer system is driven by convection processes. The offshore convection cells result in groundwater discharge to the Gulf. 5.4.2. DISPERSIVITY A second simulation was performed for cross-section two in which the dispersivity values were increased an order of magnitude. The purpose of the simulation was to provide a means of comparison between the two- and three-dimensional simulations. The streamlines and density distribution obtained from the simulation can be found in Figure 5.26. A plot of the fluxes through the surface of the model, and 113 magnified views of the plot at the shorelines of Lake Ponchartrain and the Gulf of Mexico are provided in Figures 5.27 and 5.28 respectively. Few differences are present between the results of this simulation and the base case scenario. Increasing the dispersivity values has caused the transition zones from freshwater to seawater to be more spread out. The effect of the increased amount of dispersion can be found at the shorelines of both Lake Ponchartrain and the Gulf of Mexico, where discharge and recharge rates have generally decreased. Wider transition zones from freshwater to seawater are present at these locations, decreasing fluid density gradients. Decreased gradients have reduced the driving force for downward flow, causing recharge rates to decline. This has resulted in the rate of groundwater discharge near the shorelines to decrease, as less water is supplied to the surface. An important difference to note between this simulation and the base case scenario is that no SGD is observed in Lake Ponchartrain. Discharge through that part of the coastal plain between the Gulf of Mexico and the lake extends further north in this simulation. As a result, the recharge zone of the shallow seawater recirculation system beneath Lake Ponchartrain is located further north and occupies the entire lake. The average recharge rate through the lake bed is 0.02 cm/d. Results of this simulation show that increasing the dispersivity values has little effect on the flow system. The main changes in the results can be found at the shorelines of Lake Ponchartrain and the Gulf of Mexico. At these locations, recharge rates have generally declined due to increased amounts of dispersion. Also, SGD is no longer predicted in Lake Ponchartrain. 114 5.5. DIMENSIONLESS PARAMETERS 5.5.1. PECLET NUMBER The mesh Peclet number (Pe), can be defined along a flow line as: P e = v | Z | v | L (5.1) D {Dm+a\v\) where L is the distance between nodes (L), |v| is the magnitude of the local velocity vector (L/T), D is the dispersion coefficient (L 2 /T) , D m is the molecular diffusion coefficient (L 2/T), and a is the dispersivity (L). The Peclet number is a measure of advective transport relative to the amount of diffusive and dispersive transport that occurs between nodes. Assuming that transport due to molecular diffusion is significantly less than that due to mechanical dispersion, equation (5.1) reduces to: Voss and Souza (1987) state that a Peclet number of 4 or less is required to achieve a concentration solution that does not oscillate in space; whereas Huyakorn and Pinder (1983) suggest a value of 2 or less. For the two-dimensional simulations, the Peclet number in the X-direction varied from 200 in the north and south regions of the model (dx = 1 km), to 2 at the shoreline (dx = 10 m). Therefore, throughout the majority of the two-dimensional models the Pe criterion listed above is violated. In order to check that the numeric solutions obtained from the two-dimensional models were not affected by the chosen level of discretization, a simulation was (5.2) 115 performed in which the nodal spacing was refined to reduce the Peclet number in the X-direction to a uniform value of 2. The model was not run to completion due to the extensive simulation time that was required. Preliminary results of this simulation showed that numerical oscillations were reduced. However, the density distribution and predicted discharge through the model surface were not significantly altered. Therefore, the chosen level of horizontal discretization was retained for further simulations. In the Z direction, the Peclet number varied from 400 to 2 in the two-dimensional simulations. The effect of the level of vertical discretization was also investigated using an additional two-dimensional simulation. For this run, the Peclet number varied from 200 to 1. Results of this additional simulation also showed that the results were not significantly altered. Thus, the chosen level of vertical discretization was not altered for further model runs. 5.5.2. COURANT NUMBER The Courant number (Cr) is a measure of the travel distance of a solute particle over the duration of one timestep. Among others, Anderson and Woessner (1992) have found that a Cr of less than one is required to achieve a solution free of numerical oscillation. The Cr may be defined as: t>-5£ (53) Ai where Vj is the groundwater velocity (L/T), subscript i is the coordinate direction, Ai is the nodal spacing (L) in the i direction, and At is the duration of the timestep (T). 116 From equation (5.3), it can be seen that the largest value of Cr, and hence the most unstable region, will occur where the distance between nodes is least. For the two-dimensional models, the finest nodal spacing occurred directly beneath the model surface in the vertical direction. Here, the vertical distance between nodes was 1 m. Concurrently, the greatest Cr in relation to the duration of the timestep will occur at late time in the simulation when At has reached its maximum value. The maximum timestep that was used in the two-dimensional simulations was 100 days. Using the finest nodal spacing and the greatest timestep value, it can be shown that the vertical groundwater velocity beneath the seabed would have had to exceed 0.01 m/d to violate the Cr criterion. Results from the two-dimensional simulations showed that the maximum vertical groundwater velocities in this region rarely exceeded 0.01 m/d. The maximum vertical flow velocity encountered was 0.05 m/d, which occurred in one node in the sensitivity simulation where the mud layer was removed from the first 15 km of the seabed. Groundwater velocities in the surrounding nodes were below 0.01 m/d, preventing significant numerical oscillation from occurring. In the horizontal directions, the Cr criterion was never violated. Therefore, the temporal discretization that was utilized for the two-dimensional simulations was deemed adequate for this analysis. 117 5.6. SUMMARY The results obtained from the series of two-dimensional simulations predict that no fresh groundwater is discharging directly into the Gulf of Mexico near the shoreline. Low hydraulic gradients that exist onshore, combined with the flat topography of the coastal plain, do not result in any SGD. Groundwater flow beneath the coastal plain is driven by a horizontally extensive seawater recirculation system. Recharge waters enter the Coastal Lowlands aquifer system at the shoreline, and flow inland before discharging through the coastal plain. In the base case scenario of cross-section one, the recharge zone extends 500 m inland from the coast. In the base case scenario of cross-section two, the recharge zone only reaches 20 m inland from the shoreline. But again, no SGD is expected in the near shore region. Offshore convection cells that are present in both cross-sectional models do produce saline SGD further offshore. The convection cells span the entire thickness of the aquifer system, with discharge flux rates ranging from 0.01 cm/d to 0.04 cm/d in the base case scenarios of cross-sections one and two. Table 5.2 provides a summary of the results obtained from all of the two-dimensional scenarios that were simulated. As seen from the table, near shore SGD is not observed in any scenario. The shoreline recharge zone extends a variable distance inland that ranges from 10 m to 8 km inland from the coast. SGD is predicted further offshore due to convection of saline groundwater in all of the scenarios. The volume that is discharged and the location of the discharge zones are largely dependent on the permeability of the hydrogeologic units. 118 h f = (-pypo)v + z Figure 5.1. Schematic diagram of the boundary conditions employed for the two-dimensional models A. Cross-section 1. B. Cross-section 2. Figure 5.2. Fluid density and streamlines for the base case scenario. Units are in kg/m 3 . Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. -0.08 -0 .06 -0.04 5 -0.02 E x 3 0.02 0.04 0.06 2 3 7 X(km) 2 3 9 50 100 150 200 250 300 350 X(km) Figure 5.3. Groundwater fluxes through the surface of the model for the base case scenario. The inset figure provides a magnified view of the region surrounding the shoreline, located at X = 238 km. Discharge and recharge are denoted by negative and positive values, respectively. Z(m) -28.5 J_ A 237.5 238.0 238.5 X(km) 1 1 1 1 1 1002.5 1005.0 1007.5 1010.0 1012.5 1015.0 1017.5 1020.0 1022.5 X(km) Figure 5.4. Magnified views of the region surrounding the shoreline for the base case scenario. A. Flux vectors located at nodal positions. B. Streamlines. Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. Units are in kg/m 3 . 122 Figure 5.5. Measured fluid density near the location of cross-section 1. Units are in kg/m 3 . Hydrogeologic units are superimposed on top. Modif ied from Will iamson and Grubb ( 2 0 0 1 ) . MW 13 MW 10 MW 8 MW 14 • MW 7 MW 11 MW 12 • MW 9 MW 6 MW 5 t N MW 1 MW 2 Terrebonne Bay MW 4 0 50 - I 100 km Figure 5.6. Locations of shallow groundwater wells that were sampled. All wells are completed within Permeable Zone A. Well Depth Artesian Density (m) (kg/m3) MW01 131.1 NA 1001.8 MW02 125.0 NA 1001.8 MW03 91.4 NA 1000.8 MW04 106.7 NA 1002.9 MW05 61.0 Y E S 1002.0 MW06 64.0 Y E S 1001.9 MW07 73.1 Y E S 1002.3 MW08 57.9 Y E S 1002.1 MW09 61.0 Y E S 1002.4 MW10 97.5 NA 1004.6 MW11 91.4 NA 1003.3 MW12 91.4 NA 1001.9 MW13 61.0 Y E S 1002.9 MW14 54.9 Y E S 1002.4 Table 5.1. Well depth and average fluid density values of groundwater wells that were sampled. The location of the wells can be found in Figure 5.6. Artesian conditions could not be confirmed in all wells due to pump completion. All wells were completed within Permeable Zone A. 125 Figure 5.7. F luid density and streamlines for the case where the hydraulic conductivity of Permeable Zones A , B , and C was increased by a factor of 6. Units are in kg/m J . Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. 0 . 2 0 . 3 0.4 I I I I I I I I I I I I 0 5 0 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0 3 5 0 X (km) Figure 5.8. Groundwater fluxes through the surface o f the model for the case where the hydraulic conductivity of Permeable Zones A , B , and C was increased by a factor of 6. The inset figure provides a magnified view of the region surrounding the shoreline, located at X = 238 km. Discharge and recharge are denoted by negative and positive values, respectively. Z(km) 1002.5 1005.0 1007.5 1010.0 1012.5 1015.0 1017.5 1020.0 1022.5 Figure 5.9. F lu id density and streamlines for the case where the mud layer was removed from the first 15 k m of seafloor. Units are in kg/m 3 . Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. 0.08 I I I I I I—U I I I I 0 50 100 150 200 250 300 350 X(km) Figure 5.10. Groundwater fluxes through the surface of the model for the case where the mud layer was removed from the first 15 km of the seafloor. The inset figure provides a magnified view of the region surrounding the shoreline, located at X = 238 km. Discharge and recharge are denoted by negative and positive values, respectively. Figure 5.11. Fluid density and streamlines for the case where Permeable Zone A was assigned a hydraulic conductivity equal to the mud layer. Units are in kg/m3. Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. -0.01 -0.005 0.005 "D £ x 0.01 0.015 0.02 0.025 0.03 0.035 50 100 150 200 250 300 350 X ( k m ) Figure 5.12. Groundwater fluxes through the surface of the model for the case where Permeable Zone A was assigned a hydraulic conductivity equal to the mud layer. The inset figure provides a magnified view of the region surrounding the shoreline, located at X = 238 km. Discharge and recharge are denoted by negative and positive values, respectively. 1002.5 1005.0 1007.5 1010.0 1012.5 1015.0 1017.5 1020.0 1022.5 Figure 5.13. Fluid density and streamlines for the case where the dispersivity was increased an order of magnitude. Units are in kg/m J . Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. -0.08 -0.06 -0.04 -0.02 E & 0 x Z> 0.02 0.04 0.06 0.08 K to 0.04 x m o.08 0.12 235 237 239 X(km) 0 50 100 150 200 250 300 350 X(km) Figure 5.14. Groundwater fluxes through the surface of the model for the case where the dispersivity was increased an order of magnitude. The inset figure provides a magnified view of the region surrounding the shoreline, located at X = 238 km. Discharge and recharge are denoted by negative and positive values, respectively. 1002.5 1005.0 1007.5 1010.0 1012.5 1015.0 1017.5 1020.0 1022,5 Figure 5.15. F luid density and streamlines for the case where a third type exit boundary was assigned to the onshore nodes. Units are in kg/m 3 . Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. Figure 5.16. Groundwater fluxes through the surface of the model for the case where a third type exit boundary was assigned to the onshore nodes. The inset figure provides a magnified view of the region surrounding the shoreline, located at X = 238 km. Discharge and recharge are denoted by negative and positive values, respectively. F i g u r e 5.17. F luid density and streamlines for the case where only Permeable Zones A and B were active. Units are in kg/m J . Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. Also note the change in scale. -0.08 -0.06 -0.04 O =3 0.04 X ( k m ) Figure 5.18. Groundwater fluxes through the surface of the model for the case where only Permeable Zones A and B were active. The inset figure provides a magnified view of the region surrounding the shoreline, located at X = 238 km. Discharge and recharge are denoted by negative and positive values, respectively. Figure 5.19. Fluid density and streamlines for the case where a specified concentration boundary was applied to the seafloor. Units are in kg/m 3. Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. 0 50 100 150 200 250 300 350 X (km) Figure 5.20. Groundwater fluxes through the surface of the model for the case where a specified concentration boundary was applied to the seafloor. The inset figure provides a magnified view of the region surrounding the shoreline, located at X = 238 km. Discharge and recharge are denoted by negative and positive values, respectively. 50 100 4^ O 150 200 X ( k m ) 250 300 350 I I I I 1012.5 1020.0 1030.0 1040.0 1050.0 1060.0 Figure 5.21. Fluid density and streamlines for the case that included brines. Units are in kg/m . Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. Also note the change in density scale. 0.08 I 1 1 1 1 L_U_I 1 L_ 0 50 100 150 200 250 300 350 X(km) Figure 5.22. Groundwater fluxes through the surface of the model for the case that included brines. The inset figure provides a magnified view of the region surrounding the shoreline, located at X = 238 km. Discharge and recharge are denoted by negative and positive values, respectively. Figure 5.23. Fluid density and streamlines for the base case scenario of cross-section 2. Units are in kg/m 3 . Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. -0.4 -0.3 -0.2 Figure 5.24. Groundwater fluxes through the surface of the model for the base case scenario of cross-section 2. Discharge and recharge are denoted by negative and positive values, respectively. A -0.4 -0.2 "O E & 0 x =5 0.2 0.4 137 139 X ( k m ) 1 4 1 143 B 0 (cm/c 0.4 Flux 0.8 227 229 X ( k m ) 231 233 Figure 5.25. Magnified view of the groundwater fluxes through the surface of the model for the base case scenario of cross-section 2 for (A) the northern shore of Lake Ponchartrain and (B) the shore of the Gulf of Mexico. Discharge and recharge are denoted by negative and positive values, respectively. 144 Figure 5.26. Fluid density and streamlines for cross-section two for the case where dispersivity was increased by an order of magnitude. Units are in kg/m'. Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. -0.4 -0.3 Figure 5.27. Groundwater fluxes through the surface of the model for cross-section 2 for the case where the dispersivity increased an order of magnitude. Discharge and recharge are denoted by negative and positive values, respectively. A -0.4 -0.2 _3 0.2 137 139 141 143 X(km) 0.2 227 229 X(km) 231 233 Figure 5.28. Magnified view of the groundwater fluxes through the surface of the model for cross-section two for the case where the dispersivity was increased an order of magnitude for (A) the northern shore of Lake Ponchartrain and (B) the shore of the Gulf of Mexico. Discharge and recharge are denoted by negative and positive values, respectively. 147 Cross Section Simulation Description Landward Distance From The Shoreline That Discharge Is Observed To Stop (m) Predicted Volume Of SGD From Nearshore Convection Cell(s) (m3/d per meter width) Predicted Volume Of SGD From Convection Cell(s) Located Further Offshore (m3/d per meter width) 1 Base Case Scenario 500 2.2 9.1 1 Hydraulic Conductivity of Permeable Zones A, B, and C Increased By a Factor of 6 1500 0.2 31 Mud Layer Removed From The First 15 km Of The Seafloor 10 4.3 9.1 Permeable Zone Assigned a Low Permeability 8000 0.0 0.7 Dispersivity Increased By An Order Of Magnitude 300 2.5 9.1 Onshore Nodes Assigned a Third Type Boundary Condition 2000 2.7 9.1 Model Truncated At The Base Of Permeable Zone B 200 1.9 5.8 Seafloor Nodes Assigned a First Type Boundary Condition 500 2.7 9.1 1 Simulation That Included Brines 300 1.0 4.3 2 Base Case Scenario 20 1 6.3 0.4 2 Dispersivity Increased By An Order Of Magnitude 20 1 6.2 0.4 Table 5.2. Summary of the results from obtained from the two-dimensional models. Note 1. Observed distance from the shoreline at the Gulf of Mexico. 6. THREE-DIMENSIONAL SIMULATION OF SUBMARINE GROUNDWATER DISCHARGE IN THE GULF OF MEXICO NEAR SOUTHEASTERN LOUISIANA 6.1. INTRODUCTION A three-dimensional simulation was performed in order to obtain an estimate of the spatial distribution of SGD within the study region. A three-dimensional model would also be needed to properly account for groundwater withdrawals by any major pumping centers in the region. The model domain encompassed the region depicted in Figure 4.1. Simulation of the Coastal Lowlands aquifer system in two dimensions showed that no freshwater SGD was predicted in the nearshore region of the Gulf of Mexico due to onshore hydraulic gradients. However, further offshore, SGD was predicted due to convection processes operating beneath the seafloor. A single model was run due to the extensive simulation time that was required (i.e. approximately 150 days). The simulation time was long owing to the relatively short timesteps that were required to maintain numerical stability for the large number of nodes and rather coarse discretization that was used. The model was marched forward in time from the specified set of initial conditions that were described in section 4.8. As with the two-dimensional models, the model was run until a dynamic equilibrium was achieved, which occurred after 2.5 million days of simulation time. The simulation was deemed complete when salinity fronts and fluxes through the model surface ceased to change temporally. Several assumptions that were incorporated into the three-dimensional simulation should be noted. Sea level was assumed constant at a mean elevation of 0 m. Therefore transient SGD due to tidal pumping was not included. The maximum density of 149 groundwater was set to seawater (i.e. 1025 kg/m3). Previous work by Kuiper (1994) and results of the two-dimensional simulations showed that this simplification did not significantly affect the model predicted discharge through the seafloor. Thermally induced groundwater flow and compaction-driven flow were also not included. Previous modelling efforts by Sharp et al (2001) and Person et al (1996) have shown that topography and solute gradients may be expected to be the dominant forces driving groundwater flow in similar settings. Internal sources and sinks of groundwater were not included in the three-dimensional model. Therefore, effects due to production wells located within the Gulf region could not be accounted for. Mesko et al (1990) compiled a list of the major production wells located in the Gulf region for the U.S. Geological Survey models. From this list, it was found that few wells are present within the model region. With the exception of the New Orleans region, the majority of production wells are located tens of kilometres from the coastline. In any case, it is not expected that the inclusion of wells would have a significant affect on the predicted SGD. The wells would likely intercept groundwater flowing towards the coastline from the onshore flow system. This may allow the seawater recirculation system to extend further inland, but would not result in large changes in predicted offshore discharge rates. Future modelling efforts will be required to elucidate the effect, if any, that groundwater withdrawals have on offshore discharge rates within the study region. 150 6.2. BOUNDARY CONDITIONS The boundary conditions that were assigned to the three-dimensional model are the same as were assigned in the base case simulations of the two-dimensional models. The northern boundary was assigned constant freshwater heads equal to the elevation of land surface, and a specified concentration equal to freshwater. The variable southern boundary, located at the edge of the continental shelf, was assigned constant freshwater heads calculated for the specified concentration of seawater and an assumed water level equal to mean sea level. Hydrostatic conditions were assumed for both boundaries. Constant freshwater heads were assigned to onshore nodes representing the water table, along with a specified concentration equal to freshwater. Offshore nodes were also assigned constant freshwater heads calculated for a concentration of seawater and a water level equal to mean sea level. A third type exit boundary was specified as the boundary condition for solute transport across the seabed (see section 5.2. for description). It was assumed that the lateral boundaries approximately represented groundwater flow lines. These boundaries were thus set to be impermeable to both flow and transport. The base of the Coastal Lowlands aquifer system was also specified to be impermeable. A plan view diagram of the assigned boundary conditions can be found in Figure 6.1. 151 6.3. RESULTS OF THE THREE-DIMENSIONAL SIMULATION Parameter values employed for the three-dimensional simulation are the same as for the two-dimensional base case simulations, with the exception of the dispersivity. In order to accommodate the coarser grid spacing that was used, the dispersivity values assigned to all hydrogeologic units were increased by an order of magnitude (i.e. cti = 50 m, a t = 5 m). Therefore, results will be compared to the two-dimensional cases that employed the same values of dispersivity. Results from the three-dimensional simulation will first be discussed along the regions of cross-sections 1 and 2. It is not expected that the predicted magnitude of groundwater fluxes will exactly match the two-dimensional results. The coarser grid spacing that was employed in this simulation will likely result in some minor differences being present. Also, groundwater flow paths along the cross-sections are not strictly two-dimensional. Differences in the predicted fluxes due to this reason are particularly evident in the offshore regions. 6.3.1. RESULTS NEAR CROSS-SECTION 1 Figure 6.2 provides a cross-sectional view of the density distribution and groundwater flow lines predicted near the location of cross-section 1. The figure provides a two-dimensional slice of a three-dimensional flow system (i.e. the component of flow into and out of the plane of the section is neglected). Therefore, the flow lines presented are rather more qualitative than those presented in Section 5. From Figure 6.2, it is seen that the dominant flow patterns that were predicted in the two-dimensional simulations were also predicted in three dimensions. The onshore flow system is present at the 152 northern extent of the domain, the extensive seawater recirculation system occupies the sediments beneath the coastal plain, and convection processes control offshore groundwater flow patterns. The resulting density distribution is similar to the two-dimensional results. Recharge waters enter the onshore flow system through the upland region and flow south towards the coastline. Predicted recharge rates are similar to the two-dimensional results, with a maximum flux rate of 0.037 cm/d (see Figure 6.3). As discussed in the previous chapter, the model-predicted flux rates match quite well with previously published results for the Gulf region (Williamson and Grubb, 2001). The onshore system terminates at the southern extent of the uplands at approximately X = 120 km. Here, a major discharge zone is located at the northern boundary of the coastal plain. The maximum discharge rate predicted within this zone reaches 0.097 cm/d (see Figure 6.3). Throughout the remainder of the coastal plain, the model predicts a low rate of groundwater discharge. The discharge is driven by the extensive seawater recirculation system, which spans a distance of approximately 120 km south to the shoreline (located at X = 238 km). Recharge waters enter the Coastal Lowlands aquifer system seaward of the shoreline before convecting northwards toward the coastal plain. The average recharge rate near the shoreline is 0.14 cm/d. The shoreline recharge rate is nearly identical to what was found in the two-dimensional model for the case with equivalent dispersivities. As was the case in two-dimensions, no SGD is predicted in the nearshore region. Groundwater discharge through the coastal plain ceases 10 km from shore. 153 The nearshore recharge zone extends only to X = 242 km. However, a greater volume of saline recharge is predicted to pass through the seabed in this simulation. Seawater passing through the seabed in the southern portion of the recharge zone flows in a southerly direction and becomes entrained in the offshore convection system. The system is composed of two cells. Discharge rates from the first convection cell average 0.041 cm/d, resulting in 16 m3/d of SGD per metre width. The predicted volume is almost double what was found in two-dimensions due to the increased volume of recharge at the shoreline. A second offshore recharge zone extends from X = 280 km to 310 km. Seawater flowing through the seabed in the recharge zone is discharged in the second convection cell. Discharge rates average 0.013 cm/d. The volume of SGD produced by this cell is 9.5 m3/d per metre width, similar to what was predicted in the two-dimensional model. 6.3.2. RESULTS NEAR CROSS-SECTION 2 Plots of the streamlines and density distribution, and fluxes through the model surface near the location of cross-section 2 are provided in Figures 6.4 and 6.5. From these figures, it can be seen that flow patterns and flux rates predicted along this section are also similar to the results of the two-dimensional simulations. An onshore flow system extends from the northern boundary to approximately X = 125 km. Recharge waters enter the upland region at an average flux rate of 0.021 cm/d. The recharged water flows southwards before discharging through a zone located at the southern extent of the onshore system from X = 100 km to 114 km. Discharge rates reach 0.15 cm/d. 154 Flow in the central portion of cross-section 2 is controlled by the seawater recirculation system. Recharge waters enter the sediments at the shoreline to the Gulf (located at X = 230 km) and flow north due to convection. The majority of the saline groundwater discharges through the patch of land located south of Lake Ponchartrain. A portion of the groundwater flows under the lake and discharges near its northern shore. The recharge zone located at the shoreline to the Gulf of Mexico continues to X = 246 km. Water entering the southern portion of this zone becomes entrained in the convection-driven flow patterns present in the offshore regions. In this zone, SGD due to convection occurs only in one area along the section. Discharge is predicted to pass through the seabed from X = 248 km to 278 km at an average rate of 0.037 cm/d, resulting in a volume of 11 m3/d per meter width. 6.4. THREE-DIMENSIONAL DISTRIBUTION OF RECHARGE AND DISCHARGE A plan view diagram of recharge and discharge zones is provided in Figure 6.6. This diagram shows the role topography has in determining the locations of groundwater recharge and discharge in the onshore flow system. North of the 25 m topographic contour, recharge areas are confined to zones with locally elevated topography. Groundwater discharge zones are predicted to occur within small valleys and local depressions that are superimposed on the uplands. The large discharge zone that was described in the two-dimensional models is found on the southern side of the 25 m contour. Here, the topographic gradient becomes significantly less. The discharge zone is found to extend across the entire model region, generally following the trend of the 25 m contour. 155 South of the 25 m topographic contour, the majority of the coastal plain is predicted to serve as a groundwater discharge zone. Small areas of recharge are interspersed due to small depressions present in the plain. Artesian conditions that were found in the majority of wells that were sampled confirmed that this region is primarily a discharge zone. However, in a "real world" setting, groundwater is not found to be actively discharging through the plain. Local confining units that are present in the subsurface serve to slow the upward vertical movement of groundwater. As can be seen from Figure 6.7, the magnitude of discharge through the coastal plain is predicted to be quite low. The groundwater that is discharged probably seeps into the numerous surface water bodies and marshy regions that are present throughout the coastal region. At the shoreline, the discharge zone gives way to groundwater recharge. Density contrasts that are present drive groundwater flow in this region, pulling seawater down into the seawater recirculation system. The recharge zone extends a number of kilometres inland, causing fresh recharge waters to become entrained as well. The recharge zone extends along the coastline of Louisiana throughout the study region. In the offshore regions of the model, alternating bands of discharge and recharge are present due to the presence of convection cells. These convection cells represent the only source of SGD within the offshore study region. The bands have an east-west trend, with groundwater flow paths generally following a north-south direction. The convection cells extend tens to over a hundred kilometres in length. The magnitude of groundwater recharge and discharge fluxes due to the cells is quite low because the driving force (groundwater density contrasts located at the shoreline) is located a significant distance away. \ 156 A box has been drawn on Figure 6.7 to compare the volume of SGD predicted by the model to Mississippi River discharge. The region contained within the box is similar to the region investigated by Krest et al (1999). The sum of all SGD within this box is a relatively insignificant 0.4 km3/yr, or 0.1% of the volume discharged by the river. This value is quite close to what was predicted using the two-dimensional model of cross-section 1. As alluded to in Chapter 5, whether or not the model is able to replicate the true convective regime occurring offshore cannot be determined. Other factors such as salt domes and the geothermal gradients make the system quite complex, but are beyond the scope of this work. Results from the tracer study may help to add some insight into this topic. 6.5. DIMENSIONLESS PARAMETERS 6.5.1. PECLET NUMBER For the three-dimensional simulation, a coarser level of horizontal discretization was required due to the regional scale of the analysis. Nodal spacing varied from 5 km in the northern and southern regions of the model, to 2 km in the central coastal region (from X = 80 km to X = 280 km). To accommodate the larger nodal spacing, the longitudinal dispersivity that was employed in the simulation was also increased to a value of 50 m. Using these values along with equation (5.2), it is found that the mesh Pe has values of 100 and 40. Therefore throughout the modelled region, the criterion (i.e. Pe = 2 or 4) was violated by at least an order of magnitude. 157 Due to the violation of the Pe criterion some amount of numerical dispersion was expected, particularly in the shoreline region where large density gradients are present. Indeed, a low amount of instability did develop at early time in the runs. The majority of the instability was found to smooth out over time as the solution progressed. This was especially true within the study region. Persistent instability was present along the shoreline in the eastern portion of the model. This region was located along the Mississippi coastline, where flow velocities were higher due to steeper topographic gradients. This area does however lie well outside of the tracer study region. Therefore, it is believed that these instabilities should not significantly affect the results in the area of interest. In any case, it would not have been possible to reduce the Pe and still simulate the entire model region with a realistic amount of dispersion. Simulation times were found to be exceedingly long when finer horizontal resolution was employed. The vertical discretization that was employed in the three-dimensional model was identical to the two-dimensional simulations. Therefore, the mesh Pe in the vertical direction is found to be an order of magnitude less than what was reported in Chapter 4 due to the order of magnitude increase in transverse dispersivity (i.e. 5 m). The values of Pe thus varied from 40 to 0.2. Elevated values of Pe were present at depth where the groundwater was uniformly saline, negating the criterion. Therefore, the level of vertical discretization that was used was deemed adequate for this analysis. 6.5.2. COURANT NUMBER Using equation (5.3), it can be shown that horizontal fluid velocities would need to exceed 20 m/d to violate the Cr criterion of 1 (assuming a minimum node spacing of 2 km and a maximum timestep of 100 days). This value is significantly higher than any 158 elemental groundwater velocity that was predicted by the model. In the vertical direction, the discretization was kept the same as for the two-dimensional models. Therefore, the same maximum vertical groundwater velocity applies here (i.e. 0.01 m/d). As was the case in two dimensions, this value was rarely exceeded. 6.6. SUMMARY Results of the three-dimensional simulation are consistent with the results that were obtained with the two-dimensional simulations. Groundwater flow within the modelled region can be described by the same three flow systems: onshore, seawater recirculation, and offshore convection. Some differences in the predicted magnitude of groundwater fluxes are present; however these can be attributed to the resolution used and the three-dimensional nature of groundwater flow paths. As was found in two-dimensions, no SGD is predicted to occur in the vicinity of the shoreline. Groundwater recharge is predicted to pass through the seabed along the coast throughout the study region. Low rates of SGD are predicted in the offshore portion of the model due to the convection cells. The cells occur in bands extending from east to west, with groundwater flowing in a north-south direction. However, relative to Mississippi River discharge, the predicted magnitude of SGD is insignificant. 159 h ( = {-pjp0h + z C = 0 • Figure 6.1. Schematic diagram of the boundary conditions employed for the three-dimensional simulation. 160 Figure 6.2. Fluid density and streamlines for the three-dimensional model near cross-section 1. Units are in kg/m 3 . Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. -0 .1 -0.05 0 0,05 0.1 0.15 0.2 0 50 100 ,150 200 250 300 350 X (km) Figure 6.3. Groundwater fluxes through the surface of the three-dimensional model near cross-section 1. Discharge and recharge are denoted by negative and positive values, respectively. Figure 6.4. Fluid density and streamlines for the three-dimensional model near cross-section 2. Units are in kg/m 3 . Note that the location and density of streamlines does not indicate the magnitude of groundwater fluxes. Figure 6.5. Groundwater fluxes through the surface of the three-dimensional model near cross-section 2. Discharge and recharge are denoted by negative and positive values, respectively. Figure 6.6. Predicted recharge and discharge areas for the three-dimensional model. Topographic contours have been superimposed on top. Contour interval in metres. 0 100 200 300 400 X(km) Figure 6.7. Predicted recharge and discharge rates for the three-dimensional model. The shoreline has been superimposed on top. The dotted line indicates the area that was used to calculate the volume of S G D for comparison to Mississippi River discharge. 7. COMPARISON OF HYDROGEOLOGIC MODEL SUBMARINE GROUNDWATER DISCHARGE ESTIMATES WITH TRACER-BASED ESTIMATES 7.1. INTRODUCTION Geochemical tracers have been used to study SGD in estuarine settings by a number of groups (e.g. Moore, 1996). Chemical species that have been used include, but are not limited to, Ra (e.g. Krest et al, 1999), Rn (Cable et al, 1996), He (e.g. Top et al, 2001), and Ba (e.g. Shaw et al, 1998). The basic premise behind tracer-based studies is to quantify the background (natural) concentration of the tracer in the water column, quantify all of its sources and sinks, and than estimate the rate of groundwater input based on a mass balance approach. In general, the tracer inventory is a balance between (Corbett, 2002): • Advective-diffusive exchange across the sediment-water interface • In situ production and decay • Horizontal water column advection • Eddy diffusion through the pycnocline • Air-sea exchange The relative importance of each process is dependent on the tracer used. A schematic diagram of these processes can be found in Figure 7.1. On a regional scale, tracer-based estimates provide a spatially averaged value of SGD. Small-scale variations are smoothed out by mixing in the water column (Burnett et al, 2001). Small-scale variability has been a concern in studies using seepage meters due 167 to the large number of measurements that must be taken to characterize the heterogeneity of the sediments and the resulting groundwater flow rates. Burnett et al (2001) state that, ideally, geochemical tracers should be greatly enriched in groundwater relative to coastal waters, and should be conservative and easy to measure. Adherence to these factors aids in ease of both measurement and analysis. Both Ra and Rn meet these criteria fairly well. Helium would be ideally suited for SGD studies, except for the difficulties encountered in its measurement. The following section will provide a brief summary of the chemistry of Ra, and Rn, relevant to their use as geochemical tracers, along with the calculations that were performed to quantify the magnitude of SGD in the study region. Results of the tracer-based studies that have been performed within the study region (both past and current) will then be discussed. Finally, a comparison will be made between the results derived through the hydrogeologic models described previously, and those obtained through the use of geochemical tracers. 7.2. GEOCHEMICAL TRACERS THAT HAVE BEEN USED TO MEASURE SUBMARINE GROUNDWATER DISCHARGE IN THE STUDY REGION 7.2.1. RADIUM (Ra) Krest eial (1999) measured 2 2 6Ra {Un = 1600 years) and 2 2 8Ra (\m = 5.75 years) in the Gulf of Mexico in the discharge zones of the Mississippi and Atchafalaya Rivers. Surface water samples were collected during three periods (November, 1993; March, 1994; June, 1994) to evaluate the flux of Ra from the river systems to the ocean and to determine if the differences in the discharge characteristics at the river mouths affected 168 the Ra fluxes. They found unexpectedly high concentrations of both Ra isotopes, causing them to reevaluate the factors controlling the Ra fluxes and consider SGD as a potential source. The extent of their study region, as well as sample locations from the November 1993 cruise can be found in Figure 7.2. The respective parents of 2 2 6Ra and 2 2 8Ra are 2 3 0Th (tU2 = 7.5 x 104 years) and 2 3 2Th (ti/2 = 1.4 x 1010 years). The parents are strongly bound to particles in the aquifer matrix. Ra produced through the decay of thorium (Th) is sorbed in freshwater aquifers. The ratio of Ra adsorbed on sediment to Ra dissolved in the surrounding water varies inversely with the ionic strength of solution (Moore et al, 1995). Thus in brackish and saline seawater, Ra is present in the dissolved phase. Production of Ra in the mixing zone between fresh and saline groundwater provides a constant source of the tracer. Groundwater flowing through the mixing zone will become enriched in Ra relative to seawater, allowing its use as a tracer of SGD. Krest et al (1999) used a mass balance approach to determine the amount of Ra being added to the system via SGD, given by equation (7.1): C = CR+CDes+CBol (7.1) where C is the dissolved concentration of Ra in the water column from the river component (CR), desorption from riverine particles (Coes), and input from the bottom sediments (Ceot) such as sediment resuspension, solute diffusion, and advection of groundwater. C , CR and Coes were measured directly using water and sediment samples, while CBot was calculated using equation (7.1). 169 7.2.2. RADON (Rn) Measurements of Rn in groundwater and the coastal water column of southeastern Louisiana were taken by researchers from East Carolina University (ECU) as a part of this project. Water column samples were taken during a total of six cruises to provide information regarding temporal variability in SGD. Three sampling periods took place during the fall of 2003 (i.e. October, November, and December) and during the spring of 2004 (i.e. March, April, and May) which corresponded to low and high discharge conditions in the Mississippi River, respectively. Groundwater samples were collected after each cruise. Offshore sample locations may be found in Figure 4.6. The locations of groundwater wells that were sampled can be found in Figure 5.6. Naturally occurring dissolved Rn gas is produced through the decay of 226Ra. Rn is typically elevated in groundwater relative to coastal waters due to its production from Ra in the aquifer matrix (Corbett, 1999). Due to its conservative (non-sorbing) nature, Rn may be transported from the subsurface through the seabed by both fresh and saline groundwater. It has a short half-life (ti/2 = 3.83 days), and thus does not accumulate within the water column. A simple box model was used to describe the source and sink terms of Rn in the coastal system (see Figure 7.1). The mass balance equation that was used to calculate estimates of SGD in the study region is given by: J In + J Ben + ^Resus +Jpiod ~ ^Out ~ ^Pyc ~ ^ Decay ~ R"/(it where J Ben = J Adv + ^ Diff (7-3) 170 In equation (7.2), the terms J i n and Jout account for the horizontal transport of tracer into and out of the box, Jprod and JDecay represent the production and decay of Rn in the water column, J B e n is the benthic flux of Rn across the sediment/water interface (i.e. sum of advection and diffusion), j R e s u s represents the flux of Rn due to resuspension events, and Jpyc is the eddy diffusive flux across the pycnocline. The advective flux of Rn due to groundwater flow through the seabed was calculated after all other sources and sinks had been accounted for using equation (7.2). 7.3. RESULTS OF THE TRACER-BASED ESTIMATES 7.3.1. RADIUM (Ra) Results of the surface water samples collected in the plumes of the Mississippi and Atchafalaya Rivers by Krest et al (1999) showed that concentrations of 2 2 6Ra and 2 2 8Ra were elevated above what was expected. Using equation (7.1), Krest et al (1999) could not explain the observed concentrations based on contributions from the ocean water and river water (both dissolved and sorbed) endmembers alone. They concluded that a deeper water source was causing the observed Ra enrichment. A number of possible deep water sources of the tracer were examined, including resuspension of bottom sediments and desorption of Ra, biological scavenging and subsequent remineralization of Ra, diffusion from pore spaces, and groundwater discharge across the seabed. Krest et al (1999) state that resuspension of bottom sediments could not be the source due to the length of time required to regenerate the observed level of Ra. In the case of 226Ra, they suggest that 600 years would need to pass between resuspension 171 events based on their calculations. Storm events, which suspend the sediments into the seawater and thus desorb Ra from them, are far more frequent. Biological scavenging of Ra during the formation of biogenic barite or opal has been seen to occur in seawater (e.g. Bishop, 1988). After formation, the mineral particles sink before coming to rest on the seabed. Krest et al (1999) found that their data suggested that some removal of Ra might have occurred in the surface waters of the Atchafalya plume. In seasons where the hypoxic zone forms, the barite or opal crystals could be remineralized, releasing the scavenged Ra into the surface water. Hypoxic conditions were present during the June 1994 sampling period. However, the observed Ra concentrations were not as high as during the other sampling periods, suggesting that biological scavenging and remineralization was not the primary source of Ra enrichment. Diffusion of Ra from pore spaces was also disregarded as the primary source for Ra. Krest et al (1999) calculated the maximum diffusive flux that could be expected based on the concentrations of the parent isotopes in the bottom sediments. From their calculations, they found that at most 8% of the observed enrichment could be from diffusion, leaving groundwater discharge as the only possible source. Based on the observed Ra enrichment and the concentration of Ra in a groundwater sample collected in Barataria Bay, Krest et al (1999) predicted that groundwater may be discharging across the seabed at an advective velocity of approximately 1 cm/day. Their calculation assumed that the concentration of Ra in the groundwater was spatially constant. If their estimate is assumed to represent an average value for the region, than it would be equivalent to an SGD rate of 40 km3/yr, or 10% of Mississippi River discharge. This calculation is based on a study area that has dimensions 172 similar to what used in the calculation for the results of cross-section 1 (i.e. 120 km north to south and 90 km east to west). 7.3.2. RADON (Rn) Diagrams of the predicted groundwater discharge rate across the study region as determined from Rn samples are provided in Figures 7.3 - 7.5 (data courtesy of C. McCoy and R. Corbett, January, 2005). At the time this thesis was prepared, only a preliminary analysis had been performed. The calculated discharge rates did not account for horizontal transport into and out of the study area (J;n and Jout), nor the flux of Rn from bottom sediments due to resuspension events (JR e Sus)- Neglecting horizontal transport appears to have influenced the results at the eastern edge of the study region, where inflow of Rn from the Mississippi River has increased the calculated groundwater discharge (see Figure 7.4). The rest of the study region appears unaffected. The amount of Rn in the water column due to sediment resuspension should be insignificant unless a resuspension event occurred near the time of sampling, due to the short half life of Rn (Corbett, 2002). Figure 7.3 depicts the average calculated discharge from Rn measurements taken during the fall 2003 cruises. As seen from this figure, little SGD is observed. Discharge rates are less than 0.1 cm/day throughout much of the study region. A peak in the groundwater discharge rate is found in the southwest corner of the study region. Here, upward fluxes approach 0.4 cm/day. The average discharge rate for the spring 2004 cruises is presented in Figure 7.4. Little seasonal variability in groundwater discharge appears to be present within the study 173 region, as the calculated amount of SGD has increased only slightly relative to the fall of 2003. The lack of seasonal variability suggests that the source of SGD in the study region may not be from groundwater recharged onshore, as the rate of groundwater discharge would be expected to fluctuate with the rate of precipitation. Estimated upward fluxes of groundwater reached just over 0.6 cm/d in the southwest corner of the study region. The average groundwater discharge rate for all of the sampling periods is provided in Figure 7.5. The plot shows that in general, upward groundwater flux rates are observed to increase with distance seaward of the shoreline, although the majority of the study region is estimated to have little SGD. The maximum calculated discharge rate in this diagram is 0.5 cm/d. Figure 7.5 will be used to compare the tracer-based estimates of SGD to those obtained from the hydrogeologic models. 7.4. COMPARISON BETWEEN METHODS 7.4.1. NUMERICAL MODELLING AND RADIUM (Ra) The results obtained by Krest et al (1999) are significantly different than what was predicted by the both the hydrogeologic model and the Rn tracer study. Their predicted SGD rate is an order of magnitude greater than what was predicted by the numerical simulations. Reasons for the difference in results are unclear. However, the most likely cause is the Ra concentration they used for the groundwater endmember in their calculations. A single groundwater sample, collected in Barataria Bay, was used to determine this value. If the modelling results are assumed correct, than groundwater in this region would be flowing inland within the seawater recirculation system, and would 1-74 not be expected to discharge offshore. Therefore, the concentration used would not provide a correct rate of discharge. If again, the model results are assumed to accurately depict the groundwater flow system within the region, then offshore groundwater flow paths pass deep within the sediments of the continental shelf. Here, measured Ra concentrations in saline formation waters have been found to be significantly higher (Kraemer and Reid, 1984). Krest et al (1999) state that it would take only a small fraction of this water (altered or diluted) to produce the observed Ra enrichment. Their calculations suggest that a discharge rate of only 0.03 cm/d would be required. It is likely that the true concentration of Ra in any groundwater that is discharged offshore lies somewhere between that measured by Krest et al (1999) and what was measured by Kraemer and Reid (1984). Therefore, as noted by Krest et al (1999), it is necessary to quantify the source of groundwater in order to accurately determine the rate of offshore discharge. Unfortunately, neither tracer study is able to confirm the deep flow paths predicted by the models. 7.4.2. NUMERICAL MODELLING AND RADON (Rn) Results of the hydrogeologic model and the current geochemical tracer study (i.e. Rn) are quite similar. Both methods predict no SGD immediately offshore of coastal Louisiana. Discharge is predicted to occur further offshore, at distances greater than 5 km from the shoreline. The predicted rate is low, with values generally less than 0.5 cm/d. Discharge zones that are predicted by the model to occur further offshore cannot be confirmed because they lie outside of the tracer study region. A magnified view of the model results in the tracer study region has been plotted in Figure 7.6. The average calculated rate of offshore groundwater discharge from the Rn 175 measurements (plotted in Figure 7.5) has been superimposed on top of the plot. As can be seen from this figure, the spatial location of discharge predicted by the two methods does not line up exactly. Results from the tracers suggest that the discharge zone occurs further offshore than what was predicted by the model. There may be several reasons for the discrepancy. First, the horizontal resolution used in the three-dimensional model may not have been fine enough to accurately predict the spatial location of SGD. At the location of cross-section 1, the three-dimensional model predicts the northern edge of the discharge zone to be 5 km from the shoreline, as compared to a distance of 18 km predicted by the two-dimensional model. The result from the two-dimensional model matches more closely with the tracer results, which predict the northern edge of the discharge zone to occur at a distance of approximately 15 km from shore near the location of cross-section 1. This result highlights one of the main drawbacks of using numerical modelling to simulate SGD on a regional scale. There is a trade-off between sufficient grid resolution and solution accuracy that is involved. Given current computing resources, it may not be possible to model the finer details of this process on a regional scale. The second reason may be the density of sampling points used in the Rn tracer study. As shown in Figure 4.6, water column samples were taken at only fifteen locations across the study region. Therefore, a significant amount of interpolation was required to generate the contours of SGD plotted in Figures 7.3 - 7.5. A more finely spaced sampling grid would be required to definitively determine locations of offshore discharge. A third potential reason for the discrepancy may be the inherent variability in the permeability of the seabed. The location of the discharge zone may correlate with a 176 region of increased seabed permeability relatively closer to shore. No measure of the hydraulic conductivity of the seabed was made, so this possibility cannot be ruled out. The three-dimensional model may have predicted the spatial location of the discharge zone accurately based on the parameter values that were provided. A closer match between the model and tracer results may have been attained if the heterogeneity of the seabed could have been incorporated. However, the overall match between the two methods is quite good. In terms of an ecological perspective, it appears that groundwater may not be expected to play a significant role in the development of hypoxic conditions within the Gulf of Mexico. Elevated concentrations of nitrogen may be present within groundwater flowing through the hydrogeologic units in the onshore regions. However, this groundwater is not predicted to discharge offshore. In locations where discharge is predicted, the recharge site appears to be offshore. Groundwater from this source should carry the nutrient signature of the Gulf. Therefore, any plans aimed at mitigating the development of hypoxic conditions should focus on riverine transport of nitrate. 177 Mississippi River oo Horizontal Transport S u r f a c e D e e p A t m o s p h e r i c E x c h a n g e P y c n o c l i n e In Situ P roduc t i on In Situ R e m o v a l Horizontal Transport 1 Diffusive Flux S e d i m e n t Resuspens ion A d v e c t i v e Flux Figure 7.1. Generalized diagram (or box model) depicting the sources and sinks of geochemical tracers within the study region. Note that not all processes are relevant to each tracer. Modified from Corbett (2002). Figure 7.2. Study region for the radium tracer research performed by Krest et al (1999). Circles denote sample locations. Modified from Krest et al (1999). Calculated discharge (cm/d) l-x. - O L O CO CM i — 0 0 0 0 0 0 0 0 180 0 50 km I I Figure 7.4. Average calculated rate of S G D using radon from the spring of 2004. 0 50 km J Figure 7.5. Average calculated rate of S G D using radon from all sampling periods. 150 200 250 300 Y(km) Figure 7.6. Comparison of predicted discharge rates for the three-dimensional model and the 2 2 2radon (averaged over all sampling periods). Units are in cm/d. 8. SUMMARY AND CONCLUSIONS Results of the hydrogeologic simulations show that the groundwater flow patterns of the Louisiana Gulf Coast region can be described by three systems. These include an onshore system, a seawater recirculation system, and an offshore convection system. In the base case simulations, the onshore flow system extends from the upland region at the northern boundary of the model approximately to the 25 m topographic contour. The distribution of recharge and discharge areas is controlled by topography of the land surface. All water entering this system is discharged at its southern extent, with no fresh SGD reaching offshore. The seawater recirculation system extends from the southern boundary of the onshore system to a few kilometres south of the shoreline. Recharge waters enter the system at the shoreline and flow inland before being discharged through the coastal plain. Discharge rates through the coastal plain are low, with an average value of about 0.01 cm/d. The offshore convection system spans the remainder of the flow domain. Convection is driven by small density gradients present near the shoreline. The convection cells are predicted to extend several kilometres into the subsurface, although their true nature cannot be confirmed due to the simplifications that were incorporated in the model. The model predicts that little groundwater is being discharged into the Gulf of Mexico within the study region, and none of the groundwater discharged originates from onshore recharge locations. Recharge waters passing through the seabed at the shoreline prevent any discharge from occurring in the nearshore region. SGD is predicted to occur 184 further offshore due to the convection system. The predicted rate of discharge is low, with an average value less than 0.5 cm/d. From the three-dimensional simulation, the total volume is predicted to be about 0.1% of the rate of freshwater discharged by the Mississippi River. A number of sensitivity simulations were performed in order to investigate whether there is any scenario under which an increased amount of SGD may be expected. Variations from the base case scenario that were investigated included increasing the permeability of the hydrogeologic units and the seabed, increasing the dispersivity, the effects of the chosen set of boundary conditions, and the effect of denser brines at depth. Results from all of the simulations were consistent in that little offshore groundwater discharge was predicted. Also, the three dominant flow systems were not significantly altered under any scenario. The predicted rate of SGD compares quite well with results of the tracer study that was performed using Rn. Some difference in the predicted locations of groundwater discharge is present. This may be due to the horizontal discretization that was used in the models, the low density of tracer sampling points, or variability in seabed permeability. Results from both the modelling and geochemical tracer aspects of this project differ from what was previously found by Krest et al (1999). They predicted groundwater discharge in the region to occur at a rate of approximately 1 cm/d. If a volume of SGD is calculated across a region similar to what was used for the model results, this equates to approximately 10% of annual Mississippi River discharge. This value decreases dramatically, however, if groundwater flow paths are assumed to pass deep within the 185 offshore sediments before being discharged. Clearly, identifying the pathway of any offshore discharge is an important component of any tracer study. Future research efforts should focus on describing groundwater flow in the offshore portion of the study region, particularly at depth. Salt domes and the associated dense brines are present in the subsurface throughout the study region. As shown in the two-dimensional simulations, the brines that are present may exert an influence on discharge patterns at the surface. Further complicating the complicating the situation is the increased fluid temperature that is found with depth. Future modelling efforts may be able to describe what influence these factors may have on offshore discharge. 186 A1. DEPOSITIONAL HISTORY OF THE STUDY A R E A A1.1. INTRODUCTION A brief summary of the depositional history of southern Louisiana and Mississippi is provided here. A more detailed description can be found in Hosman (1996) and Coleman and Roberts (1988) and the references contained within. The sediments making up the Coastal Lowlands aquifer system were deposited during the Cenozoic. Much of the information that is available today has been developed through interpretation of drill logs. A1.2. STRUCTURAL FEATURES The major structural features that controlled deposition in the Gulf Coastal plain are the Gulf Coast geosyncline and the Mississippi embayment (see Figure A 1.1). The formation of these features is associated with the Ouachita orogeny that lead to the deformation of the Paleozoic surface prior to Mesozoic deposition. The Gulf Coast geosyncline is the result of folding, which formed a catch basin for sediments. Further downwarping and faulting occurred as a result of the weight of accumulated sediments. Associated downwarping and rifting created the Mississippi embayment, resulting in a southward plunging trough that opens into the geosyncline. During the Mesozoic, deposition resulted in vast amounts of sediment accumulating in the geosyncline. It was filled and deepened during the Triassic and Jurassic. Lesser, but still substantial, amounts of sediments were deposited in the Mississippi embayment during this period as well. These sediments formed the floor upon which Cenozoic deposition would eventually occur. 187 Although the Gulf Coast geosyncline and the Mississippi embayment were the dominant factors controlling depositional patterns, some smaller features superimposed upon the landscape also played a role. One of these, the Sabine uplift (Figure Al.l), is located in northwest Louisiana and extends into Texas, where it occupies over 13 000 square kilometres. It is underlain by a high area in the basement rocks, of which little is known. The Sabine uplift was submerged by transgressing seas during the early Tertiary, but has remained above water ever since. The Monroe uplift is located to the east of the Sabine uplift (Figure Al.l). It is found mostly in northeast Louisiana, but also is present in neighbouring corners of Arkansas and Mississippi. It was slowly submerged during the late Cretaceous - early Tertiary, when limestone beds up to several hundred metres thick formed. During the Palaeocene, the rate of submergence increased allowing limestone formation to give way to typical deeper water marine deposition of clays. The Monroe uplift is adjacent to the Mississippi embayment to the east, and was effective in constricting sediment movement. The Mississippi embayment was also constricted on its eastern flank by the Jackson dome, located in west - central Mississippi (Figure Al.l). It was formed by an igneous intrusion that began at the end of the Cretaceous and continued at least until the Oligocene. Depositional patterns around the Jackson dome were similar to those occurring near the Monroe uplift during the late Cretaceous - early Tertiary, where hundreds of feet of limestone were also formed. Toward the end of the Palaeocene, marine clays were deposited as the transgressing sea deepened. Deposition occurred through the Eocene, until the land began to re-emerge and has been exposed since. 188 The Wiggins uplift is located in southern Mississippi and extends to the southwest corner of Alabama. It bifurcates in the southeast corner of Mississippi. The shorter segment arcs southwest towards the southeast tip of Louisiana, where it is referred to as the Hancock arch or the Harris ridge (Figure Al.l). The Terrebonne embayment, located in south-central Louisiana, is essentially an extension of the Mississippi embayment (Figure Al.l). It has a northwest - southwest trend, and was effective in trapping large volumes of sediments during depositional periods. A1.3. DEPOSITIONAL PATTERNS Sediments composing the Coastal Lowlands system were deposited during the Cenozoic. The depositional period began at the end of the Cretaceous, when the Midway Sea began to transgress. Thus began a period of successive transgressions and regressions spanning the Palaeocene and Eocene. The initial transgression lasted throughout the Palaeocene, where the sea reached into southern Missouri and Kentucky. Marine deposition continued into the early Eocene, when the sea began to retreat. Smaller and less extensive transgressions continued to occur, depositing dense marine clays. But the sea did not occupy the Mississippi embayment after the Eocene. Marine deposition was restricted to the Gulf Coastal Plain for the remainder of the Tertiary. Slow net regression continued throughout the Tertiary, until the sea had retreated as far back as the Gulf Coast geosyncline. During this period, deltaic (typically lobate shaped) deposits were left interspersed with beach and bar deposits. Organic deposits were formed in the shallow marine environment, later altering to peat and lignite deposits. Inland, fluvio-lacustrine deposits formed great thicknesses of sand and clay with 189 limited amounts of gravel. The paleo-drainage network that was present acted to rework and transport sediment throughout the basin. During the Quaternary, deltaic and marine deposition patterns continued. But, inland, depositional patterns were altered as a result of glaciation on the northern part of North America. Coarse sands and gravels were transported to the south and deposited by meltwater streams during the interglacials, forming terrace deposits. The result of the cyclic depositional patterns that persisted through the Cenozoic has been the formation of a deep sedimentary basin of alternating beds of sand and clay. The sand beds range from thin discontinuous beds interbedded with clay to massive sands hundreds of metres thick. Most of the beds have been found to be interconnected to some degree. The thickness of the sand beds and interspersed clays change abruptly throughout the basin, and individual beds are often not traceable over a regional extent. This has lead to differing interpretations of stratigraphic correlations between regions. Table Al. l gives some examples of the nomenclature that has been used in the model region. A1.4. STRATIGRAPHIC UNITS Sediments of the Jackson group form the base of the confining unit that bounds the Coastal Lowlands aquifer system. They were deposited in a marine environment, when the sea occupied the Mississippi embayment for the last time. In Louisiana and Mississippi, the basal unit of the group is known as Moody's Branch Marl. It is composed of sandy marl, with a fining upward sequence as the sediments become more calcareous. Thin interbeds of marl and impure limestone are scattered throughout the unit. The upper bed of the Jackson group is the Yazoo Formation. It is undifferentiated in Louisiana, while in Mississippi it has been divided into the Shubuta, Pachuta Marl, Cocoa 190 Sand, and North Twistwood Creek members. It contains mostly marl and limestone, along with small amounts of calcareous clay and sand. The Vicksburg group overlies the Jackson group in both states. The two groups have been differentiated based on faunal assemblages rather than changing lithology. In Louisiana, it is undifferentiated. In Mississippi, it is composed of the Tatum Limestone, Paynes Hammock, Chickasawhay, Bucatunna, Byram, Glendon, Marianna Mint Spring, Forest Hill, and Red Bluff formations. Due to their similar lithologies, the Vicksburg and Jackson groups were grouped together to form a regional confining unit that bounds the Coastal Lowlands aquifer system. The deepest water conducting sediments of the Coastal Lowlands in the model region is the Catahoula sandstone. It is the uppermost unit of the Vicksburg group. It contains primarily tuffaceous clays, sands and sandstones that were associated with volcanic activity occurring to the west. Interbedding within the unit is irregular, and individual strata are commonly discontinuous. In Mississippi, a clay unit is present above the sandstone, known as the Hattiesburg clay. The unit is up to 150 metres thick, but is not areally extensive. Miocene sediments of the two states are referred to using different nonmenclature. In Louisiana, they are grouped in the Fleming formation, which is divided into the Blounts Creek, Castor Creek, Williamson Creek, Dough Hills, Carnahan Bayou, and Lena members. The members of the formation are composed of interbedded sands, silts and clays that show a general coarsening upwards. In Mississippi, the sediments are grouped in the Pascagoula formation. It has been divided into the Fort Adams and Homochitto members. It is composed principally of clays with interbedded sands and 191 lesser amounts of sandstone. Some coarse sand and gravel beds are located near the center. The overlying Pliocene unit in Mississippi is the Graham Ferry formation, which contains deltaic sediments. Above this is the Citronelle formation, which is also recognized in southeastern Louisiana. The units are composed of nonmarine sands and clays. In southwestern Louisiana, sediments are grouped in the Foley formation. It has two units, known as the Mamou and Steep Gully members. It contains interbedded sands and clays, with the base of the formation being coarser than the top. Sediments lying above the Pliocene units have not been differentiated. Those that were deposited during the Pleistocene show features dominated by erosional and sedimentary cycles associated with glaciation. Stepped terraces formed by successive glacial meltwater events are present. The lower terraces are commonly gravel, or sand and gravely and grade to finer deposits of silt and clay. The terrace deposits merge with deltaic sediments towards the south. During the Holocene, sedimentary deposition was of two types: alluvial and coastal. Sand and gravel beds are present in the upper flood plain, where the energy of the Mississippi River was greater due to steeper topography. The river's gradient was reduced as it approached the Gulf of Mexico and flowed over the flat coastal plain. Silt and clay deposition took place in the lower flood plain. Figure A 1.2 provides a geologic map of Louisiana and Mississippi. Only those stratigraphic units which are part of the Coastal Lowlands aquifer system, or the Vicksburg-Jackson confining unit have been included. 192 A1.5. OTHER FEATURES Growth faulting has occurred throughout the Gulf Coast basin, as sediments were accumulating. All of the faulting is associated with the subsidence of the Gulf Coast geosyncline. The faults generally parallel the shoreline. Some fault zones are still active (e.g. the Baton Rouge fault zone). Another feature present in the Gulf Coast basin are numerous salt domes. The largest concentration is present in a band that extends along the coast from southeast Texas to southeast Louisiana in what is called the Gulf Coast Salt Basin. To the north of the model region two other clusters of salt domes are present. These are called the North Louisiana Salt Basin and the Mississippi Salt Basin. The locations of the salt basins can be seen in Figure A1.3. A list of all known salt domes within the study region, along with maps showing their distribution and which hydrogeologic unit they penetrate can be found in Beckman and Williamson (1990). All of the salt domes originate from the Louann salt bed. The age of this bed is currently not agreed upon, but it is generally thought to have formed sometime between the Permian and Jurassic. The salt domes began forming during the late Jurassic to early Cretaceous due to density differences between the overburden and salt bed. They formed in a series of pulses of isostatic adjustments rather than at a smooth, steady rate. They are generally one to three miles in diameter, and are surrounded by caprock (typically anhydrite). Faulting around the domes is common, as surrounding sediments sunk around the less dense salt. 193 A1.6. SUMMARY Deposition of sediments composing the Coastal Lowlands aquifer system took place during the Cenozoic. The main structural controls on deposition were the Gulf Coast geosyncline and the Mississippi embayment, although other smaller scale features played a role. Depositional environments of the region varied from marine and beach deposits early on to fluvial deltaic and terrace deposits associated with glaciation. The result has been the formation of a thick sedimentary basin composed primarily of interbedded sands, silts, and clays, with lesser amounts of limestone, lignite and gravels. Growth faults due to subsidence of the sediments are present in bands that generally parallel the shoreline. Also, numerous salt domes have formed throughout the region due to density differences between the accumulated sediments and the underlying Louann salt bed. 194 Figure Al . l . Structural features that controlled depositional patterns in the study area. 195 Erathem r stem eries Louisiana Mississippi Erathem co GO Cent ra l a n d southern Quaternary Holocene Alluvium and terrace deposits Alluvium, terrace, and loess deposits Quaternary Pleistocene Pliocene Southwestern Southeastern Citronelle Formation Graham Ferry Formation Pliocene Foley Formation Mamou Member Steep Gully Member Citronelle Formation Cenozoic Miocene Fleming Formation Blounts Creek Member Castor Creek Member Williamson Creek Member Dough Hills Member Carnahan Bayou Member Lena Member Pascagoula Formation Fort Adams Member Homochito Member Hattiesburg Clay Tertiary Oligocene Catahoula Sandstone Vicksburg Formation Vicksburg Group Catahoula Sandstone Tatum Limestone Member Paynes Hammock Formation Chickasaway Limestone Bucatunna Formation Byram Formation Glendon Formation Marianna Formation and Mint Spring Formation Forest Hill and Red Bluff Formations Eocene Jackson Group Yazoo Formation Moodys Branch Marl Jackson Group Yazoo Formation Shubuta Member Pachuta Member C o c o a Sand Member North Twistwood Creek Member Moodys Branch Marl Table Al . l . Stratigraphic units present in Louisiana and Mississippi. 196 Quaternary Alluvium (Qal) Figure A1.2. Surficial geologic map of Louisiana and Mississippi. Only those units that are part of the Coastal Lowlands aquifer system and underlying Vicksburg-Jackson confining unit have been included. Modified from Hosman (1996). 197 Figure A 1 . 3 . Location of salt basins present within the study area. 198 A2. BENCHMARK PROBLEMS A2.1. INTRODUCTION The density-dependent flow component of FRAC3DVS had not been verified against the standard test cases prior to its use for this project. Therefore, seven benchmark problems were simulated in order to verify that the code was able to accurately simulate ground water flow of variable density. The problems that were simulated include: • Box problem • Henry problem • Elder problem • HYDROlogic COde INtercomparison (HYDROCOIN), Level 1, Case 5 • Salt lake problem • Salt pool problem (low density case) • Rotational flow of three immiscible fluids (symmetric case) The first four problems listed have commonly been used to verify previous codes, while, the final three problems have been developed in recent years to provide more rigorous testing of density-dependent simulators. In the following paragraphs, a brief description of each problem and results obtained will be provided. A2.2. BOX PROBLEM This problem serves as a test of velocity consistency under hydrostatic conditions and a sharp density transition (Diersch and Kolditz, 2002). The issue of velocity consistency was first proposed by Voss and Souza (1987). They found that the discontinuity of discretized pressure gradients at element boundaries could give rise to 199 artificial velocities. By allowing both the density and local gravity vector components to vary over an element in a consistent formulation, they were able to remove these numerical artefacts from the solution. Recently, Taylor et al (2001) benchmarked a number of codes using this problem and concluded that none of the codes could match the hydrostatic conditions. The flow domain consists of a two dimensional rectangular box. Dimensions of the box were 20 m (horizontal) by 40 m (vertical), although any length would suffice. Nodal spacing was set to 0.625 m in both directions. Initially, the lower half of the box is filled with fluid of seawater density, while the overlying water is specified to be fresh. Hydrostatic heads were specified throughout, representing a stable configuration. The results of the problem should show that the density interface remains in the cell that it was initially set in. Two variations of the problem were simulated. In the first case, all boundaries were set to be impermeable. Longitudinal dispersivity was set equal to cell length, while transverse dispersivity was set to zero. In the second case, a horizontal head gradient is applied, and water is allowed to enter and leave the box. In both simulations, the interface was found to remain in the same cell throughout, with no artificial vertical velocities being generated. Results were checked against a reference simulation that did not incorporate density coupling to ensure that any solute spreading was due solely to dispersion. 200 A2.3. HENRY PROBLEM This problem was first introduced by Henry in 1964, when he developed a semi-analytical solution for a problem of confined groundwater flowing toward a seawater boundary. In order to simplify his analysis, he dropped higher order terms that were deemed insignificant. As noted by Kolditz et al (1998), he also chose to use a large value for the diffusion coefficient instead of incorporating dispersion (velocity-dependent dispersion would have precluded the formulation of his semi-analytic solution). The result is that no code has been able to match his results exactly (Guo and Langevin, 2002). Instead of striving to obtain an exact match to the semi-analytical solution, previous authors have found it sufficient to compare results between codes. Further complicating the solution, the test case was modified by Segol et al (1975) to include a fresh water outlet overlying the saltwater boundary. The usefulness of this problem as a test case has been called into question by some authors (Simpson and Clement, 2003; Weatherill et al, 2004). Simpson and Clement (2003) simulated this problem with and without density coupling. They found little difference in the results, stating that the solution to the problem was determined primarily by boundary forcing, rather than by density effects. Bearing'this in mind, the problem is still one that has been included in benchmarking tests by most codes, and was also included here. The results are depicted in Figure A2.1, along with those obtained by Frind (1982). As can be seen, the results of the two codes are quite similar. 201 A2.4. ELDER PROBLEM This standard verification problem, originally termed "the short heater problem" was created experimentally by Elder (1967) for heat flow. It was later transformed by Diersch (1981) and Voss and Souza (1987) to be applied to brine flow. It consists of a closed box, with inflow/outflow points present at the upper corners. On the surface of the box's centre a layer of brine is simulated by applying a concentration boundary. A fresh water boundary is present along the base of the domain. Initially, the box is filled with fresh water. Brine enters the domain with time by diffusion (the sole mechanism of hydrodynamic dispersion considered), and leads to the formation of two convection cells. This problem is deemed a good benchmark by Simpson and Clement (2003) because the flow patterns are totally dependent on density coupling. Figure A2.2 shows the results obtained with FRAC3DVS, along with those from SEA WAT (Guo and Langevin, 2002), SUTRA (Voss and Souza, 1987), and Elder (1967). Again, the match is seen to be quite close. Recent studies have shown that the convective flow patterns that develop are dependent upon the level of discretization that is used (e.g. Diersch and Kolditz, 2002). Where coarse meshes have been used to simulate the problem, a central downwelling pattern is formed. As the mesh is refined, the pattern switches to a central upwelling. Continued mesh refinement leads to the reappearance of the central downwelling pattern. As no analytical solution exists, the correct pattern is of yet unknown. Both convective patterns were obtained, with the results found in Figure A2.3. The coarse grid used 45 nodes in the horizontal direction (dx = 13.3 m), and 28 in the vertical direction (dz = 5.36 202 m); while the fine grid had 200 nodes in the horizontal direction (dx = 3 m) and 75 in the vertical direction (dz = 2 m). A2.5. HYDROCOIN PROBLEM, LEVEL 1, CASE 5 The goal of level 1 of the HYDROCOIN project was to evaluate the accuracy of selected ground water codes Konikow et al (1996). Case 5 of the project depicts a simple scenario of ground water flow over a salt dome, located along the middle third of the base of the domain. Flow is induced by applying a head gradient across the surface of the domain (directed from left to right in Figure A2.4). Brine is released into the flow field through dispersion, and is transported towards the right side of the domain where it may be discharged through the surface boundary. All other boundaries are specified to be impermeable. The results obtained for this problem have been found to take two forms by other authors, depending on the value of the solute diffusion coefficient that is used. Oldenburg and Pruess (1995) found that when the diffusion coefficient was set to values of 1 x 10"7 m2/s or less, the calculated concentration distribution took on a swept forward appearance, as seen in Figure A2.4 (A). When they set the diffusion coefficient to 2 x 10"6m2/s or greater, the flow field was found to be re-circulating (Figure A2.4 (B)). Results obtained using FRAC3DVS were similar to those of Oldenburg and Pruess (1995). The swept forward result was obtained using a diffusion coefficient of 0.5 x 10~7 m2/s, while the re-circulating flow pattern was found to occur with a diffusion coefficient of 3 x 10"6 m2/s. 203 A2.6. SALT LAKE PROBLEM The salt lake problem was first introduced by Wooding et al (1997). They used a Hele-Shaw cell to simulate the idealized situation of evaporation occurring through a salt lake. Their experimental data showed that the process leads to the formation of dense brines forming beneath the lake. With time, the brines form fingers which sink due to their greater density. Numerical models have been constructed by Wooding et al (1997) (stream function based model), Simmons et al (1999) (SUTRA), Mazzia et al (2001) (mixed hybrid finite element model) and Diersch and Kolditz (2002) (FEFLOW) to simulate the experimental results. All four models were able to reproduce the general features that were found in the experiment, although none were able to reproduce the results exactly. Also, differences between the results of each of the models were found. The salinity distributions obtained with FRAC3DVS are plotted in Figure A2.5, along with results from Diersch and Kolditz (2002) for comparison. It is seen that both results are quite similar. A2.7. SALT POOL PROBLEM This problem was developed based on experimental results obtained by Oswald (1998) and described in Johannsen et al (2002). It represents a three-dimensional saltwater upconing process induced through pumping. The box is initially filled with fresh water overlying saltwater (1% and 10% salt mass fraction cases), separated by a thin transition zone. Water is pumped out of one corner of the box at a constant rate, while in the opposite corner fresh water is allowed to recharge the box. During the 204 experiment, detailed measurements of the outflow concentration were obtained, allowing for comparison to numerical simulations. Attempts to reproduce the experimental results have been made by Johanssen et al (2002), Diersch and Kolditz (2002), and Oswald and Kinzelbach (2004) with varying degrees of success. The main difficulty appears to be the level of discretization required to adequately capture the upconing process, particularly with the high density cases. Johanssen et al (2002) performed a mesh convergence study. They found that a sufficient level of accuracy could be achieved for the low density case with as few as 4096 nodes; but, for the high density case a minimum of about 17 million nodes were required. For this reason, only the low density case was simulated. Results are depicted in Figure A2.6, along with those from Diersch and Kolditz (2002). As can be seen from the figure, the results obtained do not match the experimental results with a great degree of accuracy. They do, however, match the FEFLOW simulation results reasonably well, particularly at later time. A2.8. ROTATIONAL FLOW OF THREE IMMISCIBLE FLUIDS Bakker et al (2004) recently developed this problem as a test case for verifying density-dependent flow codes. The problem consists of a two dimensional box which is filled with three fluids of varying density. Adjacent fluids are separated by abrupt, diagonal interfaces, with the fluid of higher density being situated on top (see Bakker et al, 2004 for problem setup). Dispersion is purposefully not taken into account so that density differences are the only driving force for fluid motion. Both a symmetric and asymmetric case (with respect to the density distribution) were developed. Only results 205 from the symmetric case will be presented here. Initial fluid densities were 1000, 1012.5, and 1025 kg/m3. An analytical equation was developed for the initial density configuration (Bakker et al, 2004). The flow field generated from the equation shows how the fluids will rotate about the center of the domain, aiming to achieve a stable configuration. Simulations performed by the authors showed how the density distribution evolved through time. They used three different models, including SWI (Bakker, 2003), SEA WAT (Guo and Langevin, 2002) and MOCDENS3D (Oude Essink, 1998) with results from all three being quite similar. The results obtained using FRAC3DVS are found in Figure A2.7. Due to the fact that the results from all four models essentially plot on top of each other, only the results from MOCDENS3D were plotted for comparison. Comparison of the two solutions shows that a good match was obtained. A2.9. SUMMARY Seven benchmark problems were simulated in order to verify the density-dependent component of FRAC3DVS. Results obtained were compared to published results from a variety of different codes. In all cases FRAC3DVS was able to produce results that compared well with those presented by other authors, providing confidence that it could solve the set of equations accurately. Therefore, the code was chosen for the simulations that were performed for this thesis. 206 Figure A2.2. Relative concentration plots for the Elder Problem for (A) 10 years, and (B) 20 years. Contour levels are 0.2 and 0.6. 208 1 5 0 -Z f m l A , * ,» / A . A W W . ? w- ? 1 * ' i » v. V 1. . hi H r »\ ft-/1 f t • \ < ^ v J " \ ' - '< v * -I X(m) "I 6 0 0 Figure A2.3 . Vector plots for the Elder Problem at 20 years for different levels of discretization. Plots (A) and (B) have coarse and fine grid spacing, respectively. Also plotted are the 0.2 and 0.6 relative concentration contours., 209 X(m) Figure A2.4. Relative concentration and vector plots for the HYDROCOIN Problem. Swept forward (A) and re-circulating patterns (B) obtained using diffusion coefficients of 0.5 x W6 nr/s and 3 x 10° mVs, respectively. 210 Figure A2.5. Salinity values obtained for the Salt Lake Problem. Left are the results from F R A C 3 D V S , right are the results from Diersch and Kolditz ( 2 0 0 2 ) using FEFLOW. Contour interval is 2 x 10"\ 2 1 1 0.06 Experiment • • Quads; Diersch and Kolditz (2002) T i m e (s) Figure A2.6. Salinity breakthrough curves obtained at the outlet for the Salt Pool Problem. Plotted are the experimental results, simulation results from Diersch and Kolditz (2002) using FEFLOW (quadrilateral elements consisting of 274 625 nodes and symmetric halfspace using triangular elements with 140 010 nodes), and simulation results from FRAC3DVS. 212 ure A2.7. Evolution of the salinity interfaces for the Rotational Flow Problem. Results are plotted for 2000 and 10 000 days. A3. MEASURED SOLUTE CONCENTRATION AND DENSITY DATA FROM THE SHALLOW GROUNDWATER WELLS Table A3.1 provides measured solute concentrations for the shallow groundwater wells. Sample collection was performed by researchers from East Carolina University in Greenville, North Carolina. Fluid density was calculated using the following equation p = p0+EC (A3.1) where p is the fluid density (M/L ), p0 is the freshwater density (M/L ), C is the solute concentration (M/L3) and E is a dimensionless constant having an approximate value of 0.7143 for salt concentrations ranging from zero to that of seawater (Guo and Langevin, 2002). 214 SITE Concentration (mg/L) Density (kg/m3) Mar-03 MW01 1005.0 1000.7 MW02 1226.1 1000.9 MW03 1139.0 1000.8 MW04 1641.5 1001.2 MW05 2164.1 1001.5 MW06 1453.9 1001.0 MW07 1447.2 1001.0 MW08 1407.0 1001.0 MW09 1306.5 1000.9 MW10 1936.3 1001.4 MW11 1494.1 1001.1 MW12 1132.3 1000.8 MW13 1587.9 1001.1 MW14 1252.9 1000.9 Jul-03 MW01 1520.9 1001.1 MW02 1594.6 1001.1 MW03 MW04 MW05 1058.6 1000.8 MW06 857.6 1000.6 MW07 1159.1 1000.8 MW08 1259.6 1000.9 MW09 1051.9 1000.8 MW10 1507.5 1001.1 MW11 1165.8 1000.8 MW12 5487.3 1003.9 MW13 1085.4 1000.8 MW14 1072.0 1000.8 Oct-03 MW01 2914.5 1002.1 MW02 2820.7 1002.0 MW03 MW04 6499.0 1004.6 MW05 3075.3 1002.2 MW06 3283.0 1002.3 MW07 3966.4 1002.8 MW08 3283.0 1002.3 MW09 4582.8 1003.3 MW10 9179.0 1006.6 MW11 6164.0 1004.4 MW12 2010.0 1001.4 MW13 5226.0 1003.7 MW14 4556.0 1003.3 Table A3.1. Measured solute concentration and fluid density values of shallow groundwater wells. 215 SITE Concentration (mg/L) Density (kg/m3) Dec-03 MW01 3216.0 1002.3 MW02 3182.5 1002.3 MW03 MW04 MW05 3530.9 1002.5 MW06 3236.1 1002.3 MW07 3939.6 1002.8 MW08 3932.9 1002.8 MW09 4288.0 1003.1 MW10 8040.0 1005.7 MW11 6003.2 1004.3 MW12 2412.0 1001.7 MW13 4757.0 1003.4 MW14 4334.9 1003.1 Mar-04 MW01 3149.0 1002.2 MW02 3169.1 1002.3 MW03 MW04 MW05 3537.6 1002.5 MW06 3316.5 1002.4 MW07 4267.9 1003.0 MW08 3973.1 1002.8 MW09 4656.5 1003.3 MW10 9447.0 1006.7 MW11 6532.5 1004.7 MW12 2412.0 1001.7 MW13 5594.5 1004.0 MW14 4334.9 1003.1 May-04 MW-01 3082.0 1002.2 MW-02 2793.9 1002.0 MW03 MW04 MW-05 3256.2 1002.3 MW-06 3644.8 1002.6 MW-07 4288.0 1003.1 MW-08 3993.2 1002.9 MW-09 4314.8 1003.1 MW-10 8710.0 1006.2 MW-11 6030.0 1004.3 MW12 MW-13 5862.5 1004.2 MW-14 4428.7 1003.2 Table A3.1 continued. Measured solute concentration and fluid density values of shallow groundwater wells. 2 1 6 REFERENCES Bakker, M. 2003. A Dupuit formulation for modeling seawater intrusion in regional aquifer systems. Water Resources Research 39 (5) SBH12. Bakker, M., Essink, G.H.P.O., and Langevin, CD. 2004. The rotating movement of three immiscible fluids - a benchmark problem. Journal of Hydrology 287, 270-278. Beckman, J. D. and Williamson, A.K. 1990. Salt dome locations in the Gulf Coastal Plain, south-central United States. U.S. Geological Survey Water Resources Investigations Report 90-4060,44 p. Biogeochemistry. 2003. Volume 66, Issues 1-2, 202 p. Bishop, J.K.B., 1988. The barite-opal-organic carbon association in oceanic particulate matter. Nature 332, 341-343. Burnett, W.C. 1999. Offshore springs and seeps are focus of working group. EOS 80, 13-15. Burnett, W. C , Cowart, J.B., and Deetae, S. 1990. Radium in the Suwanee River and estuary: Spring and river input to the Gulf of Mexico. Biogeochemistry 10, 237-255. Burnett, W, C , Taniguchi, M., and Oberdorfer, J.A. 2001. Measurement and significance of the direct discharge of groundwater into the coastal zone. Journal of Sea Research 46, 109-116. Burnett, W.C. and Turner, J. 2001. LOICZ group investigates groundwater discharge in Australia. LOICZ Newsletter, 18, 1-4. Burnett, W.C, Chanton, J., Christoff, J., Kontar, E., Krupa, S., Lambert, M., Moore, W., O'Rourke, D., Paulsen, R, Smith, C , Smith, L., and Taniguchi, M. 2002. Assessing methodologies for measuring groundwater discharge to the ocean. EOS 83, 117-123. Burnett, W.C, Bokuniewicz, H., Huettel, M., Moore, W.S., and Taniguchi, M. 2003. Groundwater and pore water inputs to the coastal zone. Biogeochemistry, 66, 3-33. Cable, J. E., Burnett, W.C, Chanton, J., and Weatherly, G.L. 1996. Estimating groundwater discharge into the northeastern Gulf of Mexico using radon-222. Earth and Planetary Science Letters 144, 591-604. 217 Cable, J. E., Burnett, W.C, Chanton, J., Corbett, R.D. and Cable, P.H. 1997. Field evaluation of seepage meters in the coastal marine environment. Estuarine, Coastal and Shelf Science 45, 367-375. Cable, J.E., Martin, J.B., Swarzenski, P.W., Lindenberg, M.K., and Steward, J. 2004. Advection within shallow porter waters of a coastal lagoon, Florida. Ground Water 42 (7), 1011-1020. Carter, R. C , Kaufman, W.J., Orlob, G.T., Todd, D.K. 1959. Helium as a ground-water tracer. Journal of Geophysical Research 64, 2433-2439. Charette, M. A., Splivallo, R, Herbold, C , Bollinger, M.S., Moore, W.S. 2003. Salt marsh submarine groundwater discharge as traced by radium isotopes. Marine Chemistry 84, 113-121. Coffey, M., Dehairs, F., Collette, O., Luther, G, Church, T., Jickells, T. 1997. The behaviour of dissolved barium in estuaries. Estuarine, Coastal and Shelf Science 45,113-121. Coleman, J.M., and Roberts, H.H. 1988. Geo-Marine Letters 8 (63), 63-108. Corbett, D.R. 2002. "Collaborative Research: A Multiple Tracer Approach to Evaluate Groundwater Discharge into River Dominated Coastal Waters". National Science Foundation research proposal. Corbett, D. R, Burnett, W.C, Cable, P.H., Clark, S.B. 1997. Radon tracing of groundwater input into Par Pond, Savannah River site. Journal of Hydrology 203, 209-227. Corbett, D. R., Chanton, J., Burnett, W.C, Dillon, K., Rutkowski, C , Fourqurean, J.W. 1999. Patterns of groundwater discharge into Florida Bay. Limnology and Oceanography 44 (4), 1045-1055. Destouni, G. and Prieto, C. 2003. On the possibility for generic modeling of submarine groundwater discharge. Biogeochemistry 66, 171-186. Dickinson, G. 1953. Geological aspects of abnormal reservoir pressures in Gulf Coast Louisiana. American Association of Petroleum Geologists Bulletin 37 (2), 410-432. Diersch, H.-J. G. 1981. Primitive variables finite element solutions of free convection flows in porous media." Journal of Applied Math and Mechanics (ZAMM) 61, 325-337. Diersch, H. G. 1996. Interactive, Graphics-Based Finite Element Simulation System FEFLOW For Modeling Groundwater Flow, Contaminant Mass and Heat 218 Transport, WASY Institute for Water Resource Planning and System Research Ltd., Berlin Germany. Diersch, H.-J. G. and Kolditz, O. 1998. Coupled groundwater flow and transport: 2. Thermohaline and 3D convection systems. Advances in Water Resources 21, 401-425. Diersch, H.-J. G. and Kolditz, O. 2002. Variable-density flow and transport in porous media: approaches and challenges. Advances in Water Resources 25, 899-944. Domenico, P.A., and Schwartz, F.W. 1997. Physical and Chemical Hydrogeology, Second Edition. John Wiley and Sons, Toronto, 506 p. Elder, J. W. 1967. Transient convection in a porous medium. Journal of Fluid Mechanics 27, 609-623. Essink, G.H.P.O. 1998. MOC3D adapted to simulate 3D density-dependent groundwater flow. Proceedings of the MODFLOW '98 Conference, October 4-8 1998, Golden , Colorado, USA, Vol. 1, p. 291-303. Evans, D.G., Nunn, J.A., and Hanor, J.S. 1991. Mechanisms driving groundwater flow near salt domes. Geophysical Research Letters 18 (5), 927-930. Frind, E. O. 1982. Simulation of long-term transient density-dependent transport in groundwater. Advances in Water Resources 5, 73-87. Gandl, L.A. Characterization of aquifers designated as potential drinking water sources in Mississippi. U.S. Geological Survey Water-Resources Investigations Open-File Report 81-550, 90 p. Gelhar, L.W., Welty, C , and Rehfeldt, K.R. 1992. A critical review of data on field-scale dispersion in aquifers. Water Resources Research 28 (7), 1955-1974. Ground Water. 2004. Special Issue, Volume 42, Number 7, 957-1102. Guo, W. and Langevin, CD. 2002. User's guide to SEA WAT: A computer program for simulation of three-dimensional variable-density ground-water flow. U.S. Geological Survey Techniques of Water-Resources Investigations 6-A7, 77 p. Henry, H. R. 1964. Effects of dispersion on salt encroachment in coastal aquifers, sea water in coastal aquifers. U.S. Geological Survey Water Supply Paper 1613-C, 70-84. Hosman, R. L. 1996. Regional stratigraphy and subsurface geology of Cenozoic deposits, gulf coastal plain, south-central United States. U.S. Geological Survey Professional Paper 1416-G, 35 p. 219 Johannes, R.E. 1980. The ecological significance of the submarine discharge of groundwater. Marine Ecology Progress Series 3, 365-373. Johannsen, K., Kinzelbach, W, Oswald, S.E., and Wittum, G. 2002. The saltpool benchmark problem - numerical simulation of saltwater upconing in a porous medium. Advances in Water Resources 25, 335-348. Jones, P. H. 1969. Hydrology of Neogene deposits in the northern Gulf of Mexico basin. Louisiana Water Resources Research Institute Bulletin GT-2, 105 p. Jones, P.H., Hendricks, E.L., and Irelan, B. 1956. Water resources of southwestern Louisiana. U.S. Geological Survey Water-Supply Paper 1364, 460 p. Kaleris, V., Lagas, G., Marczinek, S., Piotrowski, JA. 2002. Modelling submarine groundwater discharge: an example from the western Baltic Sea. Journal of Hydrology 265, 76-99. Kohout, F. A. 1966. Submarine springs: A neglected phenomenon of coastal hydrology. Hydrology 26, 391-413. Kohout, F. A. and Kolipinski, M.C. 1967. Biological zonation related to ground water discharge along the shore of Biscayne Bay, Miami, Florida. In: Estuaries. G Lauff (ed.). American Association Advancement Science, Washington, DC, Publication No. 83, pp. 488-499. Kolditz, O., Ratke, R., Diersch, H.-J.G., and Zielke, W. 1998. Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models. Advances in Water Resources 21 (1), 27-46. Konikow, L. F., Campbell, P.J., Sanford, W.E. 1996. Modelling brine transport in a porous medium: a re-evaluation of the HYDROCOIN Level 1, Case 5 problem. In: Calibration and Reliability in Groundwater Modelling (Proceedings of the ModelCARE 96 Conference). Kovar, K. and van der Heijde, P. (eds). IAHS Publ. 237, p. 363-372. Kraemer, T.F., and Reid, D.F. 1984. The occurrence and behaviour of radium in saline formation water of the U.S. Gulf Coast region. Isotope Geoscience 2, 153-174. Krest, J. M., Moore, W.S., Rama. 1999.226Ra and 2 2 8Ra in the mixing zones of the Mississippi and Atchafalaya Rivers: indicators of groundwater input. Marine Chemistry 64, 129-152. Kuiper, L.K. 1994. Nonlinear regression flow model, gulf coast aquifer systems, south-central United States. U.S. Geological Survey Water-Resources Investigations Report 93-4020. 171 p. 220 Langevin, C. D. 2001. Simulation of ground-water discharge to Biscayne Bay, Southeastern Florida, U.S. Geological Survey Water-Resources Investigations Report 00-4251, 127 p. Langevin, C. D. 2003. Simulation of submarine ground water discharge to a marine estuary: Biscayne Bay, Florida. Ground Water 41 (6), 758-771. Li, L., Barry, D.A., Stagnitti, F., and Parlange, J.-Y. 1999. Submarine groundwater discharge and associated chemical input to a coastal sea. Water Resources Research 35 (11), 3253-3259. Martin, A. Jr. and Whiteman, CD. Jr. 1999. Hydrology of the Coastal Lowlands aquifer system in parts of Alabama, Florida, Louisiana, and Mississippi. U.S. Geological Survey Professional Paper 1416-H, 51 p. Martin, J.B., Cable, J.E., Swarzenski, P.W., and Lindenberg, M.K. 2004. Enhanced submarine ground water discharge from mixing of pore water and estuarine water. Ground Water 42 (7), 1000-1010. Mazzia, A., Bergamaschi, L., and Putti, M. 2001. On the reliability of numerical solutions of brine transport in groundwater: analysis of infiltration from a salt lake. Transport in Porous Media 43 (1), 65-86. McLaren, R.G. 2003. GridBuilder 5.5 User's Guide: A pre-processor for 2D, triangular element, finite element programs. Waterloo Centre for Groundwater Research, University of Waterloo, Waterloo, Canada, 90 p. Mesko, T.O., Williams, T.A., Ackerman, D.J., and Williamson, A.K. 1990. Groundwater-pumpage from the gulf coast aquifer systems, 1960-85, south-central United States. U.S. Geological Survey Water-Resources Investigations Report 89-4180, 182 p. Meyer, R. R. and Turcan, A.N.J. 1955. Geology and ground-water resources of the Baton Rouge area, Louisiana. U.S. Geological Survey Water Supply Paper 1296, 137 p. Milliman, J. D. 1991. In: Ocean Margin Processes in Global Change. J R. Mantoura, J.-M. Martin and R. Wollast (eds.). John Wiley and Sons, New York, p. 69-90. Mississippi River/Gulf of Mexico Watershed Nutrient Task Force. 2001. Action Plan for reducing, mitigating, and controlling hypoxia in the northern Gulf of Mexico. Washington, DC: Office of Wetlands, Oceans and Watersheds, US Environmental Protection Agency. In: Turner, R. E., Rabalais, N.N., Swenson, E.M., Kasprzak, M., and Romaire, T. 2005. Summer hypoxia in the northern Gulf of Mexico and its prediction from 1978 to 1995. Marine Environmental Research 59, 65-77. 221 Moore, W. S. 1996. Large groundwater inputs to coastal waters revealed by ZZ0Ra enrichments. Nature 380, 612-614. Moore, W. S. 1997. High fluxes of radium and barium from the mouth of the Ganges-Brahmaputra River during low river discharge suggest a large groundwater source. Earth and Planetary Science Letters 150, 141-150. Moore, W.S., Astwood, H., and Lindstrom, C., 1995. Radium isotopes in the Amazon mixing zone. Geochimica et Cosmochimica Acta 59, 4285-4298. Nyman, D.J. and Fayard, L.D. 1978. Ground-water resources of Tangipahoa and St. Tammany Parishes, southeastern Louisiana. Louisiana Department of Transportation and Development, Office of Public Works water Resources Technical Report 15, 76 p. Oberdorfer, J. A. 2003. Hydrogeologic modeling of submarine groundwater discharge: comparison to other quantitative methods. Biogeochemistry 66, 159-169. Oberdorfer, J. A., Valentino, M.A., and Smith, S.V. 1990. Groundwater contribution to the nutrient budget of Tomales Bay, California. Biogeochemistry 10, 199-216. Oldenburg, C. M. and Pruess, K. 1995. Dispersive transport dynamics in a strongly coupled groundwater-brine flow system. Water Resources Research 31 2, 289-302. Oswald, S. E. and Kinzelbach, W. 2004. Three-dimensional physical benchmark experiments to test variable-density flow models. Journal of Hydrology 290, 22-42. Paytan, A., Boehm, A.B., and Shellenbarger, G.G. 2004. Bacterial contamination and submarine groundwater discharge - A possible link. Environmental Chemistry 1, 29-30. Person, M., Raffensperger, J.P., Ge, S., and Garven, G. 1996. Basin-scale hydrogeologic modeling. Reviews of Geophysics 34 (1), 61-87. Pettijohn, R.A. 1996. Geochemistry of ground water in the Gulf Coast aquifer systems, south-central United States. U.S. Geological Survey Water-Resources Investigations Report 96-4107, 158 p. Price, R. M., Top, Z., Happell, J.D., and Swart, P.K. (2003). Use of tritium and helium to define groundwater flow conditions in Everglades National Park. Water Resources Research 39 (9), 1267, doi:10.1029/2002WR001929, 2003. 222 Prudic, D.E. 1991. Estimates of hydraulic conductivity from aquifer-test analyses and specific-capacity data, gulf coast regional aquifer systems, south-central United States. U.S. Geological Survey Water-Resources Investigations Report 90-4121, 38 p. Purkl, S. and Eisenhauer, A. 2004. Determination of radium isotopes and 2 2 2Rn in a groundwater affected coastal area of the Baltic Sea and the underlying sub-sea floor aquifer. Marine Chemistry 87, 137-149. Raloff, J. 2004. Dead waters: Massive oxygen-starved zones are developing along the world's coasts. Science News Online 165 (23), 360. Retrieved October 19, 2004, from http://www.sciencenews.org/articles/20040605/ bob9.asp. Rama and Moore, W.S. 1996. Using the radium quartet for evaluating groundwater input and water exchange in salt marshes. Geochimica et Cosmochimica Acta 60 (23), 4645-4652. Ranganathan, V., and Hanor, J.S. 1988. Density-driven groundwater flow near salt domes. Chemical Geology 74, 173-188. Rollo, J.R. 1960. Ground water in Louisiana. Louisiana Department of Public Works Water Resources Bulletin 1, 84 p. Rollo, J.R. 1966. Ground \-water resources of the greater New Orleans area, Louisiana. Louisiana Department of Public Works Water Resources Bulletin 9, 69 p. Schincariol, R.A., Schwartz, F.W., and Mendoza, CA. 1994. On the generation of instabilities in variable density flow. Water Resources Research 30 (4), 913-927. Schotting, R. J., Moser, H, and Hassanizadeh, S.M. 1999. High-concentration-gradient dispersion in porous media: experiments, analysis and approximations. Advances in Water Resources 22 (7), 665-680. Segol, G., Pinder, G.F., Gray, W.G. 1975. A Galerkin-fmite element technique for calculating the transient position of the saltwater front. Water Resources Research 11,343-347. Sharp Jr., J. M., Fenstemaker, T.R, Simmons, C.T., McKenna, T. E., and Dickinson, J.K. 2001. Potential salinity-driven free convection in a shale-rich sedimentary basin: Example from the Gulf of Mexico basin in south Texas. AAPG Bulletin 85 (12), 2089-2110. Shaw, R. D. and Prepas, EE. 1989. Anomalous, short-term influx of water into seepage meters. Limnology and Oceanography 34, 1343-1351. 223 Shaw, R. D. and Prepas, E.E. 1990. Groundwater-lake interactions. I. Accuracy of seepage meter estimations of lake seepage. Journal of Hydrology 119, 105-120. Shaw, T. J., Moore, W.S., Kloepfer, J., and Sochaski, M.A. 1998. The flux of barium to the coastal waters of the southeastern USA: The importance of submarine groundwater discharge. Geochimica et Cosmochimica Acta 62 (18), 3047-3054. Simmons, G. M. and Netherton, J. 1986. "Groundwater discharge in a deep coral reef habitat: Evidence for a new biogeochemical cycle?" In: Proceedings of the American Academy of Underwater Sciences Sixth Annual Scientific Diving Symposium, Tallahassee, Fl. Simmons, C. T., Narayan, K.A., Wooding, R.A. 1999. On a test case for density-dependent groundwater flow and solute transport models: The salt lake problem. Water Resources Research 35 (12), 3607-3620. Simpson, M. J. and Clement, T.P. 2003. Theoretical analysis of the worthiness of Henry and Elder problems as benchmarks of density-dependent groundwater flow models. Advances in Water Resources 26, 17-31. Smith, L., and Zawadzki, W. 2003. A hydrogeologic model of submarine groundwater discharge: Florida intercomparison experiment. Biogeochemistry 66, 95-110. Taniguchi, M. and Fukuo, Y 1993. Continuous measurements of ground-water seepage using an automatic seepage meter. Ground Water 31, 675-679. Taniguchi, M, Burnett, W.C., Cable, J.E., and Turner, J.V. 2002. Investigation of submarine groundwater discharge. Hydrological Processes 16, 2115-2129. Taylor, A., Hulme, P., Hughes, A., and Foot, S. 2001. Benchmarking of variable model codes against Henry's problem. In: First International Conference on saltwater intrusion in coastal aquifers - modeling, monitoring, and management. Essaouira, Morrocco, April 23-25, 2001. Therien, R., Sudicky, E.A., McClaren, R.G. 2003. FRAC3DVS: An efficient simulator for three-dimensional saturated-unsaturated groundwater flow and density-dependent, chain-decay solute transport in porous, discretely-fractured porous, or dual-porosity formations, Groundwater Simulations Group, 158 p. Top, Z., Brand, L., Corbett, D.R., Burnett, W.C., and Chanton, J. 2001. Helium and radon as tracers of groundwater input into Florida Bay. Journal of Coastal Research 17, 859-868. Turcan, A.N., Jr., Wesselman, J.B., and Kilburn, C. 1966. Interstate correlation of aquifers, southwestern Louisiana and southeastern Texas. U.S. Geological Survey Professional Paper 550-D, p D231-D236. 224 Turner, R. E., Rabalais, N.N., Swenson, E.M., Kasprzak, M., and Romaire, T. 2005. Summer hypoxia in the northern Gulf of Mexico and its prediction from 1978 to 1995. Marine Environmental Research 59, 65-77. Voss, C. I. 1984. Saturated-Unsaturated TRAnsport. U.S. Geological Survey Water-Resources Investigations Report 84-4369. Voss, C. I. and Souza, W.R. 1987. Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater-saltwater transition zone. Water Resources Research 23 (10), 1851-1866. Watson, S. J., Barry, D.A., Schotting, R.J., and Hassanizadeh, S.M. 2002. Validation of classical density-dependent solute transport theory for stable, high-concentration-gradient brine displacements in coarse and medium sands. Advances in Water Resources 25, 611-635. Weatherill, D., Simmons, C.T., Voss, C.I., and Robinson, N.I. 2004. Testing density-dependent groundwater models: two-dimensional steady state unstable convection in infinite, finite and inclined porous layers. Advances in Water Resources 27, 547-562. Webster, I. T., Hancock, G.J., and Murray, A.S. 1994. Use of radium isotopes to examine pore water exchange in an estuary. Limnology and Oceanography 39 (8), 1917-1927. Weiss, J. S. 1992. Geohydrologic units of the Coastal Lowlands aquifer system, South-Central United States. U.S. Geological Survey Professional Paper 1416-C, 32 p. Weiss, J. S. and Williamson, A.K. 1985. Subdivision of thick sedimentary units into layers for simulation of ground-water flow. Ground Water 23 (6), 767-774. Whiteman, C. D. J. 1979. Saltwater encroachment in the "600-foot" and "1500-foot" sands of the Baton Rouge area, Louisiana, 1966-78, including a discussion of saltwater in other sands. Louisiana Department of Transportation and Development, Office of Public Works Water Resources Technical Report No. 19, 49 p. Williamson, A. K., Grubb, H.F. and Weiss, J.S. 1990. Ground-water flow in the Gulf Coast aquifer systems, south-central United States - A preliminary analysis. U.S. Geological Survey Water-Resources Investigations Report 89-4071, 124 p. Williamson, A. K. and Grubb, H.F. 2001. Ground-water flow in the Gulf Coast aquifer systems, south-central United States. U.S. Geological Survey Professional Paper 1416-F, 173 p. 225 Wilson, A. M. 2003. The occurrence and chemical implications of geothermal convection of seawater in continental shelves. Geophysical Research Letters 30 (21), 2127-2130. Wilson, T. A. and Hosman, R.L. 1988. Geophysical well-log data base for the Gulf Coast aquifer systems, south-central United States. U.S. Geological Survey Open-File Report 87-677, 213 p. Wooding, R.A., Tyler, S.W., White, I., Anderson, P.A. 1997. Convection in groundwater below an evaporating salt lake. 2. Evolution of fingers or plumes. Water Resources Research 33 (6), 1219-1228. Zekster, I. S. and Dzhamalov, R.G. 1981. Groundwater discharge into the Pacific Ocean. Hydrological Sciences Bulletin 26 (3), 271-279. Zekster, I. S., Ivanov, V.A., and Meskheteli, A.V. 1973. The problem of direct groundwater discharge to the seas. Journal of Hydrology 20 (1), 1-36. 226 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052556/manifest

Comment

Related Items