“EXPLORING THE NATURE OF FREE CONVECTION IN A SABKHA WITH
ELECTRICAL IMAGING AND HYDROLOGICAL MODELING”
By
Brian Patrick Eustice

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements

for the degree of
MASTER OF SCIENCE
Geological Sciences
2011

ABSTRACT
“EXPLORING THE NATURE OF FREE CONVECTION IN A SABKHA WITH
ELECTRICAL IMAGING AND HYDROLOGICAL MODELING”
By
Brian Patrick Eustice
Free convection is an important process that has applications in many disciplines of earth
science and is a factor in a variety of environmental issues. Although the study of free
convection has been historically limited to theory, laboratory experiments, numerical modeling,
and secondary field evidence, recent studies have used electrical resistivity geophysical methods
to capture direct evidence of this phenomenon in a natural field setting. However, there is still
very little known about the transient characteristics of free convection in a field setting.
This thesis uses geophysics and hydrological modeling to investigate the occurrence,
initiating mechanisms, and transient nature of free convection in a coastal sabkha in the United
Arab Emirates. The electrical resistivity data shows that free convection in the sabkha is an
episodic process that is likely caused by the dissolution of the halite crust by precipitation and its
subsequent infiltration to the aquifer. Hydrological modeling of the system supports
precipitation as the initiating mechanism and allows for the characterization of the growth and
decay of the convective fingers.

ACKNOWLEDGMENTS

First and foremost, I would like to show sincere gratitude to my advisor, Remke Van Dam, for
his guidance, support, criticalness, and understanding which helped me produce a thesis that I am
proud of. I am genuinely grateful to the rest of committee, David W. Hyndman and Warren W.
Wood, for their guidance, encouragement, and assistance throughout the entire project. I would
like to give special thanks to Wade Kress and Dave Clark of the USGS, without whom the field
work and continued data collection would not have been possible. I would also like to give
special thanks to Craig Simmons for his modeling insights and advice and to Steve Hamilton and
David Weed of the Kellogg Biological Station for their analyses of my water samples. I am
indebted to Kaya Diker, Anthony Kendall, Mine Dogan, and Bobby Chrisman for helping in the
construction of my field equipment. Special thanks to Kaya for help with creating command
files and coding that saved me weeks of time and to Anthony for his thought-provoking
discussions and ideas - I am truly grateful to you both. I would like to acknowledge the rest of
the hydrogeology lab: Emily Luscz, Sherry Martin, Blaze Budd, Briana Jasinski, Lon Cooper,
Chelsea Mack, Ryan Nagelkirk, Erin Haacker, and Yuteng Ma for providing friendship, laughter,
and a valuable distraction when needed. I owe my deepest gratitude to my parents, Pat and
LeAnn Eustice, for their continuous love and support over the last six and a half years while I
studied “rocks.” Thanks everyone for believing in me.

iii

TABLE OF CONTENTS

LIST OF TABLES .........................................................................................................................vi
LIST OF FIGURES ......................................................................................................................vii
LIST OF SYMBOLS .....................................................................................................................ix
INTRODUCTION ..........................................................................................................................1
Free Convection ..................................................................................................................1
Previous Studies ..................................................................................................................2
Introduction .............................................................................................................2
Okavango Delta ......................................................................................................2
Padre Island National Seashore ..............................................................................3
The United Arab Emirates ......................................................................................4
Objective and Approach .....................................................................................................4
STUDY AREA ...............................................................................................................................5
Introduction .........................................................................................................................5
Groundwater Properties ......................................................................................................6
Climate ................................................................................................................................8
Characterization of the Capillary Zone ...............................................................................9
Methods ...................................................................................................................9
Physical and Electrical Properties .........................................................................10
Hydraulic Properties .............................................................................................11
CONCEPTUAL MODELS FOR CREATING A DENSITY INVERSION ................................15
Introduction .......................................................................................................................15
Precipitation ......................................................................................................................17
Evapoconcentration ...........................................................................................................18
Upwelling ..........................................................................................................................18
GEOPHYSICS ..............................................................................................................................19
Methods .............................................................................................................................19
Results ...............................................................................................................................20
1D Profile Analysis ...........................................................................................................22
2008 Data Analysis ...........................................................................................................22
2009 and 2010 Data Analysis ...........................................................................................22
Discussion .........................................................................................................................23
MODELING .................................................................................................................................26
Introduction .......................................................................................................................26
Governing Equations ........................................................................................................26
Conceptual Models ...........................................................................................................28
iv

Precipitation Induced ............................................................................................28
Upwelling Induced ................................................................................................29
Simulation Results ............................................................................................................30
Precipitation Induced ............................................................................................30
Upwelling Induced ................................................................................................34
DISCUSSION ...............................................................................................................................34
CONCLUSIONS ..........................................................................................................................36
APPENDIX A ...............................................................................................................................38
Preprocessing ....................................................................................................................38
Testing for and Removing Artifacts .....................................................................38
Inversion Settings ..............................................................................................................46
APPENDIX B ...............................................................................................................................47
Error Analysis ...................................................................................................................47
APPENDIX C ...............................................................................................................................49
2010 Sabkha Stress Test Experiment ................................................................................49
REFERENCES .............................................................................................................................52

v

LIST OF TABLES

Table 1. Water budget, adapted from Sanford and Wood (2001). ................................................. 8
Table 2. Measured heads and times from a single-ring falling-head infiltration experiments. ... 12
Table 3. Layer thickness and hydraulic conductivity to estimate travel time through the capillary
zone. .............................................................................................................................................. 15
Table 4. Ion chromatography analysis results in ppm for 5 field samples collected from depths
between 1.6 and 2.5 m (sampling sites are shown in Figure 20). ................................................. 18
Table 5. Parameters, values, and sources for COMSOL modeling. ............................................. 30

vi

LIST OF FIGURES

Figure 1. Location of the study area. Regional map (a), geophysical research site, water well and
weather station locations (b), and the layout of geophysical survey lines (c) (adapted from Van
Dam et al., 2009). For interpretation of the references to color in this and all other figures, the
reader is referred to the electronic version of this thesis. ............................................................... 6
Figure 2. Annual precipitation (mm) based on hydrological year (Oct. 1- Sept. 30) from 19832010 from WMO Weather Station 41217 near Abu Dhabi. ........................................................... 9
Figure 3. Cross section of trench wall (tape measure scale in cm) combined with a plot of
apparent resistivity and hydraulic conductivity vs. depth with the 6 layers identified. ................ 10
Figure 4. Picture of single-ring falling-head infiltrometer and graphs of cumulative infiltration
with time for depths of 0, 12, 22.5, 30, 55, and 84 cm. Linear regression lines and equations
given for steady-state conditions. ................................................................................................. 14
Figure 5. Illustration of dipole-dipole array, showing electrode locations, spacings, and image
points. ............................................................................................................................................ 20
Figure 6. Inverted resistivity models for data from lines A, B, and C collected in 2008, 2009, and
2010 (a). 1D profiles of spatially averaged resistivity with depth for lines A (b), and B (c). The
area used to generate the 1D profiles is marked with black dashed lines. See Fig. 1 for locations.
....................................................................................................................................................... 21
Figure 7. Graph of air temperature, daily precipitation, and data collection periods from 20072011............................................................................................................................................... 25
Figure 8. Simulation of groundwater convection from 50mm of precipitation at time t=0, t=2,
t=5, t=10, t=15, t=20, t=25, and t=30 weeks................................................................................. 32
Figure 9. The number of fingers in the model over time as smaller fingers coalesce into fewer
and larger fingers. ......................................................................................................................... 33
Figure 10. Convective velocity (Uc) and concentration (c) for a single saline finger (x=21.5)
from model (Figure 8) over time. ................................................................................................. 33
Figure 11. Simulation of groundwater convection caused by upwelling of less dense water at
times t=0, t=15, t=30, t=31, t=32, and t=33 years. ....................................................................... 34
Figure 12. Comparison between 2008 field data for line A (a) and COMSOL model output at
time t=15weeks (b). ...................................................................................................................... 35

vii

