U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2009–5028

Application of SEAWAT to Select Variable-Density and 

Viscosity Problems

Application of SEAWAT to Select VariableDensity and Viscosity Problems

By Alyssa M. Dausman, Christian D. Langevin, Danny T. Thorne, Jr., 

and Michael C. Sukop

Scientific Investigations Report 2009–5028

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

KEN SALAZAR, Secretary

U.S. Geological Survey

Marcia K. McNutt, Director

U.S. Geological Survey, Reston, Virginia: 2010

For more information on the USGS—the Federal source for science about the Earth, its natural and living resources, 

natural hazards, and the environment, visit http://www.usgs.gov or call 1-888-ASK-USGS

For an overview of USGS information products, including maps, imagery, and publications, 

visit http://www.usgs.gov/pubprod

To order this and other USGS information products, visit http://store.usgs.gov

Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this report is in the public domain, permission must be secured from the individual copyright owners to 

reproduce any copyrighted materials contained within this report.

Suggested citation:

Dausman, A.M., Langevin, C.D., Thorne Jr., D.T., and Sukop, M.C., 2010, Application of SEAWAT to Select VariableDensity and Viscosity Problems: U.S. Geological Survey, Scientific Investigations Report 2009-5028, 31 p.
iii

Contents

Abstract...........................................................................................................................................................1

Introduction.....................................................................................................................................................1

Purpose and Scope ..............................................................................................................................3

Acknowledgments................................................................................................................................3

Benchmark Problems....................................................................................................................................3

Problem 1: One-Dimensional Flow with Linearly Varying Viscosity Field....................................3

Problem 2: Two-Dimensional Oil Convection in Aluminum Foam .................................................6

Problem 3: Horton-Rogers-Lapwood (HRL) Convection...............................................................11

Problem 4: Double-Diffusive Finger Convection............................................................................14

Problem 5: Elder Problem ..................................................................................................................18

Problem 6: Henry-Hilleke Problem...................................................................................................21

Summary and Conclusions.........................................................................................................................25

References Cited..........................................................................................................................................26

Appendix........................................................................................................................................................29

Figures

1. Diagram showing map view of model design in SEAWAT_V4 ..............................................5

2. Plot showing results of water level values with distance in each of the 

SEAWAT_V4 models and the analytical solutions...................................................................5

 3. Schematic diagram of problem simulating temperature-dependent viscosity..................6

 4. Graph showing change in viscosity with temperature ..........................................................8

 5. Schematics showing streamline results for both variable viscosity and 

constant viscosity cases.. .........................................................................................................10

 6. Schematics showing temperature isotherms for both variable viscosity and 

constant viscosity cases ............................................................................................................. 10

 7. Schematic showing design of the Horton-Rogers-Lapwood benchmark problem 

where the length of the box is 2 times (or 4 times) the height ............................................11

 8. Schematic showing design of the Horton-Rogers-Lapwood problem in 

SEAWAT_V4.................................................................................................................................11

 9. Schematics showing temperature results for the Horton-Rogers-Lapwood 

convection problem simulated with SEAWAT_V4.................................................................13

10. Schematic showing velocity vectors for the Horton-Rogers-Lapwood problem 

where Ra* = 40............................................................................................................................13

11. Diagram showing Hele-Shaw cell set up ...............................................................................15

12. Sectional diagram showing initial conditions in the Hele-Shaw cell ................................15

13. Laboratory photographs showing results from Pringle and others (2002)........................16

14. Laboratory photographs showing results from SEAWAT_V4..............................................16

15. Schematic representation of Elder’s original experiment as simulated in 

SEAWAT_V4.................................................................................................................................18

16. Schematics showing results of SEAWAT_V4 compared to the original laboratory 

results from Elder (1967), SUTRA, and an older version of SEAWAT ................................20
iv

17. Graph showing the Nusselt number versus time showing the heat flux across 

the bottom boundary in SEAWAT_V4 ......................................................................................21

18. Diagram showing setup of the Henry-Hilleke problem showing the boundary 

conditions which include freshwater input, seawater boundary, and 

temperature boundaries............................................................................................................22

19. Plots of modeling results showing percent seawater concentrations from 

SUTRA-MS, Henry and Hilleke numerical solution, HST3D results, 

and SEAWAT_V4 ........................................................................................................................24

20. Plots of modeling results showing temperatures from SUTRA-MS, Henry and 

Hilleke numerical solution, HST3D results, and SEAWAT_V4.............................................24

21. Plots of modeling results showing magnitude and direction for velocity vectors 

for SUTRA-MS and SEAWAT_V4..............................................................................................25

Tables

1. Identification of variables............................................................................................................4

2. Input parameters for problem 2 simulating temperature-dependant viscosity 

of oil in aluminum foam................................................................................................................7

3. Input parameters for the Horton-Rogers-Lapwood problem simulated 

in SEAWAT_V4.............................................................................................................................12

4. SEAWAT_V4 input parameters for the Hele-Shaw simulation............................................17

5. Input parameters for the Elder problem simulated in SEAWAT_V4 ...................................19

6. Input parameters for the Henry-Hilleke problem simulated in SEAWAT_V4....................23
v

Conversion Factors

Inch/Pound to SI

Multiply By To obtain

foot (ft) 0.3048 meter 

foot per day (ft/d) 0.3048 meter per day 

cubic foot (ft3) 0.02832 cubic meter

pound per foot per day (lb/ft/d) 1.488 kilogram per meter per day

SI to Inch/Pound

Multiply By To obtain

meter (m) 3.281 foot

meter per second (m/s) 3.281 foot per second

meter per day (m/d) 3.281 foot per day

meter per second squared (m/s2) 3.281 foot per second squared

meter per day squared (m/d2) 3.281 foot per day squared

square meter (m2) 10.76 square foot

square meter per day (m2/d) 10.76 square foot per day

cubic meter per day (m3/d) 35.31 cubic foot per day

kilogram per second (kg/s) 2.205 pound per second

kilogram per day (kg/d) 2.205 pound per day

kilogram per meter per second (kg/m s) 0.672 pound per foot per second

kilogram per meter per day (kg/m d) 0.672 pound per foot per day

kilogram per cubic meter (kg/m3) 0.06242 pound per cubic foot

joule (J) 0.0000002 kilowatt hour

joule per kilogram per degrees Celsius (J/kg °C) 0.4535 joule per pound per degrees Celsius

watt per meter per degrees Celsius (W/m °C) 0.3048 watt per foot per degrees Celsius

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as follows:

°F=(1.8×°C)+32

Acronyms

MODFLOW MODular three-dimensional finite-difference ground-water FLOW model

MT3DMS Modular 3-D Multi-Species Transport Model

COND Hydraulic conductance

SUTRA Saturated-Unsaturated Fluid-Density-Dependent Ground-Water Flow with 

Transport of Thermal Energy or a Single Solute Species

SUTRA-MS SUTRA-simulate heat and multiple-solute transport

PCG2 Preconditioned Conjugate-Gradient 2

PAO PolyAlphaOlefin

NaCl Sodium Chloride

HST3D Heat and Solute Transport in Three-Dimensional Ground-Water Flow Systems

TVD Total Variation Diminishing

GCG Generalized Conjugate Gradient

SSOR Slice-Successive Over Relaxation
vi
Abstract

SEAWAT is a combined version of MODFLOW and MT3DMS, designed to simulate three-dimensional, variable-density, 

saturated groundwater flow. The most recent version of the SEAWAT program, SEAWAT Version 4 (or SEAWAT_V4), supports 

equations of state for fluid density and viscosity. In SEAWAT_V4, fluid density can be calculated as a function of one or more 

MT3DMS species, and optionally, fluid pressure. Fluid viscosity is calculated as a function of one or more MT3DMS species, 

and the program also includes additional functions for representing the dependence of fluid viscosity on temperature. 

This report documents testing of and experimentation with SEAWAT_V4 with six previously published problems that include 

various combinations of density-dependent flow due to temperature variations and/or concentration variations of one or more 

species. Some of the problems also include variations in viscosity that result from temperature differences in water and oil. 

Comparisons between the results of SEAWAT_V4 and other published results are generally consistent with one another, with 

minor differences considered acceptable.

Introduction

Groundwater systems are often envisioned as constant-density systems, where groundwater flow is mathematically 

straightforward. However, groundwater systems that have a shallow water table with high evapotranspiration rates, inputs of 

a highly concentrated contaminant, or saltwater intruding into the fresh groundwater at the coast will contain increased solute 

concentrations. The increased concentrations often result in spatial variations of fluid density that affect groundwater flow. Fluid 

density of seawater (1,025 kg/m3) is about 2.5 percent higher than that of freshwater (1,000 kg/m3); this difference in density 

will typically result in differing flow fields (Guo and Langevin, 2002). Therefore, the mathematical representation of variabledensity groundwater systems is more complicated than constant-density systems because the solute concentration affects the 

density and the density affects the flow.

Mathematically representing groundwater systems affected by variable-density flow can be important for research related 

to coastal systems such as saltwater intrusion or submarine groundwater discharge. Mathematical models of coastal aquifers 

were created as early as 1888 to predict the location of the saltwater-freshwater interface (Ghyben, 1888; Herzberg, 1901). The 

Ghyben-Herzberg model uses simple hydrostatics to calculate the thickness of a static freshwater lens over a static saltwater 

wedge. This model assumes a sharp interface. Hubbert (1940) modified the Ghyben-Herzberg model to take into account nonstatic freshwater where it flows toward the ocean along a sharp interface. Glover (1959) extended Hubbert’s model to include 

a zone between the shoreline and the saltwater-freshwater interface with known discharge at the coast. Cooper (1959) created 

a steady-state mathematical model that accounts for the transition zone from freshwater to seawater. Cooper’s model represents the freshwater discharging at the coast over saltwater that intrudes deep within an aquifer, and rises and returns with the 

freshwater at the coast (simulating a brackish-water zone and a convection cell in the saltwater). Later, Henry (1964) provided 

a semianalytical solution for a saltwater-freshwater interface in a dynamic coastal system where freshwater is flowing toward 

an ocean boundary. Numerous variations of the Henry problem have been solved since 1964, including a temperature-salinity 

version developed both numerically and in a laboratory setting by Henry and Hilleke (1972).

Application of SEAWAT to Select Variable-Density and 

Viscosity Problems

By Alyssa M. Dausman1, Christian D. Langevin1, Danny T. Thorne, Jr.2, and Michael C. Sukop3

1

U.S. Geological Survey, Florida :DWHU Science Center, Ft. Lauderdale

2

Georgetown College, Georgetown, Kentucky

3

Florida International University, Miami, Florida
2 Application of SEAWAT to Select Variable-Density and Viscosity Problems

A tacit assumption in most groundwater models is that isothermal conditions prevail and that temperature does not affect 

the density of the water. Therefore, the effects of viscosity variations due to changes in temperature (and often salinity), 

are also assumed to be so small they do not influence flow patterns. These assumptions are typically used, even in variabledensity models designed to represent saltwater intrusion, wastewater injection, and convection beneath salt lakes. Research 

has shown, however, that temperature and salinity can affect the flow field (Henry and Hilleke, 1972), particularly in deeper 

aquifer systems. In some deep aquifers, geothermal heating from below can cause the formation of an unstable density profile 

(Elder, 1967). 

The density of water is inversely related to temperature, often causing convective flow when warm, less dense water is 

overlain by cool, more dense water. Like fluid density, viscosity is also a function of temperature. Thus, groundwater flow 

patterns may change in response to a changing temperature distribution because of the effects on density and viscosity. The 

density of freshwater, 1,000 kg/m3, can be decreased to approximately 999.6 kg/m3 by a 10  ̊C increase in temperature. This is 

