diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 04d3a7757..31a2b50a0 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -6,6 +6,15 @@ All notable changes to this project will be documented in this file. The format is based on `Keep a Changelog`_, and this project adheres to `Semantic Versioning`_. +[Unreleased] +------------ + +Added +~~~~~ + +- Added :class:`imod.mf6.Viscosity` package to specify the viscosity of the + groundwater flow model. + [1.0.0rc7] - 2025-10-28 ----------------------- diff --git a/docs/api/mf6.rst b/docs/api/mf6.rst index 4721aef41..1bdddf888 100644 --- a/docs/api/mf6.rst +++ b/docs/api/mf6.rst @@ -396,6 +396,16 @@ Flow Packages UnsaturatedZoneFlow.copy UnsaturatedZoneFlow.is_empty UnsaturatedZoneFlow.get_regrid_methods + Viscosity + Viscosity.write + Viscosity.from_file + Viscosity.to_netcdf + Viscosity.copy + Viscosity.is_empty + Viscosity.get_regrid_methods + Viscosity.mask + Viscosity.regrid_like + Viscosity.clip_box Well Well.cleanup Well.from_imod5_data diff --git a/imod/mf6/__init__.py b/imod/mf6/__init__.py index 48ac96d53..51306e43c 100644 --- a/imod/mf6/__init__.py +++ b/imod/mf6/__init__.py @@ -57,5 +57,6 @@ from imod.mf6.timedis import TimeDiscretization from imod.mf6.uzf import UnsaturatedZoneFlow from imod.mf6.validation_settings import ValidationSettings +from imod.mf6.vsc import Viscosity from imod.mf6.wel import LayeredWell, Well from imod.mf6.write_context import WriteContext diff --git a/imod/mf6/buy.py b/imod/mf6/buy.py index 4ee39dba6..a77709a0d 100644 --- a/imod/mf6/buy.py +++ b/imod/mf6/buy.py @@ -1,26 +1,13 @@ from typing import Optional, Sequence import numpy as np -import xarray as xr from imod.logging import init_log_decorator from imod.mf6.package import Package +from imod.mf6.utilities.dataset import assign_index from imod.schemata import DTypeSchema -def assign_index(arg): - if isinstance(arg, xr.DataArray): - arg = arg.values - elif not isinstance(arg, (np.ndarray, list, tuple)): - raise TypeError("should be a tuple, list, or numpy array") - - arr = np.array(arg) - if arr.ndim != 1: - raise ValueError("should be 1D") - - return xr.DataArray(arr, dims=("index",)) - - class Buoyancy(Package): """ Buoyancy package. This package must be included when performing variable diff --git a/imod/mf6/utilities/dataset.py b/imod/mf6/utilities/dataset.py index d7afea9c8..adcb85403 100644 --- a/imod/mf6/utilities/dataset.py +++ b/imod/mf6/utilities/dataset.py @@ -53,3 +53,16 @@ def assign_datetime_coords( da["time"], unit=time_unit ) return da.assign_coords(time=time) + + +def assign_index(arg): + if isinstance(arg, xr.DataArray): + arg = arg.values + elif not isinstance(arg, (np.ndarray, list, tuple)): + raise TypeError("should be a tuple, list, or numpy array") + + arr = np.array(arg) + if arr.ndim != 1: + raise ValueError("should be 1D") + + return xr.DataArray(arr, dims=("index",)) diff --git a/imod/mf6/vsc.py b/imod/mf6/vsc.py new file mode 100644 index 000000000..df2e4a96f --- /dev/null +++ b/imod/mf6/vsc.py @@ -0,0 +1,209 @@ +from typing import Literal, Optional, Sequence, TypeAlias + +import numpy as np + +from imod.mf6.package import Package +from imod.mf6.utilities.dataset import assign_index +from imod.schemata import DTypeSchema + +ThermalFormulationOption: TypeAlias = Literal["linear", "nonlinear"] + + +class Viscosity(Package): + """ + The Viscosity Package in MODFLOW 6 is used to account for the effects of + solute concentration or temperature on fluid viscosity and thereby their + effects on hydraulic conductivity and stress-package conductance. If the + Viscosity package is used, the Groundwater Transport process must also be + used. In addition, the flow and transport models must be part of the same + simulation. The Viscosity package will adjust the conductances in the model + based on the solute concentrations. + + Parameters + ---------- + + reference_viscosity: float + Fluid reference viscosity used in the equation of state. + viscosity_concentration_slope: sequence of floats + Slope of the (linear) viscosity concentration line used in the viscosity + equation of state. This value will be used when ``thermal_formulation`` + is equal to ``"linear"`` (the default) in the OPTIONS block. When + ``thermal_formulation`` is set to ``"nonlinear"``, a value for DVISCDC + must be specified though it is not used. + reference_concentration: sequence of floats + Reference concentration used in the viscosity equation of state. + modelname: sequence of strings + Name of the GroundwaterTransport (GWT) model used for the + concentrations. + species: sequence of str + Name of the species used to calculate a viscosity value. + temperature_species_name: str + Name of the species to be interpreted as temperature. This species is + used to calculate the temperature-dependent viscosity, using all + thermal_ arguments. + thermal_formulation: str, optional + The thermal formulation to use for the temperature-dependent viscosity. + thermal_a2: float, optional + Is an empirical parameter specified by the user for calculating + viscosity using a nonlinear formulation. If thermal_a2 is not specified, + a default value of 10.0 is assigned (Voss, 1984). + thermal_a3: float, optional + Is an empirical parameter specified by the user for calculating + viscosity using a nonlinear formulation. If thermal_a3 is not specified, + a default value of 248.37 is assigned (Voss, 1984). + thermal_a4: float, optional + Is an empirical parameter specified by the user for calculating + viscosity using a nonlinear formulation. If thermal_a4 is not specified, + a default value of 133.15 is assigned (Voss, 1984). + viscosityfile: str, optional + Name of the binary output file to write viscosity information. + validate: bool, optional + Flag to indicate whether the package should be validated upon + initialization. This raises a ValidationError if package input is + provided in the wrong manner. Defaults to True. + + Examples + -------- + + The viscosity input for a single species called "salinity", which is + simulated by a GWT model called "gwt-1" are specified as follows: + + >>> vsc = imod.mf6.Viscosity( + ... reference_viscosity=8.904E-04, + ... viscosity_concentration_slope=[1.92e-6], + ... reference_concentration=[0.0], + ... modelname=["gwt-1"], + ... species=["salinity"], + ... ) + + Multiple species can be specified by presenting multiple values with an + associated species coordinate. Two species, called "c1" and "c2", simulated + by the GWT models "gwt-1" and "gwt-2" are specified as: + + >>> coords = {"species": ["c1", "c2"]} + >>> vsc = imod.mf6.Viscosity( + ... reference_viscosity=8.904E-04, + ... viscosity_concentration_slope=[1.92e-6, 3.4e-6], + ... reference_concentration=[0.0, 0.0], + ... modelname=["gwt-1", "gwt-2], + ... species=["c1", "c2"], + ... ) + + You can also specify thermal properties, even with a nonlinear thermal + formulation. + + >>> coords = {"species": ["salinity", "temperature"]} + >>> vsc = imod.mf6.Viscosity( + ... reference_viscosity=8.904E-04, + ... viscosity_concentration_slope=[1.92e-6, 0.0], + ... reference_concentration=[0.0, 25.0], + ... modelname=["gwt-1", "gwt-2"], + ... species=["salinity", "temperature"], + ... temperature_species_name="temperature", + ... thermal_formulation="nonlinear", + ... thermal_a2=10.0, + ... thermal_a3=248.37, + ... thermal_a4=133.15, + ... ) + + """ + + _pkg_id = "vsc" + _template = Package._initialize_template(_pkg_id) + + _init_schemata = { + "reference_viscosity": [DTypeSchema(np.floating)], + "viscosity_concentration_slope": [DTypeSchema(np.floating)], + "reference_concentration": [DTypeSchema(np.floating)], + } + _write_schemata = {} + + def __init__( + self, + reference_viscosity: float, + viscosity_concentration_slope: Sequence[float], + reference_concentration: Sequence[float], + modelname: Sequence[str], + species: Sequence[str], + temperature_species_name: Optional[str] = None, + thermal_formulation: ThermalFormulationOption = "linear", + thermal_a2: float = 10.0, + thermal_a3: float = 248.37, + thermal_a4: float = 133.15, + viscosityfile: Optional[str] = None, + validate: bool = True, + ): + dict_dataset = { + "reference_viscosity": reference_viscosity, + # Assign a shared index: this also forces equal lengths + "viscosity_concentration_slope": assign_index( + viscosity_concentration_slope + ), + "reference_concentration": assign_index(reference_concentration), + "modelname": assign_index(modelname), + "species": assign_index(species), + "temperature_species_name": temperature_species_name, + "thermal_formulation": thermal_formulation, + "thermal_a2": thermal_a2, + "thermal_a3": thermal_a3, + "thermal_a4": thermal_a4, + "viscosityfile": viscosityfile, + } + super().__init__(dict_dataset) + self._validate_init_schemata(validate) + + def _render(self, directory, pkgname, globaltimes, binary): + ds = self.dataset + packagedata = [] + + for i, (a, b, c, d) in enumerate( + zip( + ds["viscosity_concentration_slope"].values, + ds["reference_concentration"].values, + ds["modelname"].values, + ds["species"].values, + ) + ): + packagedata.append((i + 1, a, b, c, d)) + + d = { + "nviscspecies": self.dataset["species"].size, + "packagedata": packagedata, + } + + for varname in [ + "temperature_species_name", + "thermal_formulation", + "thermal_a2", + "thermal_a3", + "thermal_a4", + "reference_viscosity", + "viscosityfile", + ]: + value = self.dataset[varname].values[()] + if self._valid(value): + d[varname] = value + + return self._template.render(d) + + def _update_transport_models(self, new_modelnames: Sequence[str]): + """ + The names of the transport models can change in some cases, for example + when partitioning. Use this function to update the names of the + transport models. + """ + transport_model_names = self._get_transport_model_names() + if not len(transport_model_names) == len(new_modelnames): + raise ValueError("the number of transport models cannot be changed.") + for modelname, new_modelname in zip(transport_model_names, new_modelnames): + if modelname not in new_modelname: + raise ValueError( + "new transport model names do not match the old ones. The new names should be equal to the old ones, with a suffix." + ) + self.dataset["modelname"] = assign_index(new_modelnames) + + def _get_transport_model_names(self) -> list[str]: + """ + Returns the names of the transport models used by this buoyancy package. + """ + return list(self.dataset["modelname"].values) diff --git a/imod/templates/mf6/gwf-vsc.j2 b/imod/templates/mf6/gwf-vsc.j2 new file mode 100644 index 000000000..3151a3ddb --- /dev/null +++ b/imod/templates/mf6/gwf-vsc.j2 @@ -0,0 +1,25 @@ +begin options +{%- if reference_viscosity is defined %} + viscref {{reference_viscosity}}{% endif %} +{%- if temperature_species_name is defined %} + temperature_species_name {{temperature_species_name}}{% endif %} +{%- if thermal_formulation is defined %} + thermal_formulation {{thermal_formulation}}{% endif %} +{%- if thermal_a2 is defined %} + thermal_a2 {{thermal_a2}}{% endif %} +{%- if thermal_a3 is defined %} + thermal_a3 {{thermal_a3}}{% endif %} +{%- if thermal_a4 is defined %} + thermal_a4 {{thermal_a4}}{% endif %} +{%- if viscosityfile is defined %} + viscosity fileout {{viscosityfile}}{% endif %} +end options + +begin dimensions + nviscspecies {{nviscspecies}} +end dimensions + +begin packagedata +{%- for iviscspec, dviscdc, cviscref, modelname, auxspeciesname in packagedata %} + {{iviscspec}} {{dviscdc}} {{cviscref}} {{modelname}} {{auxspeciesname}}{% endfor %} +end packagedata \ No newline at end of file diff --git a/imod/tests/conftest.py b/imod/tests/conftest.py index 2e6625d6e..cdc3eb56c 100644 --- a/imod/tests/conftest.py +++ b/imod/tests/conftest.py @@ -40,6 +40,7 @@ circle_model_evt, circle_model_transport, circle_model_transport_multispecies_variable_density, + circle_model_transport_vsc, circle_partitioned, circle_result, circle_result__offset_origins, diff --git a/imod/tests/fixtures/mf6_circle_fixture.py b/imod/tests/fixtures/mf6_circle_fixture.py index a524b14d0..f14a86eb9 100644 --- a/imod/tests/fixtures/mf6_circle_fixture.py +++ b/imod/tests/fixtures/mf6_circle_fixture.py @@ -337,8 +337,86 @@ def circle_model_transport(): # Define the maximum concentration as the initial conditions, also output # options for the transport model, and assign the transport model to the # simulation as well. + transport_model["ic"] = imod.mf6.InitialConditions(start=max_concentration) + transport_model["oc"] = imod.mf6.OutputControl( + save_concentration="last", save_budget="last" + ) + + simulation["transport"] = transport_model + simulation["transport_solver"] = imod.mf6.Solution( + modelnames=["transport"], + print_option="summary", + outer_dvclose=1.0e-4, + outer_maximum=500, + under_relaxation=None, + inner_dvclose=1.0e-4, + inner_rclose=0.001, + inner_maximum=100, + linear_acceleration="bicgstab", + scaling_method=None, + reordering_method=None, + relaxation_factor=0.97, + ) + simtimes = pd.date_range(start="2000-01-01", end="2001-01-01", freq="W") + simulation.create_time_discretization(additional_times=simtimes) + return simulation + + +@pytest.fixture(scope="function") +def circle_model_transport_vsc(): + al = 0.001 + porosity = 0.3 max_concentration = 35.0 min_concentration = 0.0 + max_density = 1025.0 + min_density = 1000.0 + + simulation = make_circle_model_flow_with_transport_data(["salinity"]) + gwf_model = simulation["GWF_1"] + + slope = (max_density - min_density) / (max_concentration - min_concentration) + gwf_model["buoyancy"] = imod.mf6.Buoyancy( + reference_density=min_density, + modelname=["transport"], + reference_concentration=[min_concentration], + density_concentration_slope=[slope], + species=["salinity"], + ) + gwf_model["vsc"] = imod.mf6.Viscosity( + reference_viscosity=8.904e-04, + viscosity_concentration_slope=[3.3e-5], + reference_concentration=[min_concentration], + modelname=["transport"], + species=["salinity"], + viscosityfile="viscosity.bin", + ) + transport_model = imod.mf6.GroundwaterTransportModel(save_flows=True) + transport_model["ssm"] = imod.mf6.SourceSinkMixing.from_flow_model( + gwf_model, "salinity", save_flows=True + ) + transport_model["disv"] = gwf_model["disv"] + + # %% + # Now we define some transport packages for simulating the physical processes + # of advection, mechanical dispersion, and molecular diffusion dispersion. This + # example is transient, and the volume available for storage is the porosity, + # in this case 0.10. + + transport_model["dsp"] = imod.mf6.Dispersion( + diffusion_coefficient=1e-4, + longitudinal_horizontal=al, + transversal_horizontal1=al * 0.1, + transversal_vertical=al * 0.01, + xt3d_off=False, + xt3d_rhs=False, + ) + transport_model["adv"] = imod.mf6.AdvectionUpstream() + transport_model["mst"] = imod.mf6.MobileStorageTransfer(porosity, save_flows=True) + + # %% + # Define the maximum concentration as the initial conditions, also output + # options for the transport model, and assign the transport model to the + # simulation as well. transport_model["ic"] = imod.mf6.InitialConditions(start=max_concentration) transport_model["oc"] = imod.mf6.OutputControl( save_concentration="last", save_budget="last" diff --git a/imod/tests/test_mf6/test_circle.py b/imod/tests/test_mf6/test_circle.py index cb3142019..f31345012 100644 --- a/imod/tests/test_mf6/test_circle.py +++ b/imod/tests/test_mf6/test_circle.py @@ -156,6 +156,29 @@ def test_simulation_write_and_run_transport(circle_model_transport, tmp_path): assert concentration.shape == (52, 2, 216) +def test_simulation_write_and_run_transport_vsc(circle_model_transport_vsc, tmp_path): + simulation = circle_model_transport_vsc + + with pytest.raises( + RuntimeError, match="Simulation circle has not been written yet." + ): + circle_model_transport_vsc.run() + + modeldir = tmp_path / "circle_transport" + simulation.write(modeldir) + simulation.run() + head = simulation.open_head() + concentration = simulation.open_concentration() + + assert isinstance(head, xu.UgridDataArray) + assert head.dims == ("time", "layer", "mesh2d_nFaces") + assert head.shape == (52, 2, 216) + + assert isinstance(concentration, xu.UgridDataArray) + assert concentration.dims == ("time", "layer", "mesh2d_nFaces") + assert concentration.shape == (52, 2, 216) + + def test_simulation_clip_and_state_at_boundary(circle_model_transport, tmp_path): # Arrange simulation = circle_model_transport diff --git a/imod/tests/test_mf6/test_mf6_buoy.py b/imod/tests/test_mf6/test_mf6_buoy.py index c44d0dedc..f5c9e27fc 100644 --- a/imod/tests/test_mf6/test_mf6_buoy.py +++ b/imod/tests/test_mf6/test_mf6_buoy.py @@ -35,8 +35,7 @@ def test_buoyancy_package_simple(): end packagedata""" ) actual = buy._render(directory, "buy", globaltimes, False) - print(actual) - print(expected) + assert actual == expected @@ -71,8 +70,7 @@ def test_buoyancy_package_full(): end packagedata""" ) actual = buy._render(directory, "buy", globaltimes, False) - print(actual) - print(expected) + assert actual == expected diff --git a/imod/tests/test_mf6/test_mf6_vsc.py b/imod/tests/test_mf6/test_mf6_vsc.py new file mode 100644 index 000000000..7a72940d2 --- /dev/null +++ b/imod/tests/test_mf6/test_mf6_vsc.py @@ -0,0 +1,140 @@ +import pathlib +import textwrap + +import numpy as np +import pytest + +import imod + + +def test_viscosity_package_simple(): + vsc = imod.mf6.Viscosity( + reference_viscosity=8.904e-04, + viscosity_concentration_slope=[1.92e-6], + reference_concentration=[0.0], + modelname=["gwt-1"], + species=["salinity"], + ) + + directory = pathlib.Path("mymodel") + globaltimes = [np.datetime64("2000-01-01")] + + expected = textwrap.dedent( + """\ + begin options + viscref 0.0008904 + thermal_formulation linear + thermal_a2 10.0 + thermal_a3 248.37 + thermal_a4 133.15 + end options + + begin dimensions + nviscspecies 1 + end dimensions + + begin packagedata + 1 1.92e-06 0.0 gwt-1 salinity + end packagedata""" + ) + actual = vsc._render(directory, "vsc", globaltimes, False) + + assert actual == expected + + +def test_viscosity_package_full(): + vsc = imod.mf6.Viscosity( + reference_viscosity=8.904e-04, + viscosity_concentration_slope=[1.92e-6, 0.0], + reference_concentration=[0.0, 25.0], + modelname=["gwt-1", "gwt-2"], + species=["salinity", "temperature"], + temperature_species_name="temperature", + thermal_formulation="nonlinear", + thermal_a2=20.0, + thermal_a3=348.37, + thermal_a4=233.15, + viscosityfile="testfile.vsc", + ) + + directory = pathlib.Path("mymodel") + globaltimes = [np.datetime64("2000-01-01")] + + expected = textwrap.dedent( + """\ + begin options + viscref 0.0008904 + temperature_species_name temperature + thermal_formulation nonlinear + thermal_a2 20.0 + thermal_a3 348.37 + thermal_a4 233.15 + viscosity fileout testfile.vsc + end options + + begin dimensions + nviscspecies 2 + end dimensions + + begin packagedata + 1 1.92e-06 0.0 gwt-1 salinity + 2 0.0 25.0 gwt-2 temperature + end packagedata""" + ) + actual = vsc._render(directory, "vsc", globaltimes, False) + + assert actual == expected + + +def test_viscosity_package_update_transport_names(): + vsc = imod.mf6.Viscosity( + reference_viscosity=8.904e-04, + viscosity_concentration_slope=[1.92e-6, 0.0], + reference_concentration=[0.0, 25.0], + modelname=["gwt-1", "gwt-2"], + species=["salinity", "temperature"], + temperature_species_name="temperature", + thermal_formulation="nonlinear", + ) + vsc._update_transport_models(["gwt-1_0", "gwt-2_0"]) + directory = pathlib.Path("mymodel") + globaltimes = [np.datetime64("2000-01-01")] + + expected = textwrap.dedent( + """\ + begin options + viscref 0.0008904 + temperature_species_name temperature + thermal_formulation nonlinear + thermal_a2 10.0 + thermal_a3 248.37 + thermal_a4 133.15 + end options + + begin dimensions + nviscspecies 2 + end dimensions + + begin packagedata + 1 1.92e-06 0.0 gwt-1_0 salinity + 2 0.0 25.0 gwt-2_0 temperature + end packagedata""" + ) + actual = vsc._render(directory, "vsc", globaltimes, False) + + assert actual == expected + + +def test_viscosity_package_update_transport_names_check(): + vsc = imod.mf6.Viscosity( + reference_viscosity=8.904e-04, + viscosity_concentration_slope=[1.92e-6, 0.0], + reference_concentration=[0.0, 25.0], + modelname=["gwt-2", "gwt-1"], + species=["salinity", "temperature"], + temperature_species_name="temperature", + thermal_formulation="nonlinear", + ) + # update the transport models, but in the wrong order + with pytest.raises(ValueError): + vsc._update_transport_models(["gwt-1_0", "gwt-2_0"])