Figure 13. Cross section of inverted resistivity for line B using inversion settings from Van Dam
et al. (2009) showing layer of low resistivity below the capillary zone (a) and 1D profile (b)
showing average resistivity (Ohm•m) with depth. See Fig. 6 for area used in averaging. .......... 39
Figure 14. Cross sections of 2-layer (a) and 3-layer (b) models. Both models have a background
resistivity of 0.18 Ohm•m, a capillary zone with a resistivity of 0.71 Ohm•m. The 3-layer model
has an intermediate layer of low resistivity (0.11 Ohm•m) just below the capillary zone. .......... 40
Figure 15. Cross sections of inverted resistivity outputs from the 2-layer model (a) and the 3layer model (b) both showing a zone of low resistivity below the capillary zone. 1D profile (c)
showing average resistivity for middle 40 m of 2- and 3-layer inverted models. See Fig. 6 for
area used in averaging. See Fig. 14 for the starting model. ......................................................... 41
Figure 16. Cross sections of inverted resistivity for 2 layer (a) and 3 layer (b) models for mesh
settings of 2 and 8 divisions between electrodes respectively. 1D profiles of spatially average
resistivity for 2 layer (c) and 3 layer (d) models. See Fig. 6 for area used in averaging. ............ 43
Figure 17. 1D profiles of the two different starting models showing each model’s influence on
the resistivity at 0 m depth. See Fig. 6 for area used in averaging. ............................................. 45
Figure 18. Pseudosection of measured apparent resistivity (a) and Inverted resistivity cross
section (b) for the 2008 data line b-b’. .......................................................................................... 47
Figure 19. Pseudosection of calculated apparent resistivity (forward calculation) with 2 bad
electrodes (a) and modeled inverted resistivity cross section (b). ................................................ 48
Figure 20. Map of the study area showing the 2010 electrode arrays, infiltration zone, and
pressure transducer (background image from Google Maps taken 11/7/2010.) ........................... 50

viii

LIST OF SYMBOLS

Symbol
q
K
dh
dl

Parameter
Water flux
Hydraulic conductivity
Change in head
Change in length

g
k
β

Acceleration of gravity
permeability
Linear expansion coefficient

∆C
H
θ
v0

Change in concentration
Aquifer thickness
Porosity

μ0

Dynamic viscosity

ρ
ρ0

Density

kg/m•s
3
kg/m

Initial density

kg/m

ρs

Infiltrating density

kg/m

ρu

Upwelling density

kg/m

∆ρ
DL

Change in density

kg/m

Coefficient of molecular diffusion

m /s

DD
u
p
D
Qm

Dispersion coefficient
Darcy velocity
Fluid pressure
Vector in direction in which gravity acts

m /s
m/s
Pa
m
3
kg/(m •s)

θs

Fluid volume fraction

c
c0

Concentration

3
kg/m

Initial Concentration

kg/m

cs

Infiltrating Concentration

kg/m

cu
τ
L
qu
Uc
α

Upwelling Concentration
Tortuosity
Aquifer thickness

kg/m
m

Upwelling flux rate
Convective velocity
Dispersivity

mm/yr
m/d
m

Kinematic viscosity

Mass source term

ix

Unit
m/d
m/d
m
m
2
m/s
m2
3
kg/m
m
2
m /s

3
3
3
3

2
2

3
3
3

INTRODUCTION

Free Convection
Free convection is the mixing of fluid driven by buoyancy forces. It occurs when
temperature and salinity related changes create a density inversion that causes the system to
become critically unstable. A density inversion is a light fluid overlain by a dense fluid. When
the contrast between the fluids is great enough, the less dense fluids rise and higher density fluids
sink at rates faster than they diffuse. Free convection appears in nature in the form of convection
cells or as descending fingers of dense fluid.
Free convection is important because it causes mixing of solutes at faster rates than
advection and diffusion, over larger spatial scales, and involving larger volumes of materials. It
applications in many disciplines including: atmospheric sciences, limnology, oceanography,
geophysics, hydrogeology, and astrophysics (Diersch and Kolditz, 2002). Movement by free
convection plays an important role in a variety of environmental areas of concern such as
seawater-groundwater interactions, contaminant migration, flow through formations at
radioactive waste disposal sites, and the salinization of groundwater in arid and semi-arid regions
(Simmons, 2005; Bauer et al., 2006).
Most free convection has been studied through theory (Rayleigh, 1916; Horton and
Rogers, 1945; Lapwood, 1948), laboratory experiments (Schneider, 1963; Elder 1967), and
numerical modeling (Wooding, 1957; Elder 1967). Identification of free convection in field
settings has been limited to secondary evidence such as: unusual tritium distributions (Wood et
al., 2002), measured density inversions in groundwater, reversals in hydraulic gradients, high
fluid fluxes, theoretical considerations (high Rayleigh number), modeling of a system showing

1

convection, and salt deficits in saline lakes (Simmons, 2005). This secondary field evidence,
while supporting that free convection was possible or that mixing was taking place at rates faster
than diffusion and dispersion, did not confirm that free convection was actually occurring.
Direct evidence in a field setting has been identified as one of the missing links in the study of
free convection (Simmons et al., 2005). In addition, the transient characteristics of free
convection in field settings are poorly known. In recent years, electrical geophysical methods
have been used to capture direct field evidence of this phenomenon (Bauer et al., 2006; Van Dam
et al., 2009; Stevens et al., 2009).

Previous Studies
Introduction. There have been three major studies using geophysics to image free
convection in the subsurface. In each study, free convection was driven by density differences
caused by variations in solute concentrations and the fingering was imaged using electrical
resistivity. Electrical resistivity measurements were employed based on the relationship between
salinity variations, fluid electrical conductivity, and bulk electrical resistivity.

Okavango Delta. Free convection has been hypothesized as a possible cause of a salt
deficit in the surface water of the Okavango Delta of Botswana (Gieske, 1996). A majority of
the water in the system is lost to evapotranspiration from islands between the delta channels.
This evapotranspiration leads to the accumulation of solutes near the surface of the freshwater
aquifer creating a density inversion (Bauer et al., 2006).
Bauer et al. (2006) used a combination of surface and borehole electrical resistivity
measurements to map the distribution of salinity under two islands in the delta. Under one island,
they were able to capture evidence of a single convection finger. Through hydrodynamic
2

modeling they determined that the instabilities required for free convection would take 100 to
150 years of evapoconcentration and that the saline finger they imaged was likely over 200 years
old (Zimmerman et al., 2006).
This was the first field evidence capturing an apparent saline finger using electrical
resistivity. Hydrodynamic modeling was consistent with their geophysical findings. However,
with the development of fingering in the Okavango Delta taking 100’s of years to occur, testing
and calibrating models that simulate the transient nature of free convection is nearly impossible.

Padre Island National Seashore. The tidal flats on Padre Island, Texas, have a low
hydraulic gradient, undergo evapoconcentration, and are occasionally inundated by saline water
from Laguna Madre making them a likely location for free convection to occur (Stevens et al.,
2009). Stevens et al. (2009) used a variety of secondary evidence such as hydraulic gradients,
sampling data showing occasional density inversions, and theoretical calculations to show the
potential for free convection at their field site.
A 3D electrical resistivity survey was used to capture evidence of saline fingering
beneath Padre Island. Stevens et al. (2009) employed a surface electrode array with a 14 X 6
grid of electrodes with 2 meter spacing.
The secondary evidence presented by Stevens et al. (2009) all suggests that free
convection is indeed taking place. However, the plumes of high salinity and density inversions
that they identified in their inverted 3D resistivity results all occur along the edges and corners of
their data plots, where because of the way ERI is measured, these features must have been
extrapolated to areas where there are little or no data points.

3

The United Arab Emirates. Wood et al. (2002) suggested that density driven
convective flow was taking place in a sabkha in the United Arab Emirates (UAE) based on
measured density inversions and a distribution of elements consistent with vertical convective
mixing. The sabkha makes an ideal natural laboratory for the study of free convection in a field
setting because the aquifer sands are very homogeneous, meaning any variation measured in
resistivity will be due to variations in solute concentration.
Van Dam et al. (2009) used electrical resistivity imaging (ERI) at the sabkha to document
the existence of fingering. They used two perpendicular 2D surface arrays with 84 electrodes
each with 1 meter spacing and one array with 84 electrodes with a 0.5 m electrode spacing.
Van Dam et al. (2009) were able to, in contrast to the other two studies, image a pattern
of fingering over a large spatial scale similar to patterns shown in laboratory and modeling
studies of free convection. This study, however, is limited in the fact that it is just a single
snapshot in time and tells us very little of the characteristics of the processes involved in free
convection.

