Modeling Open Excavations in MODFLOW: High K vs. Original K with Drain Package

[Groundwater-modeler]

Hello,

When modeling an open excavation (e.g., a 50 m deep quarry or a 3 m deep trench) in MODFLOW with drainage at the bottom (using the Drain package), is it more conceptually correct to:

1. Keep the original hydraulic conductivity (K) of the surrounding material (e.g., clay or rock) in the excavation cells, or
2. Assign a very high (K) (e.g., to simulate an air-filled void) in the excavation cells?
   The excavation has no internal resistance to flow, with water entering freely through the walls and bottom (via fractures or pores) and being removed by drainage/pumping. How does this choice affect inflow and groundwater pressure in underlying layers?

To add some further content I have noticed that it matters when modeling, for instance, an excavation in clay where the clay is discretized vertically into multiple layers if I represent the cell containing the excavation with an arbitrarily high k-value even if the conductance is very high vs leaving the cells k-value untouched.
This puzzles me a bit since I would think the conductance would handle this on its own…

Sincerely
David

[emorway-usgs]

Hello @Groundwater-modeler, I'm not experienced when it comes to simulating an expanding mine in a MODFLOW grid but looking at Chapter Six of the Supplemental Technical Information Guide (distributed with MODFLOW), you might want to check out Chapter 6 (1st paragraph shown below). Sounds like the time-varying hydraulic conductivity and time-varying storage packages were built with your use case in mind (specifically the 50 m deep quarry?).

For the 2nd of your two examples, you might consider representing the 3m trench with the SFR package.

[Groundwater-modeler]

Hi @emorway-usgs,

Thanks for taking the time to answer my question. While I agree with the use of TVK for simulating progressive mining, it doesn't fully capture the concept I'm trying to convey. I also doubt that SFR will assist me with the trench excavation in this case since the excavation bottom sometimes is below and sometimes above the groundwater table.

I'm struggling to formulate the issue clearly enough to communicate my message effectively.
What I'm trying to achieve is to use a BC (or some other relevant concept) to simulate excavation without necessarily having to represent the excavation geometry in the grid. This is mainly because the excavation depth in a typical project often covers a large area (railroad, road construction, etc.). The excavation depth can range from "negligible" to moderate to large, meaning that sometimes the depth of excavation is just half a meter and sometimes it might be 10 meters. The excavations can be anywhere in this range, and the soil type varies as well.

The challenge is to capture the system's behavior in a conceptually reasonable way. My approach involves implementing drains at whatever depth the excavation bottom is for each location. However, this approach presents a few conceptual difficulties, such as:

1. Depending on the excavation depth, sometimes the drain will be placed in the first layer and sometimes in a deeper layer.
2. There appears to be a difference in how the head propagates out from the excavation cells depending on the hydraulic conductivity of the cell where the excavation drain resides.
3. Achieving a geometrically accurate excavation model would require several layers in the upper domain of the model. This would lead to a large number of cells in the model overall (not necessarily a problem, but it will influence runtime). Additionally, it would cause the geological units to stretch across multiple layers rather than being defined as layers themselves. For instance, if I want to examine drawdown in a soil type with relatively higher conductivity under a layer of clay, it's easier if that material unit is defined as a layer. However, if I need to increase vertical discretization in the model's upper domain, the unit will have to be placed in multiple layers.

I have set up two identical test models (2D - cross-section) with 8 layers. The first three layers are clay (low conductivity) from 0 to -7 m, followed by a more conductive material from -7 to -10 m, and then a couple of layers representing bedrock from -10 to -50 m.
In the first model, the excavation pit is represented only by a drain in the top layer (0 to -3 m) with a conductance of 1000 m²/s (arbitrarily high to ensure complete drainage). In the second model, everything is identical except that the excavation is represented with high conductivity rather than the original material.

There is a difference in how the heads propagate out from the pit in these two cases. I suppose it has something to do with how MODFLOW calculates the harmonic mean of heads between cells, but I don't quite understand it.

Best regards,
David

