Hi Alessandro, Hi Zeno!
According to my experiences, you generally have to distinguish betweem the numerical approaches of the BHE, Al-Khoury and Eskilson/Claesson. But before that, it is crucial to regard different model settings, f.e. such as using a full confined vs. Richards' approach, so each model will have a unique time-stepping behavior. Thus it is only possible to have a 'rule of thumb' - mostly based on user experience - for the stepping.
For the Eskilson/Claesson approach (analytical BHE calculation), my recommendation is:
You may use AB/TR schemes when the flow setup is easy (steady flow, confined or phreatic approach) or switch to FE/BE when dealing with more complicated flow settings directly. Do test runs and test for constraining the time-step by each the growth (typically values between 1.1 and 2) and the max. time-step size (the default 10 may be a good starting point, but I had models put down to 0.5 days). This setup should converge already with the first iteration, so it is acceptable to reduce the number of iterations from 12 to 1. In order to compensate, choose a smaller error criterion of 1e-5.
Before you run the model, create a nodal selection on the BHE or group of BHE's and activate the charting mode for the BC's. When running the model, check the amount of input power over time, vs. the period heat budget of the BHE in the chart. These values should be almost identical. If not, try reducing the time-steps by constraining the step size, or by using a smaller error criterion, or even by using the stricter L1 or max. error norm. When you apply an input of temperature to the BHE's, I would suggest to compare different BHE heat period budgets for different time-stepping sizes. When you run with a smaller time-step and the value does not change anymore, the convergence is fine. This is the identical strategy to running models with different mesh resolution and check for mesh convergence, only the resolution is a temporal one.
For the Al-Khoury (full numerical approach),
I would suggest to directly use a FE/BE scheme. A fine vertical resolution of the model is also crucial. At least 1m for the slice distances, according to my experience. It is not ok to run this with a single iteration, so the number of iterations should be set as default 12. I usually apply a max. error norm, this already ensures a small stepping. Smaller error criterion and constraining time-steps by growth or size are optional, but I would recommend to use a small growth of 1.1 to 1.2 and also a 1e-5 error criterion. Also do test-runs just like with the analytical approach. Since Al-Khoury is used for short period simulations, a max time-step size is hard to estimate, it may be needed to put that even in minutes.
Upwinding schemes: I try to avoid using them if possible, since you cannot control the amount of artificial dispersion which comes into the model. This will affect the physical behaviour of the model, and I would like to keep control of that. Therefore, it is a good idea to control things by the longitudinal/transverse dispersion material parameters instead and improve discretization even before that. When none of these measures help, the upwinding can be a last resort.
I hope these suggestions provide some help in your setups. However, they do only reflect my personal working experiences, so it would be really nice if others share their experiences here as well.
Best regards,
Bastian