In this paper we report on Large Eddy Simulations of natural convection in water pools with evaporation across the free surface and at the hard turbulence regime. The free surface is approximated as a free-slip top boundary. The loss of water is estimated via a dynamic and inhomogeneous evaporation model. Also, the descent of the free surface is accounted for by regularly reducing the computational domain and applying a remeshing procedure. We present results for 4 different Rayleigh numbers, ranging from \(\boldsymbol{Ra = 1.35 \times {10^8}}\) to \(\boldsymbol{Ra{ = 10^{10}}}\). Our simulations predict a slow decrease of the free-surface temperature and evaporation rate over time. This may be attributed to the descent of the free surface due to evaporation which tends to reduce the intensity of turbulent motions. On the other hand, the flow structure remains the same throughout the duration of the simulations. More specifically, the flow is organized in a large scale circulation aligned in a diagonal plane with smaller convective rolls near the corners of the domain. Also, the absence of a top hydrodynamic boundary layer enhances turbulent mixing and convective heat transfer near the free surface. This enhancement is manifested by a shift of the profile of the mean surface temperature towards the upper part of the domain, with the shift becoming more pronounced as the turbulence intensity increases. Herein we also provide results for the Nusselt number \(\boldsymbol{Nu}\) and present a new \(\boldsymbol{Nu - Ra}\) scaling for convection in pools and cavities that covers a large range of turbulence intensities.