Objective and Approach
The objective of this thesis is to investigate the occurrence and transient nature of free
convection in the UAE coastal sabkha. It will be testing the hypothesis suggested in Van Dam et
al. (2009) that free convection in the sabkha is episodic by nature and caused by precipitationinduced dissolution of the halite crust and subsequent infiltration to the aquifer.
The first part of the thesis is an introduction to the study area and is accompanied by a
description of the setting, climate, and a detailed characterization of the physical, hydraulic, and
electrical properties of the capillary zone. It will then compare the data from Van Dam et al.

4

(2009), which captured evidence of free convection in 2008 to data collected during the
following two years at the same site to test if free convection in the sabkha is episodic in nature
rather than continuous based on proposed models for creating the density inversion. The next
section will employ a hydrologic model based on measured field values to characterize the
transient nature of free convection. The hydrologic model will be used for testing the hypothesis
of the mechanism for the onset of convective fingering, the speeds at which the fingers grow, and
the rates at which they decay.

STUDY AREA

Introduction
A sabkha, which means “salt flat” in Arabic, is a supratidal environment where evaporite
minerals precipitate displacively within the sediments in the capillary zone (Boggs, 2001). The
study area is a sabkha located about 40 km southwest of the city of Abu Dhabi along the coast of
the Arabian Gulf in the United Arab Emirates (Figure 1). The field site consists of ~11.5 m of
sand from sand dunes reworked by sea level rise underlain by Miocene carbonates. The sands
are essentially homogeneous (Wood et al. 2002), except for in the top meter of the sabkha, which
is characterized by authigenic evaporites. The topography of the sabkha is relatively flat with a
gradient of 1 m per 10 km toward the coast (Sanford and Wood, 2001).

5

Figure 1. Location of the study area. Regional map (a), geophysical research site, water well and
weather station locations (b), and the layout of geophysical survey lines (c) (adapted from Van
Dam et al., 2009). For interpretation of the references to color in this and all other figures, the
reader is referred to the electronic version of this thesis.
Groundwater Properties
The region around the study area contains a series of groundwater wells that extend to
various depths throughout the sabkha aquifer and into the formation below (Figure 1).
Additionally, surface water samples were taken where ponding occurred after a rainfall event
(Wood et al., 2002). The water in the sabkha was found to have an average solute concentration
3

of 286,415 mg/L and a density of 1.1858 g/cm (Wood et al., 2002). The Tertiary formation
below has a TDS (total dissolved solids) of 101,453 mg/L and a density of 1.0753 g/cm
6

3

(densities are for 20°C) (Van Dam et al., 2009). Water and solute budgets for the sabkha have
been calculated by Wood et al. (2002). They determined that the majority of the solutes in the
sabkha aquifer have come from the slow upwelling (around 4 mm/yr) of saline water from the
Tertiary formation below. Due to the low hydraulic gradient the solutes have a residence time in
the system of ~26,000 years. Water in the system, however, has a residence time of only ~50
years. The bulk of the water flux in and out is from evaporation and precipitation. Surface water
3

samples had an average density of 1.2675 g/cm and ranged in solute concentration from
292,487-536,969 mg/L. The high concentration and density of the surface water is due to
dissolution of the dry halite crust (Wood et al., 2002).
The lateral water flux through the sabkha is calculated using Darcy’s law:
q=-K(dh/dl)

(1)

where q is the flux of water (m/d), the negative sign signifies flow from higher to lower
potential, K is the hydraulic conductivity (m/d), and dh/dl is the hydraulic gradient (m/m). Using
a falling-head permeameter, the hydraulic conductivity was measured in the sabkha sands to be
6.57 m/d (compared to 1 m/d measured by Sanford and Wood, 2001) and the hydraulic gradient
is 0.0002 (Sanford and Wood, 2001). Therefore, the lateral flux ranges from 0.073 to 0.48 m/yr.
The seepage velocity is found by dividing the lateral flux by the porosity (0.38) of the aquifer,
which gives a velocity of 1.26 m/yr (~3.5 mm/day).

7

Climate
The UAE is a hot arid sub-tropical desert. The mean minimum and maximum annual
temperatures are 21.7 and 33°C respectively. Long term precipitation values (1983-2010) from a
WMO weather station near the study area indicate that the average annual precipitation is around
32 mm/yr, mostly from November to April, with a large standard deviation since it is not
uncommon to see multiple years with no rain (Figure 2). The pan evaporation rate of fresh water
in the area is 3500mm/yr (Bottemley, 1996); however, evaporation rates from the surface of the
sabkha were measured by Sanford and Wood (2001) and found to be less than 3% of pan
evaporation (Table 1). Transpiration is not occurring because there is no vegetation on the
sabkha.
Table 1. Water budget, adapted from Sanford and Wood (2001).
Budget component
Evaporation
Precipitation (recharge)
Lateral flux in
Lateral flux out
Upward leakage
Water table depth

Flux rate (mm/yr)
69
32
80
80
4
1m

Ranges
50-88
0-132
73-480
73-480
4-5
0.5-1.3

8

Sources
Sanford and Wood (2001)
WMO Weather station 41217
Sanford and Wood (2001)
Sanford and Wood (2001)
Sanford and Wood (2001)
Wood et al. (2002); field measurement

180
160

Precipitation (mm)

140
120

100
80
60
40
20
83-84
84-85
85-86
86-87
87-88
88-89
89-90
90-91
91-92
92-93
93-94
94-95
95-96
96-97
97-98
98-99
99-00
00-01
01-02
02-03
03-04
04-05
05-06
06-07
07-08
08-09
09-10

0

Hydrological Year
Figure 2. Annual precipitation (mm) based on hydrological year (Oct. 1- Sept. 30) from 19832010 from WMO Weather Station 41217 near Abu Dhabi.

Characterization of the Capillary Zone
Methods. In 2010, a detailed characterization of the capillary zone was performed. A
trench was dug from the sabkha surface down to the water table (1.3 m) and a detailed
description of the sediment layers was made. The electrical properties of the entire vertical
column were measured in 5 cm increments using a horizontal Wenner array with an a-spacing of
5 cm (Figure 3), and the saturated hydraulic conductivities in the trench were measured for
several depths (Figure 4).

9

Figure 3. Cross section of trench wall (tape measure scale in cm) combined with a plot of
apparent resistivity and hydraulic conductivity vs. depth with the 6 layers identified.

Physical and Electrical Properties. The sediments are classified into 6 distinct layers
(Figure 3). The top 12 cm, layer 1, are characterized by a dry, friable, sand-size crust containing
crystals of gypsum 5-10 mm in diameter. This layer has the highest measured electrical
resistivity in the system ranging from 40-70 Ohm•m. Layer 2 extends from 12 to 45 cm and is
made up of a mix of sand- and clay-size grains and contains gypsum crystals 15-20 mm in
diameter. A red iron oxide band is found at 42-43 cm depth. The electrical resistivity of this
layer ranges from 0.5-1.3 Ohm•m. Layer 3, between 45 and 62 cm depth is comprised of 1-5
10

mm thick horizontal layers of blue clay. Gypsum crystals are found throughout this layer, but are
decreasing in abundance with depth. This clay layer has a very low electrical resistivity of 0.180.36 Ohm•m. The 4th layer is entirely made up of coarse-grained sand-size gypsum crystals,
ranges from 62 to 71 cm depth, and has an electrical resistivity ranging from 0.21-0.38 Ohm•m.
Below that is layer 5, 71-82 cm, which is made up of similar sized gypsum crystals but is set in a
clay matrix containing large, 15-50 mm diameter, soft anhydrite nodules. This layer has an
electrical resistivity between 0.32-0.48 Ohm•m. The material from 82 cm down is the
homogeneous fine grain sand that makes up layer six and the rest of the aquifer, the capillary
portion of which has an electrical resistivity ranging from 0.32-0.48 Ohm•m.
The top of the capillary zone in this system has a diurnal fluctuation, but remains at or
near the surface (Wood, 2002). High resistivities were measured in the top 15 cm where the
sediment is mostly dry during the day. Below that, resistivity gradually drops until around 60 cm
as the sediments become more and more saturated. From 70 cm and down the electrical
resistivity remains fairly stable at values similar to the aquifer, suggesting that this section of the
capillary zone is saturated.