only about a 0.5-percent reduction in density and can often be ignored provided that the Rayleigh Number is adequately small; 

however, when temperature changes occur over a wider range, the decrease in density can have a significant effect on groundwater flow (Henry and Hilleke, 1972). Therefore, the modeling of some groundwater systems could require simulating the 

effects of temperature on variable-density flow.

A number of mathematical approaches have been used to simulate variable-density groundwater flow, such as finite-elements 

in SUTRA (Voss, 1984), analytical elements by Strack (1995), and finite difference approaches such as earlier versions of 

SEAWAT (Guo and Langevin, 2002; Langevin and others, 2003), MOCDENSE (Sanford and Konikow, 1985), and HST3D 

(Kipp, 1987). SUTRA was modified to SUTRA-MS by Hughes and Sanford (2004) to simulate density-dependent flow resulting from multiple species such as concentration and temperature. Also, HST3D was modified by Kipp (1997) to include both 

heat and solute transport effects on variable-density flow.

SEAWAT (Langevin and others, 2003), a combined version of MODFLOW (Harbaugh and others, 2000) and MT3DMS 

(Zheng, 1990; Zheng and Wang, 1999; Zheng, 2006), is designed to simulate three-dimensional, variable-density groundwater 

flow. SEAWAT was originally designed to solve variable-density groundwater flow and transport problems with the assumption that fluid density was a function of only a single solute constituent (Guo and Bennett, 1998; Guo and Langevin, 2002; 

Langevin and others, 2003). Therefore, these previous SEAWAT versions could not be used to simulate simultaneous solute 

and heat transport. The latest version of SEAWAT, SEAWAT Version 4 (SEAWAT_V4), was designed in a general fashion so 

that fluid density and viscosity can be calculated as a function of one or more species, and heat can be represented as one of 

the species (Langevin and others, 2007). The code can be used to simulate different types of fluids, including oil (as a single 

phase) or other fluids with variable viscosity, because the effects of viscosity on flow are included in the simulations. A distinct 

diffusion coefficient for each species can also be included in SEAWAT_V4. 

For a numerical model to be reliable, the code must be rigorously tested to ensure that it accurately represents physical 

processes. Voss and Souza (1987) suggested that new modeling codes be tested with benchmark problems. A previous version 

of SEAWAT (Langevin and others, 2003) has been tested with most of the benchmark problems that are traditionally used to 

test saltwater intrusion programs. 

This study was conducted to document benchmark testing of the new features in SEAWAT_V4 (Langevin and others, 

2007), including the code’s ability to simulate a fluid other than water (for example, oil as a single phase). In the testing of this 

new code, the effects of (1) viscosity on convection, (2) heat/energy transport into and out of modeling systems, and 

(3) multiple diffusion coefficients were explored. The study also compared the results of Darcy flow simulations for which 

friction is ignored along a model boundary (referred to as a slip boundary) with other simulations that account for the resistive 

drag along a model boundary. 

Testing of SEAWAT_V4 was performed using six benchmark problems:

• One-dimensional flow through a linearly varying viscosity field (Thorne and others, 2006);

• Two-dimensional oil convection in aluminum foam (Guo and Zhao, 2005);

• Horton-Rogers-Lapwood (HRL) convection (Horton and Rogers, 1945; Lapwood, 1948);

• Double-diffusive finger convection (Pringle and others, 2002);

• Original Elder problem (Elder, 1967);

• Henry-Hilleke problem (Henry and Hilleke, 1972).

These benchmark problems are well-defined, and have been represented with laboratory experiments, analytical solutions, 

or with other modeling programs. Fluid dynamics and dimensionless numbers have been used in coordination with the benchmark problems to enable comparison between SEAWAT_V4 and other modeling codes, laboratory experiments, and systems 

of different dimensions.
Benchmark Problems 3

Purpose and Scope

This report documents the use of six benchmark problems to test the SEAWAT_V4 computer program (Langevin and 

others, 2007). Important concepts underlying the physical processes are briefly described. The effect of density and viscosity 

variations (resulting from multiple species) on groundwater flow is tested as is the ability of the code to simulate heat conduction and convection. For a complete description of the SEAWAT_V4 and related computer programs, readers are directed to 

the following sources of information: Guo and Bennett (1998), Guo and Langevin (2002), Langevin and others (2003), and 

Langevin and others (2007).

Acknowledgments

We would like to thank Weixing Guo and Alden Provost for their technical reviews of this report. Their comments and 

suggestions improved its content immensely.

Benchmark Problems

A series of benchmark problems were simulated with SEAWAT_V4. The benchmark problems use different equations and 

parameters. The governing equations for the simulations can be found in Guo and Langevin (2002), Langevin and others (2003), 

and Langevin and others (2007). The dimensionless numbers used to compare the systems are described in the appendix. All the 

variables represented in the equations are defined in table 1. 

Problem 1: One-Dimensional Flow with Linearly Varying Viscosity Field

This benchmark problem consists of one-dimensional groundwater flow between two constant-head boundaries. The purpose 

of the problem is to test the implementation of viscosity effects in SEAWAT_V4. The problem is tested for the premise that 

simulations with a linearly varying viscosity field and a constant value for reference hydraulic conductivity will give the same 

result as a simulation with a linearly varying reference hydraulic conductivity field and a constant fluid viscosity. Thorne and 

others (2006) provide an extensive discussion on the equations used for this benchmark problem. 

In SEAWAT (Guo and Langevin, 2002) and SEAWAT_V4 (Langevin and others, 2007), Darcy’s law is implemented as: 

Q C= ⋅ OND h∆ , where COND is the conductance between two adjacent nodes and ∆h is the head difference between these 

two nodes. In SEAWAT_V4, the conductance COND is defined as:

COND KWH

L

f = 

 (1)

using “equivalent freshwater” hydraulic conductivity (Guo and Langevin, 2002), such that K values used in equation 1 represent 

an aquifer that is saturated with the reference fluid at the reference temperature (normally assumed to be freshwater at 25 degrees 

Celsius). The ratio 



f is the variation in dynamic viscosity μ = μ(C,T) from a reference viscosity μf

. In SEAWAT_V4, the horizontal conductance between two cells can be averaged using a harmonic or logarithmic mean. The vertical conductance across 

two layers with adjacent cells is calculated using the harmonic average.

Conductance between cells in a model can vary by changing either the hydraulic conductivity or the viscosity ratio (



f

). 

Equation 1 reveals that the effects of changing the hydraulic conductivity and the viscosity ratio are inversely related. Therefore, 

if the viscosity ratio in a cell is increased and hydraulic conductivity in that same cell is proportionally decreased, the conductance 

between that cell and any adjacent cell will be the same as if no increase/decrease were made in the parameters.

An analytical equation can easily be derived for steady-state, confined, one-dimensional flow through a linearly varying 

hydraulic conductivity field. Two prescribed heads can arbitrarily be defined as: at x = xmin, h = h0; at x = xmax, h = h1. A linearly 

varying hydraulic conductivity field is represented by arbitrarily setting K(x) = x (for all x > 0) throughout the domain. Thus, 

Darcy’s law is q=-x(dh/dx), where q is flux and dh/dx is the hydraulic gradient. The analytical solution is obtained by setting the 

derivative of q with respect to x to zero (dq/dx = 0; that is, flux is uniform along the domain) and integrating between xmin and 

xmax subject to the boundary conditions:

h h h

x x = x x h −− [ − ] + 1 0

0 ln( ) ln( )

ln( ) ln( )

max min

min

(2)

μ

μ

μ

μ

μ

μ
4 Application of SEAWAT to Select Variable-Density and Viscosity Problems

Table 1. Identification of variables. 

L=Length, T=Time, M=Mass, D=Degrees of Temperature

Symbols Dimensions Definition

q LT -1 Darcy flux

Qin L3

T -1 Volumetric flux

H L Height of a defined system or cell

L L Length of a defined system or cell

W L Width of a defined system or cell

h L Head

x L Distance in x direction

K LT -1 Hydraulic conductivity

k L2 Permeability

Q C= ⋅ OND h∆

L L Longitudinal dispersivity

COND

T L Transverse dispersivity

Sy - Specific yield

S - Storage

θ - Porosity

D* L2

T -1 Bulk thermal diffusivity without porosity in 

the denominator

D** L2

T -1 Bulk thermal diffusivity with porosity

∆h ML-3 Density of water

COND KWH

L

f = 



s ML-3 Density of solids





f

b ML-3 Bulk density





f D-1 Thermal expansion coefficient h h h

x x = x x h −− [ − ] + 1 0

0 ln( ) ln( ) ln( ) ln( )

max min

min

M-1L3 Solutal expansion coefficient

v L2

T -1 Kinematic viscosity

v0 L2

T -1 Reference kinematic viscosity

μ ML-1T -1 Dynamic viscosity

μ0 ML-1T -1 Reference dynamic viscosity

T D Temperature

Tref D Reference temperature for reference 

dynamic/kinematic viscosity

cPfluid L2

T--2D-1 Heat capacity of fluid

cPSolid L2

T--2D-1 Heat capacity of solid

kTBulk LMT-3D-1 Bulk thermal conductivity

kTfluid LMT-3D-1 Thermal conductivity of fluid

kTsolid LMT-3D-1 Thermal conductivity of solid

g LT -2 Acceleration due to gravity

Ra - Rayleigh number

Da - Darcy number

Pr - Prandtl number

Nu - Nusselt number

Qh DT-1 Heat Flux
Benchmark Problems 5

The Viscosity (VSC) Package in SEAWAT_V4 can be used to include the effects of viscosity variations on groundwater 

flow. Viscosity can be specified on a cell-by-cell basis, or it can be calculated as a function of one or more MT3DMS species 

if the simulation includes transport. It should be noted that temperature may be one of the species simulated in the model. For 

the test problem described herein, the value of viscosity is specified as 1/x in the cases where hydraulic conductivity is kept constant; 

this should give the same results as a problem where the hydraulic conductivity is equal to x and the viscosity is held constant.

Four simulations were performed with SEAWAT_V4. The first three simulations use a constant reference hydraulic conductivity 

and a prescribed viscosity gradient in (a) the x-direction (horizontal) with 91 columns, 1 row, and 1 layer, (b) the y-direction 

(horizontal) with 1 column, 91 rows, and 1 layer, and (c) the z-direction (vertical) with 1 column, 1 row, and 91 layers. The onedimensional simulations were performed in all three layer, row, and column directions to ensure that the implementation was 

coded correctly in each direction. A fourth simulation, with a constant viscosity and a prescribed hydraulic conductivity gradient 

in the x-direction with 91 columns, 1 row, and 1 layer was also performed. Head results from these simulations should match the 

analytical solution exactly, with the possible exception of the third case (c). The analytical solution is based on a linearly varying 

hydraulic conductivity field. To exactly represent a linearly varying hydraulic conductivity field in SEAWAT_V4 that directly 

matches the analytical solution, the logarithmic mean interblock transmissivity method must be used to calculate internodal 

conductances. While logarithmic averaging was used in the first two simulations, harmonic averaging is the only option for calculating internodal conductances in the vertical direction (used in the third case). Results from the third simulation (c), however, 

should be very close to the analytical solution.

Figure 1 shows the grid with prescribed water levels on each end, 0 ft and 1 ft (0.3048 m), respectively, and the 91 columns 

(or rows, or layers), with each cell being 1 ft3 (0.02832 m3). K(x)=1 ft/d (0.3048 m/d) in the three cases with a prescribed viscosity 

gradient where the dynamic viscosity (μ) in each cell is 1/x, and the reference viscosity is 11.94 lb/ft/d (86.4 kg/m/d). In the fourth 

case, viscosity is constant at 11.94 lb/ft/d (86.4 kg/m/d) and the hydraulic conductivity in each cell varies with distance, where 

