65. Viscosity
In hydrogeologic settings where variations in viscosity significantly impact groundwater flows, the viscosity package can be activated to scale intercell conductance in response to temperature or concentration variations. Adjustment of intercell conductance may, in turn, have important implications for solute transport or geothermal model applications. For example, a migrating contaminant plume may be delayed, sped-up, or warped in settings where significant variation in groundwater temperature persists within the study site.
65.1. Example description
Application of the viscosity (VSC) package is demonstrated with a simple two-dimensional (2D) test problem. The 2D demonstration model is adapted from test case 3 in (Zheng & Wang, 1999) titled “Two-dimensional transport in a uniform flow field.” Specified heads drive groundwater flow through a confined system from left to right. No flow boundaries are specified along the other model edges. Water injected into the middle of the model domain enters with a specified concentration of 1,000 \(mg/L\). In the original problem, the simulated development of a contaminant plume was compared to an established analytical solution (Wilson & Miller, 1978) to confirm that the advective and dispersive spread of the contaminant was accurately simulated.
For this demonstration, the original groundwater flow (GWF) and groundwater solute transport (GWT) models remain unchanged and are referred to as the “baseline” GWF and GWT models hereafter (Figure 65.1). However, an additional 3 models are added to the MODFLOW 6 simulation to explore the importance of viscosity effects on the contaminant plume. This results in five models in a single MODFLOW 6 simulation. The three additional models include one each of a GWF, GWT, and groundwater energy transport (GWE) model. The only difference between the baseline GWF model and the “viscosity GWF” model (Figure 65.1) is the activation of the VSC package that accounts for the groundwater temperatures calculated by the GWE model included in the simulation.

