Peter Schätzl Grundwassermodellierer

You're welcome!

'Scale' means scale of the transport phenomenon you're dealing with. In case of a contaminant plume, this is pretty clear... When simulating concentrations of widely distributed substances, though, it gets fuzzy.

There's no 'unconfined' mode in 2D vertical models. So to simulate unconfined conditions in such a model, you need to do a fully unsaturated simulation.

This looks a bit like BCs created on the basis of a top or bottom constraint on hydraulic head (defined in the free surface settings). This, however, I'd not expect in a fully confined model... Could you check whether it's really confined?

Indeed the way of thinking is a bit different in between MODFLOW and FEFLOW... There's no unique answer to your question, but in a general system of aquifers and aquitards I'd put the aquifer head bcs to both top and bottom of the aquifer. If there's no clear separation, but rather a system of layers of slightly differing hydraulic conductivities, I'd typically tend to use the condition of the higher conductive layer a slice belongs to. And - yes - I'd also put a BC to the bottom layer.
In all cases you should keep in mind where your BC data come from to judge where to best put them. In an ideal world, vertical discretization would be fine enough to render the influence of the slices separating the geological units negligible - but the world (or the model) is never ideal as we all know.

'slice faces' are kind of a legacy mode, according to what was used in version 6.x (and earlier), thus I'd recommend to use 'faces'. For the modelling result, however, it should not make any difference. 'Thickness' is the fracture aperture, used to calculate the storage properties of the fracture. There's a detailed explanantion of the different parameters (in different cases) in the FEFLOW book.

You're welcome. You just have to keep in mind that this increases k in the unsat zone in all directions.

Hard to judge... you could try - in addition to the zero pumping phase - to add a ramp to the pumping/injection time series. However, I'd typically also link this error more with spatial discretization than temporal behaviour.

I suppose that by reducing the vertical conductivity the residual conductivity of the potentially dry first layer will be lower by a factor of 10, and recharge may not be able to infiltrate. This leads to a pressure build-up at the surface, with periodic break-throughs once saturation is high enough to increase k to a degree that lets the water drain down. If you are using the default residual water level of 1 mm for the phreatic layer(s), I suggest to increase this reasonably to increase residual saturation, hence residual k in dry layers.

The actual process of applying parameters to discrete features depends on the FEFLOW version you are using (6.x vs. 7.x). The same applies for the question of what to select in order to generate 2D discrete features: in version 6.x you'd use 'slice faces' (faces, not edges, and on a slice, not in a layer!), in 7.x you should use 'faces' rather than the legacy 'slice faces'.
It will not be exactly the same as using a separate layer, for two reasons:
1. The layer will separate the underlying and overlaying layers from each other, the discrete feature will not, as the nodes are shared.
2. The spatial discretization in vertical direction will be different, with more nodes vertically with the separate layers.
It highly depends on the parameterization of the different hydrogeological features whether the differences are significant or not. The more the flow is in horzonatal direction, and the higher the conductivity of the fracture compared to the layers above and below, the closer the results will be. Also stability will highly depend on the parameters - no general answer to that one.