  The atomic-to-molecular (HI-to-H$_2$) transition in photodissociation regions
(PDRs) has been investigated over the last several decades through analytic and
numerical modeling. However, classical PDR models typically assume uniform
density gas, ignoring the turbulent nature of the interstellar medium.
Recently, Bialy et al. (2017b, 2019) have presented a theoretical framework for
studying the HI-to-H$_2$ in a realistic turbulent medium with a non-homogeneous
density structure. Here we extend these turbulent-chemical models to explore
the possibility of tracing the magnetic field direction in turbulent PDRs using
the Gradients Technique. We utilize both subsonic and supersonic
magnetohydrodynamic numerical simulations for chemical HI/H$_2$ balance
calculations. We confirm that the density fluctuations induced by turbulence
can disperse the distribution of H$_2$ and HI fraction. We find that the energy
spectrum of moment maps gets shallower when the sonic Mach number MS increases.
We explore the ability in magnetic field tracing of gradients of higher-order
velocity centroids and compare their performance with that of traditional
velocity centroid gradients (VCGs) and with intensity gradients (IGs). We find
that the velocity gradients of the second-order centroids (VC$_2$Gs) are more
accurate than VCGs and IGs in probing the magnetic field orientation.
