I want to implement wells that screen several layers in the model domain. However, when setting up the wel package in Flopy, only one layer can be entered in the stress period data of each well, as shown in the Flopy example below:
# Create the well package
pumping_rate = -500.0
wel_sp = [[0, nrow / 2 - 1, ncol / 2 - 1, pumping_rate]]
stress_period_data = {0: wel_sp}
wel = flopy.modflow.ModflowWel(mf, stress_period_data=stress_period_data)
One solution to this problem is the implementation of a second well in the layer below and halving the pumping rate:
# Create the well package
pumping_rate = -500.0
wel_sp1 = [[0, nrow / 2 - 1, ncol / 2 - 1, pumping_rate/2]]
wel_sp2 = [[1, nrow / 2 - 1, ncol / 2 - 1, pumping_rate/2]]
stress_period_data = {0: wel_sp1, 1: wel_sp2}
wel = flopy.modflow.ModflowWel(mf, stress_period_data=stress_period_data)
This, however, is not justifiable, as soon as the layers differ in transmissivity. Is there a (built-in) way of spanning a well over several layers in Flopy?
I don't see anything wrong with your solution, except perhaps having to half the pumping rate; however, that may be necessary for your specific application.
That being said, there appears to be an implementation issue.
You mention multiple layers, but not multiple stress periods.
You are, however, implementing two distinct stress periods (well_sp1 and well_sp2).
It seems that you intend to define a single stress period with multiple cells affected by a well.
According to Flopy's documentation, it should look like this:
wel_sp = [[0, nrow / 2 - 1, ncol / 2 - 1, pumping_rate], [1, nrow / 2 - 1, ncol / 2 - 1, pumping_rate]]
I believe this is the best solution for you. If you have any doubts, you could create a multi-layer model in, say, ModelMuse and load it into Flopy to inspect the stress period data and wells to see if they match your implementation.
I hope this helps!