
Only putting global radiation at the surface will very likely overestimate the heat inflow. On the one hand, you'd probably need to consider albedo. On the other hand, there will also be radiation back into the atmosphere, plus conductive heat flow, ...

If there is only outflow on the bottom, you do not need any masstransport BC. FEFLOW will allow a free advective outflow of mass.

There's currently no way to directly add to the model. You can create an additional mesh for the neighboring area (making sure that at the common border with the existing mesh the nodes are identical), and then import both meshes (2D) into one (for example, from shapefiles). All the model setup, however, you have to repeat from there.

The PETSc AMG solver has to be selected when installing FEFLOW  just run the installer again in 'Modify' mode and select the correspondig package for installation. You may also wish to try the SAMG version 2019 to compare (you can switch in Tools  Global Settings  Parallel Computing). I personally have had trouble with SAMG 2020 when 7.5 came out, since then I am using SAMG 2019 (and haven't tried 2020 again yet).

You might want to have a look at https://www.youtube.com/watch?v=7VcGqRPvxdo for an introduction to unstructured meshes in FEFLOW. Note that handling models with unstructured meshes in general is more complicated than using the layered option.

Hi Zeb,
You can set a fixed min and max for the temperature in the properties (context menu) of Temperature in the View Components panel.
Cheers,
Peter

Stefan,
I wouldn't bee too worried about jumps in layer thickness  as long as the solver runs without any hiccups in the internal iteration such jumps should be OK. The RWD you use leads a a pseudosaturation for each partially dry element. Depending on the overall thickness of the element, the RWD leads to a corresponding residual saturation and therefore residual hydraulic conducitivity (BTW, I wish FEFLOW would use such rather than RWD). So whether 9 cm of RWD (or 4.5cm) are acceptable highly depends on the overall thickness of the layers and the sensitivity of your required model outcomes on RWD/residual conductivity. Comparing to the height of a tunnel, though, also 9 cm seem to be relatively small and you might not have to worry about the change you need to make in RWD. Why not run the model with 10 and 30 cm RWD and check the difference?

In a 'phreatic' model, FEFLOW reduces conductivity with saturation of each element. With thick layers (in you twolayer model), saturation for the large elements will be relatively large, thus conductivity reduced only moderately (but on the entire layer!!). The more you refine the model, the more likely it will be that you have dry layers, where minimum conducitivity is determined by residual water depth. Thus with thinner layers the conductivity contrasts in the model become larger, and model convergence will be harder to obtain  especially with groundwater recharge being applied to the uppermost 'dry' layer with low conductivity. However, at the same time you will also approximate the water level better, getting away from 'randomly' scaling down conductivity in the entire aquifer depending on the position of the water table within the layer.
So what you are looking for is a good balance of vertical discretization and a suitable residual water depth (as long as you use 'phreatic' mode) that allows recharge to pass through the layer(s) above the water table. This behaviour of the model shows one of the most important difficulties in groundwater modelling: Often the role of the UZ zone is not negligible, but vertical discretization does not allow for real unsaturated flow. Some ways to deal with it may be:
 Find a good balance of discretization and residual water depth
 Apply recharge in the aquifer rather than on top (only possible by using the source/sink parameter, requiring a change of the quantity by taking into account the element thickness)
 Use a simplified unsaturated model rather than 'phreatic' (but this requiring more experience on suitable parameters)

Ah  now I fully understand. Indeed this is not possible with the current design of the OpenLoop plugin.

Indeed the main difference is the variability of the transmissivity, which requires the socalled outer iteration in each time step, iterating over the step until the solution converges. In a confined simulation, this outer iteration is not necessary, the equation system can be directly solved (or with an iterative solver). So the time step length you end up with depends on the (dimensionless) error criterion you have specified plus the actual variation of the solution compared to the criterion. So if your criterion is appropriate for the application, the solution varies a lot  which could point towards an oscillation. Oscillations can occur, for example, when wells with a high abstraction rate are running in relatively lowpermeable environments, high recharge rates are applied to phreatic elements with low relative conductivity, etc.