K(x) = x. The Preconditioned Conjugate-Gradient 2 (PCG2) solver from MODFLOW2000 was used to solve for head in each 

cell. The head change criterion was set to 1 x 10-8 ft (3.048 x 10-9 m). 

Results indicate that the four SEAWAT_V4 simulations are in good agreement with the analytical solution (fig. 2). 

Although conductance is calculated using the harmonic mean in the vertical direction (as opposed to the logarithmic mean in the 

horizontal direction), differences between the models are minimal in the variable-viscosity cases. These results indicate that the 

effects of viscosity on groundwater flow have been implemented correctly in SEAWAT_V4. 

Figure 1. Map view of model design in SEAWAT_V4. h(x) is water level in feet, x is distance from bottom left 

corner of model grid. 

X=10.5 feet X=100.5 feet

h(xMIN)=0 feet h(xMAX)=1 foot

Figure 2. Results of water level values with distance in each of the SEAWAT_V4 models and 

the analytical solutions.

Analytical Solution

Variable permeability

Variable viscosity-multiple rows

Variable viscosity-multiple columns

Variable viscosity-multiple layers

EXPLANATION

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0 10 20 30 40 50 60 70 80 90 100

X, IN FEET (OR Y, OR Z, DEPENDING ON CASE)

HEAD, IN FEET
6 Application of SEAWAT to Select Variable-Density and Viscosity Problems

Problem 2: Two-Dimensional Oil Convection in Aluminum Foam

The second problem is a Lattice Boltzmann method simulation of oil convection in aluminum foam and is adopted as a 

potential new benchmark. For this problem, convection is driven by temperature differences between two vertical walls on 

opposing sides of a square cross section. Each simulation is run with and without the effects of viscosity variations. Because the 

simulated fluid is oil, there is a relatively strong dependence of viscosity on temperature (as compared to water). The purpose 

of this benchmark problem is to test the variable viscosity mechanism implemented in SEAWAT_V4 and to test the capability 

of SEAWAT_V4 to simulate thermally driven convection in oil. Nield and Bejan (1999) provide an extensive review of thermal 

convection when modeling fluids such as water.

It is important to note that this current benchmark problem using oil is new, and while it is described by Guo and Zhao 

(2005), it has not been verified against other modeling codes or laboratory tests. Thus, the results reported by Guo and Zhao 

(2005) cannot be considered definitive, and the status of this problem as a benchmark can only be considered preliminary. While 

there have been some studies on temperature-dependent viscosity of water and oil in porous media (Horne, 1975; Weber, 1975; 

Gary and others, 1982; Pasa and Titaud, 2005; Afify, 2007), there is not a specific benchmark problem simulating natural convection and variations in viscosity of single phase oil from changes in temperature in a porous medium, except for the work by 

Guo and Zhao (2005). More tests will need to be done to establish Guo and Zhao’s work as a standard benchmark problem. 

Guo and Zhao (2005) used Lattice Boltzmann methods to simulate convection of PolyAlphaOlefin (PAO) oil in aluminum 

foam. The researchers conducted a series of simulations with different Darcy and Rayleigh numbers (eqs. A7 and A9 in appendix), and also tested the effects of temperature-dependent viscosity. Schematic representation of the problem is shown in 

figure 3 and the associated model parameters, grid dimensions, and solution schemes are given in table 2.

Figure 3. Problem simulating temperature-dependent viscosity (modified 

from Guo and Zhao, 2005).

Gravity

Porous medium

T = 36 C

o

T = 6 C

o

No flow

No flow

H = 1 m

No flo

w

No flo

w
Benchmark Problems 7

Table 2. Input parameters for problem 2 simulating temperature-dependant viscosity of oil in aluminum foam.

Input Parameter Value Units Comments

Number of columns 130 - Assigned

Number of rows 1 - Assigned

Number of layers 128 - Assigned

∆x(DELR) 0.0078125 m Assigned

∆y(DELC) 0.0078125 m Assigned

∆z(DELZ) 0.0078125 m Assigned

H 1 m Length and height of “active” model domain

K, Da = 10-4 127 m/s Hydraulic conductivity

K, Da = 10-2 12,658 m/s Hydraulic conductivity

k, Da = 10-4 1.0 × 10-4 m2 Permeability

k, Da = 10-2 1.0 × 10-2 m2 Permeability Q C= ⋅ OND h∆

L 0 m Longitudinal dispersivity COND

T 0 m Transverse dispersivity

Sy 0.01 - Specific yield

S 0.01 - Storage

∆h 0.58 - Porosity

D* 4.7184 × 10-5 m2

/s Thermal diffusivity

COND KWH

L

f = 



water 768.5 kg/m3 Density of water at 21 o

C (reference temperature) 





f , Ra = 1 x 105 -9.55 × 10-5 kg/(m3 o

C) Density change with temperature





f

, Ra = 1 x 106 -9.55 × 10-4 kg/(m3 o

C) Density change with temperature h h h

x x = x x h −− [ − ] + 1 0

0 ln( ) ln( ) ln( ) ln( )

max min

min

, Ra = 1 x 107 -9.55 × 10-3 kg/(m3 o

C) Density change with temperature ( ) . . T T = ⋅ − 0 168 1 0868

, Ra = 1 x 105 1.24 × 10-7 1/ o

C Thermal expansion coefficient  T T ( ) = ⋅ +





 





0.00002394 10

248 37

133 15

.

.

, Ra = 1 x 106 1.24 × 10-6 1/ o

C Thermal expansion coefficient

∂

∂



T

, Ra = 1 x 107 1.24 × 10-5 1/ o

C Thermal expansion coefficient

v0 7.74 × 10-6 m2

/s Reference kinematic viscosity

μ0 0.00595 kg/(m s) Reference dynamic viscosity

Tref 21 o

C Reference temperature for reference dynamic 

viscosity

cPfluid 1,971.35 m2

/(s2 oC) Heat capacity of fluid

kTfluid 0.1424 m kg/(s3 o

C) Thermal conductivity of fluid

kTsolid 170 m kg/(s3 o

C) Thermal conductivity of solid

g 9.8 m/s2 Acceleration due to gravity

Cool temperature boundary 6 oC Assigned

Hot temperature boundary 36 oC Assigned

Matrix solution technique for flow PCG2 - Assigned

Head convergence value 0.001 m Assigned

Flow convergence value 1 kg/s Assigned

Advection term TVD - Assigned

Dispersion and source terms Implicit finite difference; 

GCG, SSOR - Assigned

Time-step length 1-10,000 s Varies between models for different Rayleigh 

numbers

Temperature convergence value 1 x 10-12 to 1 x 10-13 - Varies between models for different Rayleigh 

numbers
8 Application of SEAWAT to Select Variable-Density and Viscosity Problems

The viscosity (μ) of a saturating fluid is often assumed to be constant if the viscosity of the fluid is only weakly dependent 

upon temperature. However, when simulating environments with large temperature differences or that contain fluids with strong 

temperature-dependent viscosities (such as oil), variations in viscosity cannot be ignored. The change in viscosity will affect fluid 

flow and consequently, heat transfer. Viscosity varies with temperature in PAO oil according to the equation (Guo and Zhao, 2005);

( ) . . T T = ⋅ − 0 168 1 0868

(3)

it is an optional equation for temperature-dependent viscosity implemented in SEAWAT_V4. Another equation from SEAWAT_V4 

that can be used for viscosity variations with temperature when simulating water is (Hughes and Sanford, 2004):

 T T ( ) = ⋅ +





 





0.00002394 10

248 37

133 15

.

.

(4)

Comparing the change in viscosity with temperature according to equations 3 and 4 shows that viscosity variations with temperature 

are much larger when simulating oil than when simulating water (fig. 4).

The simulation uses a finite difference grid of 130 columns by 128 layers by 1 row to discretize the cross section that is 1 m 

high by 1.015625 m long (fig. 3). Insulating (that is, zero-temperature gradient) no-flow boundaries are assigned across the top 

and bottom of the box. There is a no-flow, constant “cool” boundary (Tc) of 6 oC along the right-hand side and a no-flow, constant “hot” temperature boundary (Th

) of 36 oC along the left-hand side. The columns containing the constant cool/hot temperature boundaries have a lowered hydraulic conductivity to minimize the convective heat flux from the temperature boundaries; by 

definition, conduction should be the sole heat transport mechanism at the boundary. Even though there is no water entering the 

model from the left and right boundaries, water flowing inside the modeled system can actually flow into and out of the cells at 

left and right boundaries (but not leave the system). SEAWAT_V4 is a finite difference code with the boundary enforced at the 

node in the center of a cell; therefore, water can flow into/out of the constant temperature cells and convect heat into the model 

domain. This flowing water will convect heat from the boundary as it enters back into the system if the hydraulic conductivity 

is too high. Therefore, the hydraulic conductivity value is decreased from 127 to 1 m/s at the temperature boundaries, allowing 

conductive heat flux to more accurately simulate similar boundary conditions to the Lattice Boltzmann simulations (where the 

boundary is actually enforced at the edge, not the center of a cell). Because of this lowered hydraulic conductivity at the boundaries, the “active” model domain is considered to be the 128 columns, 128 layers, and 1 row where the hydraulic conductivity is 

uniform at 127 m/s (in a 1 m by 1 m cross section). All dimensionless numbers are calculated based on a 1 m by 1 m domain and 

a hydraulic conductivity of 127 m/s.

Figure 4. Change in viscosity with temperature according to equation 3 (oil, in brown) and 

equation 4 (water, in blue).

μ

μ

0.0001

0.001

0.01

0.1

0 5 10 15 20 25 30 35 40

DYNAMIC VISCOSITY, IN KILOGRAM

PER METER PER SECOND

TEMPERATURE, IN DEGREES CELSIUS
Benchmark Problems 9

The Lattice Boltzmann simulations conducted by Guo and Zhao (2005) simulated oil in a homogeneous isotropic porous 

medium and used a generalized version of the Navier-Stokes equation that accounts for all fluid forces, including the resistive 

drag along a boundary (a no-slip boundary). Therefore, they can simulate flow with larger Darcy numbers (eq. A7 in appendix). 

The SEAWAT_V4 simulation results presented here are only valid for conditions with relatively small Darcy numbers (less 

than 10-3) because of the slip boundary condition. Table 2 lists the calculated hydraulic conductivities using two different Darcy 

numbers. The hydraulic conductivity for a Darcy number of 10-2 is extremely high (2,658 m/s), suggesting that explicit representation of boundary friction (such as that provided by the Lattice Boltzmann method) would be required. Therefore, to compare 

SEAWAT_V4 simulation results to those by Guo and Zhao (2005), a Darcy number of 10-4 is used in all the simulations so 

Darcy’s law remains valid.

Guo and Zhao (2005) performed a series of Lattice Boltzmann simulations with and without the effects of viscosity variations 

and using a range of Rayleigh numbers between 103 and 107 (eq. A9 in appendix). They also presented Darcy-Rayleigh, or Ra*, 

numbers for the same benchmark problems (eq. A10 in appendix). Guo and Zhao (2005) likely used both equations A9 and A10 

(in appendix) because they simulate flow in porous media with low to high hydraulic conductivities; therefore, as hydraulic 

conductivity is increased, the simulations moved toward simulating free fluid flow. Thermal diffusivity, used in the Rayleigh 

number, is calculated using equation A11 (in appendix), resulting in a Prandtl number (eq. A15 in appendix) of 0.164.

A total of six simulations were performed with SEAWAT_V4 for comparison with simulations from Guo and Zhao (2005). 

The Darcy number for these simulations was 10-4. Three different Rayleigh numbers were explored: 105, 106, and 107. The 

parameter 

∂

∂



