  We present a detailed phenomenological study of forward hadron neutral pion
production in deep inelastic scattering, with both the direct and the resolved
contributions calculated to NLO accuracy. A comparison of the theoretical
predictions for the various distributions with the H1 data and a study of
stability of the QCD predictions under changes of scales is the focus of this
study. We obtain a very good overall description of the recent H1 data with the
choice of scale Q^2 + E^2_{perp}, in contrast to the (Q^2 + E^2_{perp})/2
required earlier when the resolved contribution was included only at LO
accuracy. We find a more modest variation of the predictions, as the scale is
changed from (Q^2 + E^2_{perp})/2 to 2(Q^2 + E^2_{perp}), as compared to the
case where the resolved contribution was included only at LO accuracy. This
variation is of the order of the rather large experimental errors.
Unfortunately, this fact prevents us from concluding that perturbation theory
gives an unambiguous prediction for forward particle production in deep
inelastic scattering. However, the overall success of perturbative QCD in
explaining the small x_{Bj} data means that perhaps a full resummation of the
BFKL ladder is not called for. We notice the need for rather large resolved
contributions to explain the data at low x_{Bj} even at somewhat larger Q^2
values.