Hydraulic Properties. Horizontal tiers were dug in the trench at a series of depths from
the surface down to 84cm. The vertical saturated hydraulic conductivity was measured at each
tier using a single-ring falling-head infiltrometer (Figure 4). The infiltrometer consisted of a
section of pipe 15 cm in diameter with a graduated scale on the side in centimeters. The
infiltrometer was carefully inserted 5 cm into the sediment, to minimize disturbance, and then
filled with 1.5 L of sabkha water. The time was recorded at every 0.5 cm of head loss (Table 2).
The data was then plotted as cumulative infiltration (cm) per time (Figure 4). Once the

11

infiltration rate is steady, the sediment is assumed to have become saturated. A trend line was
added to the linear (steady state) section of the plots. The slope of the trend line represents the
field saturated hydraulic conductivity of the sediments in cm/s which was then converted to m/d
(Hillel, 1982). This procedure may lead to overestimates of vertical saturated hydraulic
conductivity due to the sediment cracking and creating preferential pathways when the
infiltrometer was inserted (although care was taken to prevent this) and because of non-vertical
flow below the base of the infiltrometer at 5 cm depth.
Table 2. Measured heads and times from a single-ring falling-head infiltration experiments.
Measured Time (s)
at depth:
Measured
Head (cm)

84 cm

55 cm

30 cm

22.5 cm

12 cm

0 cm

8.5
8
7.5
7
6.5
6
5.5
5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0

13
63
128
191
258
314
388
447
515
584
652
726
778
862
934
999
1047
-

0
85920
103980
170580
250979
276959
-

18
49
91
134
175
219
271
324
379
428
493
535
614
680
754
804
-

16
30
43
59
77
97
114
134
153
171
195
219
244
259
-

14
17
21
24
28
32
36
40
44
49
52
60
-

16
20
23
26
28
30
33
35
37
39
40
-



Tests at the sabkha surface (0 cm) and in the loose, dry material below (12 cm) revealed a
very high hydraulic conductivity of 204.42 m/d and 99.71 m/d, respectively. These high

12

conductivities are at least in part due to the dissolution of the soluble minerals
encountered by the infiltrating fluid.


The zone from ~0.15 m to 0.8 m depth is divided into smaller zones with distinctly
different hydraulic conductivity. The upper part consists of two zones with conductivities
of 21.25 m/d (22.5 cm) and 7.26 m/d (30 cm).



The horizontally-layered clay from 45-62 cm appeared a strong barrier to water flow. The
estimated hydraulic conductivity of 0.01 m/day implies that infiltrating rainwater would
take up to 25 days to pass this unit. If this layer is laterally continuous (it was observed in
a boring ~100 m northeast of the trench (Figure 20)) it would strongly impede
infiltration.



The hydraulic conductivity measured in the sabkha sand (84 cm) is 6.57 m/d.

13

Figure 4. Picture of single-ring falling-head infiltrometer and graphs of cumulative infiltration
with time for depths of 0, 12, 22.5, 30, 55, and 84 cm. Linear regression lines and equations
given for steady-state conditions.

Using the sediment layer thicknesses and the measured hydraulic conductivities, it is
possible to estimate the length of time it would take infiltrating water to travel from the surface
of the sabkha down to the water table (Table 3). Values in the table are approximate because the
depths of hydraulic conductivity measurements do not all correspond exactly with sedimentary
layer boundaries (Figure 3). It is estimated that a molecule of infiltrating water will take around

14

593 hours (~25 days) to reach the water table. This value is dominated almost entirely by the
very low hydraulic conductivity of the clay layer from 45 to 62 cm depth.
Table 3. Layer thickness and hydraulic conductivity to estimate travel time through the capillary
zone.
Measurement Depth
(cm)
0
12
22.5
30
55
84
Total

Approximate Layer
Thickness (cm)
12
10.5
7.5
15
17
68
130

Measured K
(m/d)
204.42
99.71
21.25
7.26
0.01
6.57

Estimated Travel Time
(hours)
0.014
0.025
0.085
0.496
590.3
2.49
593.4

CONCEPTUAL MODELS FOR CREATING A DENSITY INVERSION

Introduction
Three major conceptual models have been proposed for the formation of a density
inversion in the UAE sabkha. The first involves a precipitation-induced dissolution of the halite
crust, which infiltrates through the capillary zone into the aquifer. The second proposes that
evapoconcentration leads to a buildup of a high salinity layer near the surface of the aquifer (Van
Dam et al., 2009). The third suggests that the slow upwelling of low-density water from the
Miocene formation below the aquifer creates the inversion.
A major factor in these conceptual models is that they must create a large enough density
inversion to cause the system to become unstable causing free convection to occur. Stability of
layered fluids is calculated using the Rayleigh number; a dimensionless value describing the

15

ratio of buoyancy driven flow to dispersive and viscous forces (Rayleigh, 1916). The Rayleigh
number is calculated using the equation:

(2)
where:
2

g is the gravitational acceleration (m/s );
2

k is the permeability (m );
β is the linear expansion coefficient for density;
∆C is the difference in concentration by mass fraction = Cmax-Cmin (kg/kg);
H is the thickness of aquifer (m);
θ is the porosity (dimensionless);
2

DL is the fluid coefficient of molecular diffusion (m /s).
2

v0 is the kinematic viscosity = µ0/ρ0 (m /s);
µ0 is the dynamic viscosity (kg/m•s);
3

ρ0 is the underlying water density (kg/m );

2

To evaluate the stability of a system, a critical Rayleigh number Rac=4π was derived for
an infinite horizontal layer with uniform porosity, no flow, and constant temperature boundaries
(Lapwood, 1948). When a system’s Rayleigh number is below Rac, it is considered stable and

16

diffusive forces control solute movement. A system with a Rayleigh number higher than Rac is
likely to be unstable as buoyancy forces dominate and free convection occurs.

Precipitation
Density inversions may be created in the sabkha through the dissolution of the halite crust
by rain water and the subsequent infiltration to the water table. Concentrations differences
measured between ponded surface water (rain water that has dissolved the halite crust) and the
5

aquifer were found as high as 250,000 mg/L (Wood et al., 2002) with a Ra=2.5x10 , a value four
orders of magnitude higher than the critical value for convection. For this model to work,
however, there needs to be enough halite in the system for the rainwater to dissolve so that its
concentration and density become high enough to create an inversion.
Ion chromatography measurements of water samples collected in 2010 show that sodium
and chloride are the dominant ions in solution, making halite the most common salt to precipitate
at the surface (Table 4). Sodium, the limiting ion for halite in the sabkha, has a concentration of
92,802 mg/L. If we take a 1 m by 1 m area of the surface, based on the evaporative flux rate
measured to be between 50-88 mm/yr by Sanford and Wood (2001), there will be 50-88 L of
water evaporated per year. The amount of halite precipitated can be estimated by taking 50 L X
92,802 mg/L = 4.64 kg Na (8.17 kg Na using 88 L of water). The equivalent mass of
precipitated Cl will be 7.16 kg (12.6 kg). Dividing the mass of precipitated NaCl by the density
3

3

of 2.16 g/cm , we end up with 5461.0 - 9611.3 cm of halite per year in the 1x1 m square, which
is equivalent to a deposition rate of 5.48-9.64 mm/yr. Using the Pitzer equation in PREEQC
(Parkhurst and Appelo, 1999), the estimated deposition rates, and the solubility of halite, there

17

would need to be between 32-57 mm/yr of precipitation to dissolve the halite at the same rate as
it is being deposited. This indicates that with an average rainfall of 32 mm/yr there should be a
surplus of dissolvable halite in the system of up to 4.5 mm/yr.
Table 4. Ion chromatography analysis results in ppm for 5 field samples collected from depths
between 1.6 and 2.5 m (sampling sites are shown in Figure 20).
Sample

Cl

Br

2-B
3-B
4-B
5-B
6-B
Average

149786
144466
145117
145748
146679
146359

130
124
126
133
126
128

NO3N
40
39
39
38
39
39

SO4
1773
1913
1947
1979
1928
1908

Ca

Mg

Na

K

3702
3638
3654
3734
3663
3678

8116
7834
7793
7944
7886
7915

94166
92110
92156
93264
92315
92802

3605
3455
3430
3512
3487
3498