T affects the calculated thermal expansion coefficient t ( β

ρ

ρ t T = ∂

∂





 



 1 ) used in the Rayleigh number. Therefore, by

increasing the change in density with the change in temperature, the Rayleigh number is also increased. Two cases were run with 

each Rayleigh number: (1) a constant viscosity case using the reference viscosity of 0.00595 kg/m/s (table 2), and (2) a variable 

viscosity case using equation 3. Solvers and time-step options for the models are presented in table 2. Models were run until 

temperatures reached an equilibrium configuration.

Results from SEAWAT_V4 are compared with results from Guo and Zhao (2005) as shown in figures 5 and 6. Differences 

between the simulation series are likely caused by differences in boundary conditions (that is, slip vs. no-slip boundaries and 

the decreased hydraulic conductivity boundary in SEAWAT_V4 to reduce convective heat flux). The difference in streamline 

patterns (fig. 5) between the series may also result from differences in visual display of the output. Quantitative values for the 

streamlines (or temperature contours) were not provided in the paper by Guo and Zhao (2005).

In general, the overall streamline and temperature patterns for results from Guo and Zhao (2005) and those using SEAWAT_V4 are 

similar, particularly in the constant viscosity cases, but there are important differences. For constant viscosity, the flow patterns 

are radially symmetric. For the SEAWAT_V4 simulations, when viscosity variations are included, the flow patterns are asymmetric (fig. 5), where the center of the convection cell shifts toward the warmer wall. This shift results from the decrease in viscosity near the warmer wall and increase in viscosity near the cooler wall. In other words, the center of the convection cell tends 

to move toward the warmer wall because the overall resistance to flow is less in the warmer parts of the domain. 

The temperature contours for the Lattice Boltzmann and SEAWAT_V4 simulations appear similar (fig. 6). Looking closely, 

the change in temperature contours between the variable and constant viscosity cases can be seen in figure 6, as well as a difference in temperature contours among the different Ra* numbers (eq. A10 in appendix). For a given Ra*, the isotherms for constant viscosity cases have a slightly higher slope near the boundaries and a steeper slope away from the boundaries than the variable viscosity cases. In figure 6, for the model where Ra* = 10, the temperature isotherms are nearly vertical or have very high 

slopes for both the constant and variable viscosity cases, revealing that conduction plays a significant role in heat transport (the 

reference viscosity used in the calculation of the Ra* is from table 2, which is the viscosity at the reference temperature, 21  ̊C, 

the average of the end-member temperatures). For the variable viscosity cases, the temperature contours are slightly less vertical, suggesting there is more heat transfer by way of convection than the similar constant viscosity case with the same Ra*; this 

pattern is most clearly seen where Ra* = 10. As Ra* increases there is an increased tendency for convection versus conduction. 

This is also indicated by the flow patterns (fig. 5), which become more elongated in the horizontal direction as Ra* increases 

from 10 to 1,000. A similar pattern is observed with the isotherms (fig. 6); the isotherms are more horizontal as Ra* increases.

This benchmark problem was used to test the capability of SEAWAT_V4 to simulate variations in the viscosity of oil from 

differences in temperature. Results from SEAWAT_V4 and Lattice Boltzmann methods reported by Guo and Zhao (2005) differ 

in some of the streamline patterns, but are similar for temperature patterns. Moreover, both models show similar responses in the 

simulation of temperature when the effects of viscosity variations are explicitly represented. The consistent response between 

both approaches indicates that viscosity variations can affect the shape of a convection cell as well as increase the relative contribution of convective to conductive flow. Laboratory tests or simulations by other codes may need to be done before this problem 

can be considered a standard benchmark problem. 
10 Application of SEAWAT to Select Variable-Density and Viscosity Problems

Figure 6. Temperature isotherms for both variable viscosity and constant viscosity cases. The gray/black isotherms are from Guo and 

Zhao (2005) and are shown adjacent to temperature isotherms from SEAWAT_V4 in color. A Darcy number of 10-4 was used for all of the 

simulations (calculated from eq. A7 in appendix). Ra is Rayleigh Number from equation A9 in appendix, and Ra* is Rayleigh Number from 

equation A10 in appendix.

Figure 5. Streamline results for both variable viscosity and constant viscosity cases. The gray/black streamlines are from Guo and 

Zhao (2005) and are shown adjacent to resulting streamlines from SEAWAT_V4 in blue. A Darcy number of 10-4 was used for all of the 

simulations (calculated from eq. A7 in appendix). Ra is Rayleigh Number from equation A9 in appendix, and Ra* is Rayleigh Number from 

equation A10 in appendix.

Ra = 105, Ra*= 10

Ra = 107, Ra*= 1,000

Ra = 106, Ra*= 100

NO VARIABLE VISCOSITY

NO VARIABLE VISCOSITY VARIABLE VISCOSITY

VARIABLE VISCOSITY

NO VARIABLE VISCOSITY VARIABLE VISCOSITY

Ra = 105, Ra*= 10

Ra = 107, Ra*= 1,000

Ra = 106, Ra*= 100

Ra = 105, Ra*= 10

Ra = 106, Ra*= 100

Ra = 107, Ra*= 1,000

NO VARIABLE VISCOSITY

NO VARIABLE VISCOSITY VARIABLE VISCOSITY

VARIABLE VISCOSITY

NO VARIABLE VISCOSITY VARIABLE VISCOSITY

Ra = 105, Ra*= 10

Ra = 106, Ra*= 100

Ra = 107, Ra*= 1,000
Benchmark Problems 11

Problem 3: Horton-Rogers-Lapwood (HRL) Convection

The HRL convection problem is based on the pioneering work of Rayleigh (1916), Horton and Rogers (1945), and 

Lapwood (1948). The HRL problem consists of an infinitely long layer of fully saturated porous media with a specified height of 

H (fig. 7). The temperature assigned to the lower boundary is fixed at a value higher than the temperature assigned for the upper 

boundary. Thus, water density decreases with depth. For the HRL problem, one of two things will occur depending on the values 

used to define the system. One possibility is that the fluid may remain stationary, and heat will conduct from the lower boundary 

to the upper boundary. Under this condition, the entire layer is characterized by a linear vertical temperature gradient. The other 

possibility is that convection cells will form. Horton and Rogers (1945) showed that convection will occur when Ra* (eq. A13 in 

appendix) is greater than 4p2 (a value of approximately 39.48).

Instead of trying to approximate an infinitely long layer, Weatherhill and others (2004) showed that the HRL problem can be 

represented as a rectangular box provided that the length of the box, L, is an integer multiple of the height of the box, H (fig. 7). 

Or, for simplification, the aspect ratio, A, defined as (Weatherhill and others, 2004):

A L

H = (5)

must be an even integer if 4p2 is to apply as the critical Ra* value. If A is an even integer, then the central vertical plane of an 

upwelling or downwelling plume will coincide with a boundary and the critical Ra* value is 4p2. If A is not an even integer, 

convective flow will still occur, but the critical Ra* number will not be 4p2 because the center zones of upwelling/downwelling 

will not be at the boundary.

The HRL convection problem was simulated with SEAWAT_V4 to test the ability of the code to simulate thermally driven 

convection. The physical aspects of the problem are shown schematically in figure 8 with the model parameters, grid dimensions, 

and solution schemes given in table 3. The model simulates freshwater flow in a domain of 1 row, 416 columns, and 54 layers. 

Figure 7. Design of the Horton-Rogers-Lapwood benchmark problem where the length of the box is 2 times 

(or 4 times) the height.

Figure 8. Design of the Horton-Rogers-Lapwood problem in SEAWAT_V4. The colors from blue to red show 

the initial concentrations for the simulations from 0 °C in blue to 100 °C in red. Height, length, and boundary 

conditions are also displayed.

H = 156 m

T = 100 °C L = 8*H = 1248 m

Constant Temperature Boundary T = 0 °C

Constant Temperature Boundary

perturbation

Constant

Head

Constant

No Flow Head

No Flow

H

L = 2*H

Hot Surface

Cool Surface

L = 4*H
12 Application of SEAWAT to Select Variable-Density and Viscosity Problems

Table 3. Input parameters for the Horton-Rogers-Lapwood problem simulated in SEAWAT_V4.

Input Parameter Value Units Comments

Number of columns 416 - Assigned

Number of rows 1 - Assigned

Number of layers 54 - Assigned

∆x(DELR) 3 m Assigned

∆y(DELC) 3 m Assigned

∆z(DELZ) 3 m Assigned

K, where Ra* = 39 0.0385 m/d Hydraulic conductivity

K, where Ra*= 40 0.03949 m/d Hydraulic conductivity

K, where Ra*= 200 0.19744 m/d Hydraulic conductivity

k, where Ra* = 39 4.54 × 10-14 m2 Permeability

k, where Ra*= 40 4.66 × 10-14 m2 Permeability

k, where Ra*= 200 2.33 × 10-13 m2 Permeability Q C= ⋅ OND h∆

L 0 m Longitudinal dispersivity

COND

T 0 m Transverse dispersivity

Sy None-steady state - Specific yield

S None-steady state - Storage

∆h 0.1 - Porosity

D* 0.308 m2

/d Thermal diffusivity

COND KWH

L

f = 



water 1,000 kg/m3 Density of water at 4 o

C

∂T

∂ -2 kg/(m3 oC) Density change with temperature

 0.002 1/ oC Thermal expansion coefficient

v0 0.0864 m2

/d Reference kinematic viscosity

μ0 86.4 kg/(m d) Reference dynamic viscosity

Tref 4 o

C Reference temperature for reference 

dynamic viscosity

g 7.32 × 1010 m/d2 Acceleration due to gravity

Cool temperature boundary 0 oC Assigned

Hot temperature boundary 100 oC Assigned

Matrix solution technique for flow PCG2 - Assigned

Head convergence value 1 × 10-10 m Assigned

Flow convergence value 1 kg/d Assigned

Advection term Implicit finite difference-central in 

space weighting - Assigned

Dispersion and source terms Implicit finite difference; GCG, 

SSOR - Assigned

Time-step length Calculated using courant of 0.01 s Assigned

Temperature convergence value 1 × 10-12 - Assigned
Benchmark Problems 13

The top and bottom boundaries were set as constant temperature boundaries with a decreased hydraulic conductivity of 

1 x 10-5 m/d. This low hydraulic conductivity value was assigned only to the top and bottom layers in order to minimize the convective heat flux across the temperature boundaries; by definition, conduction should be the sole heat transport mechanism at the 

boundary. Because of the decreased hydraulic conductivity in the top and bottom layers, H is calculated as the distance between 

the bottom of layer 1 and the top of layer 54 (156 m). Constant-head boundary cells, with an arbitrary value, were placed at the 

upper left and upper right corners (fig. 8) to ensure convergence of the flow and transport equations and facilitate a unique solution 

to the benchmark problem. The left and right hand sides of the model are insulated no-flow boundaries. The model only simulates 

heat transfer; solute transport is not included in this simulation. Although temperatures of zero and 100 oC are imposed at the 

boundaries, the processes of evaporation and condensation (i.e. freezing or boiling) are not considered.

Three simulations were performed with SEAWAT_V4 using hydraulic conductivity values leading to Ra* values of 39, 40, and 

200 (table 3). These simulations were identical except that the initial temperatures for the simulation for which Ra* = 200 were set 

as the temperatures resulting from the simulation for which Ra* = 40. For the simulation for which Ra* = 40, the temperature at a 

single node was slightly perturbed to initiate convection in the model. 

Results show that for the first simulation (Ra*= 39) there is no convection, and all heat transport is solely through conduction 

(fig. 9A). As the Rayleigh number increases (fig. 9B and 9C), convection cells form. Convection can clearly be seen in the 