Figure 65.1 The MODFLOW 6 simulation setup includes two GWF models (blue rectangles), two GWT models (green rectangles), and one GWE model (orange rectangle). Solid arrows highlight the exchange of information among the models included in the simulation. Dashed arrows show the generation of output used to highlight the effects of viscosity on a solute transport model.
Moreover, for this demonstration, two MODFLOW 6 simulations are run that vary the magnitude of the temperature gradient between the upper and lower model boundaries. In the first scenario, a 30\(^{\circ}C\) difference is simulated between the upper and lower boundary (Table 65.1). A 90\(^{\circ}C\) gradient exists in the second MODFLOW 6 simulation. The initial temperatures in either simulation reflect a linear gradient between the upper and lower boundary temperatures (Table 65.1). Water entering the left side of the model domain enters at a temperature consistent with the initial temperature setup (i.e., varies linearly between the specified temperatures of the upper and lower boundaries).
Scenario |
Scenario Name |
Parameter |
Value |
---|---|---|---|
1 |
30degGrad |
temp upper (\(^{\circ}C\)) |
4.0 |
temp lower (\(^{\circ}C\)) |
34.0 |
||
2 |
90degGrad |
temp upper (\(^{\circ}C\)) |
4.0 |
temp lower (\(^{\circ}C\)) |
94.0 |
Inputs to the VSC package, including parameters for the nonlinear viscosity formulation, are listed in table Table 65.2. Within MODFLOW 6, specifically the VSC package (within the GWF model), pointers to the simulated groundwater temperature enable the VSC package to modify intercell conductance calculated by the node-property flow (NPF) package. Altered intercell conductance in turn affects the flow solution used by the viscosity-affected GWT model (Figure 65.1) that would otherwise be the same as the baseline GWT model.
Parameter |
Value |
---|---|
Number of layers |
1 |
Number of rows |
31 |
Number of columns |
46 |
Column width (\(m\)) |
10.0 |
Row width (\(m\)) |
10.0 |
Layer thickness (\(m\)) |
10.0 |
Top of the model (\(m\)) |
0.0 |
Porosity |
0.3 |
Simulation time (\(days\)) |
365 |
Horizontal hydraulic conductivity (\(m/d\)) |
1.0 |
Volumetric injection rate (\(m^3/d\)) |
1.0 |
Concentration of injected water (\(mg/L\)) |
1000.0 |
Longitudinal dispersivity (\(m\)) |
10.0 |
Ratio of transverse to longitudinal dispersivity |
0.3 |
Reference viscosity (\(\frac{kg}{m \cdot sec}\)) |
8.904e-4 |
Viscosity representation (\(-\)) |
nonlinear |
User-specified parameter in nonlinear formulation (\(-\)) |
10.0 |
User-specified parameter in nonlinear formulation (\(-\)) |
248.37 |
User-specified parameter in nonlinear formulation (\(-\)) |
133.16 |
Species ID for use in a VSC package (\(-\)) |
0 |
not used for the nonlinear formulation, but must be specified |
0.0 |
Reference temperature for viscosity (\(^{\circ}C\)) |
20.0 |
Groundwater temperature (\(^{\circ}C\)) |
4.0 |
Advection scheme (\(-\)) |
TVD |
Starting concentration (\(mg/L\)) |
0.0 |
Heat capacity of aquifer material (\(\frac{J}{kg \cdot ^{\circ} C}\)) |
800.0 |
Heat capacity of water (\(\frac{J}{kg \cdot ^{\circ} C}\)) |
4180.0 |
Density of water (\(kg/m^3\)) |
1000.0 |
Density of dry solid aquifer material (\(kg/m^3\)) |
1200.0 |
Latent heat of vaporization (\(kJ/kg\)) |
2500.0 |
Thermal conductivity of the aquifer material (\(\frac{W}{m \cdot ^{\circ} C}\)) |
0.25 |
Thermal conductivity of water (\(\frac{W}{m \cdot ^{\circ} C}\)) |
0.58 |
Concentration of boundary inflow (\(mg/L\)) |
0.0 |
65.2. Example Results
For the scenario with a 30\(^{\circ}C\) temperature gradient, output from the baseline and viscosity-affected GWT models is collected, differenced, and displayed in figure
Figure 65.2. Warm colors show where higher
concentrations are simulated by the baseline GWT model. An alternative way to think about this result is that ignoring the effects of viscosity can lead to an over-prediction of concentrations in portions of the migrating plume. Conversely, wherever the differenced concentration field is negative, meaning the concentrations in the viscosity-affected model run are greater than in the baseline model run, the calculated concentrations in the no viscosity model may be considered under predicted. Interestingly, the plume is asymmetric along the direction of flow. Below y coordinate values of approximately 130 \(m\), concentrations in the baseline GWT simulation appear to be over-predicted. For y coordinate values greater than 130 \(m\), there is a mixture of under- and over-prediction of concentrations within the baseline model. The leading edge of the plume (left side) shows a general over-prediction of the concentration in the baseline GWT model while the down-gradient region of the plume is under-predicted.
In the 90\(^{\circ}C\) temperature gradient scenario, results are markedly different relative to the 30\(^{\circ}C\) temperature gradient scenario (Figure 65.3. First, the relative position of the under- and over-predicted concentrations is reversed in the 90\(^{\circ}C\) temperature gradient scenario compared to the 30\(^{\circ}C\) temperature gradient scenario. That is, the region of under-predicted concentrations is to the left of the over-predicted region in the 90\(^{\circ}C\) temperature gradient scenario. Second, the down-gradient portion of the plume bends to the right relative to the direction of flow. Third, the extent and magnitude of the spread of the plume is greater in the 90\(^{\circ}C\) temperature gradient scenario owing to the higher temperatures that reduce viscosity and lead to higher conductance in the associated GWF model.
In summary, settings with stark contrasts in groundwater temperature may also have significant variations in the flow field as a result of viscosity effects. Thus, activation of the viscosity package may appropriately adjust the intercell conductance values in response to viscosity and reduce compansatory adjustment of the hydraulic conductivity field during a calibration routine, for example.

Figure 65.2 The difference between a flow-transport model setup that ignores the effects of viscosity (baseline simulation) on a migrating plume and a model setup that simulates a 30\(^{\circ}C\) temperature gradient in the groundwater to demonstrate the effect of viscosity on the migrating plume. Differences are calculated after 365 days of simulation time. The difference is calculated by subtracting the 2D concentration field calculated by the model that accounts for viscosity from the model run that ignores viscosity.

Figure 65.3 Similar to figure Figure 65.2, the difference between two concentration fields calculated by two closely related model runs. In one model setup, the effect of viscosity is ignored whereas the other model setup simulates a 90\(^{\circ}C\) temperature gradient within the model domain. As with the 30\(^{\circ}C\) temperature gradient, differences are calculated after 365 days of simulation time.
65.3. References Cited
Wilson, J. L., & Miller, P. J. (1978). Two-dimensional plume in uniform ground-water flow. Journal of the Hydraulics Division, 104(4), 503–514. https://doi.org/10.1061/JYCEAJ.0004975
Zheng, C., & Wang, P. P. (1999). MT3DMS—a modular three-dimensional multi-species transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; documentation and user’s guide.
65.4. Jupyter Notebook
The Jupyter notebook used to create the MODFLOW 6 input files for this example and post-process the results is: