Thermostatically controlled loads (TCLs) can effectively support network operation through their intrinsic flexibility and play a pivotal role in delivering cost effective decarbonization. This paper proposes a scalable distributed solution for the operation of large populations of TCLs providing frequency response and performing energy arbitrage. Each TCL is described as a price-responsive rational agent that participates in an integrated energy/frequency response market and schedules its operation in order to minimize its energy costs and maximize the revenues from frequency response provision. A mean field game formulation is used to implement a compact description of the interactions between typical power system characteristics and TCLs flexibility properties. In order to accommodate the heterogeneity of the thermostatic loads into the mean field equations, the whole population of TCLs is clustered into smaller subsets of devices with similar properties, using k-means clustering techniques. This framework is applied to a multi-area power system to study the impact of network congestions and of spatial variation of flexible resources in grids with large penetration of renewable generation sources. Numerical simulations on relevant case studies allow to explicitly quantify the effect of these factors on the value of TCLs flexibility and on the overall efficiency of the power system.