velocity field as shown in figure 10 for the Ra* = 40 simulation. Upwelling brings the warmer water closer to the top boundary; downwelling occurs and brings warmer water closer to the bottom boundary (fig. 9). According to Weatherhill and others 

(2004), the wavelength, l, of the convection cells is 2H and the number of convection cells should equal L/H (the aspect ratio). 

Thus, for the values used here, there should be eight convection cells, which is the number observed in figure 9B and 9C and 

in figure 10.

Figure 9. Temperature results for the Horton-Rogers-Lapwood convection problem simulated with SEAWAT_V4. 

(A) Simulation with Ra* = 39, (B) simulation with Ra* = 40, and (C) simulation with Ra* = 200. The colors from blue 

to red are the same as for the initial concentrations for the simulations from 0 °C in blue to 100 °C in red.

Figure 10. Velocity vectors for the Horton-Rogers-Lapwood problem where Ra* = 40.

A

B

C
14 Application of SEAWAT to Select Variable-Density and Viscosity Problems

Ataie-Ashtiani and Aghayi (2006) analyzed the suitability of the HRL convection problem for testing variable-density 

modeling codes. Although the analysis of the critical Ra* number may be useful for benchmarking a code, the wavelength test 

(l = 2H) does not apply if calculated Ra* numbers are significantly higher than 4p2

. The significantly larger Ra* numbers can 

cause more convection cells to develop in the model (greater than the aspect ratio, A). Similar results to Ataie-Ashtiani and 

Aghayi (2006) were also observed with SEAWAT_V4 simulations of the HRL problem; however, the research here revealed that 

initial temperatures affected the number of convection cells. Results from our simulation with an Ra* number of 200 revealed 

that if initial temperatures were stable such as in figure 9A, the number of convection cells would be greater than the number 

expected, 8, which is what Weatherhill and others (2004) observed. However, if the initial temperatures were assigned from the 

results of the Ra* = 40 simulation (fig. 9B), then eight convection cells would result (as shown in figs. 9C and 10). The conflicting steady-state results for the simulation with an Ra* = 200 likely result from numerical errors, as Ataie-Ashtiani stated, “It 

seems that the steady state solution of the Infinite Horizontal Box (IHB) problem is very sensitive to numerical errors (including 

truncation and round off errors). So even when steady state solutions of the IHB problem converge to a specific solution... still 

it cannot be assured that this stationary solution is unique. Such a phenomenon is observed in some other numerical problems 

which were based on formation of convective instabilities in porous media.”

The results shown indicate that SEAWAT_V4 can correctly simulate heat conduction and convection according to results from 

an analytical solution (Nield and Bejan, 1999). However, when simulating density-dependent flow systems with high Ra* numbers, 

it is important to note that initial temperatures can potentially alter the results. Regardless, results from this benchmark problem 

lend further confidence in the application of SEAWAT_V4 to problems that involve heat conduction and convection.

Problem 4: Double-Diffusive Finger Convection

Double-diffusive finger convection is an interesting mixing process because density-driven fingering can occur from an initially 

stable density configuration (a less-dense fluid on top of a more-dense fluid). Consider a case where two miscible solutions with 

different viscosities are separated by a sharp interface. The solute in the overlying fluid has a molecular diffusion coefficient that 

is less than the molecular diffusion coefficient of the solute in the underlying fluid. Under these conditions, double-diffusive finger 

convection can occur in response to small perturbations along the interface, even though the overlying fluid is less dense than the 

fluid on the bottom. As explained by Hughes and others (2005), “a parcel of fluid perturbed downward across the interface takes 

on solute mass from the surrounding fluid faster than it diffuses solute mass so the parcel continues to fall.” 

Pringle and others (2002) documented the process of double-diffusive finger convection in the laboratory using a Hele-Shaw 

cell (figs. 11 and 12). Prior to convection, the cell consisted of a less dense sucrose solution overlying a more dense sodium-chloride 

(NaCl) solution. The Hele-Shaw cell was tilted to an angle of 25 degrees from the horizontal plane (figs. 11 and 12). A conservative dye was mixed with the NaCl solution to allow for visual observation of the finger formation and convection. The dye was 

assumed to have a minimal effect on the density of the NaCl solution. A sequence of laboratory photographs clearly documents the 

formation and convection of fingers (fig. 13). Actual times, t, are related to the dimensionless times, t*, listed in figure 13 by the 

following equation:

(6)

where t is elapsed time in seconds and H is the height of the cell (0.1625 m), and D* is the diffusion coefficient. By the end of 

the experiment, the dye is almost entirely mixed throughout the cell.

This double-diffusive finger convection observed in the laboratory experiment was simulated using a modified version of 

the SUTRA program, SUTRA-MS (Hughes and Sanford, 2004; Hughes and others, 2005). A similar approach is used here to test 

SEAWAT_V4 and its ability to simulate finger convection with multiple diffusion coefficients and variable viscosity. The input 

parameters, grid dimensions, and solution schemes for the numerical simulation are shown in table 4. These parameters correspond 

with those listed by Pringle and others (2002). The simulation domain is the same size as the Hele-Shaw cell (0.2541 m long by 

0.1625 m high). The model domain is surrounded on all four sides by no-fluid-flux and zero diffusive flux conditions (fig. 12). 

Constant-head cells with an arbitrary head value of 0.1625 m were placed at the top left and top right corners of the model. Random 

perturbations were applied to the initial concentrations along the interface following the procedure of Hughes and Sanford (2004).

The laboratory photographs shown in figure 13 can be directly compared to the dye concentrations from the SEAWAT_V4 

simulation shown in figure 14. The similarity of results indicates that SEAWAT_V4 is capable of representing the complex process of double-diffusive finger convection. This is an additional demonstration of the correct implementation of multiple diffusion 

coefficients, variable viscosity, and variable density flow in the SEAWAT_V4 code.
Benchmark Problems 15

12 bit CCD camera

2033 x 2048 pixel array

with 4096 gray level range

Light

box

Hele-Shaw cell

Concrete support

structure

Figure 11. Hele-Shaw cell set up by Pringle and others (2002). Note that the angle of 

the Hele-Shaw cell from the vertical is the reason that gravity in the simulation is set 

to 4.14 meters per second squared.

Figure 12. Initial conditions in the Hele-Shaw cell.

Sucrose

LENGTH = 0.2541 METER

HEIGHT = 0.1625 METER

CONSTANT HEAD = 0.1625 METER

Dye

and

NaCl
16 Application of SEAWAT to Select Variable-Density and Viscosity Problems

C/Co

0 1

A B C

D E F

G H I

J K L

EXPLANATION

C/Co

0 1

A B C

D E F

G H I

J K L

EXPLANATION

Figure 13. Results from Pringle and others (2002) at (A) t* = 4.06 x 10-5, (B) t* = 1.29 x 10-4, (C) t* = 3.96 x 10-4, 

(D) t* = 3.35 x 10-4, (E) t* = 4.35 x 10-4, (F) t* = 5.36 x 10-4, (G) t* = 6.03 x 10-4, (H) t* = 7.37 x 10-4, (I) t* = 8.04 x 10-4, 

(J) t* = 1.04 x 10-3, (K) t* = 1.78 x 10-3, and (L) t* = 3.19 x 10-3. Color sequence is relative concentrations of the 

dye (representative of the sodium-chloride solution).

Figure 14. Results from SEAWAT_V4 at (A) t* = 4.06 x 10-5, (B) t* = 1.29 x 10-4, (C) t* = 3.96 x 10-4, 

(D) t* = 3.35 x 10-4, (E) t* = 4.35 x 10-4, (F) t* = 5.36 x 10-4, (G) t* = 6.03 x 10-4, (H) t* = 7.37 x 10-4, 

(I) t* = 8.04 x 10-4, (J) t* = 1.04 x 10-3, (K) t* = 1.78 x 10-3, and (L) t* = 3.19 x 10-3. Color sequence is 

relative concentrations of the dye (representative of the sodium-chloride solution).
Benchmark Problems 17

Table 4. SEAWAT_V4 input parameters for the Hele-Shaw simulation.

Input Parameter Value Units Comments

Number of columns 1,024 - Assigned

Number of rows 1 - Assigned

Number of layers 656 - Assigned

∆x(DELR) 2.48 × 10-4 m Assigned

∆y(DELC) 2.48 × 10-4 m Assigned

∆z(DELZ) 2.48 × 10-4 m Assigned

K 1.0784 × 10-2 m/s Hydraulic conductivity

D* NaCl 1.477 × 10-9 m2

/s Molecular diffusion

D* Sucrose 4.878 × 10-10 m2

/s Molecular diffusion

D* Dye 5.670 × 10-10 m2

/s Molecular diffusion

Q C= ⋅ OND h∆

NaCl 0.0 kgNaCl/kgfluid Minimum concentration

CONDNaCl 0.03463 kgNaCl/kgfluid Maximum concentration

∆hSucrose 0 kgSucrose/kgfluid Minimum concentration

COND KWH

L

f = 



Sucrose 0.05234 kgSucrose/kgfluid Maximum concentration





f Dye 0.0 kgDye/kgfluid Minimum concentration





f Dye 1 kgDye/kgfluid Maximum concentration

h h h

x x = x x h −− [ − ] + 1 0

0 ln( ) ln( ) ln( ) ln( )

max min

min

L 0 m Longitudinal dispersivity

( ) . . T T = ⋅ − 0 168 1 0868

T 0 m Transverse dispersivity

 T T ( ) = ⋅ +





 





0.00002394 10

248 37

133 15

.

.

1.0 - Porosity

∂

∂



Twater 998 kg/m3 Density of freshwater

t NaCl 689 kgfluid

2

/(kgNaCl m3) Density change with concentration

β ρ

ρ t T = ∂

∂





 



 1

Sucrose 371 kgfluid

2

/(kgSucrose m3) Density change with concentration

A L

H =

Dye 0 kgfluid

2

/(kgDye m3) Density change with concentration

μ 0.001 kg/(m s) Reference dynamic viscosity

∂

∂



T

NaCl 1.59 × 10-3 kgfluidkg/(kgNaCl m s) Viscosity change with concentration

∂

∂



T

Sucrose 2.75 × 10-3 kgfluidkg/(kgSucrose m s) Viscosity change with concentration

∂

∂



T Dye 0 kgfluidkg/(kgDye m s) Viscosity change with concentration

g 4.14 m/s2 Acceleration due to gravity

Matrix solution technique for flow PCG2 - Assigned

Head convergence value 1.0 × 10-8 m Assigned

Flow convergence value 1 kg/s Assigned

Length of stress period 57,600 s Steady-state model, one stress period

Advection term Implicit finite difference-central in 

space weighting - Assigned

Dispersion and source terms Implicit finite difference; GCG - Assigned

Maximum transport step length 10 s Assigned

Concentration convergence value 1 × 10-10 - Assigned

∂μ

∂μ

∂μ

∂C

∂C

∂C
18 Application of SEAWAT to Select Variable-Density and Viscosity Problems

Problem 5: Elder Problem

The classic laboratory experiment reported by Elder (1967) has been widely used to test numerical codes (that is, Voss and 

Souza, 1987; Guo and Langevin, 2002). The original Elder problem consists of two-dimensional, thermally driven convection 

caused by a heat source centrally located at the bottom of the domain (fig. 15). As the water is heated, it rises causing convection 

within the system.

To test variable-density groundwater flow and solute transport codes, Voss and Souza (1987) modified the original Elder 

problem. Instead of heating water from below, convection in the modified version is driven by a constant-concentration solute boundary at the top that causes the formation of dense brine. The modified Elder problem is based on a set of parameters 

(table 5) that result in an identical Ra* number as used in the original experiment (from eq. A13 in appendix). 