Evapoconcentration
The second conceptual model for the development of density inversions in this system is
evapoconcentration. As the water is being evaporated from sabkha, its salts are being left
behind, which in turn increase the concentration of the fluids near the top of the aquifer. This
progressive increase of concentration near the top of the aquifer, if occurring at a faster rate than
diffusive mixing, could eventually cause the system to reach a critical concentration gradient and
induce convection (Nield et al., 2008).

Upwelling
The third conceptual model for the formation of density inversions is vertical transport of
less dense fluid into the aquifer from the formation below. Water with a concentration difference
of approximately 185,000 mg/L is upwelling from the Tertiary formation below the aquifer at a

18

rate of ~4 mm/yr. The Rayleigh number for the sabkha and the underlying formation is 1.92x10

5

2

and is well above the Rac of 4π ≈39.5 and considered unstable enough for convection to occur.

GEOPHYSICS

Methods
Electrical resistivity imaging is a proven method for detecting variations in groundwater
salinity. Electric current is easily conducted through the pore water by electrolytic conduction
(Reynolds, 1997). Conductivity of the pore water, the electrolyte, is related to its total dissolved
solids (TDS) concentration. The more dissolved ions in solution the higher the conductivity.
Since electrical resistivity is the inverse of conductivity, the higher the conductivity the lower the
resistivity. The differences in TDS for the precipitation and upwelling conceptual models have
5

been measured at ~10 mg/L. These large variations in salinity would be detected as variations
in electrical resistivity.
The field site was outfitted with two perpendicular electrical resistivity lines, A and B,
with 84 electrodes each, at spacing of 1m. Another 84 electrode line, line C, was set parallel to
line B with an electrode spacing of 0.5 m (Figure 1C). Data were collected using two electrode
configurations, dipole-dipole and pole-pole arrays. Pole-pole surveys, which are relatively fast
but have limited sensitivity to lateral variations in resistivity, were used for vertical soundings.
These soundings confirmed the hydrostratigraphic model for the field site (Van Dam et al.,
2009).

19

For the purpose of this study, the dipole-dipole array (Figure 5) was chosen because of its
high sensitivity to lateral variations in resistivity (Reynolds 1997). A disadvantage of this type
of array is that it has a lower signal to noise ratio, requiring longer measurement times. The
survey was set up using Advanced Geosciences, Inc. (AGI) SS Command Creator. The
command file had a maximum n length of 8 m and a maximum dipole length of 5 m resulting in
2580 image points along each line.

Figure 5. Illustration of dipole-dipole array, showing electrode locations, spacings, and image
points.

Results
Inversions of the electrical resistivity data were performed using Advanced Geosciences
EarthImager2D (Figure 6) Initial inversion settings for the software were tested and optimized
using simplified models of the subsurface resistivity distribution. Initial settings and
preprocessing steps are listed in Appendix A. 1D resistivity profiles were created from the 2D
inversion results of lines A and B by taking the average resistivity of the middle 40m plotted
versus depth (Figure 6b and 6c).

20

Figure 6. Inverted resistivity models for data from lines A, B, and C collected in 2008, 2009, and 2010 (a). 1D profiles of spatially
averaged resistivity with depth for lines A (b), and B (c). The area used to generate the 1D profiles is marked with black dashed lines.
See Fig. 1 for locations.

21

1D Profile Analysis
The 1D profiles illustrate that the average resistivity does not exhibit much change from
year to year. This suggests that the average salinity profile stays fairly constant. The shallow
values have a maximum resistivity of around 0.9-1.0 Ohm•m, which are significantly smaller
than those measured in the soil pit of around 70 Ohm•m. This inconsistency is because the
smallest dipole separation of 1 m in the surface resistivity survey does not have the resolution to
image such shallow features.

2008 Data Analysis
The three inverted dipole-dipole resistivity lines collected in 2008 (Figure 6a) each
display similar features. The top 75-80 cm shows the highest resistivity with values around 1.0
Ohm•m. This layer of high resistivity is representative of the capillary zone. Below the capillary
zone the resistivity sharply drops and vertical structures of alternating resistivity of around 0.1
and 0.3 Ohm•m start to appear. Most of the vertical structures of low resistivity extend to a
depth of around 7 m; however a few reach close to the bottom of the profile. The shallow
vertical structures of low resistivity have a wavelength of 6.5-7 m while the deeper structures
display an irregular pattern with a wavelength between 10-15 m where present.

2009 and 2010 Data Analysis
The inverted electrical resistivity profiles from the 2009 and 2010 field campaigns are all
very similar (Figure 6a). Change in resistivity occurs vertically with depth and there is very little
lateral variation. Like the 2008 data, the layer from the surface to 80 cm has the highest
resistivity and represents the capillary zone. Below that, the resistivity decreases to a value of

22

around 0.25 Ohm•m for ~1 m. The next 5-6 m have the lowest resistivity with an average value
of around 0.19 Ohm•m. From about 7 m depth to the bottom of the profile, the resistivity stays
around 0.25 Ohm•m.

Discussion
The 2008 inverted resistivity models shows significant lateral variations in resistivity
(Figure 6a), which is consistent with unstable fingering. The vertical low resistivity features
appear to be high density fingers that are descending through the aquifer while the high
resistivity features are ascending lighter fluids. There are no obvious vertical features and very
little lateral variation in resistivity in the 2009 and 2010 data (Figure 6a). This lack of lateral
variation in resistivity suggests concentrations in the groundwater are only changing vertically
between strata and convective fingering is not occurring. Analysis of the geophysical data thus
suggests that free convection is not occurring continuously in this system. The episodic nature of
the convection along with the 1D and 2D profiles give some insight into which conceptual model
is creating the density inversion.
Evapoconcentration in an arid region would likely be a continuous process of building up
high concentration brines at the surface of the water table. If evapoconcentration is what triggers
free convection, it is likely that there would be a layer of low resistivity just below the capillary
zone present in the 2008 data or building back up by the 2010 data; however, this is not the case.
Evapoconcentration was initially proposed as a conceptual model because the results published
by Van Dam et al. (2009) displayed a layer of low resistivity below the capillary zone (Figure
13). The layer of low resistivity seen in their data was shown to be an artifact of the inversion
process and does not actually exist (Appendix A). Additionally, the presence of a halite crust

23

and the high measurements of TDS in ponded rain water suggest that the bulk of the salts are
being precipitated out of solution when it is evaporated instead of being left behind in the
aquifer. The combination of the lack of a low resistivity layer and the presence of a halite crust
suggest that evapoconcentration is not a likely mechanism for creating a density inversion.
Upwelling from the formation below the aquifer can be seen in the 1D and 2D profiles as
the gradual increase in resistivity from around 4 m to 12 m depth (Figure 6a,b, and c). This
gradient is present during all three years even though free convection is not, suggesting that
although it may contribute to the density gradient it is likely not the mechanism that initiates
convection.
Dissolution of the halite crust by rainwater that becomes saturated with solutes and
infiltrates to the water table is another possible mechanism for the creation of a density inversion.
Figure 7 shows daily temperature, daily precipitation in mm, and geophysical data collection
periods since 2007. Two major precipitation events of 50 mm and 64 mm took place 18 weeks
and 6 weeks, respectively, before the first data collection period in March 2008. Precipitation
was insignificant in the year prior to the 2009 field campaign. There was a total of 18 mm of
rain with the maximum event being 4 mm. Three large precipitation events of around 42 mm
each occurred 47, 36, and 7 weeks before the third data collection period in 2010.
A conservative halite solute budget based on evaporative flux as described earlier shows
that the two rainfall events before the 2008 data collection and the rainfall event 47 weeks before
the 2010 data collection would be able to dissolve enough halite to make the system unstable.
The other two large events before the 2010 data collection would likely not have enough halite to
dissolve in order to reach a density inversion.

24

Based on the amount of dissolvable halite in the system and the timing of the
precipitation events in relation to the data collection campaigns it is likely that free convection in
the sabkha is precipitation-induced. A rainfall of 50 mm has the potential to become completely
saturated with dissolved halite. Completely saturated rainwater at 27.4 ºC (average air
temperature) would have a sodium concentration of ~140,000 mg/L, which, when compared to
sodium concentrations of 92,802 mg/L measured in the aquifer, illustrates the potential for a
density inversion. Assuming that the rainfall event 18 weeks prior to the 2008 data collection is
responsible for fingering and taking into account the ~25 days for infiltration through the
capillary zone (Table 3), it can be calculated that the fingers had around 15 weeks to develop.

