U.S. Department of the Interior

U.S. Geological Survey

Techniques and Methods Book 6, Chapter A22

SEAWAT Version 4: A Computer Program for Simulation

of Multi-Species Solute and Heat Transport
Blank page
SEAWAT Version 4: A Computer Program 

for Simulation of Multi-Species Solute and 

Heat Transport

Christian D. Langevin1, Daniel T. Thorne, Jr.2, Alyssa M. Dausman1, Michael C. Sukop3, 

and Weixing Guo4

Techniques and Methods Book 6, Chapter A22

U.S. Department of the Interior

U.S. Geological Survey

1 Florida Integrated Science Center, U.S. Geological Survey, Fort Lauderdale, Florida.

2 Department of Mathematics, Physics, and Computer Science, Georgetown College, Georgetown, Kentucky.

3 Department of Earth Sciences, Florida International University, Miami, Florida.

4 Schlumberger Water Services, Fort Myers, Florida.
U.S. Department of the Interior

DIRK KEMPTHORNE, Secretary

U.S. Geological Survey

Mark D. Myers, Director

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

For product and ordering information:

World Wide Web: http://www.usgs.gov/pubprod

Telephone: 1-888-ASK-USGS

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

natural hazards, and the environment:

World Wide Web: http://www.usgs.gov

Telephone: 1-888-ASK-USGS

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:

Langevin, C.D., Thorne, D.T., Jr., Dausman, A.M., Sukop, M.C., and Guo, Weixing, 2008, SEAWAT Version 

4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport: U.S. Geological Survey 

Techniques and Methods Book 6, Chapter A22, 39 p.
iii

Preface

This report describes an enhanced version of the SEAWAT computer program that can be obtained on the 

Internet as follows: http://water.usgs.gov/software/ground_water.html/

The performance of the program has been tested in a variety of applications. Future applications, 

however, might reveal errors that were not detected in the test simulations. Users are requested to send 

notification of any errors found in this User Guide or the model program to the contact listed on the SEAWAT 

website. Updates might be made to both the report and to the model program. Users can check for updates at 

the Internet address shown above.

SEAWAT Version 4 is a replacement for SEAWAT-2000. SEAWAT-2000 users are encouraged to use 

this new version, even if the features documented in this report are not required for a particular application. SEAWAT Version 4 is backward compatible with datasets created for SEAWAT-2000 and, therefore, 

SEAWAT-2000 is no longer supported or maintained.
iv
Contents

Preface ...........................................................................................................................................................iii

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

Introduction.....................................................................................................................................................2

Mathematical Formulation ...........................................................................................................................3

Governing Equations for Flow and Transport...................................................................................3

Fluid Density...........................................................................................................................................4

Fluid Viscosity........................................................................................................................................6

Implementation of Fluid Viscosity Effects on Ground-Water Flow........................................................7

Use of MT3DMS to Simulate Heat Transport............................................................................................7

Thermal Equilibration between Fluid and Solid ...............................................................................8

Heat Conduction..................................................................................................................................10

Specified Heat Flux.............................................................................................................................11

Additional Simulation Options in SEAWAT_V4........................................................................................12

Flow-Field Update Control.................................................................................................................12

Options for Constant-Head Boundaries..........................................................................................12

Input Instructions and Evaluation of Temperature Output....................................................................15

MT3DMS Dispersion (DSP) Input File .............................................................................................15

Variable-Density Flow (VDF) Process Input File ............................................................................16

Viscosity (VSC) Package Input File..................................................................................................20

Time-Variant Constant-Head (CHD) Package Input File...............................................................22

Evaluating Temperature Output........................................................................................................23

Example Problem .........................................................................................................................................23

Case 1....................................................................................................................................................25

Case 2....................................................................................................................................................27

Case 3....................................................................................................................................................28

Case 4....................................................................................................................................................30

Case 5....................................................................................................................................................32

Case 6....................................................................................................................................................34

Case 7....................................................................................................................................................36

References Cited..........................................................................................................................................38

Figures

1. Diagrams showing density options used at constant-head boundaries to calculate 

boundary reference head..........................................................................................................14

2-8. Graphs showing transient movement of the:

2. Freshwater-seawater transition zone for Case 1 ................................................................... 26

3. Freshwater-seawater transition zone for Case 2 ................................................................... 28

4. Salinity and temperature transition zones for Case 3 ........................................................... 29

5. Salinity and temperature transition zones for Case 4 ........................................................... 31

6. Salinity and temperature transition zones for Case 5 ...........................................................33

7. Salinity and temperature transition zones for Case 6 ........................................................... 35

8. Salinity and temperature transition zones for Case 7 ........................................................... 37 

v
Tables

1. Values used in equations that relate dynamic viscosity to temperature in degrees Celsius.............6

2. Values of specific heat capacity and density for selected substances and materials........................9

3. Values of thermal conductivity for selected substances and materials ...............................................10

4. Description of the options available for the Time-Variant Constant Head (CHD) Package. .............23

5. Input values for the example problem ..........................................................................................................24

6. Description of the cases simulated for the example problem.................................................................25

Conversion Factors

Multiply By To obtain

meter (m) 3.281 foot (ft) 

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

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

square meter per day (m2/d) 10.7636 square foot per day (ft2/d)

square meter per second (m2/s) 10.7636 square foot per second (ft2/s)

Celsius (°C) °F = (1.8 × °C) + 32 Fahrenheit (°F)

Acronyms and Abbreviations

BAS6 Basic PCG2 Pre-Conditioned Conjugate Gradient

BCF Block Centered Flow RCT Reactions

CHD Time Variation Constant Head SSM Source and Sink Mixing

DSP Dispersion TDS Total Dissolved Solids

GCG Generalized Conjugate Gradient VDF Variable Density Flow

GWF Ground Water Flow VSC Viscosity

HUF2 Hydrogeologic Unit Flow g/L gram per liter

IMT Integrated MT3DMS Transport J/kgoK Joule per kilogram degrees Kelvin

LPF Layer Property Flow

vi
SEAWAT Version 4: A Computer Program for 

Simulation of Multi-Species Solute and 

Heat Transport

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

and Weixing Guo4

Abstract

The SEAWAT program is a coupled version of MODFLOW and MT3DMS designed to simulate threedimensional, variable-density, saturated ground-water flow. Flexible equations were added to the program to allow 

fluid density to be calculated as a function of one or more MT3DMS species. Fluid density may also be calculated 

as a function of fluid pressure. The effect of fluid viscosity variations on ground-water flow was included as an 

option. Fluid viscosity can be calculated as a function of one or more MT3DMS species, and the program includes 

additional functions for representing the dependence on temperature. Although MT3DMS and SEAWAT are not 

explicitly designed to simulate heat transport, temperature can be simulated as one of the species by entering 

appropriate transport coefficients. For example, the process of heat conduction is mathematically analogous to 

Fickian diffusion. Heat conduction can be represented in SEAWAT by assigning a thermal diffusivity for the 

temperature species (instead of a molecular diffusion coefficient for a solute species). Heat exchange with the 

solid matrix can be treated in a similar manner by using the mathematically equivalent process of solute sorption. 

By combining flexible equations for fluid density and viscosity with multi-species transport, SEAWAT Version 

4 represents variable-density ground-water flow coupled with multi-species solute and heat transport. SEAWAT 

Version 4 is based on MODFLOW-2000 and MT3DMS and retains all of the functionality of SEAWAT-2000.

SEAWAT Version 4 also supports new simulation options for coupling flow and transport, and for representing 

constant-head boundaries. In previous versions of SEAWAT, the flow equation was solved for every transport 

timestep, regardless of whether or not there was a large change in fluid density. A new option was implemented 

in SEAWAT Version 4 that allows users to control how often the flow field is updated. New options were also 

implemented for representing constant-head boundaries with the Time-Variant Constant-Head (CHD) Package. 

These options allow for increased flexibility when using CHD flow boundaries with the zero-dispersive flux solute 

boundaries implemented by MT3DMS at constant-head cells.

This report contains revised input instructions for the MT3DMS Dispersion (DSP) Package, Variable-Density 

Flow (VDF) Package, Viscosity (VSC) Package, and CHD Package. The report concludes with seven cases of an 

example problem designed to highlight many of the new features.

1 Florida Integrated Science Center, U.S. Geological Survey, Fort Lauderdale, Florida.

2 Department of Mathematics, Physics, and Computer Science, Georgetown College, Georgetown, Kentucky.

3 Department of Earth Sciences, Florida International University, Miami, Florida.

4 Schlumberger Water Services, Fort Myers, Florida.
2 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Introduction

SEAWAT is a coupled version of MODFLOW (Harbaugh and others, 2000) and MT3DMS (Zheng and Wang, 

1999; Zheng, 2006) designed to simulate three-dimensional, variable-density ground-water flow and multi-species 

transport. The Variable-Density Flow (VDF) Process in SEAWAT is based on the constant-density Ground-Water 

Flow (GWF) Process of MODFLOW-2000. The VDF Process uses the familiar and well-established MODFLOW 

methodology to solve the variable-density ground-water flow equation (Langevin and others, 2003). The MT3DMS 

part of SEAWAT, referred to as the Integrated MT3DMS Transport (IMT) Process, solves the solute transport 

equation.

Previous versions of SEAWAT (Guo and Bennett, 1998; Guo and Langevin, 2002; Langevin and others, 2003) 

were developed for isothermal conditions, and thus, did not solve simultaneous solute and heat transport with 

the combined effects of concentration and temperature on variable-density flow. This restriction was due to the 

implementation of the equation of state for fluid density, which could be a function of only solute concentration. 

Furthermore, previous versions of SEAWAT did not represent the effects of fluid viscosity variations, which are most 

commonly caused by temperature variations. The updated version of SEAWAT described herein, called SEAWAT 

Version 4, or SEAWAT_V4, retains all of the capabilities of SEAWAT-2000, but is also able to simulate solute and 

heat transport simultaneously, as well as the effect that fluid viscosity variations have on resistance to ground-water 

flow.

Heat transport and solute transport contain many similarities (Anderson, 2005). Their mathematical 

representation is nearly identical when the terms describing heat transport are formulated in equivalent solute 

expressions. SEAWAT_V4 leverages these similarities by using MT3DMS to simulate heat transport. With this 

approach, temperature is one of the simulated species. Thorne and others (2006) provide a description of the 

approach and show how to assign appropriate thermal parameters for the temperature species.

Although SEAWAT_V4 was designed primarily to accommodate simultaneous solute and heat transport, the 

features that allow for this capability have been implemented in a general fashion that may be useful for other 

applications where temperature variations are not important. Other features have also been added to improve overall 

functionality. The new features are:

•	 The equation of state for fluid density can be a function of one or more MT3DMS species and, optionally, a 

pressure term.

•	 The effect of fluid viscosity on the internodal conductance terms can be included in a simulation by using the 

