For large-eddy simulation with a finite-difference scheme, a simple stochastic subgrid-scale (SGS) model is introduced which describes the effects of random scs motions on the resolved (filtered) scales of incompressible turbulent motions. The model extends the Smagorinsky-Lilly model by adding realizable random stresses and fluxes which are constructed as quadratic expressions of Gaussian random velocity and temperature fields. The random components reduce the correlations between stresses and strain rates to in between 0.16 and 0.5, in agreement with observations. The random stresses (fluxes) also induce random accelerations (temperature changes) with a k(4) power spectrum. Such random sources backscatter energy (variance) from scs motions to resolved scale motions when temporally correlated with finite timescales. The timescales are different for momentum and heat flux. The analysis of the model provides an upper estimate of the magnitude of backscatter which is close to previous predictions. The analysis identifies the influence of the quasi-normal assumption and of numerical filters and determines the variance of the pressure fluctuations induced by the random accelerations at grid scales. Backscatter increases the scs turbulent Prandtl number to a degree depending strongly on the numerical filter. Tests of the model in large-eddy simulation of isotropic turbulence show energy decay rates in close agreement with expected rates when the stochastic scs model is included. Backscatter cannot be simulated with reduced diffusivities or filter widths.