On a connected, oriented, smooth Riemannian manifold without boundary we consider a real scalar field whose dynamics is ruled by $E$, a second order elliptic partial differential operator of metric type. Using the functional formalism and working within the framework of algebraic quantum field theory and of the principle of general local covariance, first we construct the algebra of locally covariant observables in terms of equivariant sections of a bundle of smooth, regular polynomial functionals over the affine space of the parametrices associated to $E$. Subsequently, adapting to the case in hand a strategy first introduced by Hollands and Wald in a Lorentzian setting, we prove the existence of Wick powers of the underlying field, extending the procedure to smooth, local and polynomial functionals and discussing in the process the regularization ambiguities of such procedure. Subsequently we endow the space of Wick powers with an algebra structure, dubbed E-product, which plays in a Riemannian setting the same r\^ole of the time ordered product for field theories on globally hyperbolic spacetimes. In particular we prove the existence of the E-product and we discuss both its properties and the renormalization ambiguities in the underlying procedure. As last step we extend the whole analysis to observables admitting derivatives of the field configurations and we discuss the quantum M\o ller operator which is used to investigate interacting models at a perturbative level.