A statistical model is proposed for the turbulent premixed combustion regime when chemical reactions occur in interacting thin layers (flamelets). The statistics of the reacting turbulent scalar field are described in terms of a reference scalar field (rsf) with an extra probability dimension. A transport equation for rsf is derived from that for the probability density function of the reacting scalar, and the presence of interacting flamelets is reckoned through a closure hypothesis on the conditionally averaged scalar dissipation. It appears that molecular diffusivity, chemistry, and flamelet interaction terms enter explicitly in the resulting rsf equation, and in this framework microscale mixing is described within a well-posed Cauchy problem. The model proposed has been applied to simulations of homogeneous mixing with and without chemical reactions. Predicted trends for scalar variance and the flatness factor quantatively follow the Direct Numerical Simulations (DNS) data. For the reacting case the rsf model results in a nearly bimodal probability density function shape. Simulations of V-shaped lean methane- and hydrogen-air flames have also been performed with the rsf equation coupled with the standard k - epsilon model. The results fit the experimental data fairly well in spite of simple flow dynamics modelling within a boundary layer approximation.