Figure 7. Graph of air temperature, daily precipitation, and data collection periods from 20072011.
25

MODELING

Introduction
The geophysics appears to demonstrate that free convection is occurring in the sabkha,
but also that it is not a continuous process. Because the geophysics only offers snapshots in
time, numerical modeling of the sabkha system was done to investigate the transient nature of
convective fingering. Modeling was performed using COMSOL Multiphysics software
(COMSOL Multiphysics, 2011).

Governing Equations
Fluid flow in a porous media is controlled in this model by applying Darcy’s law:

(3)
where:
u is the Darcy velocity or specific discharge vector (m/s);
p is the fluid pressure (Pa);
μ is the dynamic viscosity (kg/m•s);
2

g is the magnitude of gravitational acceleration (m/s );
∇D is a unit vector in the direction over which the gravity acts (m).
2

κ is the permeability (m ); calculated by

(4)
K is the hydraulic conductivity;
26

3

ρ is its density (kg/m ); calculated by
ρ=ρ0+βc

(5)

3

ρ0 is the initial groundwater density (kg/m );
β is the linear expansion coefficient;

Darcy’s law is then inserted into a continuity equation:

(6)
where:
θ is the porosity;
3

Qm is a mass source term (kg/(m ·s)).

This results in:

(7)

Since the model assumes that solute transport is taking place in a fully saturated medium
and that the solutes are conservative (thus ignoring adsorption and desorption terms) transport
can be modeled using a simplified advection-dispersion equation:

(8)
where:

27

θs is the fluid volume fraction (porosity);
3

c is the salt concentration (kg/m );
u is the Darcy velocity (m/s);
τ is the fluid phase tortuosity factor τ =θ

1/3

;

DD is the dispersion term;
2

DL is the fluid diffusion coefficient (m /d).

Conceptual Models
Precipitation Induced. The sabkha system was modeled in 2D using the hypothesis that
free convection is precipitation-induced. For simplicity, the starting conditions of the model
assume that a rainfall has already become saturated from the dissolution of the halite crust and
has infiltrated through the capillary zone to the top of the water table. The model domain is
divided into two layers. The bottom layer is 10m thick and represents the sandy aquifer. The top
layer’s thickness is dependent on the amount of precipitation that occurred divided by the
porosity (0.38) so that it represents a completely saturated layer that is an extension of the
aquifer. For example, a rainfall event of 0.10 m would result in a top layer with a thickness of
0.26 m. For the Darcy’s law interface, the sides and bottom are considered to be no flow
boundaries. Zimmerman et al. (2006) have modeled the effects of hydraulic gradients on the
formation of convective fingering. Their results demonstrate that low gradients (~0.0001) have
almost no effect on finger development for a site where the concentration difference and a
Rayleigh number are each an order of magnitude smaller than that of the sabkha. Lateral flow in
the sabkha based on a hydraulic gradient of 0.0002 is 1.26 m/yr and is insignificant compared to
28

the velocity of vertical density driven flow and will therefore not be considered in the model.
Additionally, the influence on the system from the upwelling of 4 mm/yr from the formation
below the aquifer is orders of magnitude smaller than that of a large precipitation event and was
therefore ignored. The top boundary is treated as an atmospheric pressure boundary and the
initial pressure of the entire domain is equal to
p0+ρ•g•L

(9)

where:
p0 is the initial pressure (Pa•s);
3

ρ is the calculated density (kg/m );
2

g is the gravitational acceleration (m/s );
L is the domain height (m).

The solute transport calculations treat all four boundaries as no flux boundaries;
therefore, no solutes may be transported across them by advection, dispersion, or diffusion. The
top layer in the domain represents the infiltrating water with the starting concentration cs and
density ρs and the bottom layer, the aquifer, begins with c0 and ρ0.

Upwelling Induced. A second model was created based on the hypothesis that
upwelling of less dense fluid from the Miocene formation below the aquifer is causing free
convection. The model consists of one domain, a 10 m by 30 m sandy aquifer. Just as in the
precipitation induced model, the top boundary is an atmospheric pressure boundary and the side

29

boundaries are no flux and no flow. The bottom boundary of the model was assigned an inward
mass flux based on the concentration and density of the upwelling water (Table 5) at a rate of 4
mm/yr.

Simulation Results
Precipitation Induced. The starting model assumes an infiltration layer 50 mm thick to
simulate a rainfall event similar to the one that coincided with the 2008 data collection. Model
parameters are based on values measured in the field in an attempt to simulate the sabkha aquifer
as accurately as possible. Table 5 lists the values used in the model and their sources.
Table 5. Parameters, values, and sources for COMSOL modeling.
Parameter

Value

Unit

Domain height (average sabkha depth) (L)

10

[m]

Density of aquifer water (ρ0)

1180

[kg/m ]

Density of infiltrating water (ρs)

1305

[kg/m ]

Density of upwelling water (ρu)

1075

[kg/m ]

Concentration of aquifer water (c0)

275.8

[kg/m ]

Concentration of infiltrating water (cs)

440

[kg/m ]

Concentration of upwelling water (cu)

101.5

[kg/m ]

Dynamic viscosity (μ)

0.001

[kg/m s]

Hydraulic conductivity (K)

6.57

[m/d]

Permeability (κ)

6.566x10-

Porosity (θ)

0.38

3
3
3
3
3
3

12

30

2

[m ]

Source
Sanford and
Wood (2001)
Wood et
al.(2002)
Wood et
al.(2002)
Wood et
al.(2005)
Wood et
al.(2002)
Wood et
al.(2002)
Wood et
al.(2005)
Assumed for
fresh water
Measured
Calculated from
K, μ, ρ0,and g
Wood et
al.(2002)

Table 5 (cont’d)
9

2

Simmons et al.
(2001)
SchulzeMakuch (2005)
Calculated
Assumed for
sea-level

Coefficient of molecular diffusion (DL)

1x10-

[m /s]

Dispersivity (α) (lateral and transverse)

0.08

[m]

Linear expansion coefficient (β)

0.8384

Pressure (p)

0

[atm]

Acceleration of gravity (g)

9.81

[m/s ]

Constant

Upwelling flux rate (qu)

4.0

[mm/yr]

Sanford and
Wood (2001)

2

The model simulation ran for 52 weeks. Fingers of high concentration quickly developed
at the top and begin descending towards the bottom of the model (Figure 8). In the first weeks,
the model developed approximately 28 small fingers. Over time, the small fingers began to
coalesce into fewer and larger fingers. At the end of 10 weeks there were 18 fingers and by 25
weeks there were 13 fingers (Figure 9).
The rate of finger descent was relatively fast at first, but slowed over time due to the
dispersion of the solutes (Figure 10). Eventually, as the concentration differences between the
fingers and the aquifer became less and less the rate of descent became very slow. After about
30 weeks, advection slowed down to the point where there is little movement of the fingers. At
this point diffusive forces are becoming dominant and the fingers begin to disappear.
Convective velocity (Uc) was calculated using the equation:

(10)
where:
Uc is the convective velocity (m/d);
K is the hydraulic conductivity (m/d);

31

θ is the porosity;
3

∆ρ is the density difference from initial density (ρ-ρ0) (kg/m );
3

ρ0 is the initial density (kg/m ).

Figure 8. Simulation of groundwater convection from 50mm of precipitation at time t=0, t=2,
t=5, t=10, t=15, t=20, t=25, and t=30 weeks.

32

Figure 9. The number of fingers in the model over time as smaller fingers coalesce into fewer
and larger fingers.

Figure 10. Convective velocity (Uc) and concentration (c) for a single saline finger (x=21.5)
from model (Figure 8) over time.

33

Upwelling Induced. The simulation results show that upwelling from below the aquifer
at a rate of 4 mm/yr has the potential to initiate free convection (Figure 11). A layer of low
concentration fluid continuously builds up at the bottom of the aquifer. The model shows that it
takes around 30 to 31 years of upwelling in order to create a layer that is unstable enough for
convection to occur.

Figure 11. Simulation of groundwater convection caused by upwelling of less dense water at
times t=0, t=15, t=30, t=31, t=32, and t=33 years.

