One of the essential issues for sustainable development is obtaining energy from the ground medium. An important aspect of this issue is the parameter of, among other things, the thermal conductivity of the soil. There is still a need to find accurate models describing the thermal conductivity of soil at various degrees of saturation. The authors are seeking a model in the framework of analytical micromechanics that describes the relationship between the shape of ellipsoidal inclusions θ of water and air in pore space and the variable degree of saturation. Using the Mori–Tanaka approximation scheme and the simulated annealing optimization method, we first attempt to find the best pair of quadratic functions θw(Sr) and θa(Sr) (for water and air inclusions, respectively). Then, we verify the possibility of applying only one distribution of θ values, identical for both water and air. These efforts are verified using thermal conductivities of selected real soil samples.