For Euclidean-signature quantum field theories with renormalizable self-interactions, Holland and Hollands have shown the operator product expansion (OPE) coefficients satisfy "flow equations'': For a (renormalized) self-coupling parameter $\lambda$, the partial derivative of any OPE coefficient with respect to $\lambda$ is given by an integral over Euclidean space of a sum of products of other OPE coefficients evaluated at $\lambda$. 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. The purpose of this paper is to generalize the Holland and Hollands results for flat Euclidean space to curved Lorentzian spacetimes in the context of the solvable "toy model'' of massive Klein-Gordon scalar field theory on globally-hyperbolic curved spacetimes, with the squared mass, $m^{2}$, viewed as the "self-interaction parameter''. There are a number of difficulties that must be overcome to carry out this program. Even in Minkowski spacetime, a serious difficulty arises from the fact that all integrals must be done over a compact region of spacetime to ensure convergence. However, there does not exist any Lorentz-invariant function of compact support, so any flow relations that involve only integration over a compact region cannot be Lorentz covariant. We show how covariant flow relations can be obtained by the addition of "counterterms'' that cancel the non-covariant dependence on the cutoff function in a manner similar to that used in the Epstein-Glaser renormalization scheme. The necessity of integration over a finite region also effectively introduces an "infrared cutoff scale'' $L$ into the flow relations, which gives rise to undesirable behavior of the OPE coefficients under scaling of the metric and $m^{2}$. (This behavior also occurs in the Euclidean case.) We show how to modify the flow relations so that the dependence on $L$ is systematically removed, thereby yielding flow relations compatible with almost homogeneous scaling of the fields. A potentially even more serious difficulty arises in curved spacetime simply due to the fact that the flow relations involve integration over a spacetime region. Such an integration will cause the OPE coefficients to depend non-locally on the spacetime metric, in violation of the requirement that the quantum fields should depend locally and covariantly on the spacetime metric. We show how this difficulty can be overcome by replacing the metric with a local polynomial approximation carried to suitable order about the OPE expansion point. We thereby obtain local and covariant flow relations for the OPE coefficients of Klein-Gordon theory in an arbitrary curved Lorentzian spacetime. As a byproduct of our analysis, we prove that 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 $F_{n}(x_{1},\dots,x_{n};z)$ to the identity OPE coefficients, $C_{\phi\cdots\phi}^{I}(x_{1},\dots,x_{n};z)$, for the elementary $n$-point products. We thereby obtain an explicit construction of any renormalization prescription for the nonlinear Wick fields in terms of the OPE coefficients $C_{\phi\cdots\phi}^{I}$. 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 in Lorentzian spacetimes beyond our "toy model'' based upon the associativity properties of the OPE coefficients. We illustrate our method by applying it to the flow relations of $\lambda\phi^{4}$-theory.