Inverse problems (IPs) appear in many science and mathematics branches, mainly in situations where quantities of interest are different from those we can measure; model parameter values must be esti-mated from the observed data. Bayesian approach to IPs or recently known as Bayesian uncertainty quantification (UQ), is becoming ever more popular to approach these kinds of problems. The posterior distribution of the unknown parameters quantifies the uncertainty of these parameters' possible values consistent with the observed data. However, explicit analytic forms are usually not available for the posterior distributions, so sampling approaches such as the Markov chain Monte Carlo are required to characterize it. These methods involve repeated forward map (FM) solutions to evaluate the likelihood function. We usually do not have an analytical or computationally precise implementation of the FM and necessarily involves a numerical approximation, leading to a numerical/approximate posterior distribution. Thus, the numerical solution of the FM will introduce some numerical error in the posterior distribution. In this thesis, a bound on the relative error (RE), in posterior expectations for some functional, is found. This is later used to justify a bound in the error of the numerical method used to approximate the FM for Bayesian UQ problems. Moreover, our contribution is a numerical method to compute the after-the-fact (i.e., a posteriori) error estimates of the numerical solution of a class of semilinear evolution partial differential equations and show the potential applicability of the result previously obtained for bounding the RE in posterior statistics for Bayesian UQ problems. Several numerical examples are presented to demonstrate the efficiency of the proposed algorithms. We also present a Bayesian approach to infer a fault displacement from geodetic data in a slow slip event (SSE). Our physical model of the slip process reduces to a multiple linear regression subject to constraints. Assuming a Gaussian model for the geodetic data and considering a multivariate truncated normal (MTN) prior distribution for the unknown fault slip, the resulting posterior distribution is also MTN. Regarding the posterior, we propose an ad hoc algorithm based on Optimal Directional Gibbs that allows us to efficiently sample from the resulting high-dimensional posterior distribution without supercomputing resources. As a by-product of our approach, we are able to estimat