Over the past years, there has been an increase in interest in the quantification and propagation of model uncertainties to quantities of interest. In the context of computational mechanics, the most common type of uncertainty under consideration is that of model parameters. These are then propagated via a numerical model to the solution in the forward problem, or inferred from measurements in the inverse problem. This numerical model, be it FEM or any other numerical method, will have its own inherent epistemic uncertainty arising from the discretization error. Neglecting this numerical uncertainty can lead to inaccurate parameter estimates. This is especially true in an inverse modeling context, where coarse discretizations are often necessary to limit computational costs of the repeated forward model calls. Simple approaches to modeling discretization error, such as additive noise or a squared exponential Gaussian process are often unable to identify regions of large discretization error for more complex geometries. In this work, our goal is capture the discretization error associated with the finite element method by modeling it as epistemic uncertainty, represented by a probability density. Doing so would account for the discretization error in a statistically consistent manner in Bayesian data assimilation, inverse modeling and optimization contexts. Multiple approaches to arriving at such a probability density are presented, which can be divided into two categories: force-perturbation and stiffness-perturbation. In the first category, the stiffness matrix is fixed and a distribution is assumed over the forcing term, which implicitly defines a probability density over the solution space. Conversely, in the second category, the forcing term is fixed and the components of the stiffness matrix are perturbed, which again produces an implicit probability density over the solution space. For both the force-perturbation and stiffness-perturbation approach, there are different ways of introducing the perturbation, each of which with its own advantages and disadvantages. These different approaches are compared in terms of convergence behavior, quality of error estimation in both the solution field and other quantities of interest, and computational efficiency.