DISCUSSION
The hypothesis that free convection is precipitation induced is supported by both the
geophysics and the modeling. A comparison between the 2008 field data and the modeling
output is shown in Figure 12. As calculated earlier, the convective fingering seen in the
geophysics data is assumed to have been developing for a period of about 15 weeks, therefore, it
is being compared to the model output at 15 weeks. The modeled finger lengths at 15 weeks

34

correlate with the features in the geophysical data. The model is not intended to be an exact
replication of field data and is illustrating 2D representations of 3D features; however, it does
show that fingers can occur, grow, and decay in this system on the timescales that we have
observed with geophysical measurements.
a)

b)

Figure 12. Comparison between 2008 field data for line A (a) and COMSOL model output at
time t=15weeks (b).

Secondary evidence suggests that upwelling may contribute to the density inversion, but
the geophysics shows that this is likely not the main cause since convection in the aquifer in not
a continuous process. Modeling has shown that upwelling is capable of creating density
inversions that, given enough time, can become unstable enough to cause free convection. The
amount of time it takes for upwelling fluids to make the system unstable enough for convection

35

make it an unlikely source due to the fact that it would most likely be erased by the faster and
more frequent precipitation induced convective mixing.
Evapoconcentration was not evident in the geophysics data and as discussed earlier is an
unlikely source because of the presence of the precipitated halite crust. Additionally, Appendix
A shows that the array typed and inversion method used is not sensitive enough to detect a thin
low-resistivity layer between the water table and the capillary zone.

CONCLUSIONS
The analysis of the electrical resistivity data from three consecutive years at a sabkha in
the UAE confirms the hypothesis that free convection is not a continuous event. The episodic
nature of the process, the system’s high Rayleigh number, the timing of rain events, and the
unlikelihood of the other conceptual models strongly suggests that the convection is
precipitation-induced.
Modeling of a density inversion based on a rain event of 50 mm and measured field
values demonstrates that free convection in this system is possible. The timing for the
development of fingers in the model coincides with the timing of a large rain event prior to the
2008 data collection. The decay rate of the fingering and lack of rainfall explains why fingers
where not apparent in subsequent data collections.

36

APPENDICES

37

APPENDIX A

Preprocessing

Testing for and Removing Artifacts
The inverted resistivity data presented in Van Dam et al. (2009) contains two potential
artifacts. The first is a layer of very low resistivity directly below the capillary zone. The second
potential artifact is the low resistivity value found at the ground surface, where the resistivity has
been measured to be much higher (Figure 3). This appendix will describe the testing process to
see if these features were artifacts, and if so, how they could be removed.
The inverted resistivity images of the 2008 data presented in Van Dam et al. (2009) show
a layer of high resistivity just below the surface with a layer of very low resistivity directly
below it near the water table (Figure 13). Initial data analysis of the 2009 data showed a similar
low-resistivity layer, but without the fingering observed in the 2008 data. This layer of low
resistivity has been hypothesized to be the accumulation of dense water though
evapoconcentration and possibly the driving force of the free convection. Alternatively, it has
been suggested that this layer is an artifact of the inversion process caused by the sudden change
from the high resistivity capillary zone into the low resistivity aquifer.

38

Figure 13. Cross section of inverted resistivity for line B using inversion settings from Van Dam
et al. (2009) showing layer of low resistivity below the capillary zone (a) and 1D profile (b)
showing average resistivity (Ohm•m) with depth. See Fig. 6 for area used in averaging.

Testing of this layer as a potential artifact was done by creating simple two and three
layer models, forward calculation of the synthetic resistivity measurement data (using the same
electrode configurations and array combinations), and subjecting them to the same inversion
methods as in Van Dam et al (2009). The two layer model contains only the highly resistive
capillary zone and the low resistivity aquifer. The three layer model has a layer of even lower

39

resistivity between the capillary zone and the aquifer to represent a zone of higher salinity due to
evapoconcentration (Figure 14).

Figure 14. Cross sections of 2-layer (a) and 3-layer (b) models. Both models have a background
resistivity of 0.18 Ohm•m, a capillary zone with a resistivity of 0.71 Ohm•m. The 3-layer model
has an intermediate layer of low resistivity (0.11 Ohm•m) just below the capillary zone.

For both the two and three layer models, the inverted data show a low resistivity layer
followed by a thin layer of high resistivity before returning to the background resistivity value
(Figure 15). The occurrence of these layers in both models suggests they are indeed artifacts
created during the inversion process. These results, however, do not rule out the possibility of
the existence of these layers in the sabkha. Therefore, the next step was to adjust the inversion
process in order to remove the artifacts from the models and then use the updated inversion
methods to reprocess the data.
40

c)

Figure 15. Cross sections of inverted resistivity outputs from the 2-layer model (a) and the 3layer model (b) both showing a zone of low resistivity below the capillary zone. 1D profile (c)
showing average resistivity for middle 40 m of 2- and 3-layer inverted models. See Fig. 6 for
area used in averaging. See Fig. 14 for the starting model.

41

Since the artifacts were likely caused by the way the inversion software handled the sharp
boundary from high to low resistivity, the first step that was taken was to refine the mesh
discretization. In the initial processing of the data (Van Dam et al., 2009), the number of mesh
divisions between electrodes was set at 2. The finer the mesh setting the more accurate the
forward models become, but computing time is longer. Since computing time was not an issue,
the number of mesh divisions was increased to 8, the highest recommended by the software.
To compare the results of the two different meshes, an average resistivity value was
calculated at 20 cm intervals over the middle 40 meters of the inverted resistivity sections.
These were then plotted in 1D profiles of depth versus resistivity (Figure 16c and d). When the
mesh is set at 2 divisions it creates a sinusoidal type curve that shows high amplitude changes in
resistivity with depth, which are reflected in the inverted sections as the layers of low and high
resistivity. With the number of mesh divisions set at 8 the 1D profile still show minor
sinusoidal curves, however their amplitude is muted and the values are much closer to the
expected background values. These results also demonstrate that the method (or array type) used
in not sensitive to the thin low-resistivity layer (0.11 Ohm•m) in Figure 14b.

42

Figure 16. Cross sections of inverted resistivity for 2 layer (a) and 3 layer (b) models for mesh
settings of 2 and 8 divisions between electrodes respectively. 1D profiles of spatially average
resistivity for 2 layer (c) and 3 layer (d) models. See Fig. 6 for area used in averaging.

The 1D profiles indicate another artifact caused by the original inversion settings that are
not consistent with the models or the field site. The inverted resistivity value at the ground
43

surface is calculated to between 0.2-0.3 Ohm•m while the model’s surface layer is set at 0.71
Ohm•m (Figure 16c and d). This error is likely caused by the inversion starting model being set
to “average apparent resistivity” with a value of 0.23 Ohm•m which is normally reserved for
borehole or noisy surface data sets. With this setting, the inversion process begins with a model
that is entirely one apparent resistivity, 0.23 Ohm•m in this case, which it then uses to create the
first forward model. During the iteration process a new forward model is created from the
previous model and it is possible that the artifact is a numerical edge effect at the surface caused
by the starting model.
To investigate if the starting model was the cause of the artifact, the model was run again
using “Pseudosection” as the starting model. This setting uses an interpolated pseudosection of
apparent resistivity based on the input data as the starting model. 1D profiles created from the
new inversion outputs show the resistivities at the surface have a value constant with the models
(Figure 17).

44

Figure 17. 1D profiles of the two different starting models showing each model’s influence on
the resistivity at 0 m depth. See Fig. 6 for area used in averaging.

This analysis of the inversion settings has shown that artifacts were created in the way the
data was initially processed. By changing the mesh size and using a more suitable starting
model, the artifacts have been minimized or removed. For processing the 2008, 2009, and 2010
data in this thesis, the number of mesh divisions will be set at 8 and the starting model will be the
interpolated pseudosection. All other settings will be the same as those used in Van Dam et al.
(2009).

45