The original Elder problem is used here to further test the ability of SEAWAT_V4 to simulate thermally driven convection. 

However, there is no analytical solution to the Elder problem; therefore, SEAWAT_V4 results are compared with (1) the original 

numerical results from Elder (1967), (2) results from a previous version of SEAWAT (Guo and Langevin, 2002), and (3) results 

from SUTRA (Voss, 1984; Voss and Souza, 1987; Prasad and Simmons, 2005).

The design of the Elder problem simulated in SEAWAT_V4 is shown in figure 15 with the model parameters, grid dimensions, 

and solution schemes given in table 5. This numerical model simulates freshwater flow with variations in temperature for 

7,300 days (or 20 years). The top boundary was set at a constant temperature of zero. A relatively low hydraulic conductivity 

value of 1 x 10-5 m/d was also assigned to the top layer to minimize the convective heat flux from the boundary. The bottom 

boundary has a constant temperature of 100 °C from columns 23 to 66 and also has a relatively low hydraulic conductivity 

value. A constant head of zero is assigned at the bottom left and bottom right corners to ensure convergence of the flow solution. 

The initial temperature is set everywhere to zero. The temperature values of zero and 100 °C are arbitrary in this case; the value 

for ∂

∂



T

was calculated such that the Ra* number, from equation A13 in the appendix, is equal to 400, the value for the original 

Elder problem. (As stated previously, the actual temperature values are not important for certain benchmark problems. What is 

pertinent is that the correct Ra* of 400 is calculated for the “active” model domain in figure 15. Different temperature values 

could have been used with different values for ∂

∂



T or the permeability modified to calculate an Ra* number equal to 400.)

Results from SEAWAT_V4 are similar to the original Elder results, SUTRA results from Voss and Souza (1987), and results

from a previous version of SEAWAT (fig. 16). There are visual differences in the isotherm (or isochlor) contours and timing of the 

maximum penetration depth of the 60-percent isotherm. The discrepancy in results is likely due to differences in grid dimensions, 

time-step length, or solution schemes. Variations in discretization and solution schemes have been shown to affect simulated results 

for the Elder problem revealing the sensitive nature of the Elder problem. This has been studied extensively in the following papers: 

Frolkovic and Schepper (2000), Diersch and Kolditz (2002), Simpson and Clement (2004), Al-Maktoumi and others (2007).

Porous medium

Warm temperature boundary

Cool temperature boundary

Constant

Head

Constant

Head

L = 600 METERS

H = 150 METERS

L = 300 METERS

Figure 15. Schematic representation of Elder’s original experiment as simulated in SEAWAT_V4. The constant 

head cells at the bottom left and right hand side are used in the numerical simulation; however, do not exist in 

Elder’s original lab set-up.
Benchmark Problems 19

Table 5. Input parameters for the Elder problem simulated in SEAWAT_V4.

Input Parameter Value Units Comments

Number of columns 88 - Assigned

Number of rows 1 - Assigned

Number of layers 52 - Assigned

∆x(DELR) 6.818 m Assigned

∆y(DELC) 6.818 m Assigned

∆z(DELZ) 3 m Assigned

K 0.411 m/d Hydraulic conductivity

k 4.85 × 10-13 m2 Permeability

Q C= ⋅ OND h∆

L 0 m Longitudinal dispersivity

COND

T 0 m Transverse dispersivity

∆h 0.1 - Porosity

D* 0.308 m2

/d Bulk thermal diffusivity

COND KWH

L

f = 



water 1,000 kg/m3 Density of water at 4 o

C





f -2 kg/(m3 oC) Density change with temperature





f

0.002 1/ oC Thermal expansion coefficient

v0 0.0864 m2

/d Reference kinematic viscosity

μ0 86.4 kg/(m d) Reference dynamic viscosity

Tref 0 o

C Reference temperature for reference 

dynamic viscosity

g 7.323 × 1010 m/d2 Acceleration due to gravity

Cool Temperature Boundary 0 oC Assigned

Hot Temperature Boundary 100 oC Assigned

Matrix solution technique for flow PCG2 - Assigned

Head convergence value 1 × 10-7 m Assigned

Flow convergence value 1,000 kg/d Assigned

Flow time-step length 60 time steps with time step multiplier of 5 d Assigned

Advection term Implicit finite difference: upstream weighting - Assigned

Dispersion and source terms Implicit finite difference; GCG, SSOR - Assigned

Transport time-step length Calculated using courant of 0.1 d Assigned

Temperature convergence value 1 × 10-10 - Assigned
20 Application of SEAWAT to Select Variable-Density and Viscosity Problems

Figure 16. Results of SEAWAT_V4 compared to the original laboratory results from Elder (1967), SUTRA 

(upside-down), and an older version of SEAWAT (upside-down). Modified from Guo and Langevin, 2002.

20

60

60

20

20

60

20

1 YEAR

10 YEARS

2 YEARS

15 YEARS

3 YEARS

20 YEARS

20

60

20

20

60

60

60

20

60

SEAWAT_V4 LINE OF RELATIVE TEMPERATURE, IN PERCENT

SEAWAT LINE OF RELATIVE SALINITY CONCENTRATION, IN PERCENT (GUO AND LANGEVIN, 2002)

SUTRA LINE OF RELATIVE SALINITY CONCENTRATION, IN PERCENT (VOSS AND SOUZA, 1987)

ELDER LINE OF RELATIVE SALINITY CONCENTRATION, IN PERCENT (VOSS AND SOUZA, 1987)

EXPLANATION

20

20

20

20

20

20

60

20

60

20

60

60
Benchmark Problems 21

Prasad and Simmons (2003; 2005) and Elder (1967) calculated the Nusselt number (eq. A16 from the appendix) as a function 

of time for their Elder simulations. The temporal variations in the Nusselt number for SEAWAT_V4, the SUTRA results from 

Prasad and Simmons (2003), and the results from Elder (1967) are shown in figure 17. The results are similar between the 

three models; however some differences are evident. These differences could be due to using different modeling codes, spatial 

discretization, temporal discretization, or solution schemes. The similarities between the graphs reveal how the simulations are 

alike in the timing and amount of the transfer of energy (or solute in Prasad and Simmons, 2003), as well as the energy transfer 

mechanism. Overall the three different models have a similar total flux compared to the expected conductive or diffusive flux. 

(For more details, see the discussion of the Nusselt number in the appendix.) 

The results from the Elder problem simulation demonstrate that SEAWAT_V4 represents complex convective flow patterns 

of thermally driven variable-density flow in a manner similar to other codes that have simulated this problem.

Problem 6: Henry-Hilleke Problem

Henry and Hilleke (1972) describe a temperature and salinity variation of the Henry (1964) problem, which has been used 

as a standard benchmark problem for testing variable-density groundwater flow and solute transport codes. Henry and Hilleke 

(1972) described a laboratory experiment and numerical modeling effort designed to represent the coastal part of a carbonate 

platform where groundwater flow is affected by salinity and temperature variations. 

The Henry and Hilleke (1972) problem has also been used as a numerical benchmark problem that is designed to simulate 

the laboratory experiment. Hughes and Sanford (2004) use this benchmark problem to test two other variable-density modeling 

codes, HST3D (Kipp, 1987; 1997) and SUTRA-MS (Hughes and Sanford, 2004), which include temperature and concentration 

effects on variable-density flow.

6

10

14

18

0 5 10 15 20

TIME, IN YEARS

8

12

16

NUSSELT NUMBER 20

Figure 17. The Nusselt number versus time showing the heat flux across the bottom boundary in 

SEAWAT_V4 in black. This is compared to the concentration flux across the top boundary (in the 

inverted Elder problem) simulated using SUTRA from Prasad and Simmons (2003) in red, and the 

heat flux across the bottom boundary from the original Elder problem in blue (Elder, 1967).
22 Application of SEAWAT to Select Variable-Density and Viscosity Problems

Figure 18. Setup of the Henry-Hilleke problem showing the boundary conditions 

which include freshwater input, seawater boundary, and temperature boundaries.

The design of the Henry and Hilleke problem is shown in figure 18. Model parameters, grid dimensions, and solution schemes 

are given in table 6. The cross section is 1.025 m by 1.025 m in the x-z vertical plane (the aspect ratio, from equation 5, equals 

one). However, the rectangular area defined by the centers of the four corner cells is 1 m by 1 m. This allows the SEAWAT_V4 

solution domain, in which nodes are at cell centers, to better correspond with the solution domain for the SUTRA-MS representation, in which nodes are at the corners of elements. SUTRA-MS is a finite-element code and allows for elements to be placed 

along boundaries. The right boundary has seawater concentrations of 35.7 kg/m3, and the left boundary has a constant influx 

of freshwater. The model has constant temperatures cells surrounding the domain (fig. 18). The right boundary has a constant 

temperature of 5 °C. The lower left corner has a constant temperature of 50 °C, the upper left corner has a constant temperature 

of 38.75 °C, and temperatures vary linearly between them. Temperatures also vary linearly along the top and bottom boundaries. Equations used to calculate the constant temperature at each cell along the boundaries are shown in figure 18. The top and 

bottom boundaries have a relatively low hydraulic conductivity value of 432 m/d to minimize the convective heat flux from the 

temperature boundaries. 

SEAWAT_V4 simulation results are similar to SUTRA-MS and HST3D (Hughes and Sanford, 2004; Thorne and others, 2006), 

as well as the original Henry and Hilleke problem (1972). Results from SUTRA-MS, HST3D, Henry-Hilleke, and SEAWAT_V4 

are shown in figures 19 to 21. Simulated concentration and temperature contours using SUTRA-MS and SEAWAT_V4 are almost 

identical (figs. 19 and 20), and are similar to those using HST3D. Resulting velocity vectors using SUTRA-MS and SEAWAT_V4 

also show similar patterns (fig. 21). 

Concentration and temperature results using SEAWAT_V4, SUTRA-MS, and HST3D are more similar to each other than 

the numerical solution from Henry-Hilleke, although the general shape of the contours is similar. The differences are likely 

because the Henry and Hilleke (1972) numerical solution is a “simplified form of the variable-density flow and transport equation” (Hughes and Sanford, 2004) that only becomes similar to the full equations used in SUTRA-MS and SEAWAT_V4 when 

the aspect ratio (eq. 5) is much smaller than one. Hughes and Sanford (2004) give a more detailed discussion of the effects of 

mesh resolution on the numerical simulation of the Henry and Hilleke problem. Their results show that the smaller the aspect 

ratio, the more similar the results are from Henry and Hilleke (1972) and SUTRA-MS. (The numerical results are not compared 

to the laboratory results because the dimensions are different between the two.)

These results indicate SEAWAT_V4 produces results similar to other codes for the Henry-Hilleke problem, and thus is as 

capable of simulating simultaneous solute and heat transport coupled with variable-density groundwater flow.

T(1,1) = 5 C°

T(0,0) = 50 °C T(1,0) = 5 C°

T(0,1) = 38.75 °C

Porous medium

L = 1 METER

H = 1 METER

Freshwater = 7.2 m /d Qin 3

Seawater Concentration

Boundary = 35.7 kg/m3

T(x,1) = 38.75 - 33.75x

T(x,0) = 50 - 45x

T(0,y) = 50 - 11.25y

T(1,y) = 5
Benchmark Problems 23

Table 6. Input parameters for the Henry-Hilleke problem simulated in SEAWAT_V4.

Input Parameter Value Units Comments

Number of columns 41 - Assigned

Number of rows 1 - Assigned

Number of layers 41 - Assigned

∆x(DELR) 0.025 m Assigned

∆y(DELC) 1 m Assigned

∆z(DELZ) 0.025 m Assigned

K 864 m/d Hydraulic conductivity

Q C= ⋅ OND h∆

