Parameterizing incident solar radiation over complex topography regions in Earth System Models (ESMs) remains a challenging task. In ESMs, downward solar radiative fluxes at the surface are typically computed using plane parallel radiative transfer schemes, which do not explicitly account for the effects of a three-dimensional topography, such as shading and reflections. To improve the representation of these processes, we introduce and test a parameterization of radiation-topography interactions tailored to the Geophysical Fluid Dynamics Laboratory (GFDL) ESM land model. The approach presented here builds on an existing correction scheme for direct, diffuse and reflected solar irradiance terms over three-dimensional terrain. Here we combine this correction with a novel hierarchical multivariate clustering algorithm which explicitly describes the spatially varying downward irradiance over mountainous terrain. Based on a high-resolution digital elevation model, this combined method first defines a set of sub–grid land units ("tiles") by clustering together sites characterized by similar terrain-radiation interactions (e.g., areas with similar slope orientation, terrain and sky view factors). Then, based on terrain parameters characteristic for each tile, correction terms are computed to account for the effects of local 3-D topography on shortwave radiation over each land unit. We develop and test this procedure based on a set of Monte Carlo ray tracing simulations approximating the true radiative transfer process over three dimensional topography. Domains located in three distinct geographic regions (Alps, Andes, and Himalaya) are included in this study to allow for independent testing of the methodology over surfaces with differing topographic features. We find that accounting for the sub–grid spatial variability of solar irradiance originating from interactions with complex topography is important as these effects lead to significant local differences with respect to the plane-parallel case, as well as with respect to grid–cell scale average topographic corrections. Finally, we quantify the importance of the topographic correction for a varying number of terrain clusters and for different radiation terms (direct, diffuse, and reflected radiative fluxes) in order to inform the application of this methodology in different ESMs with varying sub-grid tile structure.