Igor Pavlovskii Dalhousie University Post-Doctoral Fellow

Your saturation data indicates that the aquifer's hydraulic conductivity is not sufficient to sustain such pumping rate without pumping well running dry (which not surprising given that 5E4 m³/day is extremely high pumping rate).
You can check how much pumping your aquifer can sustain by setting a head constraint equal to elevation for your well nodes. This would limit pumping to the rate at which wells remain in the saturated zone (the rate value can checked by looking into water budget for the well nodes). In terms of setup, I would remove well nodes at the top of the aquifer - these have very limited available head are much more likely to go dry.
I would also check some analytical solutions for partially penetrating wells - this should give you a ballpark estimate of pumping rates feasible in this setup.

Are you importing raster directly into FeFlow?
It is generally much more effective to sample elevations at points locations externally
See [b]Reply # 2[/b] at [i]https://52.136.242.250/index.php?topic=2127.0[/i]

Right-click on Head BC in the "Data" section => select "add parameter" => select "max flow rate constraint"
This would add a new parameter that you can assign to nodes with Head BC.
To prevent outflow in a specific node assign max flow rate constraint value of 0.

http://www.feflow.info/html/help72/feflow/09_Parameters/Boundary_Conditions/Flow/flow_boundary_conditions.html
Note different sign convention for BCs and constraints

FeFlow uses "head" to refer specifically to freshwater head. So, in case of saline waters  it increases with depth even in conditions without vertical flow (provided you have non-zero density ratio).

Equation linking freshwater and saltwater heads can be found at
http://www.feflow.info/html/help73/feflow/09_Parameters/Boundary_Conditions/Flow/hydraulic_head_bc_saltwaterhead.html

1) I see no point in using Cauchy boundary in this case. Since you are using infinitely high transfer rate, you can simply use fixed head boundary instead.
2) The easiest option to make boundary conditions one way is to impose flow rate constraint.

Flow rate constraint for fluid transfer BCs (option for fixed head BC is equivalent):
[url=http://www.feflow.info/html/help/default.htm?turl=HTMLDocuments%2Freference%2Fparameters%2Fboundaryconditions%2Fflow%2Fbcc_flowrate_head_transfer.htm]http://www.feflow.info/html/help/default.htm?turl=HTMLDocuments%2Freference%2Fparameters%2Fboundaryconditions%2Fflow%2Fbcc_flowrate_head_transfer.htm[/url]

I assume you want to setup Initial Conditions for the run (not boundary conditions):
just set up hydraulic head at all nodes (not just specific slice) to whatever [b]absolute[/b] elevation your groundwater table should be in the model.
Based on your post I think -20 m. I suspect your model top is actually at 0 m and bottom is at -200 m. You can change these through Edit => 3D layer configuration

[quote author=Jimmy Lin link=topic=22183.msg29509#msg29509 date=1620707863]
Why the recharge in gravel only can work when area below water table?
[/quote]
I was referring specifically to [b]negative[/b] recharge (i.e. water being extracted from the system)
Simply drying the surface above gravel will not force water to flow upward (neither in reality, nor in properly converging model)
In contrast, in fine sediments capillary fringe can extend metres above water table. In such cases, when drying surface is within such fringe it can source water from groundwater making negative recharge values a feasible option.

Convergence in the non-discretised case may actually be a numerical artifact (see attached figure).
Negative recharge values may work in gravel only if corresponding area is consistently below water table. Otherwise, water cannot simply go upwards beyond capillary fringe (that is only few cm in gravel). You may want to consider other BCs located at the bottom of the model to simulate area-averaged pumping.

Further stabilisation can be done by
- applying CVFEM (tickbox in Problem Class => Unsaturated Flow menu)
- switching to Modified Van Genuchten and reducing delta parameter value as described at
[url=http://www.feflow.info/uploads/media/Frind.pdf]www.feflow.info/uploads/media/Frind.pdf[/url]

Porosity and thermal conductivity are elemental properties. So, I assume you have exported elemental properties to excel and edited it there, but preserved the element numbers from the model.
The standard workflow is the following:
1) right-click on the data tab => select "add map"
2) right-click on the imported map "defaut" property  => select "define coordinate fields" => select appropriate fields from drop-down menu
3) right-click on the imported map => select "link to parameters" => use graphical user interface to link fields in the excel and model parameters (remember to select field with element number in the topology section)
4) select all elements of the model => right-click on the link in the "linked attributes" menu of the imported map => select "activate link for data assignment" => click on ✓ in the assignment tab (the same one uses to assign properties manually)

NB: FeFlow needs Excel to be installed on the computer to be able to import .xlsx files

The exported element centroids should have had both element and layer numbers
You need to preserve these fields during evaluation and export from Leapfrog.
Given that you have mentioned .xlsx, I think you may have deleted these during one of intermediate steps