L 0 m Longitudinal dispersivity

COND

T 0 m Transverse dispersivity

∆h 0.35 - Porosity

D*m 2.0571 m2

/d Molecular diffusivity

D* 20.571 m2

/d Thermal diffusivity

COND KWH

L

f = 



f 1,000 kg/m3 Density of freshwater





f

s 1,025 kg/m3 Density of seawater





f 0.7 - Density change with concentration

h h h

x x = x x h −− [ − ] + 1 0

0 ln( ) ln( ) ln( ) ln( )

max min

min

-0.375 kg/(m3 oC) Density change with temperature

v0 0.0864 m2

/d Reference kinematic viscosity

μ0 86.4 kg/(m d) Reference dynamic viscosity

Tref 0 o

C Reference temperature for reference 

dynamic viscosity

Qin 7.2 m3

/d

Total flow in along left hand boundary. 

Split into 41 nodes (0.1756 m3/d in 

each node)

Matrix solution technique for 

flow PCG2 - Assigned

Head convergence value 1 × 10-7

m Assigned

Flow convergence value 1 × 10-7 kg/d Assigned

Length of stress period 0.2 d One stress period, one flow time step

Advection term TVD - Assigned

Dispersion and source terms Implicit finite difference; 

GCG, SSOR - Assigned

Transport time-step length 6.9444 × 10-4 d Assigned

Concentration/temperature convergence value 1 × 10-6 - Assigned
24 Application of SEAWAT to Select Variable-Density and Viscosity Problems

A

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

B

DISTANCE, IN METERS

HEIGHT, IN METERS

1.0

0.8

0.6

0.4

0.2

0

0.05

0.20

0.50

0.80

0.95

0.05

0.20

0.50

0.80

0.95

EXPLANATION

EXPLANATION

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

DISTANCE, IN METERS

HEIGHT, IN METERS

A

B

1.00

0.75

0.50

0.25

0

9.50

11.75

16.25

27.50

38.75

43.25

EXPLANATION

EXPLANATION

9.50

11.75

16.25

27.50

38.75

43.25

Figure 19. Modeling results showing percent seawater concentrations from (A) SUTRA-MS in color with the solid black 

lines, Henry and Hilleke numerical solution shown as the 50 percent contour and the solid red line, and HST3D results 

shown as the dashed black line (Hughes and Sanford, 2004), and (B) SEAWAT_V4 in color with the solid black lines. 

Figure 20. Modeling results showing temperatures in degrees Celsius from (A) SUTRA-MS in color with the solid 

black lines, Henry and Hilleke numerical solution shown as the 27.5 °C contour in the solid red line, and HST3D results 

shown as the dashed black line (Hughes and Sanford, 2004), and (B) SEAWAT_V4 in color with solid black lines.
Summary and Conclusions 25

Summary and Conclusions

Numerical models of groundwater systems are often used as planning tools for improving water supply and management, 

as well as understanding groundwater flow processes. Simulating variable-density groundwater systems such as coastal 

aquifers, which include saltwater and freshwater, wastewater disposal sites, salt lakes, or deep aquifers affected by geothermal 

heating, requires the use of a numerical modeling code that solves the variable-density flow equation. This equation needs to 

include not only the effects of concentration on density, but also temperature. 

Depending on the temperature or concentration of the water, the effects of these variables on viscosity may need to be 

included. SEAWAT_V4 supports equations of state for fluid density and viscosity. The code can be used to simultaneously 

simulate salinity and temperature effects on variable-density flow in water, as well as in oil (as a single phase). Variations in 

viscosity from changes in temperature or concentration have also been included, as well as the ability to implement distinct 

diffusion coefficients to accommodate multiple species. Because thermal diffusivity can be an order of magnitude higher 

than molecular diffusion, the capacity to include multiple diffusion coefficients is pertinent when simulating salinity and 

temperature concurrently.

This report documents the first thorough verification of SEAWAT_V4 by comparing results of simulations with published 

numerical modeling, analytical, or laboratory problems for six benchmark problems. Results from the benchmark problems 

reveal that SEAWAT_V4 reasonably simulates the following processes: (1) density-dependent flow due to changes in temperature and/or concentration of one or more species, (2) temperature or concentration-dependent fluid viscosity, and (3) molecular 

diffusion of multiple species (including solute and heat). Results from the benchmark problems also reveal that SEAWAT_V4 

reasonably simulates convective versus conductive transport of heat/energy through a porous medium. 

In conclusion, SEAWAT_V4, has been expanded to simulate single phase liquids with concentrations from fresh to brine. 

Included in the new version is the ability to simulate temperature and concentration effects on density and viscosity. This 

expanded version of the code could be used to aid in solving various field problems such as water supply and waste disposal. 

Figure 21. Modeling results showing magnitude and direction for velocity vectors 

for (A) SUTRA-MS (Hughes and Sanford, 2004), and (B) SEAWAT_V4.

A B

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

DISTANCE, IN METERS

HEIGHT, IN METERS

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9
26 Application of SEAWAT to Select Variable-Density and Viscosity Problems

References Cited

Afify, A.A., 2007, Effects of temperature-dependent viscosity with Soret and Dufour numbers on non-Darcy MHD free convective 

heat and mass transfer past a vertical surface embedded in a porous medium: Transport in Porous Media, v. 66, p. 391-7401.

Al-Maktoumi, A., Lockington, D.A., and Volker, R.E., 2007, SEAWAT 2000: Modeling unstable flow and sensitivity to discritization 

levels and numerical schemes: Hydrogeology Journal; Earth and Environmental Science, DOI: 10.1007/s10040-007-0164-2.

Ataie-Ashtiana, B., and Aghayi, M. M., 2006, A note on benchmarking of numerical models for density dependent flow in 

porous media: Advances in Water Resources, v. 29, p. 1918-1923.

Bird, R.B., Stewart, W.E., and Lightfoot, E.N., 2006, Transport Phenomena (2nd Edition): New York, John Wiley & Sons. 

Cooper, H. H., Jr., 1959, A hypothesis concerning the dynamic balance of fresh water and salt water in a coastal aquifer: Journal 

of Geophysical Research. v. 64, p. 461-467.

Diersch, H.J.G., and Kolditz, O., 2002, Variable-density flow and transport in porous media: Approaches and Challenges: 

Advances in Water Resources, v. 25; p. 899-944.

Elder, J.W., 1967, Transient convection in a porous medium: Journal of Fluid Mechanics, v. 27, no. 3, p. 609-623.

Frolkovič, P., and De Schepper, H., 2000, Numerical modelling of convection dominated transport coupled with density driven 

flow in porous media: Advances in Water Resources, v. 24, no.1, p. 63-72.

Gary, J., Kassory, D.R., Tadjeran, H. and Zebib, A., 1982, The effect of significant viscosity variation on convective heat transport 

in water-saturated porous media: Journal of Fluid Mechanics, v. 117, p. 233-249.

Ghyben, W.B., 1888, Nota in verband met de voorgenomen putboring nabij Amsterdam, Tijdschrift van Let Koninklijk Inst. Van Ing.

Glover, R.E., 1959, The pattern of fresh-water flow in a coastal aquifer: Journal of Geophysical Research. v. 64, p. 457-460.

Guo, Weixing, and Bennett, G.D., 1998, Simulation of saline/fresh water flows using MODFLOW, in E. Poeter and others, 

MODFLOW ’98 Conference, Golden, Colorado, 1998, Proceedings: Golden, Colorado, v. 1, p. 267-274. 

Guo, Weixing, and Langevin, C.D., 2002, User’s guide to SEAWAT: A computer program for simulation of three-dimensional 

variable-density groundwater flow: U.S. Geological Survey Techniques of Water-Resources Investigations, book 6, chap. A7, 77 p. 

Guo, Zhaoli, and Zhao, T.S., 2005, Lattice Boltzmann simulation of natural convection with temperature-dependent viscosity in 

a porous cavity: Progress in Computational Fluid Dynamics v. 5, nos. 1/2, p. 110-117.

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, M.G., 2000, MODFLOW-2000, the U.S. Geological Survey modular 

ground-water model—User guide to modularization concepts and the ground-water flow process: U.S. Geological Survey 

Open-File Report 00-92, 121 p. 

Henry, H.R., 1964, Effects of dispersion on salt encroachment in coastal aquifers, in H.H Cooper, Sea Water in Coastal Aquifers: 

U.S. Geological Survey Water-Supply Paper 1613-C, p. C71-C84.

Henry, H.R., Hilleke, J.B., 1972, Exploration of multiphase fluid flow in a saline aquifer system affected by geothermal heating: 

University of Alabama Bureau of Engineering Research Contract No. 14-08-0001-12681, submitted to U.S. Geological Survey, 

Washington D.C., 105 p.

Herzberg, A., 1901, Die Wasserversorgung einiger nordseebader: J. Gasbeleucht, Wasserversorg, v. 44, p. 815-819.

Holzbecher, E.O., 1998, Modeling density-driven flow in porous media: Principles, numerics, software: New York, SpringerVerlag, 286 p.

Horne, R.N., 1975, Transient effects in geothermal convective systems, Phd. Thesis, University of Auckland.

Horton, C.W., and Rogers Jr., F.T., 1945, Convection currents in a porous medium: Journal of Applied Physics, v. 16, p. 360-370. 
References Cited 27

Hubbert, M.K., 1940, The theory of ground-water motion: Journal of Geology, v. 48, p. 785-944.

Hughes, J.D., and Sanford, W.E., 2004, SUTRA-MS, a version of SUTRA modified to simulate heat and multiple-solute transport: 

U.S. Geological Survey Open-File Report 2004-1207, 141 p.

Hughes, J.D., Sanford, W.E., and Vacher, H.L., 2005, Numerical simulation of double-diffusive finger convection: Water 

Resources Research, v. 41, 16 p.

Kipp, K.L., Jr., 1987, HST3D—A computer code for simulation of heat and solute transport in three-dimensional ground-water 

flow systems: U.S. Geological Survey Water-Resources Investigations Report 86-4095, 517 p.

Kipp, K.L., Jr., 1997, HST3D—Guide to the revised heat and solute transport simulator: HST3D—Version 2: U.S. Geological 

Survey Water-Resources Investigations Report 97-4157, 149 p.

Langevin, C.D., Thorne, D.T., Jr., Dausman, A.M., Sukop, M.C., and Guo, Weixing, 2007, SEAWAT Version 4: A computer 

program for simulation of multi-species solute and heat transport: U.S. Geological Survey Techniques and Methods, book 6, 

chap. A22, 39 p.

Langevin, C.D., Shoemaker, W.B., and Guo, Weixing, 2003, MODFLOW-2000, the U.S. Geological Survey modular ground-water 

model—Documentation of the SEAWAT-2000 version with the variable-density flow process (VDF) and the integrated MT3DMS 

transport process (IMT): U.S. Geological Survey Open-File Report 03-426, 43 p.

Lapwood, E.R., 1948, Convection of a fluid in a porous medium: Proceedings. Cambridge Phil. Society. A 225, p. 508-521.

Nield, D.A., and Bejan, A.,1999, Convection in porous media (2nd edition): New York, Springer-Verlag.

Nusselt, W. 1944. Technische Thermodynamik II (Theorie der Warmekraftmaschined), Sammlung Goschen, Bd. 1151, 144 p.

Pasa, G., and Titaud, O., 2005, A class of viscosity profiles for oil displacement in porous media or Hele-Shaw cell: Transport in 

Porous Media, v. 3, no. 3, p. 269-286.

Prasad, A., and Simmons, C.T., 2003, Unstable density-driven flow in heterogeneous porous media: A stochastic study of the 

Elder [1967b] “short heater” problem: Water Resources Research v. 39, No. 1.

Prasad, A., and Simmons, C.T., 2005, Using quantitative indicators to evaluate results from variable-density ground-water flow 

models: Hydrogeology Journal v. 13, p. 905-914.

Pringle, S.E., Glass, R.J., and Cooper, C.A., 2002, Double-diffusive finger convection in a Hele-Shaw cell—An experiment exploring 

the evolution of concentration fields, length scales and mass transfer: Transport in Porous Media, v. 47, no. 2, p. 195-214.

Rayleigh, Lord, 1916, On convection currents in a horizontal layer of fluid when the higher temperature is on the under side. 

Philos Mag., v. XXXII, p. 529-546.

Simpson, M.J., and Clement, T.P., 2004, Improving the worthiness of the Henry problem as a benchmark for density-dependent 

ground-water flow models: Water Resources Research v. 40, No. 1.

Sanford, W.E., and Konikow, L.F., 1985, A two-constituent solute-transport model for ground water having variable density: 

U.S. Geological Survey Water-Resources Investigations Report 85-4279, 88 p.

Strack, O.D.L., 1995, A Dupuit-Forchheimer model for three-dimensional flow with variable density: Water Resources Investigation, 

v. 12, p. 3007-3017.

Thorne, D., Langevin, C.D., and Sukop, M.C., 2006, Addition of simultaneous heat and solute transport and variable fluid viscosity 

to SEAWAT: Computer and Geosciences v. 32, p. 1758-1768.

Voss, C.I., 1984, A finite-element simulation model for saturated-unsaturated, fluid-density-dependent ground-water flow with 

energy transport or chemically-reactive single-species solute transport: U.S. Geological Survey Water-Resources Investigations 

Report 84-4369, 409 p.

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, v. 23, no. 10, p. 1851-1866.
28 Application of SEAWAT to Select Variable-Density and Viscosity Problems

Weatherhill, D., Simmons, C.T., Voss, C.E., and Robinson, N.I., 2004, Testing density-dependent ground-water models: twodimensional steady state unstable convection in infinite, finite and inclined porous layers: Advances in Water Resources, v. 27, 

p. 547-562.

Weber, J.E., 1975, The boundary-layer regime for convection in a vertical porous layer: International Journal of Heat and Mass 

Transfer, v. 18, April 1975, p. 569-573.

Zheng, C., 1990, MT3D: A modular three-dimensional transport model for simulation of advection, dispersion and chemical 

reactions of contaminants in ground-water systems:. Report to the U.S. Environmental Protection Agency, Ada, Oklahoma, 

USA.

Zheng, Chunmiao, 2006, MT3DMS v5.2 Supplemental user’s guide: Technical Report to the U.S. Army Engineer Research and 

Development Center, Department of Geological Sciences, University of Alabama, 24 p.

Zheng, Chunmiao, and Wang, P.P., 1999, MT3DMS—A modular three-dimensional multispecies transport model for simulation 

of advection, dispersion and chemical reactions of contaminants in ground-water systems; documentation and user’s guide: 

U.S. Army Corps of Engineers Contract Report SERDP-99-1. 
Appendix 29

Appendix

Benchmark problems from different systems created with different codes can be compared based on the calculation of 

dimensionless numbers such as the Nusselt or Rayleigh. Numbers, many of which are used to analyze flow, viscosity, and heat 

transfer through the liquid in a porous medium system are discussed herein because they are used to evaluate and compare the 

benchmark problems. Two systems that are described by the same dimensionless equations (in particular, the same values of the 

dimensionless numbers) and that have the same initial and boundary conditions (including geometric similarity) are mathematically equivalent (Bird and others, 2006). Because the equations for temperature and concentration effects on a flow regime are 

analogous (Thorne and others, 2006; Langevin and others, 2007), if dimensionless numbers and dimensionless initial and boundary 

conditions are equal, a comparison between the two systems can be made. All variables used in this report are defined in table 1.

Darcy Number

The Darcy number (Da) is a dimensionless number that is a measure of the validity of Darcy’s law (Nield and Bejan, 1999) 

where k is the permeability, and H is a characteristic length, such as the height of the system being modeled:

(A7)

The velocity calculated by Darcy’s law:

(A8)

is valid when assuming flow is through an area of porous media. Darcy’s law is not valid for general free-fluid or open channel flow.

The Darcy number is important when analyzing the effects of model boundaries on simulations of groundwater flow 

through a porous medium using Darcy’s law. A slip boundary condition is assumed when simulating flow based on Darcy’s law, 

basically meaning that flowing water is not slowed by the presence of a model boundary, but retains its momentum and slips 

along the boundary by ignoring friction at the model boundary. A particle flowing through a porous medium is impeded by that 

medium, therefore decreasing the momentum. The model boundary frictional effect on the momentum of flow is considered 

minimal (almost zero) in comparison to the porous medium effect on the reduction in momentum; therefore the model boundary effect is ignored (or considered a slip boundary). However, when simulating free fluid flow or a porous medium with high 

hydraulic conductivities, the frictional effects of the model boundary on momentum can no longer be ignored. Therefore, a no-slip 

boundary condition has to be implemented and Darcy’s law cannot be used to simulate flow in this environment.

According to Nield and Bejan (1999), Darcy’s law holds true for Darcy numbers less than 10-3. When Darcy numbers are 

greater than 10-3, the system is close to an open-channel or free fluid flow system; thus it violates the slip boundary condition 

assumption in normal applications of Darcy’s law. SEAWAT_V4 solves a variable-density flow form of Darcy’s law (Langevin 

and others, 2007); therefore all benchmark problems have to have Darcy numbers less than 10-3. 

This research used the Darcy number in the benchmark problem which tests SEAWAT_V4 against another modeling 

method, Lattice Boltzmann, which solves two-dimensional oil convection in aluminum foam (Guo and Zhao, 2005). The Lattice 

Boltzmann method is capable of solving free fluid flow; therefore, it can solve problems with Darcy numbers greater than 10-3. 

However, for our purpose in comparing the two flow systems in porous media, the Darcy numbers were limited to less than 10-3.

Rayleigh Number

The Rayleigh number is a dimensionless number which defines the dominance of convection as opposed to conduction of 

heat (or solute) transfer. High Rayleigh number systems are characterized by convective flow, and low Rayleigh number systems 

have stable density profiles and heat is transferred by conduction. The Critical Rayleigh number is the value where the change 

from convection to conduction of the heat through the water will occur. All else being equal, systems with equal Rayleigh numbers should have similar solute or heat transport patterns.
30 Application of SEAWAT to Select Variable-Density and Viscosity Problems

The equation used to calculate the Rayleigh number for free thermal convection (when porosity = 1) is:

(A9)

(Rayleigh, 1916; Holzbecher, 1998; Guo and Zhao, 2005). However, another Rayleigh number, herein referred to as Ra*, is often 

used for porous media. , where Ra* is the Darcy number ( ) times the Rayleigh number (equation A9), Ra*

is also called the Darcy-Rayleigh number (the Darcy number is discussed in the previous section). Essentially, the permeability is

included in the equation for Ra* to account for the resistance of fluid flow from the porous medium. Therefore, the equation is:

(A10)

(Lapwood, 1948; Nield and Bejan, 1999), where thermal diffusivity (a physical property that controls the rate at which heat is 

transferred by conduction) is calculated using the following equation (Holzbecher, 1998):

(A11)

The bulk thermal conductivity (kTbulk) is calculated using:

(A12)

Another version of the Rayleigh number for a porous medium uses the following equation (Prasad and Simmons, 2003; Weatherhill 

and others, 2004):

(A13)

Equation A13 uses porosity in the denominator which is different from Ra* in equation A10 without porosity. This is because 

thermal diffusivity is calculated using the following equation (Thorne and others, 2006):

(A14)

Porosity cancels out in the denominator for Ra*; therefore equations A10 and A13 are both correct and equal for a porous 

medium, depending on how thermal diffusivity is calculated.

Thermal Rayleigh numbers are analogous to solutal Rayleigh numbers because they are described by the same dimensionless 

equations. The coefficients used in calculation of the dimentionless Rayleigh numbers are analogous, such as:

• The change in temperature term, , in equations A9, A10 or A13 may be substituted with the change in concentration

term, , for solute transport problems.

• The thermal and solutal expansion coefficients have a similar overall effect on the Rayleigh number. The thermal expansion 

coefficient is: . is negative when simulating temperature, because as temperature increases, density 

decreases. When simulating concentrations and calculating the solutal expansion coefficient, is now , and as concentration increases, density increases (the sign is now positive). However, the sign cancels out with the value for in 

equations A9, A10, and A13 (or ).

• The molecular diffusion coefficient can be used instead of the thermal diffusion coefficient and it has a similar effect in 

the solute transport process.

θρ
Appendix 31

Therefore, the dimensionless Rayleigh number in systems with variable-density flow resulting from changes in temperature can 

be compared to systems with flow affected by changes in concentration.

The different types of Rayleigh numbers presented are used in some of the benchmark problems to compare the results 

from SEAWAT_V4 to previous research. The following benchmark problems use at least one of the different types of Rayleigh 

numbers discussed:

• Two-dimensional oil convection in aluminum foam (Guo and Zhao, 2005);

• Horton-Rogers-Lapwood (HRL) convection (Horton and Rogers, 1945; Lapwood, 1948); and

• Original Elder problem (Elder, 1967).

Prandtl Number

The Prandtl number is a ratio which aids in determining the relative thickness of the momentum boundary layer to the 

thermal boundary layer (Nield and Bejan, 1999; Guo and Zhao, 2005):

(A15)

According to Nield and Bejan, low Prandtl numbers (much less than one) indicate that hydrodynamic wall effects have little 

influence on heat transfer from the wall because the diffusive/thermal boundary layer is greater than the momentum (viscous) 

boundary layer.

The Prandtl number is used on the benchmark problem two-dimensional oil convection in aluminum foam. Guo and Zhao 

(2005) compute Prandtl and Darcy numbers, showing mathematically “that for a given Prandtl number, the non-linear drag force 

[along a model boundary] becomes negligibly small only for a small Da, or in the Darcy Regime”, and therefore Darcy’s law 

applies only for small Darcy numbers (as discussed previously). The Prandtl number is identical for both the model created with 

the Lattice Boltzmann method and the SEAWAT_V4 model so the two systems are comparable.

Nusselt Number

The Nusselt number is a dimensionless number which enables comparison of the actual heat (or solute) transfer across a 

boundary in a model layer to the heat (or solute) transfer through the entire system if heat transfer occurred solely through conduction 

(or diffusion in the case of solute transport). The Nusselt number, defined as

(A16)

is the ratio of the total heat flux through a layer to the steady state diffusive flux across the model (Nusselt, 1944; Prasad and 

Simmons, 2005). Here, represents a volumetric temperature gradient across a layer. Systems where all heat 

transfer is solely by conduction have a Nusselt number of one. Systems which are not stable with respect to their density profiles 

and heat transfer occurs by free convection in addition to heat conduction, would have Nusselt numbers greater than one. The 

end result is a dimensionless number representing the energy (or mass) transferred into the system. 

The Nusselt number changes with time and is typically analyzed by graphing the Nusselt number versus time for a 

system/model. The graphs for different systems are compared to one another to analyze the timing, type of heat/solute transfer 

(convection or conduction), and the amount of heat transfer. Benchmark problems with similar graphs of the Nusselt number 

with time are considered analogous. The Nusselt number is used to compare the Elder problem from SEAWAT_V4 to two other 

benchmark problems, the original Elder problem (Elder, 1967), and a more recent simulation by Prassad and Simmons (2003).

Pr

Dausman, A.M., and others—Application of SEAWAT to Select Variable-Density and Viscosity Problems—Scientific Investigations Report 2009-5028

Printed on recycled paper