Earth System Models (ESMs) traditionally operate at large horizontal resolutions, on the order of 100 km, which can obscure the effects of smaller scale heterogeneity. When examining land surface states and fluxes in ESMs, one common approach to mitigate this issue is to divide the sub-grid land surface into distinct homogeneous clusters and then resolve the water, energy, and biogeochemical processes on each cluster or tile. The literature, as well as work in the Coupling of Land and Atmospheric Subgrid Parameterizations (CLASP) project, indicates that surface heterogeneity has important implications for atmospheric processes as well. Previous work using large-eddy simulation (LES) shows that spatial variability in surface heating can produce significant secondary circulations closely related to the type and scale of heterogeneity that are not captured by single column models. This presentation aims to address this persistent weakness by using a clustering or tiling approach, similar to that used with land surface processes, for the atmosphere. To accomplish this task, we run the Cloud Layers Unified By Binomials (CLUBB) single column model, a sub-grid turbulence and cloud parameterization scheme, over a 100 km box centered at the Southern Great Plains site in Oklahoma for a variety of surface and atmospheric conditions. The model is run independently over multiple surface clusters defined by surface sensible heat fluxes in the given domain. Results indicate that significant differences exist for some cases between the single column and multicolumn cases for liquid water path (LWP) as well as the turbulent kinetic energy (TKE) budget, and that results converge on consistent results with a fairly low number of clusters (i.e., atmospheric columns). We follow this up with a connected multicolumn setup where each column is dynamically connected with the other columns throughout the run to qualitatively capture the circulations observed in the LES output. The existing results show promise for capturing the effects of subgrid scale surface flux heterogeneity on the lower atmosphere in ESMs with the application of a multicolumn CLUBB setup.