Given that land surface processes most often depend nonlinearly on soil moisture and that coarse-grained soil moisture fields are often used in modern large-scale applications, it becomes necessary to account in some way for the subgrid variability of physiographic properties and their interactions with subgrid soil moisture dynamics in order to derive robust estimates of grid-scale land surface fluxes. Albertson and Montaldo <2003> present a conservation equation for subgrid soil moisture spatial variability, which includes production and destruction terms related to the various interactions between moisture fields and flux fields (i.e., infiltration, drainage, evapotranspiration, and horizontal redistribution). Here we present a rigorous closure of the conservation equations of the first and second moments of the soil moisture field (omitting topographic redistribution in this effort). Closed terms include correlation coefficients between soil moisture and physiographic properties. The approach presented here provides accounting for the effects of dynamic interaction between soil moisture spatial variability and the underlying soil texture and vegetation density fields. Closure is achieved in terms of the spatial mean soil moisture evolution, with closure coefficients being functions of the correlation coefficient between the underlying soil texture and vegetation density patterns (i.e., the physiographic setting). Making use of these relationships, the proposed approach uses known subgrid soil texture and vegetation density fields to accurately simulate the temporal trajectory of subgrid variance of soil moisture, which supports improved estimates of grid-scale land surface fluxes. The approach is tested synthetically (closure model versus distributed model) over a wide range of physiographic characteristics and meteorological conditions, with results demonstrating that the closure model based approach reproduces the temporal dynamics of the soil moisture variability and provides improved grid-scale land surface flux estimates. In fact, the second-order closure was shown to reduce the errors in integrated land surface flux estimates, as compared to those derived from zero-order approximation estimates, by an order of magnitude in certain conditions.