Liquid cooling systems have better heat dissipation capabilities than air based ones, and are expected to become a standard choice in future data centers, due to the ever increasing power density and heat rejection needs of the compute infrastructure. A convenient side-effect of implementing liquid cooling is that it facilitates the efficient recovery of the heat waste. However designing and managing these heat recovery infrastructures benefit from having control-oriented models that can accurately describe how different operating conditions of the to-be-cooled heat sources will affect the thermal status of the coolant. The aim of this manuscript is to derive control-oriented models of liquid immersion cooling systems, i.e., systems where the compute infrastructure is immersed in a vessel filled with dielectric fluid. More specifically we derive, starting from physical interpretations, a general lumped-parameters gray box dynamical model that has-as inputs-the electrical consumption of the heat sources and the working point of the heat recovery system, and has-as outputs-the temperature distribution of the coolant in the most relevant points of the system. Beyond proposing this modelling methodology we also validate the generalization capabilities of the obtainable models. In specific, we test the achievable statistical performances in a field case, plus compare with the ones of classical black box system identification strategies. We thus report that in the considered field case our gray box model reached a fit index of 91.08% when simulating test sets, while the best black box model we have been able to identify reached (on the same test sets) fit indexes of only 72.56%.