In this paper, we formulate the Load Flow (LF) problem in radial electricity distribution networks as an unconstrained Riemannian optimization problem, consisting of two manifolds, and we consider alternative retractions and initialization options. Our contribution is a novel LF solution method, which we show belongs to the family of Riemannian approximate Newton methods guaranteeing monotonic descent and local superlinear convergence rate. To the best of our knowledge, this is the first exact LF solution method employing Riemannian optimization. Extensive numerical comparisons on several test networks illustrate that the proposed method outperforms other Riemannian optimization methods (Gradient Descent, Newtons), and achieves comparable performance with the traditional Newton-Raphson method, albeit besting it by a guarantee to convergence. We also consider an approximate LF solution obtained by the first iteration of the proposed method, and we show that it significantly outperforms other approximants in the LF literature. Lastly, we derive an interesting comparison with the well-known Backward-Forward Sweep method.