Hydraulic fracture closure plays an important role during the shut-in period after a stimulation operation and during production from a reservoir. Fracture closure characteristics define the productivity of the created fractures. Fast fracture closure can prevent proppant from settling uniformly whereas gradual fracture closure can help in making the induced fracture system more productive. Conventional fracture closure models either use over-simplified analytical estimates of the stress dependent permeability of the reservoir or explicitly calculate the fracture width using empirical relationships. Fluid leak-off, non-uniform proppant placement in the fracture and heterogeneity in the reservoir further complicate the process. In this work, we show that the incorporation of a contact model in a fully implicit geomechanical fracture simulator is very important in determining the final propped width of the fracture. This simulator implicitly couples the reservoir solid deformation, reservoir fluid flow, fracture fluid flow and wellbore flow. A Barton-Bandis type relation describes the non-linear relationship between the magnitude of a normal contact force and the fracture aperture. We implicitly couple the Barton-Bandis contact force with the governing equations of displacement and use the Newton-Raphson method to solve the coupled non-linear system of equations. We discuss two applications in this work. We first apply our model to simulate a diagnostic fracture injection test. Pressure data in the fracture is recorded and the response is analyzed to validate the fracture closure stress. For the second application, we simulate production from a single partially propped fracture. The fracture conductivity as a function of stress is controlled by the contact model and proppant distribution. Simulation results highlight the importance of the implicit incorporation of fracture closure models in geomechanical fracturing simulators to accurately simulate fracturing, shut-in and production applications.