The singular properties of quantum fields have posed an intransigent obstruction to formulating a mathematically well-defined theory of interacting quantum fields at nonzero coupling. To date, a systematic renormalization of the unavoidable "ultraviolet" divergences produced by pointwise products of quantum fields has only been achieved order-by-order in perturbation theory. In the coincidence limit, the behavior of products of quantum fields is characterized by the coefficients of their operator product expansion (OPE). For Euclidean quantum field theories, Holland and Hollands have shown the OPE coefficients satisfy "flow equations": For interaction parameter λ, the partial derivative of any OPE coefficient with respect to λ is given by an integral over Euclidean space of a sum of products of other OPE coefficients. These Euclidean flow equations were proven to hold order-by-order in perturbation theory, but they are well defined non-perturbatively and, thus, provide a possible route towards giving a non-perturbative construction of the interacting field theory. In this thesis, we generalize these results for flat Euclidean space to curved Lorentzian spacetimes in the context of the solvable "toy model" of massive Klein-Gordon scalar field theory, with the squared mass viewed as the "self-interaction parameter". Even in Minkowski spacetime, a serious difficulty arises from the fact that all integrals must be taken over a compact spacetime region to ensure convergence but any integration cutoff necessarily breaks Lorentz covariance. We show how covariant flow relations can be obtained by adding compensating "counterterms" in a manner similar to that of the Epstein-Glaser renormalization scheme. We also show how to eliminate dependence on the "infrared-cutoff scale" L, thereby yielding flow relations compatible with almost homogeneous scaling of the fields. In curved spacetime, the spacetime integration will cause the OPE coefficients to depend non-locally on the spacetime metric, in violation of the requirement that quantum fields should depend locally and covariantly on the metric. We show how this potentially serious difficulty can be overcome by replacing the metric with a suitable local polynomial approximation about the OPE expansion point. We thereby obtain local and covariant flow relations for the OPE coefficients of Klein-Gordon theory in curved Lorentzian spacetimes. As a byproduct of our analysis, we prove the field redefinition freedom in the Wick fields (i.e. monomials of the scalar field and its covariant derivatives) can be characterized by the freedom to add a smooth, covariant, and symmetric function to the identity OPE coefficients for the elementary n-point products. We thereby obtain an explicit construction of any renormalization prescription for the nonlinear Wick fields in terms of these elementary OPE coefficients. The ambiguities inherent in our procedure for modifying the flow relations are shown to be in precise correspondence with the field redefinition freedom of the Klein-Gordon OPE coefficients. In an appendix, we develop an algorithm for constructing local and covariant flow relations beyond our "toy model" based on the associativity properties of OPE coefficients. We illustrate our method by applying it to the flow relations of an interacting scalar field with a quartic potential.