Which set of numerical settings are generally likely to produce the most accurate results is ambitious to assess. Each model is unique with respect to the configuration/combination of the question being investigated, the problem class, the model conceptualization, the spatiotemporal discretization and so on. Therefore I try to give a rather general answer.
The difference between the [b]divergent form[/b] and the [b]convective form[/b] of the transport equation is mainly given by the convective term. The convective term of the divergent transport equation involves a divergent expression of the velocity field and the species/heat. In contrast, the convective term of the convective transport equation takes a gradient relationship into account. Both forms are physically equivalent. However, in their discretized forms these two equations lead to different formulations of Boundary Conditions (BC), which may or may not result into distinct budgets. In a practical sense, a Mass-flux BC (Neumann) and a Mass nodal sink/source BC take advective and dispersive transport into account if you solve the system by the divergent transport equation. In contrast, if you adopt the same BC’s, but solve the system by the convective form you only take the dispersive transport component into account. Under certain circumstances the divergent form may also lead to instabilities at "free-outflow" boundaries.
In some situations where the advective driver overwhelms the diffusive/dispersive driver the objective of a bounded numerical solution cannot always be satisfied. [b]Upwinding[/b] techniques may help to get a bounded solution by stabilizing the oscillating numerical behavior by non-physical dispersion. However, Upwinding techniques have the disadvantage with respect to the accuracy, because stability does not imply accuracy. The different Upwinding options differ in the context of how numerical dispersion is introduced. For example, the Streamline Upwind method stabilizes the solution where advection is dominating. The stabilization spatially correlates with the advective trajectory, while stabilization in transversal direction is neglected, thus crosswind damping is neglected. In contrast, the Shock Capturing method takes the gradients of the species/heat into account rather than the trajectory of groundwater flow. Accordingly, Shock Capturing works in both the longitudinal and transversal direction respectively. In a practical sense, I would try to solve the problem without Upwinding. Instead, I would try to improve the model (e.g. mesh). If all model improvements do result in the solution I expect, I would possibly intend to use the “most harmless” Upwinding.
Regarding the options for the predictor-corrector scheme FEFLOW provides the possibility to impose additional time-step constraints. Additional time-step constraints are provided by a [b]growth factor between subsequent time-steps[/b] and a [b]maximum time-step size[/b]. This option is required if the automatic predictor-corrector returns too large time-steps with respect to the required temporal accuracy. This situation is rarely the case. However, in some applications additional time-step constraints may become useful. A possible application is density-dependent flow. Another application are geothermal simulations involving Borehole Heat Exchangers (BHE’s) where the inlet temperature at a specific time-step is calculated from the outlet temperature of the previous time-step (e.g. temperature difference or power).
I hope this rather general answer helps.
