Hello,
I’m new at using OGS and I’ve been trying to do some models with HT processes.
Right now i’m trying to reproduce the first series of models from http://dx.doi.org/10.1016/j.jvolgeores.2013.02.008
Models’ description :
 2D cross section : rectangular box (dimensions : 24 x 6 km)
 Domain separated in 2 units : upper unit is impervious sediments and lower unit is granit (more permeable).
 Thermal boundary conditions correspond to a fixed temperature (10 °C) at the surface and a fixed heat flow at the base of the model (90 mW/m2 at 6 km of depth). Lateral boundaries are thermally insulating.
 A fixed atmospheric pressure (1e+05 Pa) is imposed at the surface while a zero flux condition is imposed on all other boundaries.
 Initial thermal field (a steadystate conductive regime) evolves towards a steadystate convective regime in a few Myrs
 The initial pressure field increases linearly with depth.
I can’t upload files as i’m a new user but here is a filesender link where you can find my project file and the associated meshes :
https://filesender.renater.fr/?s=download&token=4703a9abfde3421e80aa225de0c7b032
So if you don’t mind I had some questions :

When running the simulation, i can’t obtain a result because it doesn’t converge.
I found that the problem comes with the function viscosity that i put in the prj file because when i try to put a constant viscosity of 1e03 it works just fine. So I’m wondering what is the issue with this function ? It is the same as in the article and i checked multiple times the expression of the derivative i put in. 
At the same time I have issues with the units, because from what i understand thermal dirichlet conditions are in kelvins ? But for example in the benchmark Constant Viscosity (sorry i can’t put a third link as i’m a new user) they put temperature in celsius ? So right now my conditions are in celsius, so do i need to change them into kelvin.

And again at the same time, the density and viscosity functions that i put in my prj file require a temperature in celsius. So maybe this is the issue ? That i need to transform this functions that require celsius into functions that require kelvins ? But i’m not really sure how to do that.

I was also wondering if dvalues is always necessary when working with functions ?

In most benchmark that i saw with HT processes, the thermal dispersivity is always set at 0, in what cases can i see it change ?

When working with HT processes, I’m not really sure what to use in the <nonlinear_solvers> , <linear_solvers> and <convergence_criterion> blocks ? For example can you comment on what i used in my prj file please ? Is it efficient ? Or how can i improve my results ? For example if i change something will it also improve the total running time ? I plan on working on 3D models later, will it be enough ?
I don’t know if it makes a lot of sense,
But thanks in advance,
RoseNelly