Inversion Settings
Minimum Voltage (mv) = 0.02
Minimum V/I (Ohm) = 2E-5
Minimum apparent resistivity (Ohm•m) = 0.01
Maximum apparent resistivity (Ohm•m) = 10000
Maximum repeat error (%) = 5
Maximum reciprocal error (%) = 5
Remove negative apparent resistivity in ERT data: No
Keep All Data (no data removal): No
Inversion Method: Smooth model inversion
Vertical axis: Positive Upward
Y Coordinate = Depth
Min electrode spacing X (m) = 0.01
Min electrode spacing Z (m) = 0.01
Forward Modeling Method: Finite element method
Forward system solver: Cholesky decomposition method
Boundary condition type: Dirichlet
Number cells or elements between two electrodes = 8
Lower-layer-thickness / Upper-layer-thickness = 1.1
Depth of Inverted Model / Depth of Pseudosection = 1.1
Max number of iteration of nonlinear inversion = 8
Stop RMS error = 3%
Minimum error reduction between two iterations = 5%
Stop at Max number of iterations: Yes
Stop when RMS is small enough: Yes
Stop when RMS can not be reduced: No
Res Data reweighting: No
Use Reciprocal Error: No
Stop when L2 norm is small enough: Yes
Initial Lagrange multiplier or roughness factor = 10
Roughness conditioner = 0.2
Starting model: pseudosection.
Minimum resistivity =
0.1 Ohm•m
Maximum resistivity = 100000.0 Ohm•m
Number of elements combined horizontally = 1
Number of elements combined vertically = 1
Vertical / Horizontal roughness ratio = 0.5
Estimated noise of resistivity data = 3%
Initial damping factor of resistivity = 10
Starting iteration of quasi Newton method = 21
IP inversion method: No IP Inversion

46

APPENDIX B

Error Analysis
The raw apparent resistivity data from the ERI survey conducted in 2008 shows evidence
of at least two bad electrodes (Figure 18a). Bad electrodes are easily identified in the raw data
by diagonal trends of anomalous data (see Figure 5). The inversion software is capable of
removing erroneous data points and does so by deleting the bad image point. Even though these
bad image points were removed, the inverted resistivity cross section from the 2008 data line bb’ (Figure 18, b), shows feature that resemble the ones created by the bad electrodes.

Figure 18. Pseudosection of measured apparent resistivity (a) and Inverted resistivity cross
section (b) for the 2008 data line b-b’.

To test if deleting multiple image points from the apparent resistivity pseudosection could
still have an effect on the inverted data, a two layer model (figure 14a) was used and 2 bad
electrodes were simulated using the data from the 2008 b-b’ line. The modeled apparent
resistivity pseudosection (Figure 19a) shows the same erroneous diagonal lines as the 2008 data
(Figure 13a). The data were then inverted using the same method as the 2008 data (Figure 19b).

47

Figure 19. Pseudosection of calculated apparent resistivity (forward calculation) with 2 bad
electrodes (a) and modeled inverted resistivity cross section (b).

Analysis of these results shows that there is indeed no effect caused by deleting multiple
image points from bad electrodes. This confirms that the features present in the 2008 resistivity
cross section of line b-b’ are the result of free convection and not an artifact creating during
inversion.

48

APPENDIX C

2010 Sabkha Stress Test Experiment
In the 2010 visit, an 11.5 by 12 m grid of 240 electrodes with 0.5 m electrode spacing and two 8electrode borehole arrays with 0.25 m electrode spacing were installed at the field site to monitor
an infiltration experiment using dipole-dipole configurations. Over the course of 15 days, ~6350
liters of saline water were extracted from the sabkha 110 m southwest of the study area. Over a
6-day period the sabkha water was then sprayed evenly over a 10 m diameter circle, centered in
the study area (Figure 20), resulting in approximately 80 mm of infiltration. The electrical
resistivity monitoring system was programmed to perform daily measurements of the study area
and then upload them for processing using a wireless USB modem. In addition to ERT
monitoring, six groundwater sampling wells were installed, four inside the infiltration area and
two outside, at depths between 1.6 and 2.5 m. Groundwater samples were collected from the
wells daily throughout the experiment and then periodically over the next three months.

49

Figure 20. Map of the study area showing the 2010 electrode arrays, infiltration zone, and
pressure transducer (background image from Google Maps taken 11/7/2010.)

50

REFERENCES

51

REFERENCES

Bauer, P., Supper, R., Zimmermann, S., and Kinzelbach, W. (2006). Geoelectrical imaging of
groundwater salinization in the Okavango Delta, Botswana, 60, 126–141.
Boggs, S. (2001). Principles of sedimentology and stratigraphy third edition. Upper Saddle
River, NJ: Prentice Hall.
COMSOL Multiphysics (Version 4.2) [Software]. (2011). Burlington, MA: COMSOL, Inc.
Diersch, H.-J. G., and Kolditz, O. (2002). Variable-density flow and transport in porous media:
approaches and challenges. Advances in Water Resources, 25(8-12), 899-944.
Elder J. W. (1967). Steady free convection in a porous medium heated from below. Journal of
Fluid Mechanics, 27, 29–48.
Elder J. W. (1967). Transient convection in a porous medium. Journal of Fluid Mechanics, 27,
609–23.
Gieske, A. (1996). Vegetation driven groundwater recharge below the Okavango Delta
(Botswana) as a solute sink mechanism — an indicative model. Botswana Journal of Earth
Sciences 3, 33–37.
Hillel D. (1982). Introduction to Soil Physics. Academic Press, Inc., Orlando, FL (1982), pp.
212-215.
Horton, C.W., and Rogers Jr., F.T. (1945). Convection currents in a porous medium. Journal of
Applied Physics. 16, 367–370.
Lapwood, E. R. (1948). Convection of a fluid in a porous medium. Proceedings of the
Cambridge Philological Society 44, 508-521.
Nield, D. A., and Bejan, A. (2006). Convection in Porous Media, 3rd ed., Springer, New York.
Nield, D. A, Simmons, C. T., Kuznetsov, A. V., and Ward, J. D. (2008). On the evolution of salt
lakes: Episodic convection beneath an evaporating salt lake. Water Resources Research,
44(2), 1-13.
Parkhurst, D. L., and Appelo, C. A. J. (1999). User’s guide to PHREEQC (version 2)--A
computer program for speciation, batch-reaction, one-dimensional transport, and inverse
geochemical calculations: U.S. Geological Survey Water-Resources Investigations Report
99-4259, 312

52

Rayleigh, L. (1916). On convection currents in a horizontal layer of fluid when the higher
temperature is on the underside. Philosophical Magazine 32 (6), 527– 546.
Sanford, W. E., and Wood, W. W. (2001). Hydrology of the coastal sabkhas of Abu Dhabi,
United Arab Emirates. Hydrogeology Journal, 9(4), 358-366.
Schneider, K.J. (1963). Investigation of the influence of free thermal convection on heat transfer
through granular material: In: Proceedings of the 11th International Congress of
Refrigeration, Oxford: Pergamon Press, 11(4), 247–253.
Schulze-Makuch, D. (2005). Longitudinal dispersivity data and implications for scaling
behavior. Ground water, 43(3), 443-456.
Simmons, C. T., Fenstemaker, T. R., and Sharp Jr., J. M. (2001). Variable-density groundwater
flow and solute transport in heterogeneous porous media: approaches, resolutions and future
challenges. Journal of Contaminant Hydrology, 52(1-4), 245-75.
Simmons, C. T. (2005). Variable density groundwater flow: From current challenges to future
possibilities. Hydrogeology Journal, 13(1), 116-119.
Stevens, J. D., Sharp Jr., J. M., Simmons, C. T., and Fenstemaker, T. R. (2009). Evidence of free
convection in groundwater: Field-based measurements beneath wind-tidal flats. Journal of
Hydrology, 375(3-4), 394-409.
Van Dam, R. L., Simmons, C. T., Hyndman, D. W., and Wood, W. W. (2009). Natural free
convection in porous media: First field documentation in groundwater. Geophysical
Research Letters, 36(11), L11403.
Wood, W. W., Sanford, W. E., and Al Habshi A. R. S. (2002). Source of solutes to the coastal
sabkha of Abu Dhabi Source of solutes to the coastal sabkha of Abu Dhabi. Geological
Society Of America Bulletin, 114(4), 259-268.
Wood, W. W., Sanford, W. E., and Frape, S. K. (2005). Chemical openness and potential for
misinterpretation of the solute environment of coastal sabkhat. Chemical Geology, 215(1-4),
361-372.
Wooding R. A. (1957). Steady state free thermal convection of liquid in a saturated permeable
medium. Journal of Fluid Mechanics, 2, 273–285.
Zimmermann, S., Bauer, P., Held, R., Kinzelbach, W., and Walther, J. (2006). Salt transport on
islands in the Okavango Delta: Numerical investigations. Advances in Water Resources,
29(1), 11-29.

53