Statistically planar turbulent partially premixed flames for different initial intensities of decaying turbulence have been simulated for global equivalence ratios =0.7 and =1.0 using three-dimensional, simplified chemistry-based direct numerical simulations (DNS). The simulation parameters are chosen such that the flames represent the thin reaction zones regime combustion. A random bimodal distribution of equivalence ratio phi is introduced in the unburned gas ahead of the flame to account for the mixture inhomogeneity. The results suggest that the probability density functions (PDFs) of the mixture fraction gradient magnitude |delta| (i.e., P(|delta|)) can be reasonably approximated using a log-normal distribution. However, this presumed PDF distribution captures only the qualitative nature of the PDF of the reaction progress variable gradient magnitude |delta c| (i.e., P(|delta c|)). It has been found that a bivariate log-normal distribution does not sufficiently capture the quantitative behavior of the joint PDF of |delta| and |delta c| (i.e., P(|delta|, |delta c|)), and the agreement with the DNS data has been found to be poor in certain regions of the flame brush, particularly toward the burned gas side of the flame brush. Moreover, the variables |delta| and |delta c| show appreciable correlation toward the burned gas side of the flame brush. These findings are corroborated further using a DNS data of a lifted jet flame to study the flame geometry dependence of these statistics.