The basic equations of Darcy flow coupled to thermal transport are:
$$
\nabla \cdot \vec{u} = 0
\\
\vec{u} = -\frac{\mathbf{K}}{\mu} (\nabla p - \rho \vec{g})
\\
\sigma \frac{\partial T}{\partial t} + \vec{u}\cdot\nabla T - k\nabla^2 T = 0
$$
where $$\vec{u}$$ is the fluid velocity, $$\mathbf{K}$$ is the permeability tensor, $$\mu$$ is fluid viscosity, $$p$$ is the pressure, $$\rho$$ is the fluid density, $$\vec{g}$$ is the gravity vector, $$T$$ is the temperature, $$\sigma$$ is the ratio of the heat capacity of the medium to the heat capacity of the fluid, and $$k$$ is the effective thermal diffusivity of the saturated medium.

We shall henceforth assume that $$\vec{g} = \vec{0}$$.  If gravity is important, it is possible to include buoyancy effects in the fluid density by introducing Boussinesq's "slightly compressible" approximation, $$\rho = \rho_0 (1 - \alpha(T - T_0))$$.  Using this assumption, taking the divergence of the second equation, and imposing the divergence-free condition leads to the following system of two equations in the unknowns $$p$$ and $$T$$:
$$
-\nabla \cdot \frac{\mathbf{K}}{\mu} \nabla p  = 0
\\
\sigma \frac{\partial T}{\partial t} + \vec{u}\cdot\nabla T - k\nabla^2 T = 0
$$
where $$\vec{u} = -\frac{\mathbf{K}}{\mu} \nabla p$$.  That is, the pressure satisfies Laplace's equation (with possibly non-constant coefficients) and the temperature satisfies the linear convection-diffusion equation.