I'm working with a Modflow-based software and I am trying to model a GWHP system and its influence on GW temperatures. The aquifer is supposed to be homogeneous
Peclet and Courant numbers will always give you a good hint. In addition to that... I always make a first discretization, with equal length layers covering the aquifer thickness. Then check stability. If needed, then refine. If, for some reason, you also want to model the unsaturated zone, set some extra layers there. You can download a few recent studies we have made on heat transport for shallow geothermal in our LinkedIn storefront (see the link).