Flow nets are one of the effective tools for illustrating groundwater flow rate and direction. Flow nets of a thermal groundwater system under different reference temperatures are examined in this paper. A 3D mathematical model describing water flow and heat transport is established. The stream function and the equivalent hydraulic head are used to describe the flow field in terms of mass flux. The partial differential governing equations are given by the mass conservation equation of the fluid phase and by the advection–diffusion transport equations of heat. They are solved with the standard Galerkin finite element method. Influence of variation in temperature on thermal groundwater flow nets is illustrated with a simplified 2D vertical steady-state mathematical model. The results show that the equipotential lines are not always orthogonal to the streamlines in a geothermal system of a confined aquifer when the heat flow from below is present and when the temperature increases, the equipotential lines are clearly curved and the skew intersection angles between the streamlines and the equipotential lines decrease. For the node of coordinate of (5000, 0) in the aquifer, when the temperature difference between the top boundary and the bottom boundary is 4 and 10 °C, the skew intersection angle is 85.54° and 78.73°, respectively. When the hydraulic gradients are low, the equipotential lines are clearly curved. However, when the hydraulic gradients are high, the bending of the equipotential lines is not obvious. When the hydraulic gradient increases, the skew intersection angles change relatively little. For the same node when the hydraulic head of the left boundary is 600, 700, 800 and 1,000 m, the skew intersection angle is 78.73°, 84.3°, 86.19° and 87.71°, respectively. When the reference temperature increases, the hydraulic heads and the skew intersection angles increase.