Adaptive mesh refinement is an important technology to increase the accuracy of numerical simulations while keeping computational costs controlled. Core to the adaptive refinement process are error estimation strategies, which guide which regions of the domain should be refined. In the context of Finite Element methods for elliptic problems, robust and efficient error estimators can be derived based on the Prager-Synge identity [1]. In this approach, H 1 approximations for the primal variable (e.g., pressure) are combined with H(div) approximations for the dual variable (e.g., velocity) to obtain guaranteed upper bounds for the approximation error in the energy norm. Another approach that has gained popularity in recent years is goal-oriented error estimation [2, 3], in which the error is estimated according to a Quantity of Interest (QoI) rather than standard norms. The idea is that the QoI, represented by a functional, captures the aspects of the solution that are more relevant for the application in question. In this work, we study these two error estimation strategies in the context of mixed Finite Element approximations. Using the Darcy problem as a model, we illustrate the core differences between the two approaches, highlighting the main strengths and challenges of each. We apply both techniques to a relevant problem of Petroleum Engineering of fluid flow in a porous medium near a horizontal well. In this application, the three-dimensional Darcy equations are coupled with the one-dimensional Hagen-Poiseuille equation that models the flow inside the wellbore. Both the Darcy velocity and the wellbore flow rate are approximated using H(div) spaces, implemented using hierarchical basis functions, as described in [4]. Numerical results show how the different approaches impact the adaptive refinement process and the final solution accuracy.