Merge pull request #231 from ArtesiaWater/zeeland

Updates based on model in Zeeland

docs/examples/12_layer_generation.ipynb

@@ -162,7 +162,7 @@
    "id": "98b6d820",
    "metadata": {},
    "source": [
-    "When we plot the model top, we now see that in the cells with an equal resolution to that of Regis (or smaller) no infomration is lost is lost."
+    "When we plot the model top, we now see that in the cells with an equal resolution to that of Regis (or smaller) no information is lost."
    ]
   },
   {

nlmod/dims/base.py

@@ -225,12 +225,18 @@ def extrapolate_ds(ds, mask=None):
             continue
         data = ds[key].data
         if ds[key].dims == dims:
-            data[mask] = data[~mask][i]
+            if np.isnan(data[mask]).sum() > 0:  # do not update if no NaNs
+                data[mask] = data[~mask, i]
         elif ds[key].dims == ("layer",) + dims:
             for lay in range(len(ds["layer"])):
-                data[lay][mask] = data[lay][~mask][i]
+                if np.isnan(data[lay, mask]).sum() > 0:  # do not update if no NaNs
+                    data[lay, mask] = data[lay, ~mask][i]
         else:
-            raise (Exception(f"Dimensions {ds[key].dims} not supported"))
+            logger.warning(
+                f"Data variable '{key}' not extrapolated because "
+                f"dimensions are not {dims}."
+            )
+            # raise (Exception(f"Dimensions {ds[key].dims} not supported"))
         # make sure to set the data (which for some reason is sometimes needed)
         ds[key].data = data
     return ds

nlmod/dims/grid.py

@@ -344,8 +344,8 @@ def refine(
         If False nan layers are kept which might be usefull if you want
         to keep some layers that exist in other models. The default is True.
     model_coordinates : bool, optional
-        When model_coordinates is True, the features supplied in refinement features are
-        allready in model-coordinates. Only used when a grid is rotated. The default is
+        When model_coordinates is True, the features supplied in refinement_features are
+        already in model-coordinates. Only used when a grid is rotated. The default is
         False.
 
     Returns
@@ -924,7 +924,7 @@ def da_to_reclist(
     else:
         if first_active_layer:
             fal = get_first_active_layer(ds)
-            cellids = np.where((mask) & (fal != fal.attrs["_FillValue"]))
+            cellids = np.where((mask.squeeze()) & (fal != fal.attrs["_FillValue"]))
             layers = col_to_list(fal, ds, cellids)
         elif only_active_cells:
             cellids = np.where((mask) & (ds["idomain"][layer] == 1))

nlmod/dims/layers.py

@@ -61,6 +61,73 @@ def calculate_thickness(ds, top="top", bot="botm"):
     return thickness
 
 
+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.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        dataset containing information about top and bottom elevations
+        of layers
+    kv : str, optional
+        name of data variable containing vertical conductivity, by default 'kv'
+    thickness : str, optional
+        name of data variable containing thickness, if this data variable does not exists
+        thickness is calculated using top and botm. By default 'thickness'
+    top : str, optional
+        name of data variable containing tops, only used to calculate thickness if not 
+        available in dataset. By default "top"
+    botm : str, optional
+        name of data variable containing bottoms, only used to calculate thickness if not 
+        available in dataset. By default "botm"
+
+    Returns
+    -------
+    c : xarray.DataArray
+        DataArray containing vertical resistance (c). NaN where layer thickness is zero
+    """
+
+    if thickness in ds:
+        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.dims['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 hasattr(c, "long_name"):
+        c.attrs["long_name"] = "resistance"
+    if hasattr(c, "standard_name"):
+        c.attrs["standard_name"] = "c"
+    if hasattr(thickness, "units"):
+        if hasattr(ds[kv], 'units'):
+            if ds[kv].units == "m/day" and thickness.units in ['m', 'mNAP']:
+                c.attrs["units"] = "day"
+            else:
+                c.attrs["units"] = ""
+        else:
+            c.attrs["units"] = ""
+
+
+    return c
+
+
 def split_layers_ds(
     ds, split_dict, layer="layer", top="top", bot="botm", return_reindexer=False
 ):

nlmod/gwf/gwf.py

@@ -335,6 +335,7 @@ def ghb(
     da_name=None,
     pname="ghb",
     auxiliary=None,
+    layer=None,
     **kwargs,
 ):
     """get general head boundary from model dataset.
@@ -359,6 +360,9 @@ def ghb(
         package name
     auxiliary : str or list of str
         name(s) of data arrays to include as auxiliary data to reclist
+    layer : int or None
+        The layer in which the boundary is added. It is added to the first active layer
+        when layer is None. The default is None.
 
     Raises
     ------
@@ -384,14 +388,15 @@ def ghb(
     mask_arr = _get_value_from_ds_datavar(ds, "cond", cond, return_da=True)
     mask = mask_arr != 0
 
+    first_active_layer = layer is None
     ghb_rec = grid.da_to_reclist(
         ds,
         mask,
         col1=bhead,
         col2=cond,
-        first_active_layer=True,
+        layer=layer,
+        first_active_layer=first_active_layer,
         only_active_cells=False,
-        layer=0,
         aux=auxiliary,
     )
 
@@ -449,6 +454,9 @@ def drn(
         this is deprecated, name of the drn files in the model dataset
     pname : str, optional
         package name
+    layer : int or None
+        The layer in which the boundary is added. It is added to the first active layer
+        when layer is None. The default is None.
 
     Returns
     -------
@@ -470,15 +478,14 @@ def drn(
     mask = mask_arr != 0
 
     first_active_layer = layer is None
-
     drn_rec = grid.da_to_reclist(
         ds,
         mask=mask,
         col1=elev,
         col2=cond,
+        layer=layer,
         first_active_layer=first_active_layer,
         only_active_cells=False,
-        layer=layer,
     )
 
     if len(drn_rec) > 0:
@@ -599,7 +606,14 @@ def sto(
 
 
 def chd(
-    ds, gwf, mask="chd_mask", head="chd_head", pname="chd", auxiliary=None, **kwargs
+    ds,
+    gwf,
+    mask="chd_mask",
+    head="chd_head",
+    pname="chd",
+    auxiliary=None,
+    layer=0,
+    **kwargs,
 ):
     """get constant head boundary at the model's edges from the model dataset.
 
@@ -619,6 +633,9 @@ def chd(
         package name
     auxiliary : str or list of str
         name(s) of data arrays to include as auxiliary data to reclist
+    layer : int or None
+        The layer in which the boundary is added. It is added to the first active layer
+        when layer is None. The default is 0.
     chd : str, optional
         deprecated, the new argument is 'mask'
 
@@ -640,7 +657,15 @@ def chd(
     mask = maskarr != 0
 
     # get the stress_period_data
-    chd_rec = grid.da_to_reclist(ds, mask, col1=head, aux=auxiliary)
+    first_active_layer = layer is None
+    chd_rec = grid.da_to_reclist(
+        ds,
+        mask,
+        col1=head,
+        layer=layer,
+        aux=auxiliary,
+        first_active_layer=first_active_layer,
+    )
 
     chd = flopy.mf6.ModflowGwfchd(
         gwf,
@@ -653,7 +678,7 @@ def chd(
     )
     if (auxiliary is not None) and (ds.transport == 1):
         logger.info("-> adding CHD to SSM sources list")
-        ssm_sources = ds.attrs["ssm_sources"]
+        ssm_sources = list(ds.attrs["ssm_sources"])
         if chd.package_name not in ssm_sources:
             ssm_sources += [chd.package_name]
             ds.attrs["ssm_sources"] = ssm_sources

nlmod/gwf/recharge.py

@@ -123,7 +123,10 @@ def ds_to_evt(gwf, ds, pname="evt", nseg=1, surface=None, depth=None, **kwargs):
 
     # get stress period data
     if ds.time.steady_state:
-        mask = ds["evaporation"] != 0
+        if "time" in ds["evaporation"].dims:
+            mask = ds["evaporation"].isel(time=0) != 0
+        else:
+            mask = ds["evaporation"] != 0
         rate = "evaporation"
     else:
         evt_name_arr, evt_unique_dic = _get_unique_series(ds, "evaporation", pname)
@@ -136,7 +139,7 @@ def ds_to_evt(gwf, ds, pname="evt", nseg=1, surface=None, depth=None, **kwargs):
         ds,
         mask,
         col1=surface,
-        col2=rate,
+        col2=ds[rate].squeeze(),
         col3=depth,
         first_active_layer=True,
         only_active_cells=False,

nlmod/gwf/wells.py

@@ -1,7 +1,11 @@
+import logging
+
 import flopy as fp
 import numpy as np
 from tqdm import tqdm
 
+logger = logging.getLogger(__name__)
+
 
 def wel_from_df(
     df,
@@ -23,7 +27,7 @@ def wel_from_df(
             cid1 = gwf.modelgrid.intersect(irow[x], irow[y], irow[top], forgive=False)
             cid2 = gwf.modelgrid.intersect(irow[x], irow[y], irow[botm], forgive=False)
         except Exception:
-            print(
+            logger.warning(
                 f"Warning! well {wnam} outside of model domain! ({irow[x]}, {irow[y]})"
             )
             continue
@@ -36,20 +40,23 @@ def wel_from_df(
             idomain_mask = gwf.modelgrid.idomain[kt : kb + 1, i, j] > 0
         # mask only active cells
         wlayers = np.arange(kt, kb + 1)[idomain_mask]
-        for k in wlayers:
-            if len(cid1) == 2:
-                wdata = [(k, icell2d), irow[Q] / len(wlayers)]
-            elif len(cid1) == 3:
-                wdata = [(k, i, j), irow[Q] / len(wlayers)]
-
-            if aux is not None:
-                if not isinstance(aux, list):
-                    aux = [aux]
-                for iaux in aux:
-                    wdata.append(irow[iaux])
-            if boundnames is not None:
-                wdata.append(irow[boundnames])
-            well_lrcd.append(wdata)
+        if len(wlayers) > 0:
+            for k in wlayers:
+                if len(cid1) == 2:
+                    wdata = [(k, icell2d), irow[Q] / len(wlayers)]
+                elif len(cid1) == 3:
+                    wdata = [(k, i, j), irow[Q] / len(wlayers)]
+
+                if aux is not None:
+                    if not isinstance(aux, list):
+                        aux = [aux]
+                    for iaux in aux:
+                        wdata.append(irow[iaux])
+                if boundnames is not None:
+                    wdata.append(irow[boundnames])
+                well_lrcd.append(wdata)
+        else:
+            logger.warning(f"Well '{wnam}' not added, no active screened layers.")
 
     wel_spd = {0: well_lrcd}
 

nlmod/gwt/gwt.py

@@ -168,6 +168,10 @@ def ssm(ds, gwt, sources=None, **kwargs):
     if build_tuples and sources is not None:
         sources = [(ipkg, "AUX", "CONCENTRATION") for ipkg in sources]
 
+    if len(sources) == 0:
+        logger.error("No sources to add to gwt model!")
+        raise ValueError("No sources to add to gwt model!")
+
     ssm = flopy.mf6.ModflowGwtssm(gwt, sources=sources, **kwargs)
     return ssm
 

nlmod/gwt/output.py

@@ -171,6 +171,8 @@ def freshwater_head(ds, hp, conc, denseref=None, drhodc=None):
     if "z" not in ds:
         if "thickness" not in ds:
             thickness = calculate_thickness(ds)
+        else:
+            thickness = ds.thickness
         z = ds["botm"] + thickness / 2.0
     else:
         z = ds["z"]

nlmod/modpath/modpath.py

@@ -436,15 +436,17 @@ def pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5):
     return pg
 
 
-def sim(mpf, pg, direction="backward", gwf=None, ref_time=None, stoptime=None):
+def sim(
+    mpf, particlegroups, direction="backward", gwf=None, ref_time=None, stoptime=None
+):
     """Create a modpath backward simulation from a particle group.
 
     Parameters
     ----------
     mpf : flopy.modpath.mp7.Modpath7
         modpath object.
-    pg : flopy.modpath.mp7particlegroup.ParticleGroupNodeTemplate
-        Particle group.
+    particlegroups : ParticleGroup or list of ParticleGroups
+        One or more particle groups.
     gwf : flopy.mf6.mfmodel.MFModel or None, optional
         Groundwater flow model. Only used if ref_time is not None. Default is
         None
@@ -484,7 +486,7 @@ def sim(mpf, pg, direction="backward", gwf=None, ref_time=None, stoptime=None):
         referencetime=ref_time,
         stoptimeoption=stoptimeoption,
         stoptime=stoptime,
-        particlegroups=pg,
+        particlegroups=particlegroups,
     )
 
     return mpsim