[emorway-usgs]

@Groundwater-modeler, regarding the following:

Depending on the excavation depth, sometimes the drain will be placed in the first layer and sometimes in a deeper layer.

I think it's possible that the input format to the DRN package will support you? Recall that you can specify new drains - and effectively turn-off other drains specified previously - by taking advantage of the PERIOD blocks. So could you do something like the following:

BEGIN PERIOD 1
1 1 1 100.0 1000.0
END PERIOD 1

BEGIN PERIOD 2
1 1 1 100.0 0.0
2 1 1 90.0 1000.0
END PERIOD

The idea being that there is a DRN in the upper most layer at a specified depth of 100 in the first stress period that is then effectively shut-off in the second stress period by resetting the conductance to 0.0. In addition, a "new" DRN is listed for layer 2 at a depth of 90.0 in the second stress period.

You may be able to play a similar game in SFR by taking advantage of the BEDK keyword available with sfrsetting. You would need to set up all the various trenches in the PACKAGEDATA block and then activate whichever one you want by varying bedk to be 0.0 (off) or >0.0 (on).

[Groundwater-modeler]

I’m not sure I follow. I dont really have any issues with stress periods as the problem is not related to progressing depth.

The problem is rather that I cant find a BC or way around the issue where an excavation has a geometri that doesnt fit the grid cells and as such cant be represented as void space with tvk. In trying to express conductance in a such that inflow and head in surrounding cells are reflected accurately.
It seems as the drain cells internal resistance is playing a problematic role here as it limits the flow through the cell regardless of the conductance of my drain.

I have not found any other useful BCs with an off Switch once head falls below the level of my bcs head.
I dont want the bc to start feeding water to my model or create short circuit effects where two ajacent cells have different heads leading to one cell feeding water to the other.

[emorway-usgs]

I'm afraid I may not have satisfactory solution for you @Groundwater-modeler.

Regarding your comment about "internal resistance playing a problematic role": the connection between the groundwater cell and the boundary condition is governed by the conductance specified in the drain package as well as the properties and distance between the cell's computational node and the position (elev) of the boundary condition. If what you mean by this comment is that you've raised the conductance in the drain package to a relatively high value and further increases result in little (if any) response in the heads, it may be because it is the cell properties - and not the specified conductance of the drain - that are controlling the interaction with the boundary condition.

As best I can tell from the description of your problem, I do think the drain package is what you want. When the calculated head in a cell drops below the drain elevation, the drain will no then 'introduce' water into the cell, as you point. It was hard to tell from your description, but I think what you want to do is set the drain elevation to be equal to the bottom of the excavation depth. If you then need to vary properties above and below the excavation depth, I'm afraid you may need to consider refining the vertical discretization a bit further (but here again I may not be understanding what you are seeking to represent in the model).

Sorry, but I think this is the best I can offer. Others may have a better grasp of you're trying to accomplish.

[Groundwater-modeler]

Thank you. I believe you've accurately captured the essence of my problem with your description.

I feared that the solution might necessitate additional vertical discretization. Although I suspect that implementing it not always a viable option. It would be ideal to have a boundary condition (perhaps a new type) that functions similarly to a drain but is more independent of the grid. Perhpas something that could be to adjust the cells' horizontal and vertical conductivity based on the volume the BC occupy within the cell.

I often find it challenging, if not nearly impossible, to implement more vertical discretization. If anyone has good suggestions, I would appreciate them. When working on soil excavation projects with moderate depth in relatively large models, discretizing the soil layers almost always results in extremely thin layers in areas where the bedrock approaches the surface. This is especially true if the soil is defined as the space between the bedrock surface and the ground surface. This in turn almost certainly lead to model instability and convergence issues (regardless if NWT is used or not).

Alternatively, one could experiment with different layering concepts, such as defining a fixed thickness for the upper layers and assigning properties to the cells based on their position relative to the surfaces. However, I feel that this method diminishes the ability to observe the effects within a specific material unit if the layers do not correspond to some type of geological material.