Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add between_layers to calculate_resistance, and set to False #386

Merged
merged 2 commits into from
Nov 18, 2024
Merged

Add between_layers to calculate_resistance, and set to False #386

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a656e69
Add between_layers to calculate_resistance, and set to False
rubencalje Nov 13, 2024
6ab5312
Implement changes requested by @dbrakenhoff
rubencalje Nov 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 33 additions & 24 deletions nlmod/dims/layers.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,10 +121,11 @@ def calculate_transmissivity(
return T


def calculate_resistance(ds, kv="kv", thickness="thickness", top="top", botm="botm"):
"""Calculate vertical resistance (c) between model layers from the vertical
conductivity (kv) and the thickness. The resistance between two layers is assigned
to the top layer. The bottom model layer gets a resistance of infinity.
def calculate_resistance(
ds, kv="kv", thickness="thickness", top="top", botm="botm", between_layers=False
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe set to None, warn users that to get the default they need to set between_layers=False, and the future default will be True?

):
"""Calculate vertical resistance (c) of model layers from the vertical conductivity
(kv) and the thickness.

Parameters
----------
Expand All @@ -142,6 +143,12 @@ def calculate_resistance(ds, kv="kv", thickness="thickness", top="top", botm="bo
botm : str, optional
name of data variable containing bottoms, only used to calculate thickness if not
available in dataset. By default "botm"
between_layers : bool, optional
If True, calculate the resistance between the layers, which MODFLOW uses to
calculate the flow. The resistance between two layers is then assigned to the
top layer, and the bottom model layer gets a resistance of infinity.
If False, calculate the resistance of the layers themselves. The defauls is
Copy link
Collaborator

@dbrakenhoff dbrakenhoff Nov 13, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

typo: default

False.

Returns
-------
Expand All @@ -152,26 +159,28 @@ def calculate_resistance(ds, kv="kv", thickness="thickness", top="top", botm="bo
thickness = ds[thickness]
else:
thickness = calculate_thickness(ds, top=top, bot=botm)

# nan where layer does not exist (thickness is 0)
thickness_nan = xr.where(thickness == 0, np.nan, thickness)
kv_nan = xr.where(thickness == 0, np.nan, ds[kv])

# backfill thickness and kv to get the right value for the layer below
thickness_bfill = thickness_nan.bfill(dim="layer")
kv_bfill = kv_nan.bfill(dim="layer")

# calculate resistance
c = xr.zeros_like(thickness)
for ilay in range(ds.sizes["layer"] - 1):
ctop = (thickness_nan.sel(layer=ds.layer[ilay]) * 0.5) / kv_nan.sel(
layer=ds.layer[ilay]
)
cbot = (thickness_bfill.sel(layer=ds.layer[ilay + 1]) * 0.5) / kv_bfill.sel(
layer=ds.layer[ilay + 1]
)
c[ilay] = ctop + cbot
c[ilay + 1] = np.inf
if between_layers:
# nan where layer does not exist (thickness is 0)
thickness_nan = xr.where(thickness == 0, np.nan, thickness)
kv_nan = xr.where(thickness == 0, np.nan, ds[kv])

# backfill thickness and kv to get the right value for the layer below
thickness_bfill = thickness_nan.bfill(dim="layer")
kv_bfill = kv_nan.bfill(dim="layer")

# calculate resistance
c = xr.zeros_like(thickness)
for ilay in range(ds.sizes["layer"] - 1):
ctop = (thickness_nan.sel(layer=ds.layer[ilay]) * 0.5) / kv_nan.sel(
layer=ds.layer[ilay]
)
cbot = (thickness_bfill.sel(layer=ds.layer[ilay + 1]) * 0.5) / kv_bfill.sel(
layer=ds.layer[ilay + 1]
)
c[ilay] = ctop + cbot
c[ilay + 1] = np.inf
else:
c = thickness / ds[kv]

if hasattr(c, "long_name"):
c.attrs["long_name"] = "resistance"
Expand Down