Skip to content

Commit b0dc733

Browse files
JoerivanEngelenJoerivanEngelen
authored
Issue #1597 add vsc package (#1633)
Fixes #1597 # Description Adds viscosity package. All optinional options which have a default value in MODFLOW 6, have these enforced here as well. This to explicitly write these values in the .vsc file. # Checklist - [x] Links to correct issue - [x] Update changelog, if changes affect users - [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737`` - [x] Unit tests were added - [ ] **If feature added**: Added/extended example - [x] **If feature added**: Added feature to API documentation
1 parent 87bbed8 commit b0dc733

File tree

12 files changed

+512
-18
lines changed

12 files changed

+512
-18
lines changed

‎docs/api/changelog.rst‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,15 @@ All notable changes to this project will be documented in this file.
66
The format is based on `Keep a Changelog`_, and this project adheres to
77
`Semantic Versioning`_.
88

9+
[Unreleased]
10+
------------
11+
12+
Added
13+
~~~~~
14+
15+
- Added :class:`imod.mf6.Viscosity` package to specify the viscosity of the
16+
groundwater flow model.
17+
918
[1.0.0rc7] - 2025-10-28
1019
-----------------------
1120

‎docs/api/mf6.rst‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -396,6 +396,16 @@ Flow Packages
396396
UnsaturatedZoneFlow.copy
397397
UnsaturatedZoneFlow.is_empty
398398
UnsaturatedZoneFlow.get_regrid_methods
399+
Viscosity
400+
Viscosity.write
401+
Viscosity.from_file
402+
Viscosity.to_netcdf
403+
Viscosity.copy
404+
Viscosity.is_empty
405+
Viscosity.get_regrid_methods
406+
Viscosity.mask
407+
Viscosity.regrid_like
408+
Viscosity.clip_box
399409
Well
400410
Well.cleanup
401411
Well.from_imod5_data

‎imod/mf6/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -57,5 +57,6 @@
5757
from imod.mf6.timedis import TimeDiscretization
5858
from imod.mf6.uzf import UnsaturatedZoneFlow
5959
from imod.mf6.validation_settings import ValidationSettings
60+
from imod.mf6.vsc import Viscosity
6061
from imod.mf6.wel import LayeredWell, Well
6162
from imod.mf6.write_context import WriteContext

‎imod/mf6/buy.py‎

Lines changed: 1 addition & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,26 +1,13 @@
11
from typing import Optional, Sequence
22

33
import numpy as np
4-
import xarray as xr
54

65
from imod.logging import init_log_decorator
76
from imod.mf6.package import Package
7+
from imod.mf6.utilities.dataset import assign_index
88
from imod.schemata import DTypeSchema
99

1010

11-
def assign_index(arg):
12-
if isinstance(arg, xr.DataArray):
13-
arg = arg.values
14-
elif not isinstance(arg, (np.ndarray, list, tuple)):
15-
raise TypeError("should be a tuple, list, or numpy array")
16-
17-
arr = np.array(arg)
18-
if arr.ndim != 1:
19-
raise ValueError("should be 1D")
20-
21-
return xr.DataArray(arr, dims=("index",))
22-
23-
2411
class Buoyancy(Package):
2512
"""
2613
Buoyancy package. This package must be included when performing variable

‎imod/mf6/utilities/dataset.py‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,3 +53,16 @@ def assign_datetime_coords(
5353
da["time"], unit=time_unit
5454
)
5555
return da.assign_coords(time=time)
56+
57+
58+
def assign_index(arg):
59+
if isinstance(arg, xr.DataArray):
60+
arg = arg.values
61+
elif not isinstance(arg, (np.ndarray, list, tuple)):
62+
raise TypeError("should be a tuple, list, or numpy array")
63+
64+
arr = np.array(arg)
65+
if arr.ndim != 1:
66+
raise ValueError("should be 1D")
67+
68+
return xr.DataArray(arr, dims=("index",))

‎imod/mf6/vsc.py‎

Lines changed: 209 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,209 @@
1+
from typing import Literal, Optional, Sequence, TypeAlias
2+
3+
import numpy as np
4+
5+
from imod.mf6.package import Package
6+
from imod.mf6.utilities.dataset import assign_index
7+
from imod.schemata import DTypeSchema
8+
9+
ThermalFormulationOption: TypeAlias = Literal["linear", "nonlinear"]
10+
11+
12+
class Viscosity(Package):
13+
"""
14+
The Viscosity Package in MODFLOW 6 is used to account for the effects of
15+
solute concentration or temperature on fluid viscosity and thereby their
16+
effects on hydraulic conductivity and stress-package conductance. If the
17+
Viscosity package is used, the Groundwater Transport process must also be
18+
used. In addition, the flow and transport models must be part of the same
19+
simulation. The Viscosity package will adjust the conductances in the model
20+
based on the solute concentrations.
21+
22+
Parameters
23+
----------
24+
25+
reference_viscosity: float
26+
Fluid reference viscosity used in the equation of state.
27+
viscosity_concentration_slope: sequence of floats
28+
Slope of the (linear) viscosity concentration line used in the viscosity
29+
equation of state. This value will be used when ``thermal_formulation``
30+
is equal to ``"linear"`` (the default) in the OPTIONS block. When
31+
``thermal_formulation`` is set to ``"nonlinear"``, a value for DVISCDC
32+
must be specified though it is not used.
33+
reference_concentration: sequence of floats
34+
Reference concentration used in the viscosity equation of state.
35+
modelname: sequence of strings
36+
Name of the GroundwaterTransport (GWT) model used for the
37+
concentrations.
38+
species: sequence of str
39+
Name of the species used to calculate a viscosity value.
40+
temperature_species_name: str
41+
Name of the species to be interpreted as temperature. This species is
42+
used to calculate the temperature-dependent viscosity, using all
43+
thermal_ arguments.
44+
thermal_formulation: str, optional
45+
The thermal formulation to use for the temperature-dependent viscosity.
46+
thermal_a2: float, optional
47+
Is an empirical parameter specified by the user for calculating
48+
viscosity using a nonlinear formulation. If thermal_a2 is not specified,
49+
a default value of 10.0 is assigned (Voss, 1984).
50+
thermal_a3: float, optional
51+
Is an empirical parameter specified by the user for calculating
52+
viscosity using a nonlinear formulation. If thermal_a3 is not specified,
53+
a default value of 248.37 is assigned (Voss, 1984).
54+
thermal_a4: float, optional
55+
Is an empirical parameter specified by the user for calculating
56+
viscosity using a nonlinear formulation. If thermal_a4 is not specified,
57+
a default value of 133.15 is assigned (Voss, 1984).
58+
viscosityfile: str, optional
59+
Name of the binary output file to write viscosity information.
60+
validate: bool, optional
61+
Flag to indicate whether the package should be validated upon
62+
initialization. This raises a ValidationError if package input is
63+
provided in the wrong manner. Defaults to True.
64+
65+
Examples
66+
--------
67+
68+
The viscosity input for a single species called "salinity", which is
69+
simulated by a GWT model called "gwt-1" are specified as follows:
70+
71+
>>> vsc = imod.mf6.Viscosity(
72+
... reference_viscosity=8.904E-04,
73+
... viscosity_concentration_slope=[1.92e-6],
74+
... reference_concentration=[0.0],
75+
... modelname=["gwt-1"],
76+
... species=["salinity"],
77+
... )
78+
79+
Multiple species can be specified by presenting multiple values with an
80+
associated species coordinate. Two species, called "c1" and "c2", simulated
81+
by the GWT models "gwt-1" and "gwt-2" are specified as:
82+
83+
>>> coords = {"species": ["c1", "c2"]}
84+
>>> vsc = imod.mf6.Viscosity(
85+
... reference_viscosity=8.904E-04,
86+
... viscosity_concentration_slope=[1.92e-6, 3.4e-6],
87+
... reference_concentration=[0.0, 0.0],
88+
... modelname=["gwt-1", "gwt-2],
89+
... species=["c1", "c2"],
90+
... )
91+
92+
You can also specify thermal properties, even with a nonlinear thermal
93+
formulation.
94+
95+
>>> coords = {"species": ["salinity", "temperature"]}
96+
>>> vsc = imod.mf6.Viscosity(
97+
... reference_viscosity=8.904E-04,
98+
... viscosity_concentration_slope=[1.92e-6, 0.0],
99+
... reference_concentration=[0.0, 25.0],
100+
... modelname=["gwt-1", "gwt-2"],
101+
... species=["salinity", "temperature"],
102+
... temperature_species_name="temperature",
103+
... thermal_formulation="nonlinear",
104+
... thermal_a2=10.0,
105+
... thermal_a3=248.37,
106+
... thermal_a4=133.15,
107+
... )
108+
109+
"""
110+
111+
_pkg_id = "vsc"
112+
_template = Package._initialize_template(_pkg_id)
113+
114+
_init_schemata = {
115+
"reference_viscosity": [DTypeSchema(np.floating)],
116+
"viscosity_concentration_slope": [DTypeSchema(np.floating)],
117+
"reference_concentration": [DTypeSchema(np.floating)],
118+
}
119+
_write_schemata = {}
120+
121+
def __init__(
122+
self,
123+
reference_viscosity: float,
124+
viscosity_concentration_slope: Sequence[float],
125+
reference_concentration: Sequence[float],
126+
modelname: Sequence[str],
127+
species: Sequence[str],
128+
temperature_species_name: Optional[str] = None,
129+
thermal_formulation: ThermalFormulationOption = "linear",
130+
thermal_a2: float = 10.0,
131+
thermal_a3: float = 248.37,
132+
thermal_a4: float = 133.15,
133+
viscosityfile: Optional[str] = None,
134+
validate: bool = True,
135+
):
136+
dict_dataset = {
137+
"reference_viscosity": reference_viscosity,
138+
# Assign a shared index: this also forces equal lengths
139+
"viscosity_concentration_slope": assign_index(
140+
viscosity_concentration_slope
141+
),
142+
"reference_concentration": assign_index(reference_concentration),
143+
"modelname": assign_index(modelname),
144+
"species": assign_index(species),
145+
"temperature_species_name": temperature_species_name,
146+
"thermal_formulation": thermal_formulation,
147+
"thermal_a2": thermal_a2,
148+
"thermal_a3": thermal_a3,
149+
"thermal_a4": thermal_a4,
150+
"viscosityfile": viscosityfile,
151+
}
152+
super().__init__(dict_dataset)
153+
self._validate_init_schemata(validate)
154+
155+
def _render(self, directory, pkgname, globaltimes, binary):
156+
ds = self.dataset
157+
packagedata = []
158+
159+
for i, (a, b, c, d) in enumerate(
160+
zip(
161+
ds["viscosity_concentration_slope"].values,
162+
ds["reference_concentration"].values,
163+
ds["modelname"].values,
164+
ds["species"].values,
165+
)
166+
):
167+
packagedata.append((i + 1, a, b, c, d))
168+
169+
d = {
170+
"nviscspecies": self.dataset["species"].size,
171+
"packagedata": packagedata,
172+
}
173+
174+
for varname in [
175+
"temperature_species_name",
176+
"thermal_formulation",
177+
"thermal_a2",
178+
"thermal_a3",
179+
"thermal_a4",
180+
"reference_viscosity",
181+
"viscosityfile",
182+
]:
183+
value = self.dataset[varname].values[()]
184+
if self._valid(value):
185+
d[varname] = value
186+
187+
return self._template.render(d)
188+
189+
def _update_transport_models(self, new_modelnames: Sequence[str]):
190+
"""
191+
The names of the transport models can change in some cases, for example
192+
when partitioning. Use this function to update the names of the
193+
transport models.
194+
"""
195+
transport_model_names = self._get_transport_model_names()
196+
if not len(transport_model_names) == len(new_modelnames):
197+
raise ValueError("the number of transport models cannot be changed.")
198+
for modelname, new_modelname in zip(transport_model_names, new_modelnames):
199+
if modelname not in new_modelname:
200+
raise ValueError(
201+
"new transport model names do not match the old ones. The new names should be equal to the old ones, with a suffix."
202+
)
203+
self.dataset["modelname"] = assign_index(new_modelnames)
204+
205+
def _get_transport_model_names(self) -> list[str]:
206+
"""
207+
Returns the names of the transport models used by this buoyancy package.
208+
"""
209+
return list(self.dataset["modelname"].values)

‎imod/templates/mf6/gwf-vsc.j2‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
begin options
2+
{%- if reference_viscosity is defined %}
3+
viscref {{reference_viscosity}}{% endif %}
4+
{%- if temperature_species_name is defined %}
5+
temperature_species_name {{temperature_species_name}}{% endif %}
6+
{%- if thermal_formulation is defined %}
7+
thermal_formulation {{thermal_formulation}}{% endif %}
8+
{%- if thermal_a2 is defined %}
9+
thermal_a2 {{thermal_a2}}{% endif %}
10+
{%- if thermal_a3 is defined %}
11+
thermal_a3 {{thermal_a3}}{% endif %}
12+
{%- if thermal_a4 is defined %}
13+
thermal_a4 {{thermal_a4}}{% endif %}
14+
{%- if viscosityfile is defined %}
15+
viscosity fileout {{viscosityfile}}{% endif %}
16+
end options
17+
18+
begin dimensions
19+
nviscspecies {{nviscspecies}}
20+
end dimensions
21+
22+
begin packagedata
23+
{%- for iviscspec, dviscdc, cviscref, modelname, auxspeciesname in packagedata %}
24+
{{iviscspec}} {{dviscdc}} {{cviscref}} {{modelname}} {{auxspeciesname}}{% endfor %}
25+
end packagedata

‎imod/tests/conftest.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@
4040
circle_model_evt,
4141
circle_model_transport,
4242
circle_model_transport_multispecies_variable_density,
43+
circle_model_transport_vsc,
4344
circle_partitioned,
4445
circle_result,
4546
circle_result__offset_origins,

0 commit comments

Comments
 (0)