Viscosity (VSC) Package and the Layer Property Flow (LPF) Package. The equation for viscosity can be a 

function of one or more MT3DMS species. There is also an option to designate the temperature species, and 

to specify one of several non-linear equations to represent the dependence of viscosity on temperature.

•	 An option has been implemented to solve the VDF equation to update the flow field only if fluid densities 

have changed by some user-specified value. This option may substantially reduce runtimes for simulations in 

which a new flow solution is not necessarily required for each transport timestep.

•	 Auxiliary variables can now be used with the Time-Variant Constant-Head (CHD) Package to designate the 

fluid density value associated with the prescribed head. This option allows for more control of boundary heads 

specified with the CHD Package.
Mathematical Formulation 3

This user guide provides a brief description of new options available in SEAWAT_V4, describes the required 

input data files, and demonstrates the options with a simple example problem. This user guide describes only these 

new features, and does not provide all of the background material necessary for designing successful SEAWAT 

simulations. Background material is available in the published user guides for SEAWAT (Guo and Langevin, 2002), 

SEAWAT-2000 (Langevin and others, 2003), MODFLOW-2000 (Harbaugh and others, 2000), and MT3DMS (Zheng 

and Wang, 1999; Zheng, 2006).

Mathematical Formulation

This section presents the governing equations and the equations of state used to represent fluid density and 

viscosity as a function of state variables. The equations presented here do not explicitly allow for dissolution and 

precipitation of minerals that may likely result from large temperature changes and the heat lost or gained from those 

reactions. The formulation presented here is for water in the liquid phase only and does not represent the transition to 

ice or vapor.

Governing Equations for Flow and Transport

The VDF Process solves the following form of the variable-density ground-water flow equation (tensors and 

vectors shown in bold):

s sq s t

C

t C

h h z S ' 0

,0

0

0

0

0 U U U T

U

U U

P

P U  w

w

w

w



w

w » 

1⁄4

o «

¬

a

 ̧

 ̧

1

·

 ̈

 ̈

©

§    K0   , (1)

where 

ρ0 is the fluid density [ML-3] at the reference concentration and reference temperature; 

μ is dynamic viscosity [ML-1T-1];

K0 is the hydraulic conductivity tensor of material saturated with the reference fluid [LT-1];

h0 is the hydraulic head [L] measured in terms of the reference fluid of a specified concentration and 

temperature (as the reference fluid is commonly freshwater). Note that h0 was referred to as hf in previous 

SEAWAT publications.

Ss,0 is the specific storage [L-1], defined as the volume of water released from storage per unit volume per unit 

decline of h0;

t is time [T];

θ is porosity [-];

C is salt concentration [ML-3]; and 

q's is a source or sink [T-1] of fluid with density ρs.

The IMT Process solves the following form of the solute transport equation:

     k

s s

k k

k k

b d C C q C

t

K C 1        ' w

w

 ̧

 ̧

1

·

 ̈

 ̈

©

§

 TD q

T

T

U

, (2)

where 
4 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

ρb is the bulk density (mass of the solids divided by the total volume) [ML-3], 

Kd

k is the distribution coefficient of species k [L3M-1], 

Ck is the concentration of species k [ML-3], 

D is the hydrodynamic dispersion coefficient tensor [L2T-1], 

q is specific discharge [LT-1], and

Cs

k is the source or sink concentration [ML-3] of species k.

In previous versions of SEAWAT, P

P0

 was approximated as one, and thus viscosity effects were neglected, 

and fluid density was treated as a simple linear function of only one solute species. In SEAWAT_V4, fluid density 

and viscosity can be calculated using concentrations from one or more solute species. These equations of state relate 

the density and viscosity terms in equation 1 to one or more of the MT3DMS species concentrations (Ck) in equation 2. As will be discussed later, equation 2 can also be used to represent heat transport. In this case, one of the 

MT3DMS solute species is used to represent temperature, and the mathematical terms describing heat transport are 

converted into forms analogous to those used in equation 2. By allowing one or more of the MT3DMS species to 

affect fluid density and viscosity, variable-density ground-water flow can be coupled with simultaneous solute and 

heat transport.

Fluid Density

The equation of state, which relates fluid density to solute concentration, temperature, and pressure, is typically 

written as a function of the volumetric expansion coefficient (Nield and Bejan, 1999). As presented in Diersch and 

Kolditz (2002), the volumetric expansion coefficients for solute concentration (βC), temperature (βT), and pressure 

(βP) are:

T P

C C ,

1  ̧

1

·  ̈

©

§

w

w U

U

E 

,

 

(3)

C P

T T ,

1  ̧

1

·  ̈

©

§

w

w U

U

E , and (4)

C T

P P ,

1  ̧

1

·  ̈

©

§

w

w U

U

E ,

 

(5)

where T is temperature [oK] and P is pressure [ML-1T-2]. Each expansion coefficient is expressed for one dependent 

variable (that is concentration, temperature or pressure) while holding the others fixed. If it also can be assumed that 

each volumetric expansion coefficient is constant over the range of expected values, a general equation of state for 

fluid density can be written as:

    = − + −+ − 0 00 0 exp  C TP (CC TT PP ) ( ) ( )  

.

 (6)

A linearized form of equation 6 is typically implemented in numerical flow and transport codes, such as SUTRA 

(Voss, 1984) and HST3D (Kipp, 1987), and is obtained using the first term of a series expansion (Bear, 1972). The 

linearized form is written as:
Mathematical Formulation 5

      U U0  U0EC C C0  U0ET T  T0  U0EP P  P0 , (7)

or by expanding equation 7 around the reference density:

      0 0 0 P P0 P

T T

T

C C

C  w

w  

w

w  

w

w  U U U U U , (8)

where:

C C U E U 0 w

w

, (9)

T T U E U 0 w

w

, and (10)

P P U E U 0 w

w

. (11)

For the transport of multiple solute species, where the concentration of each species can affect fluid density, Hughes 

and Sanford (2004) implemented the following equation of state (with the exception of the pressure term) in 

SUTRA-MS:

      0 0

1

0 0 P P

P

T T

T

C C

C

NS

k

k k

k  w

w  

w

w  

w

w  ¦ 

U U U U U , (12)

where NS is the number of solute species that affect fluid density. A difference between SEAWAT and SUTRA-MS 

is that SEAWAT is formulated in terms of a reference head, whereas SUTRA-MS is formulated in terms of pressure. 

For this reason, the third term on the right-hand side of equation 12 is reformulated here using the height of a water 

column l [L] of density, ρ0. This variable, which can be thought of as the pressure head in terms of the reference 

density, is related to pressure by:

P = l ρ0 g, (13)

where l = h0 - z. With some simple rearrangement, the final form of the equation of state used and implemented in 

the program is:

      0 0

1

0 0 A A

A  w

w  

w

w  

w

w ¦ 

U U U U U T T

T

C C

C

NS

k

k k

k , (14)

where wA

wU can be calculated from the volumetric expansion coefficient for pressure, βP, using:

U gE P

U 2 0 ww

A . (15)
6 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Fluid Viscosity

Dynamic viscosity is considered here to be a function of only temperature and solute concentration, which is 

the typical approach (Diersch and Kolditz, 2002), and neglects the weak dependence of viscosity on fluid pressure. 

This approach is also used in other ground-water modeling codes such as HST3D (Kipp, 1987) and SUTRA-MS 

(Hughes and Sanford, 2004). A general equation for representing fluid viscosity as a function of concentration and 

temperature is: 

    0

1

0 0 T T

T

C C

C

NS

k

k k

k  w

w  

w

w ¦ 

P P P P . (16)

This linear equation is one of the options available in SEAWAT_V4. 

Over many temperature ranges, the linear approximation does not adequately represent the effect of temperature 

on dynamic viscosity. For this reason, an alternative equation for dynamic viscosity has been implemented:

  ¦   

 w

w 

NS

k

k k

T k C C

C

T

1

0

P P P . (17)

When using equation 17, users have several options for representing dynamic viscosity as a function of temperature. 

These options include:

   ̧

 ̧

1

·

 ̈

 ̈

©

§

  4

3

1 2

T A

A

PT T A A , (18)

  >  @ 5

1 2 3 4

A

PT T A  A  A T  A , and (19)

  2

1

A

PT T A T . (20)

Values for the A constants are specified by the user. These three equations are based on viscosity and temperature 

relations encountered in the literature. For example, SUTRA (Voss, 1984) uses the following viscosity and temperature relation:

   ̧

1

·  ̈

©

§

   133.15

248.37

7 239.4 10 10 T P T x 

,

 (21)

where T is expressed in degrees Celsius. To use this same relation in SEAWAT_V4, the user would select equation 

18 and specify values of 239.4 x 10-7, 10, 248.37, and 133.15 for A1, A2, A3, and A4, respectively. Values for A can be 

found in table 1.

Table 1. Values used in equations that relate dynamic viscosity to temperature in degrees Celsius.

[--, not applicable for this equation]

A1 A2 A3 A4 A5

Equation

(in report) Reference

2.394 x 10-5 10 248.37 133.15 -- 18 Voss (1984)

1 x 10-3 1 1.5512 x 10-2 -20 -1.572 19 Pawlowski (1991)

0.168 -1.0868 -- -- -- 20 Guo and Zhao (2005)*

*Relation is for oil viscosity as a function of temperature (between 5 and 170 degrees Celsius). Relation is not applicable to water.
Use of MT3DMS to Simulate Heat Transport 7

Implementation of Fluid Viscosity Effects on Ground-Water Flow

The variable-density form of Darcy’s Law used in previous versions of SEAWAT was modified to include the 

effects of viscosity variations. The equations for including the viscosity terms in the finite-difference representation 

are given in Thorne and others (2006). The following are some general considerations for including viscosity effects 

in a SEAWAT simulation:

•	 The effects of viscosity variations are activated by using the Viscosity (VSC) Package;

•	 The effects of these viscosity variations are implemented only for the Layer-Property Flow (LPF) Package; the 

Block-Centered Flow (BCF) and Hydrogeologic Unit Flow (HUF2) Packages cannot be used with the VSC 

Package;

•	 Flows to or from boundaries (for example, river or general-head boundaries) are not adjusted for the effects of 

viscosity variations; and

•	 The VSC Package was designed specifically for the VDF Process, which means that the “VDF” file type must 

be included in the name file.

The LPF Package was modified to include the effects of viscosity variations for each of the different horizontal 

conductance weighting options. The LPF Package includes three weighting options for calculating the horizontal 

internodal conductances: harmonic, logarithmic, and a mixed method that uses the arithmetic mean for saturated 

thickness and the logarithmic mean for hydraulic conductivity. The LPF package supports only harmonic weighting 

for internodal conductances in the vertical direction. 

Including the effects of viscosity, variations on ground-water flow required two programming modifications. 

First, the hydraulic conductivity of each cell is multiplied by P

P0 (where μ represents the cell viscosity) before the 

internodal conductance is calculated using one of the weighting options. Second, the LPF Package was modified 

to ensure that the internodal conductances are recalculated each time the flow equation is solved; in a standard 

MODFLOW simulation, internodal conductances for a confined layer are calculated once at the start of a simulation 

and remain constant for the entire simulation. Thorne and others (2006) verified the viscosity implementation 

through a match with results from an equivalent analytical solution.

Use of MT3DMS to Simulate Heat Transport

The standard stand-alone version of MT3DMS can be used in its present form to approximate heat transport 

in porous media. For example, Martin and others (2001) used MT3DMS to simulate heat transport on Grand 

Cayman Island. A key assumption in using the standard stand-alone version of MT3DMS to represent heat 

transport is that temperature does not affect fluid density and viscosity, and thus, it is not necessary to couple 

flow and transport. Another limitation prior to MT3DMS version 5.2 was that simultaneous solute and heat 

transport during the same simulation could not be performed because only a single molecular diffusion coefficient 

(representative of all species) could be entered. Similarly, previous versions of SEAWAT could be used to 

approximate heat transport, provided that: (1) temperature was the only species included in the simulation; and 

(2) the effect of temperature variation on fluid density was the only feedback on ground-water flow (that is 

viscosity did not affect ground-water flow).
8 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

In SEAWAT_V4, both heat and multi-species solute transport can be performed in the same simulation, and 

as shown in the previous section, concentration and temperature affect ground-water flow through the density and 

viscosity terms. Simulation of simultaneous solute and heat transport is performed by simply treating temperature 

as a species in MT3DMS. The ability to enter different molecular diffusion coefficients for each species was 

implemented in MT3DMS version 5.2. This MT3DMS version has been implemented in SEAWAT_V4, thereby 

eliminating an important limitation to representing simultaneous solute and heat transport.

Analogous forms of the solute and heat transport equations highlight similarities between the two processes. 

These similarities are used to calculate the correct input values required by MT3DMS when representing heat 

transport. The following is a general form of the solute transport equation solved by MT3DMS:

    k

s s

k k k

m

k k

b d D C C q C

t

K C 1     ' »

1⁄4

o «

¬

a  ̧ 1

·  ̈

©

§   

w

w

 ̧

 ̧

1

·

 ̈

 ̈

©

§

 q q Į

T

T T

T

U

,

 

(22)

where Dm

k

 is the molecular diffusion coefficient [L2T-1] for species k, and α is the dispersivity tensor [L]. The 

following is a form of the heat transport equation, which was manipulated by Thorne and others (2006) to highlight 

the similarity with the solute transport equation:

    s s

Pfluid

Tbulk

Pfluid

s Psolid T T q T

c

k

t

T

c

c ' 1 1  

»

»

1⁄4

o

«

«

¬

a  ̧

 ̧

1

·

 ̈

 ̈

©

§  

w

w

 ̧

 ̧

1

·

 ̈

 ̈

©

§   q q Į

TU T

T T

U

U

T

T

, (23)

where 

ρs is the density of the solid (mass of the solid divided by the volume of the solid) [ML-3],

cPsolid is the specific heat capacity of the solid [L2 T-2 oK-1],

cPfluid is the specific heat capacity of the fluid [L2 T-2 oK-1],

kTbulk is the bulk thermal conductivity of the aquifer material [ML3 T-2 oK-1], and 

Ts is the source temperature [oK]. 

Note that ρb, ρs, and θ are related by: ρb = ρs (1-θ ).

Specific relations are required when using MT3DMS to represent thermal equilibration between the fluid 

and solid, heat conduction, and a specified heat flux. There are no special considerations for representing advective 

heat transport (convection). 

Thermal Equilibration between Fluid and Solid

The analogy between solute and heat transport is evident in equations 22 and 23. In equation 22, the factor 

preceding the time derivative on the left-hand side of the equation represents solute sorption. A similar term is 

shown in equation 23 and represents an equilibration of temperature between the fluid and solid. For solute and heat 

transport, this term acts as a “retardation coefficient” and causes a concentration or temperature front to move slower 

than the average linear flow velocity. To represent thermal equilibrium between the fluid and solid using MT3DMS, 

the thermal distribution factor, Kd_temp [M-1L3], can be calculated for the temperature species using the following 

equation:

Pfluid

Psolid

tempd c

c K

U _ .

 (24)
Use of MT3DMS to Simulate Heat Transport 9

This expression contains the fluid density, which means that the user must assume a fluid density value prior to 

running the simulation. This is a limitation with the present implementation in SEAWAT_V4. It is unlikely, however, 

that model results will be highly sensitive to errors in this term, so use of a representative fluid density in the calculation of Kd_temp should provide satisfactory results for most simulations. Another limitation is that the specific heat 

capacity of water, and perhaps other substances, may not be constant over large temperature ranges. For example, 

Holzbecher (1998) reports that although the specific heat capacity for water remains constant at a value of about 

4,200 J/kg oK for temperatures less than 100 oC, the value increases by about 40 percent over the range from 100 oC 

to 300 oC. Table 2 provides a list of specific heat capacity values for a variety of substances.

The linear sorption isotherm should be selected to represent the thermal exchange between the water and 

the solid. By using the linear isotherm, the thermal exchange between the fluid and solid is constant for different 

temperatures (that is, the exchange is not a function of temperature). For solute sorption (which is used here to 

represent thermal exchange), other isotherms are sometimes required to express the relation between the sorbed and 

dissolved solute.

Table 2. Values of specific heat capacity (cP ) and density for selected 

substances and materials.

[All values are at 25 °C unless otherwise noted. The rock densities listed here 

represent ρs . J/kgoK, joules per kilogram degrees Kelvin; kg/m3, kilogram per cubic 

meter]

Substance cP (J/kg°K)* ρ (kg/m3)

Freshwater (liquid) 4.186 x 103 1.00 x 103

Freshwater (ice) 2.108 x 103 (at 0 °C) 0.92 x 103

Seawater 3.93 x 103 (at 2.2 °C) 1.025 x 103

Air 1.005 x 103 1.205 (at 20 °C)

Granite 7.9 x 102 **2.75 x 103

Marble 8.8 x 102 **2.80 x 103

Sandstone 9.2 x 102 **2.80 x 103

Limestone 8.4 x 102 **2.70 x 103

Dolostone 9.2 x 102 **2.90 x 103

Clay 9.2 x 102 *1.8 x 103 - 2.6 x 103

Quartz 7.1 x 102 (at 0 °C) **2.64 x 103

Calcite 8.4 x 102 (average over 0–100 °C) **2.71 x 103

*Engineering Tool Box (2007).

**Cote and Konrad (2005).
10 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Heat Conduction

The solute transport process of molecular diffusion is mathematically similar to the heat transport process 

of thermal conduction. In both processes, transport occurs in response to a gradient (Fick’s law for molecular 

diffusion and Fourier’s law for heat conduction). The relation between equations 22 and 23 indicates that the 

molecular diffusion coefficient entered for the temperature species (Dm_temp) is:

Pfluid

Tbulk

m temp c

k D

TU _ .

 (25)

This thermal conduction term, Dm_temp [L2T-1], is referred to here as the bulk thermal diffusivity. Hughes and 

Sanford (2004) present two equations for calculating the bulk thermal conductivity from the solid and fluid thermal 

conductivities (kTfluid and kTsolid , respectively). A simple arithmetic mean that is weighted by volume can be calculated using:

  Tbulk Tfluid Tsolid k Tk  1T k . (26)

Evidence indicates that the geometric mean may also be appropriate for calculating the bulk thermal conductivity 

(Ingebritsen and Sanford, 1998). The bulk thermal conductivity can be calculated using the following geometric 

mean equation:

T  T  1

Tbulk Tfluid Tsolid k k k . (27)

Whether using the arithmetic mean or the geometric mean, the value for Dm_temp entered into MT3DMS (as a 

molecular diffusion coefficient) to represent heat conduction could be calculated using equation 25 and either 

equation 26 or 27. Table 3 provides a list of thermal conductivity values for a variety of substances. Note that many 

of the thermal conductivity values in table 3 are for dry rock; bulk thermal conductivity values for saturated rock 

may be higher.

Table 3. Values of thermal conductivity (kT )for selected substances and materials. 

[All values are at 25 degrees Celsius, except for freshwater (ice), which is at 0 degrees Celsius]

Substance kT (W/[m oK]) Reference

Freshwater (liquid) 0.58 Engineering Tool Box (2007)

Freshwater (ice) 2.18 Engineering Tool Box (2007)

Air 0.024 Engineering Tool Box (2007)

Granite 1.73 - 3.98 Marble Institute of America (2007)

Quartzite 5.38 Marble Institute of America (2007)

Marble 2.07 - 2.94 Marble Institute of America (2007)

Sandstone 1.83 - 2.90 Marble Institute of America (2007)

Limestone 1.26 - 2.15 Marble Institute of America (2007)

Dolostone 3.8 Cote and Konrad (2005)

Clay 1 Ingebritsen and Sanford (1998)

Quartz 3 Engineering Tool Box (2007)

Calcite 3.59 Horai (1971)
Use of MT3DMS to Simulate Heat Transport 11

In MT3DMS versions prior to 5.2, the molecular diffusion coefficient was entered as a single value for each 

model layer. When using molecular diffusion to represent heat conduction, it may be necessary to enter variations 

of Dm within a model layer to characterize heterogeneity patterns. In version 5.2, MT3DMS was modified to allow 

Dm to be entered as a three-dimensional array. Thus, heat conduction can be represented using a heterogeneous 

distribution for thermal conductivity, porosity, and fluid density. Although it is possible that the thermal conductivity 

can be anisotropic—for example, where strong grain alignments characterize the aquifer fabric—this is not expected 

to be an important factor in most cases and only isotropic thermal conductivities can be entered into the model.

There are several limitations with using the present implementation of the molecular diffusion process in 

MT3DMS to represent heat conduction:

•	 The approach requires an isotropic representation of thermal conductivity;

•	 A representative fluid density must be used to calculate Dm_temp , which is fixed for the entire simulation 

regardless of whether or not the fluid density changes; and

•	 Over large temperature ranges, thermal conductivity (and therefore Dm_temp) is a function of temperature 

(Holzbecher, 1998).

These limitations are not considered to be overly restrictive; however, users are encouraged to evaluate their effects. 

For example, to quantify the limitation of having a fluid density term in equations 24 and 25, the user could recalculate the three-dimensional Kd and Dm input arrays using fluid densities calculated from a previous simulation.

Specified Heat Flux

A heat flux expressed as energy per unit time, E [ML 2T-3], may be specified for any model cell by using 

the “mass-loading source” option (ITYPE = 15) in the MT3DMS Source and Sink Mixing (SSM) Package. The 

appropriate “solute equivalent,” denoted here as MT  [oKL3T-1] to enter in the SSM Package can be calculated using 

the following relation:

Pfluid

T

c

E M

U

  

. 

 (28)

This calculated input value for MT  can be entered on a cell-by-cell basis in the SSM input file. 
12 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Additional Simulation Options in SEAWAT_V4

Two new simulation options have been added in SEAWAT_V4. The first option allows the user to control 

whether or not the flow field is recalculated based on a user-specified tolerance for the maximum change in fluid 

density at any cell. The second option allows for more flexibility in specifying heads with the CHD Package. 

Flow-Field Update Control

In previous versions of SEAWAT, the flow and transport equations are solved for each transport timestep. For 

most problems that require the use of a variable-density code, such as SEAWAT, this is a conservative and reasonable 

approach. For some problems, however, it may not be necessary to solve the flow equation during each transport 

timestep. This is the case for standard applications of MT3DMS, where it is assumed that the dilute solute has no 

affect on fluid density. 

To reduce computer runtimes, a new option was added to SEAWAT_V4 that allows the user to control how 

often the flow field is updated. Because this option cannot be used with the iterative coupling option in SEAWAT 

(NSWTCPL is assigned a value greater than 1), the flow-field update control is activated by assigning a -1 to 

NSWTCPL in the VDF Package input file. The user also specifies a value for the maximum change in fluid density 

(DNSCRIT). If the maximum change in fluid density (between the present transport timestep and the last transport 

timestep for which the flow field was updated) is less than DNSCRIT, then SEAWAT_V4 does not update the flow 

field. Once the maximum change in fluid density exceeds DNSCRIT, then SEAWAT_V4 updates the flow field by 

solving the flow equation. Sanford and Konikow (1985) implemented a similar option in the MOCDENSE computer 

code. Under the present implementation, SEAWAT_V4 will always solve for the flow field:

•	 If the density change since the last flow field update exceeds DNSCRIT at one or more cells,

•	 For the first transport timestep of the flow timestep, or

•	 For the last transport timestep of a flow timestep.

Users can evaluate the loss of accuracy that will result from using the flow-field update control. For a seawater 

intrusion problem, where the fluid density may vary between freshwater (1,000 kg/m3) and seawater (1,025 kg/m3), 

users may want to assess the compromise between shorter computer runtimes and loss of accuracy for different 

percentages of the density range. To update the flow field if the density changes by more than 10 percent of the 

freshwater-seawater range, users would specify a value of 2.5 kg/m3 for DNSCRIT.

As a cautionary note, parameter estimation programs may require extremely precise sensitivities for accurate 

identification of an optimal parameter set. Because the flow-field update control can reduce accuracy, in perhaps a 

non-linear manner, use of the flow-field update control is not recommended when precise sensitivities are required 

by a parameter estimation program.

Options for Constant-Head Boundaries

Assignment and proper use of constant-head boundaries in MODFLOW is relatively straightforward. There are 

some widely recognized problems in that a constant-head boundary can provide or remove unrealistic quantities of 

water, but in general, measured heads from the field can be assigned directly to constant-head boundaries without 

problem. For variable-density systems, assignment of heads to model boundaries is less straightforward. The 

complication arises from the dependence of heads on fluid density. This complication can be avoided, provided the 

user recognizes that head values (whether observed in the field or simulated by a model) have an associated density 

to which that head is referenced.
Additional Simulation Options in SEAWAT_V4 13

SEAWAT-2000 reads and writes heads in terms of the density at the point where the head is measured (or 

simulated). These head values are called “point-water heads,” “actual heads,” “heads in terms of the native aquifer 

water,” or just, “heads.” After the point-water heads are read, SEAWAT converts them into heads in terms of the 

reference fluid density (ρ0). Because ρ0 is normally set to the density of freshwater, these reference heads are 

commonly referred to as equivalent freshwater heads. For heat and solute transport simulations or multi-species 

simulations where fluid density is dependent upon more than one species, the term “equivalent freshwater head” 

does not provide enough information and requires an additional qualifier, such as “at the reference temperature.” The 

term “reference head,” therefore, more clearly states the intended meaning. SEAWAT solves the flow equation in 

terms of the reference head; however, SEAWAT writes point-water heads to the output files. These point-water heads 

are referenced to the simulated density at that location and time.

The dependence of head on fluid density has important implications for assigning boundary heads, particularly 

when the density of the boundary head changes during the simulation. This commonly occurs when using the 

MT3DMS constant-head boundary condition (ITYPE = 1 in the SSM Package) to represent the ocean or other 

coastal water body. When this solute boundary type is applied, the concentration at the cell is calculated as part of 

the simulation, subject to the condition that there is no dispersive flux across the cell boundary. Additionally, any 

fluid entering the aquifer from the boundary cell is assigned a concentration equal to that specified by the user in the 

SSM input file. When a constant-concentration boundary condition is used (ITYPE = -1), the following issues do not 

apply. 

As described in the SEAWAT-2000 documentation (Langevin and others, 2003), simulation results may vary 

slightly depending on whether the Basic (BAS6) or the CHD Package is used to assign constant-head boundaries. 

If the BAS6 Package is used to specify constant-head cells, the initial head is converted to a reference head using the 

initial fluid density array, and this reference head value remains constant during the simulation. If the CHD Package 

is used, the reference head value assigned to the boundary cell is updated prior to each transport timestep using the 

fluid density from the previous transport timestep. Thus, if the density changes in a constant-head cell, SEAWAT 

may produce different results depending on the method used to simulate the constant-head cell.

Previous implementations of the BAS6 and CHD methods for simulating constant-head cells in SEAWAT did 

not provide a robust approach for assigning boundary head values. Consider the model section shown in figure 1A 

where a constant-head boundary is assigned to the last column of cells to represent the ocean. The solute boundary 

assigned in the SSM input file (using ITYPE = 1) specifies that any flow into the model from the boundary is 

assigned seawater salinities. When flow is reversed, however, from the aquifer into the boundary, the salinity 

value calculated at the boundary cell may not necessarily equal the salinity value specified in the SSM input file. 

If seawater salinities were assigned everywhere as the initial condition for this simulation, the BAS6 approach 

would provide a reasonable representation of the head boundary; the reference head value at the boundary would be 

calculated using seawater densities (the initial densities), and the reference head value would be held constant for the 

duration of the simulation. 

The BAS6 approach has two potential limitations: (1) the constant head (ocean level) cannot be changed during 

the simulation, and (2) this condition may not be appropriate for estuaries or boundaries in which the discharge can 

change the salinity of the boundary feature. In such cases, it is necessary to reference the specified point-water head 

to the calculated density at the boundary cell, rather than the initial density. Both of these limitations are avoided 

by using the CHD Package to simulate constant-head cells. With the CHD Package, the boundary head can change 

during the simulation, and the head value used for each transport timestep is referenced to the calculated density 

value at that cell. Although this method was appropriate for shallow estuaries or other features where ground-water 

discharge could change the surface water salinity, it was not appropriate for deep ocean boundaries, where it can 

reasonably be assumed that ocean salinity is not affected by fresh ground-water discharge. 
14 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

bnd

P

 g

P

h

5z

g d  z 

Freshwater

Seawater

Transition zone

OPTION

0 (default) 1 2

User-specified head













Datum

A

B

CONSTANT-HEAD OCEAN BOUNDARY

z z z

P

5  g

X

Z

EXPLANATION

bnd

P



g

h

is pressure

is density

is gravity

is head

z is elevation

Seawater

Freshwater

Transition zone

is density of the boundary

Figure 1. Density options used at constant-head boundaries to calculate boundary reference 

head: (A) a typical model section where the density at the boundary cells is calculated based on 

the flow system, and (B) different representations of the ocean boundary. 
Input Instructions and Evaluation of Temperature Output 15

To increase flexibility in assigning boundary heads, three additional options were added for the CHD Package in 

SEAWAT_V4. As described later, the options are activated using auxiliary variables in the CHD Package. Figure 1B 

shows three of the four options available for representing boundary heads with the CHD Package. The options differ 

only in the assumed vertical density profile that is used to convert the assigned boundary head to the reference head 

(referenced to ρ0) used in the solution to the flow equation. Option 0 is the default option, and was first implemented 

in SEAWAT-2000. With Option 0, the input head represents the point-water head in terms of the calculated density at 

the cell. With Option 1, the user also enters (within the CHD Package) the density (ρbnd) for which the input head is 

referenced. This option is appropriate for most types of ocean boundaries, where ρbnd would be assigned equal to the 

density of seawater. For Option 2, the density value used in the conversion of the input head to reference head is the 

average density of the overlying active cells. This option may be useful for inland settings with vertically stratified 

density distributions. With Option 2, the input head value can be conceptualized as environmental head (Lusczynski, 

1961). Environmental heads are unique in that there is no vertical flow if the vertical gradient of the environmental 

head is zero. Thus, if the same input head value is assigned to all of the constant-head cells in a vertical column, this 

boundary type will result in hydrostatic conditions at that column. With Option 3, the final option, the user enters the 

reference head (freshwater head) directly and there is no adjustment in the program for density. 

An important consequence of using either Option 1 or 2 is that the output heads written by SEAWAT_V4 may 

not equal the constant-head values that were entered as part of the CHD Package. This outcome is produced because 

SEAWAT writes point-water heads in terms of the calculated densities. Thus, extreme caution is necessary when 

using output heads from one model run as input to another when the CHD Package is used with either Option 1 or 2.

Input Instructions and Evaluation of Temperature Output

Multi-diffusive solute transport, whereby a unique diffusion coefficient can be entered for each MT3DMS 

species, is activated through the input file for the MT3DMS Dispersion (DSP) Package. The effect of multiple 

species on fluid density is invoked through the input file for the VDF Process. The effect of fluid viscosity on 

ground-water flow is activated through a new package in SEAWAT_V4, called the Viscosity (VSC) Package. The 

VSC Package works only with the LPF Package and is not compatible with the BCF or HUF2 Packages. Moreover, 

the VSC Package is compatible only with the VDF Process, and will not function for constant-density simulations 

(the VDF option is not included in the name file). Details on the input structure for the DSP, VDF, VSC, and CHD 

input files and an evaluation of simulated temperatures are presented herein.

MT3DMS Dispersion (DSP) Input File

In MT3DMS version 5.2, users can enter different diffusion coefficients for each species. With this option, 

the diffusion coefficient is entered as a three-dimensional array with one value for each model cell. In addition, 

users can enter molecular diffusion coefficients for solute species and a thermal conductivity parameter (eq. 25) for 

the temperature species. The new multi-diffusion option is described by Zheng (2006) and consists of entering the 

keyword “$MultiDiffusion” (case insensitive and without quotation marks) to the first non-commented line of the 

DSP input file. If this keyword is detected, MT3DMS (and SEAWAT_V4) uses the RARRAY utility (described in 

Zheng and Wang, 1999) to read a two-dimensional diffusion coefficient array for each layer. This process is repeated 

for each mobile solute component.
16 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Variable-Density Flow (VDF) Process Input File

Input for the VDF Process is read from the file that is specified in the name file with “VDF” as the file type. All 

single-valued variables are free format if the option “FREE” is specified in the BAS6 Package input file; otherwise, 

the variables have 10-character fields. New input variables or variables that have new meaning in SEAWAT_V4 are 

underlined.

FOR EACH SIMULATION

0. [#Text]

Item 0 is optional and the symbol # must be in column 1. Item 0 can be repeated multiple times, but it cannot be used 

in the file after the first data record.

1. MT3DRHOFLG MFNADVFD NSWTCPL IWTABLE

2. DENSEMIN DENSEMAX

Read Item 3 if NSWTCPL > 1 or NSWTCPL = -1

3. DNSCRIT 

Read Item 4 if MT3DRHOFLG ≥ 0

4. DENSEREF DRHODC(1)

Read Item 4a through 4c if MT3DRHOFLG = -1

4a. DENSEREF DRHODPRHD PRHDREF

4b. NSRHOEOS

4c. MTRHOSPEC(NSRHOEOS) DRHODC(NSRHOEOS) CRHOREF(NSRHOEOS)

Item 4c is read as a table, with one row representing each MT3DMS species to be included in the equation 

of state for fluid density. There is no requirement that this table be sorted by MTRHOSPEC.

5. FIRSTDT

FOR EACH STRESS PERIOD

Read Items 6 and 7 only if MT3DRHOFLG = 0

6. INDENSE

Read Item 7 only if INDENSE is greater than zero

7. [DENSE(NCOL,NROW)] – U2DREL

Item 7 is read for each layer in the grid.
Input Instructions and Evaluation of Temperature Output 17

Explanation of Variables Read from the VDF Process Input File

Text—is a character variable (199 characters) that starts in column 2. Any characters can be included in Text. The 

“#” character must be in column 1. Except for the name file, lines beginning with # are restricted to these first lines 

of the file. Text is printed when the file is read.

MT3DRHOFLG—is the MT3DMS species number that will be used in the equation of state to compute fluid 

density. This input variable was formerly referred to as MTDNCONC (Langevin and others, 2003).

If MT3DRHOFLG = 0, fluid density is specified using items 6 and 7, and flow will be uncoupled with 

transport if the IMT Process is active.

If MT3DRHOFLG > 0, fluid density is calculated using the MT3DMS species number that corresponds 

with MT3DRHOFLG. A value for MT3DRHOFLG greater than zero indicates that flow will be coupled 

with transport.

If MT3DRHOFLG = -1, fluid density is calculated using one or more MT3DMS species. Items 4a, 4b, and 

4c will be read instead of item 4.

The dependence of fluid density on pressure head can only be activated when MT3DRHOFLG = -1. A value for 

MT3DRHOFLG of -1 indicates that flow will be coupled with transport.

MFNADVFD—is a flag that determines the method for calculating the internodal density values used to conserve 

fluid mass.

If MFNADVFD = 2, internodal conductance values used to conserve fluid mass are calculated using a 

central-in-space algorithm.

If MFNADVFD ≠ 2, internodal conductance values used to conserve fluid mass are calculated using an 

upstream-weighted algorithm.

NSWTCPL—is a flag used to determine the flow and transport coupling procedure. 

If NSWTCPL = 0 or 1, flow and transport will be explicitly coupled using a one-timestep lag. The explicit 

coupling option is normally much faster than the iterative option and is recommended for most applications.

If NSWTCPL > 1, NSWTCPL is the maximum number of non-linear coupling iterations for the flow and 

transport solutions. SEAWAT-2000 will stop execution after NSWTCPL iterations if convergence between 

flow and transport has not occurred.

If NSWTCPL = -1, the flow solution will be recalculated only for:

The first transport step of the simulation, or

The last transport step of the MODFLOW timestep, or

The maximum density change at a cell is greater than DNSCRIT.

IWTABLE—is a flag used to activate the variable-density water-table corrections (Guo and Langevin, 2002, eq. 82).

If IWTABLE = 0, the water-table correction will not be applied.

If IWTABLE > 0, the water-table correction will be applied.

DENSEMIN—is the minimum fluid density. If the resulting density value calculated with the equation of state is less 

than DENSEMIN, the density value is set to DENSEMIN.
18 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

If DENSEMIN = 0, the computed fluid density is not limited by DENSEMIN (this is the option to use for 

most simulations).

If DENSEMIN > 0, a computed fluid density less than DENSEMIN is automatically reset to DENSEMIN.

DENSEMAX—is the maximum fluid density. If the resulting density value calculated with the equation of state is 

greater than DENSEMAX, the density value is set to DENSEMAX. 

If DENSEMAX = 0, the computed fluid density is not limited by DENSEMAX (this is the option to use for 

most simulations).

If DENSEMAX > 0, a computed fluid density larger than DENSEMAX is automatically reset to 

DENSEMAX.

DNSCRIT—is a user-specified density value. If NSWTCPL is greater than 1, DNSCRIT is the convergence criterion, in units of fluid density, for convergence between flow and transport. If the maximum fluid density difference 

between two consecutive implicit coupling iterations is not less than DNSCRIT, the program will continue to iterate 

on the flow and transport equations, or will terminate if NSWTCPL is reached. If NSWTCPL is -1, DNSCRIT is 

the maximum density threshold, in units of fluid density. If the fluid density change (between the present transport 

timestep and the last flow solution) at one or more cells is greater than DNSCRIT, then SEAWAT_V4 will update the 

flow field (by solving the flow equation with the updated density field).

DENSEREF—is the fluid density at the reference concentration, temperature, and pressure. For most simulations, 

DENSEREF is specified as the density of freshwater at 25 °C and at a reference pressure of zero.

DRHODC—formerly referred to as DENSESLP (Langevin and others, 2003), is the slope of the linear equation of 

state that relates fluid density to solute concentration. In SEAWAT_V4, separate values for DRHODC can be entered 

for as many MT3DMS species as desired. If DRHODC is not specified for a species, then that species does not affect 

fluid density. Any measurement unit can be used for solute concentration, provided DENSEREF and DRHODC are 

set properly. DRHODC can be approximated by the user by dividing the density difference over the range of endmember fluids by the difference in concentration between the end-member fluids. 

For many applications, the end-member fluids are freshwater and seawater. Therefore, if the measurement 

units for both density and solute concentration are either kilograms per cubic meter or pounds per cubic feet, the 

value for DRHODC is about 0.7143. Users may also prefer to use relative concentrations, such that freshwater has a 

concentration of zero, and seawater has a concentration of 1.0. If relative concentrations are used with density values 

in units of kilograms per cubic meter, the appropriate DRHODC value is about 25.

SEAWAT_V4 offers the flexibility to use a variety of different solute concentrations to calculate fluid density, 

provided DRHODC is specified correctly. For example, users can utilize MT3DMS to simulate the concentration 

of total dissolved solids (TDS), chloride concentration, or the relative seawater concentration on a scale from 0 

to 1. Use of TDS as the solute species for the fluid density calculation is straightforward, and is described in Guo 

and Langevin (2002, p. 18). If chloride is used for the fluid density calculation, a constant ratio between chloride 

concentration and TDS normally has to be assumed. For seawater and mixtures of seawater with freshwater, the ratio 

between chloride concentration and TDS is about 55 percent (19 g/L of chloride per 35 g/L of TDS). In some cases, 

particularly at relatively low TDS concentrations, this constant ratio may no longer hold. In these cases, caution 

is necessary in defining the relation between concentration and fluid density. If the linear relation implemented 

in SEAWAT_V4 is not appropriate, users may need to add a new equation of state and recompile the program. 

Regardless of the approach, simple hand calculations of fluid density are suggested as a way of ensuring that the 

correct input values are being used.
Input Instructions and Evaluation of Temperature Output 19

DRHODPRHD—is the slope of the linear equation of state that relates fluid density to the height of the pressure 

head (in terms of the reference density). Note that DRHODPRHD can be calculated from the volumetric expansion coefficient for pressure using equation 15. If the simulation is formulated in terms of kilograms and meters, 

DRHODPRHD has an approximate value of 4.46 × 10-3 kg/m4. A value of zero, which is typically used for most 

problems, inactivates the dependence of fluid density on pressure.

PRHDREF—is the reference pressure head. This value should normally be set to zero.

NSRHOEOS—is the number of MT3DMS species to be used in the equation of state for fluid density. This value is 

read only if MT3DRHOFLG = -1.

MTRHOSPEC—is the MT3DMS species number corresponding to the adjacent DRHODC and CRHOREF.

CRHOREF—is the reference concentration (C0) for species, MTRHOSPEC. For most simulations, CRHOREF 

should be specified as zero. If MT3DRHOFLG > 0, CRHOREF is assumed to equal zero (as was done in previous 

versions of SEAWAT).

FIRSTDT—is the length of the first transport timestep used to start the simulation if both of the following two conditions are met:

1. The IMT Process is active, and 

2. Transport timesteps are calculated as a function of the user-specified Courant number (the MT3DMS 

 input variable, PERCEL, is greater than zero)

INDENSE—is a flag. INDENSE is read only if MT3DRHOFLG is equal to zero. 

If INDENSE < 0, values for the DENSE array will be reused from the previous stress period. If it is the first 

stress period, values for the DENSE array will be set to DENSEREF.

If INDENSE = 0, values for the DENSE array will be set to DENSEREF.

If INDENSE ≥ 1, values for the DENSE array will be read from item 7.

If INDENSE = 2, values read for the DENSE array are assumed to represent solute concentration, and will 

be converted to density values using the equation of state.

DENSE—is the fluid density array read for each layer using the MODFLOW-2000 U2DREL array reader. The 

DENSE array is read only if MT3DRHOFLG is equal to zero. The DENSE array may also be entered in terms of 

solute concentration, or any other units, if INDENSE is set to 2 and the constants used in the density equation of 

state are specified appropriately.
20 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Viscosity (VSC) Package Input File

Input for the VSC Package is read from the file that is specified in the name file with “VSC” as the file type. 

The VDF Process must be active in order to use the VSC Package. In addition, the VSC Package works only with 

the LPF Package and is not compatible with the BCF or HUF2 Packages. All variables in the VSC input file are read 

using free format.

FOR EACH SIMULATION

0. [#Text]

Item 0 is optional and the symbol # must be in column 1. Item 0 can be repeated multiple times. 

1. MT3DMUFLG

2. VISCMIN VISCMAX

Read item 3 if MT3DMUFLG ≥ 0

3. VISCREF DMUDC(1) CMUREF(1)

Read item 3a through 3c if MT3DMUFLG = -1

3a. VISCREF

3b. NSMUEOS MUTEMPOPT

3c. MTMUSPEC(NSMUEOS) DMUDC(NSMUEOS) CMUREF(NSMUEOS)

Item 3c is read NSMUEOS times, once for each MT3DMS species to be included in the viscosity equation.

3d. MTMUTEMPSPEC AMUCOEFF(MUNCOEFF)

Item 3d is read if MUTEMPOPT > 0. The size of the AMUCOEFF coefficient array depends on the selected 

equation for viscosity.

FOR EACH STRESS PERIOD

Read items 4 and 5 only if MT3DMUFLG = 0

4. INVISC

Read item 5 only if INVISC is greater than zero

5. [VISC(NCOL,NROW)] – U2DREL

Item 5 is read for each layer in the grid.

Explanation of Variables Read from the VSC Package Input File

Text—is a character variable (199 characters) that starts in column 2. Any characters can be included in Text. The 

“#” character must be in column 1. Except for the name file, lines beginning with # are restricted to these first lines 

of the file. Text is printed when the file is read.

MT3DMUFLG—is the MT3DMS species number that will be used in the equation to compute fluid viscosity.
Input Instructions and Evaluation of Temperature Output 21

If MT3DMUFLG ≥ 0, fluid density is calculated using the MT3DMS species number that corresponds with 

MT3DMUFLG.

If MT3DMUFLG = -1, fluid viscosity is calculated using one or more MT3DMS species.

VISCMIN—is the minimum fluid viscosity. If the resulting viscosity value calculated with the equation is less than 

VISCMIN, the viscosity value is set to VISCMIN.

If VISCMIN = 0, the computed fluid viscosity is not limited by VISCMIN (this is the option to use for most 

simulations).

If VISCMIN > 0, a computed fluid viscosity less than VISCMIN is automatically reset to VISCMIN.

VISCMAX—is the maximum fluid viscosity. If the resulting viscosity value calculated with the equation is greater 

than VISCMAX, the viscosity value is set to VISCMAX. 

If VISCMAX = 0, the computed fluid viscosity is not limited by VISCMAX (this is the option to use for 

most simulations).

If VISCMAX > 0, a computed fluid viscosity larger than VISCMAX is automatically reset to VISCMAX.

VISCREF—is the fluid viscosity at the reference concentration and reference temperature. For most simulations, 

VISCREF is specified as the viscosity of freshwater.

DMUDC—is the slope of the linear equation that relates fluid viscosity to solute concentration.

NSMUEOS—is the number of MT3DMS species to be used in the linear equation for fluid viscosity (this 

number does not include the temperature species if the nonlinear option is being used). This value is read only 

if MT3DMUFLG = -1. A value of zero indicates that none of the MT3DMS species have a linear effect on fluid 

viscosity (the nonlinear temperature dependence may still be activated); nothing should be entered for item 3c in this 

case.

MUTEMPOPT—is a flag that specifies the option for including the effect of temperature on fluid viscosity.

If MUTEMPOPT = 0, the effect of temperature on fluid viscosity is not included or is a simple linear 

relation that is specified in item 3c.

If MUTEMPOPT = 1, fluid viscosity is calculated using equation 18. The size of the AMUCOEFF array in 

item 3e is 4 (MUNCOEFF = 4).

If MUTEMPOPT = 2, fluid viscosity is calculated using equation 19. The size of the AMUCOEFF array in 

item 3e is 5 (MUNCOEFF = 5).

If MUTEMPOPT = 3, fluid viscosity is calculated using equation 20. The size of the AMUCOEFF array in 

item 3e is 2 (MUNCOEFF = 2).

If NSMUEOS and MUTEMPOPT are both set to zero, all fluid viscosities are set to VISCREF.

MTMUSPEC—is the MT3DMS species number corresponding to the adjacent DMUDC and CMUREF.

CMUREF—is the reference concentration for species, MTMUSPEC.

DMUDC—is the slope of the linear equation that relates fluid viscosity to solute concentration.

MTMUTEMPSPEC—is the MT3DMS species number that corresponds to temperature. This value must be between 

1 and NCOMP and should not be listed in MTMUSPEC of item 3c.
22 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

AMUCOEFF—is the coefficient array of size MUNCOEFF. AMUCOEFF is “A” in equations 18, 19, and 20.

MUNCOEFF—is the size of the AMUCOEFF array.

INVISC—is a flag. INVISC is read only if MT3DMUFLG is equal to zero. 

If INVISC < 0, values for the VISC array will be reused from the previous stress period. If it is the first 

stress period, values for the VISC array will be set to VISCREF.

If INVISC = 0, values for the VISC array will be set to VISCREF.

If INVISC ≥ 1, values for the VISC array will be read from item 5.

If INVISC = 2, values read for the VISC array are assumed to represent solute concentration, and will be 

converted to viscosity values.

VISC—is the fluid viscosity array read for each layer using the MODFLOW-2000 U2DREL array reader. The VISC 

array is read only if MT3DMUFLG is equal to zero. The VISC array may also be entered in terms of solute concentration (or any other units) if INVISC is set to 2, and the simple linear expression in item 3 can be used to represent 

the relation to viscosity.

Time-Variant Constant-Head (CHD) Package Input File

The new options for specifying the density value assumed to correspond to the head values entered in the 

CHD Package are activated using auxiliary variables (Harbaugh and others, 2000). If the auxiliary variable 

“CHDDENSOPT” (in all uppercase) is detected in the CHD input file, then the program will read the CHD density 

option as an additional field. If all the values assigned to CHDDENSOPT are either 0 or 2, then no additional 

auxiliary variable is required. If, however, one or more CHD cells are assigned a value of 1 for CHDDENSOPT, then 

SEAWAT_V4 requires a user-specified density value for those cells. This density value must be entered through an 

additional auxiliary variable with the name “CHDDENS” (again in all uppercase). All three CHDDENSOPT options 

may be used within a single simulation. If option 1 is used for one or more CHD cells, the CHDDENS value will 

be required for all the CHD cells. The CHDDENS value will only be used, however, for those CHD cells with a 

CHDDENSOPT of 1. Table 4 provides a list of these options.

An example of a simple CHD input file is shown below:

 3 AUX CHDDENSOPT AUX CHDDEN

 3 0

 1 1 1 0.0000000 0.0000000 0 1025.000

 2 1 1 0.0000000 0.0000000 1 1025.000

 3 1 1 0.0000000 0.0000000 2 1025.000

This CHD input file will result in the assignment of CHD cells in layers 1, 2, and 3 of row 1 and column 1. CHD 

option 0 will be used for the CHD cell in layer 1. This option means that the assigned head value will be converted 

to reference head using the simulated density at that cell. CHD option 1 will be used for the CHD cell in layer 2. 

This option means that the assigned head value will be converted to a reference head using a density value of 1,025 

kg/m3 (entered in auxiliary variable “CHDDEN”). CHD option 2, used for the CHD cell in layer 3, indicates that 

the assigned head value is an environmental head that will be converted to a reference head using an average density 

value of the overlying CHD cells. 
Example Problem 23

Evaluating Temperature Output

Evaluating simulated temperatures is straightforward, and can be done using the same types of post-processing 

programs that are used to evaluate concentrations. Calculation of thermal energies, however, can be confusing 

because the “mass” values reported in the listing and mass budget files do not directly correspond with temperature 

or heat energy. The “mass” values reported by SEAWAT_V4 are calculated by multiplying temperature by a 

volumetric flux and the length of the transport timestep. Thus, for the temperature species, a reported mass value

might have units of m3 oC (or 

m C

d

3

° for a mass flux). These “mass” values (MT) can be converted into thermal energy

(E), using the following simple relation: E = MT ρcPfluid .

Example Problem

An example problem was developed to illustrate the new model features presented in this report. The problem 

consists of a two-dimensional cross section of a confined coastal aquifer initially saturated with relatively cold 

seawater at a temperature of 5 oC. Warmer freshwater with a temperature of 25 oC is injected into the coastal aquifer 

along the left boundary to represent flow from inland areas (fig. 1). The warmer freshwater flows to the right, where 

it discharges into a vertical ocean boundary. The ocean boundary is represented with hydrostatic conditions based 

on a fluid density calculated from seawater salinities at 5 oC. No-flow conditions are assigned to the top and bottom 

boundaries. A complete list of the input values used for the problem is given in table 5. This problem is a simplified 

representation of what might occur in a coastal carbonate platform.

Table 4. Description of the options available for the Time-Variant

Constant Head (CHD) Package.

CHDDENSOPT CHDDEN Description

0 Not used

Default option. Input head is converted to 

reference head (h0) using calculated fluid 

density value at the cell.

1 Required

Input head is converted to reference head

(h0) using the density value specified in the 

CHDDEN auxiliary variable.

2 Not used

Input head is an environmental head and 

is converted to a reference head (h0) using 

the average of the calculated fluid density 

values of the overlying active cells.

3 Not used Input head is the reference head (h0). 

No conversion is made.
24 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Table 5. Input values for the example problem.

[m, meter; m/d, meter per day; J/(kg oC), joules per kilogram degrees Celsius; kg/m3, kilogram per cubic meter; W/(m oC), watts per meter degrees Celsius; 

m2/d, square meter per day; m3/kg, cubic meter per kilogram; kg/(m3 oC), kilogram per cubic meter degrees Celsius; kg/m4, kilogram per meter to the 4th power]

Input parameter Value Units Comments

Grid properties

Number of columns (NCOL) 200 – Assigned

Number of rows (NROW) 1 – Assigned

Number of layers (NLAY) 100 – Assigned

Δx (DELR) 10 m Assigned

Δy (DELC) 1 m Assigned

Δz (DZ) 10 m Assigned

Aquifer properties

Horizontal hydraulic conductivity (Kh0) 10 m/d

Vertical hydraulic conductivity (Kv0) 0.1 m/d Assigned in terms of the reference fluid (at reference 

concentration and reference temperature)

Specific storage (Ss) 1.00 × 10-5 m

Porosity (θ) 0.35 – Assigned

Longitudinal dispersivity (αL) 1 m Assigned

Transverse dispersivity (αT) 0.1 m Assigned

Heat capacity of the solid (CPsolid) 835 J/(kg oC) Approximate value for calcite

Density of the solid (ρs) 2,710 kg/m3 Approximate value for calcite

Bulk density (ρb) 1,761.5 kg/m3 Calculated from θ and ρs, ρb= ρs(1- θ)

Thermal conductivity of solid (kTsolid) 3.59 W/(m oC) Approximate value for calcite

Bulk thermal conductivity (kTbulk) 2.547 W/(m oC) Volumetric weighted average (35% water, 65% calcite)

Diffusion coefficient (Dm_salinity) 1.00 × 10-10 m2/d Assigned

Bulk thermal diffusivity (Dm_temp) 0.150309621 m2/d Calculated using eq. 25

Fluid properties

Heat capacity of the fluid (CPfluid) 4,183 J/(kg oC) Approximate value for water

Thermal conductivity of water (kTfluid) 0.61 W/(m oC) Approximate value for water

Distribution coefficient for salinity (Kd_salinity) 0 m3/kg Assigned

Distribution coefficient for temperature (Kd_temp) 2.00 × 10-4 m3

/kg Calculated using eq. 24

Reference density (ρ0) 1,000 kg/m3 Assigned

δρ/δC 0.7 – Approximate change in density over change in salinity

Reference concentration for density (Cρ0) 0 kg/m3 Assigned

δρ/δT -0.375 kg/(m3

 oC) Approximate change in density over change in temperature

Reference temperature (T0) 25 oC Assigned

Reference viscosity 0.001 kg/(m s) Assigned

δμ/δC 1.92 × 10-6 m2

/d Approximate change in viscosity over change in salinity

Reference concentration for viscosity (Cμ0) 0 kg/m3 Assigned

Constants used for viscosity-temperature relation Based on eq. 18

A1 2.394 × 10-5 – From Hughes and Sanford (2004)

A2 10 – From Hughes and Sanford (2004)

A3 248.37 – From Hughes and Sanford (2004)

A4 133.15 – From Hughes and Sanford (2004)

δρ/δl 4.46 × 10-3 kg/m4 Calculated from water compressibility and eq. 15

Boundary conditions

Qinflow 10 m2

/d Assigned

Tinflow 25 oC Assigned

Cinflow 0 kg/m3 Assigned

hocean 0 m Head (or stage) of the ocean—not the freshwater head

Cocean 35 kg/m3 Assigned

Tocean 5 oC Assigned
Example Problem 25

To demonstrate the process of developing a numerical simulation of this problem, seven different cases are 

presented. For all of the cases, salinity is represented using species 1, and temperature is represented using species 2. 

For each case, a new process or feature is implemented in the model. A description of the seven cases is given in 

table 6.

Flow and transport were solved using identical solution procedures in all seven cases. Flow was solved using 

the Pre-Conditioned Conjugate Gradient (PCG2) Package. The head and flow convergence parameters were set at 

1 × 10-5 m and 1 kg/d (mass per time), respectively. The advective term of the solute transport equation was solved 

using the Total-Variation Diminishing (TVD) option; the remaining terms (dispersion, source/sink mixing, reaction) 

were solved using the iterative Generalized Conjugate Gradient (GCG) solver with SSOR pre-conditioner. Each case 

was run in transient mode for 200,000 days (about 550 years). A Courant number of 1.0 was used to calculate the 

length of the transport timesteps. A length of 0.01 day was specified for the first transport timestep. 

Table 6. Description of the cases simulated for the example problem.

Case

Salinity affects Temperature affects Heat

conduction

Thermal 

equilibration 

with solid

Boundary 

heat 

conduction

Pressure 

affects 

Density Viscosity Density Viscosity density

1 x

2 x x

3 x x x

4 x x x x

5 x x x x x

6 x x x x x x x

7 x x x x x x x x

Case 1

Case 1 is the simplest case and can be simulated with previous versions of SEAWAT. It is a variable-density 

simulation in which the fluid density is a function only of salinity. Although temperature is represented as species 2, 

it is represented as a conservative tracer in which thermal equilibration with the solid as well as heat conduction are 

not represented. 

The name file for case 1 is as follows:

#SEAWAT LISTING FILE

LIST 7 coastal_flow.lst

#

#

#SEAWAT VDF INPUT FILE

VDF 80 coastal_flow.vdf

#

#

#MODFLOW INPUT FILES

BAS6 1 coastal_flow.bas

DIS 29 coastal_flow.dis

LPF 11 coastal_flow.lpf

WEL 12 coastal_flow.wel

OC 22 coastal_flow.oc

PCG 19 coastal_flow.pcg

#

#
26 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

#MT3DMS INPUT FILES

BTN 81 coastal_flow.btn

ADV 82 coastal_flow.adv

DSP 83 coastal_flow.dsp

SSM 84 coastal_flow.ssm

GCG 86 coastal_flow.gcg

#

#

#OUTPUT FILES

DATA(BINARY) 50 coastal_flow.cbb

DATA(BINARY) 30 coastal_flow.hds

All of the files listed in the name file are either output files created during the simulation, or standard MODFLOW or 

MT3DMS input files. The only exception to this is the VDF input file, which is specific to SEAWAT. The contents of 

the VDF input file (coastal_flow.vdf) are as follows:

 1 1 0 0 MT3DRHOFLG MFNADVFD NSWTCPL IWTABLE

 0.0000 0.0000 DENSEMIN DENSEMAX

 1000.0000 0.7000 DENSEREF DRHODC

1.0000e-02 FIRSTDT

Three important values in this input file are MT3DRHOFLG = 1, DENSEREF = 1,000.0 kg/m3, and DRHODC = 

0.7. An MT3DRHOFLG value of 1 indicates that fluid density will be calculated using species 1 (salinity) in the 

equation of state. DENSEREF is the reference density (ρ0

), and DRHODC is the slope of the linear relation between 

density and salinity (δρ/δC). With these values, density will vary linearly between 1,000.0 kg/m3

 for freshwater and 

1,024.5 kg/m3 for seawater.

The transient movement of the transition zone between freshwater and seawater is shown in figure 2 for five 

different simulation times. Temperature is not shown because the results are identical to the salinity results. Model 

results show that the transition zone, which is initially vertical, tilts toward the upper-right corner during the 

simulation due to the effects of fluid density. At the end of the 200,000-day simulation, the interface slopes gently 

upward toward the outflow boundary.

Figure 2. Transient movement of the freshwater-seawater transition zone for Case 1. Gray shading represents 

salinity values between 1 and 99 percent of seawater. The black lines represent a salinity value equal to 

50-percent seawater.
Example Problem 27

Case 2

The effect of temperature on fluid density is activated for Case 2 (table 6). A linear expression is used to relate 

fluid density to temperature. The slope of that linear relation ( wT

wU ) is about -0.375 kg

m C3n . The negative sign indicates 

the fluid density increases as temperature decreases. As described earlier, more than one MT3DMS species can 

affect fluid density by entering -1 for MT3DRHOFLG in the VDF input file:

 -1 1 0 0 MT3DRHOFLG MFNADVFD NSWTCPL IWTABLE

 0.0000 0.0000 DENSEMIN DENSEMAX

 1000.0000 0.0000 0. DENSEREF DRHODPRHD PRHDREF

 2 NSRHOEOS

 1 0.7000 0. MTRHOSPEC(1) DRHODC(1) CRHOREF(1)

 2 -0.3750 25. MTRHOSPEC(2) DRHODC(2) CRHOREF(2)

1.0000e-02 FIRSTDT

By entering a -1 for MT3DRHOFLG and a value of 2 for NSRHOEOS, the contributions of species 1 and 2 to the 

density equation of state are both included (as shown in lines 5 and 6 of the VDF input file).

Results from the Case 2 simulation are shown in figure 3. As before, temperature results are not shown because 

they correspond directly with salinities. The 50-percent seawater line for the five different times from Case 1 are also 

shown in figure 3 (as dashed lines) to illustrate the effect of including the temperature dependence on fluid density. 

By adding the temperature effect on fluid density, the transition zone moves more quickly at the top of the model 

than at the center or bottom of the model. The final shape of the transition zone is also sloped more steeply than for 

Case 1.

The Case 2 simulation could also have been performed using only a single species. Because no transport 

processes operate differently on salinity than on temperature, it is possible to reformulate the density equation of 

state to include the effects of both temperature and salinity. For example, temperature is a simple linear function of 

salinity for Case 2, where:

 	 4 20

7

inflow ocean

inflow ocean

inflow ocean

T T

T CT T C

C C

¥  ́ 
   
 
  ¦ μ § ¶ 
 . (29)

By substituting this equation into the equation of state for fluid density, the resulting DENSEREF and DRHODC 

values would be 1,000.0 kg/m3 and 0.914, respectively. If these values had been used for the Case 1 simulation, 

the results would be identical to the results shown for Case 2. This is true only because salinity and temperature 

are directly related. Heat conduction is included in the next case, and as a consequence, temperature is not directly 

related to salinity.
28 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Case 3

In Case 3, heat conduction is included in the simulation by specifying a value for bulk thermal diffusivity 

(Dm_temp) within the DSP input file. A bulk thermal diffusivity value of 0.15 m2/d (1.74 × 10-6 m2/s) was calculated 

using equation 25; therefore, the following DSP input file is used for Case 3:

$MultiDiffusion

 0 1.000(10E12.4) -1 DISPERSIVITY Layer 1

 0 1.000(10E12.4) -1 DISPERSIVITY Layer 2

(layers 3 through 99 not shown)

 0 1.000(10E12.4) -1 DISPERSIVITY Layer 100

 0 0.100(10E12.4) -1 HZ DISPERSIVITY RATIO

 0 0.100(10E12.4) -1 VT DISPERSIVITY RATIO

 0 1.0E-10 -1 MOL DIFF SP 1 LAYER 001

 0 1.0E-10 -1 MOL DIFF SP 1 LAYER 002

(layers 3 through 99 not shown)

 0 1.0E-10 -1 MOL DIFF SP 1 LAYER 100

 0 0.15 -1 MOL DIFF SP 2 LAYER 001

 0 0.15 -1 MOL DIFF SP 2 LAYER 002

(layers 3 through 99 not shown)

 0 0.15 -1 MOL DIFF SP 2 LAYER 100

For Case 3, the evolution of the salinity and temperature fields are no longer identical; the temperature transition 

zone is much more diffuse than the salinity transition zone (fig. 4). 

Figure 3. Transient movement of the freshwater-seawater transition zone for Case 2. Gray shading represents 

salinity values between 1 and 99 percent of seawater. The solid black lines represent a salinity value equal to 

50-percent seawater for Case 2. Dashed lines represent the 50-percent seawater line for Case 1.
Example Problem 29

Figure 4. Transient movement of the (a) salinity and (b) temperature transition zones for Case 3. Gray shading 

represents the 1 to 99 percent salinity or temperature values. The solid black lines represent the salinity or 

temperature value equal to 50 percent for Case 3. For comparison, dashed lines represent the 50-percent line for 

Case 2.
30 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Case 4

In Case 4, thermal equilibration between the fluid and the solid matrix is included in the simulation. This feature 

is activated using the MT3DMS Reactions (RCT) Package, which has been added to the name file. This “reaction” 

process is activated only for the temperature species (species 2) by entering 2.0 × 10-4 m3/kg for Kd_temp. The 

following RCT input file is used for Case 4:

 1 0 2 0

 0 1760.000(10E12.4) -1 Blk Dnsty Lay 1

 0 1760.000(10E12.4) -1 Blk Dnsty Lay 2

(layers 3 through 99 not shown)

 0 1760.000(10E12.4) -1 Blk Dnsty Lay 100

 00.0000e+00(10E12.4) -1 Kd Species 1 Lay 1

 00.0000e+00(10E12.4) -1 Kd Species 1 Lay 2

(layers 3 through 99 not shown)

 00.0000e+00(10E12.4) -1 Kd Species 1 Lay 100

 02.0000e-04(10E12.4) -1 Kd Species 2 Lay 1

 02.0000e-04(10E12.4) -1 Kd Species 2 Lay 2

(layers 3 through 99 not shown)

 02.0000e-04(10E12.4) -1 Kd Species 2 Lay 100

 00.0000e+00(10E12.4) -1 SP2 Sp 1 not used Lay 1

 00.0000e+00(10E12.4) -1 SP2 Sp 1 not used Lay 2

(layers 3 through 99 not shown)

 00.0000e+00(10E12.4) -1 SP2 Sp 1 not used Lay 100

 00.0000e+00(10E12.4) -1 SP2 Sp 2 not used Lay 1

 00.0000e+00(10E12.4) -1 SP2 Sp 2 not used Lay 2

(layers 3 through 99 not shown)

 00.0000e+00(10E12.4) -1 SP2 Sp 2 not used Lay 100

Simulated salinities and temperatures for Case 4 are shown in figure 5. Simulated salinities are similar to the 

salinities from Case 3; however, the effect of including thermal equilibration with the solid substantially affects 

the transient response of the simulated temperatures. Inspection of figure 5 indicates that the temperature front in 

Case 4 is moving about half as fast as the temperature front in Case 3. This is consistent with a thermal retardation 

coefficient of about 2.0, which can be calculated directly from the thermal properties listed in table 5 and the 

pre-factor of the time derivative on the left-hand side of equation 23.
Example Problem 31

Figure 5. Transient movement of the (a) salinity and (b) temperature transition zones for Case 4. Gray shading 

represents the 1 to 99 percent salinity or temperature values. The solid black lines represent the salinity or 

temperature value equal to 50 percent for Case 4. For comparison, dashed lines represent the 50-percent line for 

Case 3.
32 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Case 5

Heat conduction at the seawater boundary is included in the Case 5 simulation by changing the MT3DMS 

boundary type for the temperature species to ITYPE = -1. By using the constant-head boundary type (ITYPE = 

1 in the SSM input file), only advective transport (no dispersive transport) is allowed across the boundary. With 

this condition, the solute concentration (or temperature) value at the constant-head boundary can change from the 

value specified in the SSM input file. Without a dispersive flux across the boundary, there is no heat conduction. By 

setting ITYPE = -1 for temperature, the assigned temperature value will remain constant at the boundary, and heat 

conduction to or from the boundary will occur. 

The SSM input file for Case 5 is as follows: 

T F F F F F

 300

 300

 1 1 200 3.500e+01 1 3.500000e+001 -1 

 2 1 200 3.500e+01 1 3.500000e+001 -1 

(layers 3 through 99 not shown)

 100 1 200 3.500e+01 1 3.500000e+001 -1 

 1 1 200 3.500e+01 -1 -1 5.000000e+000 

 2 1 200 3.500e+01 -1 -1 5.000000e+000 

(layers 3 through 99 not shown)

 100 1 200 3.500e+01 -1 -1 5.000000e+000 

 1 1 1 0.000e+00 2 0.000000e+000 2.500000e+001 

 2 1 1 0.000e+00 2 0.000000e+000 2.500000e+001 

(layers 3 through 99 not shown)

 100 1 1 0.000e+00 2 0.000000e+000 2.500000e+001 

The first line indicates which flow packages are active. The “T” indicates that the well package is active. The second 

line specifies the maximum number of sources and sinks present in the simulation. The third line indicates the 

number of sources or sinks that are specified in the subsequent data record. The subsequent 300 lines contain layer, 

row, column, a dummy value, ITYPE, and the concentration value for each species. A negative concentration value 

indicates that this boundary type does not apply for this species. Thus, the first 100 lines indicate the right boundary 

has a zero dispersive flux (ITYPE = 1) for salinity (species 1) and that any solute entering from the boundary will 

have a concentration of 35 kg/m3. The next 100 lines with an ITYPE = -1 indicate that temperature will be held 

constant at 5 oC. The remaining 100 lines specify the salinity and temperature values for the inflow along the left 

boundary.

For this particular problem, the use of a constant temperature boundary does not have a substantial effect on 

the results (fig. 6). Simulated salinities for Case 5 are the same as those for Case 4. Including the boundary, heat 

conduction affects simulated temperatures slightly, but only near the vertical ocean boundary.
Example Problem 33

Figure 6. Transient movement of the (a) salinity and (b) temperature transition zones for Case 5. Gray shading 

represents the 1 to 99 percent salinity or temperature values. The solid black lines represent the salinity or 

temperature value equal to 50 percent for Case 5. For comparison, dashed lines represent the 50-percent line for 

Case 4.
34 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Case 6

For Case 6, the effects of fluid viscosity variation on salinity and temperature have been included in the 

simulation. Viscosity effects are included by activating the VSC Package in the name file. The input file for the VSC 

Package is as follows:

#Viscosity Package

-1 MT3DMUFLG

0. 0. VISCMIN VISCMAX

8.904e-4 VISCREF

1 1 NSMUEOS MUTEMPOPT

1 1.923e-6 0. MTMUSPEC(1) DMUDC(1) CMUREF(1)

2 239.4e-7 10. 248.37 133.15 MTMUTEMPSPEC A1 A2 A3 A4

The first line is a comment. The second line indicates that viscosity is a function of one or more MT3DMS species. 

The third line contains the viscosity limiters, which are not active because of the zero value. The fourth line is the 

reference viscosity. The fifth line contains two important flags. The first indicates that one of the species has a 

linear effect on fluid density, and that information for the species will be read on the following line. The second flag 

contains information about the effects of temperature on fluid viscosity. A value of 1 indicates that fluid viscosity 

is calculated using equation 18, and that additional information regarding the temperature/viscosity relation will be 

read later. The sixth line contains the MT3DMS species number that has a linear effect on fluid density, the slope of 

this relation, and the reference concentration this relation is based upon. The last line indicates that species 2 is the 

temperature species and provides constants used in equation 18.

Results from the Case 6 simulation are shown in figure 7, and suggest that viscosity variations have a minimal 

effect on simulated salinities and temperatures for this particular problem. There are some minor differences between 

Cases 5 and 6, particularly at later times, but the overall results are similar.
Example Problem 35

Figure 7. Transient movement of the (a) salinity and (b) temperature transition zones for Case 6. Gray shading 

represents the 1 to 99 percent salinity or temperature values. The solid black lines represent the salinity or 

temperature value equal to 50 percent for Case 6. For comparison, dashed lines represent the 50-percent line for 

Case 5.
36 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

Case 7

For the last case, the effects of pressure (the weight of the overlying column of fluid) on fluid density 

are considered. The pressure effects are included through the VDF input file (shown below) by entering 

4.46 × 10-3 kg m-4 for the change in density with the change in pressure head (DRHODPRHD), and a value of zero 

for the reference pressure head (PRHDREF). 

 -1 1 -1 0 MT3DRHOFLG MFNADVFD NSWTCPL IWTABLE

 0.0000 0.0000 DENSEMIN DENSEMAX

 0.0000 DNSCRIT

 1000.0000 4.46E-03 0. DENSEREF DRHODPRHD PRHDREF

 2 NSRHOEOS

 1 0.7000 0. MTRHOSPEC(1) DRHODC(1) CRHOREF(1)

 2 -0.3750 25. MTRHOSPEC(2) DRHODC(2) CRHOREF(2)

1.0000e-02 FIRSTDT

Results from the Case 7 simulation are shown in figure 8. The compressibility of water due to pressure appears 

to have a minor effect on the simulated salinities and temperatures, particularly at later times. This minor effect is 

largely due to the deep aquifer system (1 km) considered here. For shallow surficial aquifers, the effect of pressure 

on fluid density is expected to be negligible.
Example Problem 37

Figure 8. Transient movement of the (a) salinity and (b) temperature transition zones for Case 7. Gray shading 

represents the 1 to 99 percent salinity or temperature values. The solid black lines represent the salinity or 

temperature value equal to 50 percent for Case 7. For comparison, dashed lines represent the 50-percent line for 

Case 6.
38 SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport

References Cited

Anderson, M.P., 2005, Heat as a ground water tracer: Ground Water, v. 43, no. 6, p. 951-968.

Bear, J., 1972, Dynamics of fluids in porous media: New York, Dover Publications, Inc., 764 p.

Cote, Jean, and Konrad, J.M., 2005, Thermal conductivity of base-course materials: Canadian Geotechnical Journal 

v. 42, p. 61-78.

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.

Engineering Tool Box, 2007, The Engineering Toolbox: Tools and Basic Information for Design, Engineering and 

Construction of Technical Applications: Accessed on May 8, 2007, at http://www.EngineeringToolBox.com

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

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

Guo, Weixing, and Langevin, C.D., 2002, User’s guide to SEAWAT: A computer program for simulation of threedimensional variable-density ground-water 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. 

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

Springer-Verlag, 286 p.

Horai, K.I., 1971, Thermal conductivity of rock-forming minerals: Journal of Geophysical Research, v. 76, no. 5, 

p. 1278-1308.

Hughes, J.D., and Sanford, W.E., 2004, SUTRA-MS a version of SUTRA modified to simulate heat and multiplesolute transport: U.S. Geological Survey Open-File Report 2004-1207, 141 p.

Ingebritsen, S.E., and Sanford, W.E., 1998, Groundwater in geologic processes: United Kingdom, Cambridge 

University Press.

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.

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.

Lusczynski, N.J., 1961, Head and flow of ground water of variable density: Journal of Geophysical Research, v. 66, 

no. 12, p. 4247-4256.

Marble Institute of America, 2007, R-value for natural stone: Accessed May 9, 2007, at http://www.marble-institute.

com/industryresources/rvalue.cfm

Martin, R.J., Bender, S.F., Gaulke, S.W., and Wallace, J., 2001, Simulation of groundwater flow and heat transport 

on Grand Cayman Island, in Seo, Poeter, Zheng, and Poeter, eds., MODFLOW 2001 and Other Modeling 

Odysseys, Conference Proceedings, p. 776-782.
References Cited 39

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

Pawlowski, J., 1991, Veranderliche stoffgroben in der ahnlichkeitstheorie, salle+sauerlander: Frankfurt, 108 p.

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.

Thorne, D.T., Jr., 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 WaterResources Investigations Report 84-4369, 409 p.

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.
Langevin, Thorne, Dausman, Sukop, and Guo 3⁄4 TM 6 A22 3⁄4 - SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport 

,Printed on recycled paper