specified flux boundaries

[mnfienen]

Unless I'm missing something, it looks like specified flux boundary conditions are only supported for inset TMR models in MF2005/NWT. However, may be simple refactoring to support MF6 as well. Like to discuss before trying to get it going.

[ntdosch]

@mnfienen and @aleaf

• would like to clean up some logic in use of file names in yml file as they are passed to Tmr
• consider consolidating where constant head and constant flux are constructed. Constant heads setup in setup_perimeter_boundary and constant fluxes are setup in setup_wel_data. Constant heads are setup in more general model building parts of the package, constant fluxes are setup within well building. Does it make sense to make these more consistent?
• will it be necessary to navigate which cell faces supply fluxes from parent to inset?
• at least for constant flux, need to navigate to Kh information based on model type. Currently an assert statement looks for UPW only.

Propose to build around existing tests to implement these changes.

[aleaf]

Hey @ntdosch and @mnfienen, I've been slowly chipping away at a refactoring of how perimeter boundaries are done, in the tmr.TmrNew class on the feature_gen_perimeter branch. These changes apply to specified head or flux boundaries, although currently only specified head is implemented. From a user interface/input standpoint, the biggest change is that perimeter boundaries are now specified within their respective package blocks, e.g.

chd:
  perimeter_boundary:
    shapefile: 'shellmound/tmr_parent/gis/irregular_boundary.shp'
    parent_head_file: 'shellmound/tmr_parent/shellmound.hds'

for specified heads; and explicit mapping of parent to inset layers is no longer needed. Instead, the parent head solution is linearly interpolated in 3D to the inset model boundary cell locations (using a method similar to scipy.interpolate.griddata), at each stress period.

The options for specifying where perimeter boundary should be applied have also been expanded. There are basically 3 options, which are implemented by TmrNew.get_inset_boundary_cells():

1. no specification of where the perimeter boundary should be applied (e.g. a shapefile) and by_layer=False: perimeter BC cells are applied to active cells that coincide with edge of the maximum areal footprint of the active model area. In places where the edge of the active area is inside of the max active footprint, no perimeter cells are applied.
2. no specification of where the perimeter boundary should be applied and by_layer=True: same as above option but the active footprint is defined by layer from the idomain array. This is probably not a great option in many cases, as it would often lead to boundary cells being included in the model interior (along layer pinch-outs, for example).
3. specification of perimeter boundary cells with a shapefile. The locations of cells can be explicitly specified this way, but they still must coincide with the edge of the active extent in each layer (modflow-setup won't put perimeter cells in the model interior). (Open) Polyline or Polygon shapefiles can be used; in either case a buffer is used to align the supplied features with the active area edge, which is determined using the Sobel edge detection filter in Scipy.

TmrNew.get_inset_boundary_values() is what gets the actual values at each stress period. Right now it has a conditional that if self.boundary_type == 'head', it gets the heads. I'm thinking we could do the same for the fluxes. The only difference would be some additional math to translate parent cell budget fluxes to volumetric quantities at each inset boundary cell. Something like:

• read the left/right and front/back fluxes from the parent cell budget output
  • top and bottom can be neglected if we assume the perimeter to vertically extend through all layers, with an absences being in cells that are inactive)
  • or 3D TMR could be supported by including the top/bottom fluxes
• convert those to a specific discharge vector (independent of area)
• rotate the specific discharge vector to the orientation of the inset model
• convert the rotated specific discharge vector back to fluxes thru the inset perimeter cell faces
• TmrNew.get_inset_boundary_cells() already returns the cell face of the boundary for each perimeter cell; use this to convert the perimeter cell face fluxes to a volumetric quantity added or removed from each perimeter cell.
• within MF6model and MFnwtModel, add the same logic (as in the CHD package) to the WEL package
• the tests for specified head vs specified flux could follow the same logic as well
  • test_mf6_tmr_shellmound.py::test_irregular_perimeter_boundary looks at the interpolation
  • some of the tests in test_tmr.py test some basic mechanical issues; some of these are a little half-baked/unsatisfying
  • the test_bcs.py::test_remove_inactive_bcs test uses the basic_model_instance fixture (3 different models) to test the setup_chd() for MODFLOW-NWT and MODFLOW 6. Would be good to do something similar for setup_wel() in the context of perimeter boundaries.

Does this make any sense? Am I missing something?

[mnfienen]

@aleaf this is great - thanks. You are farther along than we saw in the develop branch.

I do wonder if we will want to add a rotated model to the basic_model_instance test to be sure all the vector rotation from parent to inset works out. Seems we will want to operate on all that in model coordinates only....something to consider.

My only other thought at the moment is whether it makes sense to invoke the CHD or WEL package automatically if these inset boundaries are requested from the parent. At least in the develop branch, CHDs are generated if the perimeter_boundary_type option contains "head" in the configuration YML file but one must invoke the WEL package separately if perimeter_boundary_type contains "flux". Perhaps your refactoring to date has already handled that.

[aleaf]

@mnfienen @ntdosch FYI, just updated the develop branch with the contents of feature_gen_perimeter. Also added a doc section for explaining some of the broader concepts in modflow-setup, and within it, the above info on perimeter boundaries. I think we should keep developing the specified flux boundaries in the feature branch but regardless, would be good to rebase so that you're working with the latest version.

[mnfienen]

@aleaf - by rebase, do you mean we should rebase the feature branch with develop?

[aleaf]

@mnfienen, yes, you want the feature branch to have the same history as develop, plus the new commits you've added.

git checkout develop
git fetch
git rebase <upstream>/develop
git checkout feature_gen_perimeter
git rebase develop

(where points to aleaf/modflow-setup)

[mnfienen]

great - that's what I thought. Thanks!