The most dynamic electromagnetic coupling between the magnetosphere and ionosphere occurs in the polar upper atmosphere. It is critical to quantify the electromagnetic energy and momentum input associated with this coupling as its impacts on the ionosphere and thermosphere system are global and major, often leading to considerable disturbances in near-Earth space environments. The current general circulation models of the upper atmosphere exhibit systematic biases that can be attributed to an inadequate representation of the Joule heating rate resulting from unaccounted stochastic fluctuations of electric fields associated with the magnetosphere-ionosphere coupling. These biases exist regardless of geomagnetic activity levels. To overcome this limitation, a new multiresolution random field modeling approach is developed, and the efficacy of the approach is demonstrated using SuperDARN data carefully curated for the study during a largely quiet 4 hours period on February 29, 2012. Regional small-scale electrostatic fields sampled at different resolutions from a probabilistic distribution of electric field variability conditioned on actual SuperDARN LOS observations exhibit considerably more localized fine-scale features in comparison to global large-scale fields modeled using the SuperDARN Assimilative Mapping procedure. The overall hemispherically integrated Joule heating rate is increased by a factor of about 1.5 due to the effect of random regional small-scale electric fields, which is close to the lower end of arbitrarily adjusted Joule heating multiplicative factor of 1.5 and 2.5 typically used in upper atmosphere general circulation models. The study represents an important step towards a data-driven ensemble modeling of magnetosphere-ionosphere-atmosphere coupling processes.