This paper presents an approach for developing Reduced Order Models (ROMs) of complex geomechanical systems from high-fidelity numerical simulations. The approach uses proper orthogonal decomposition of the state variables for the model reduction, combined with long short-term memory neural networks to update the simulations in response to changing boundary conditions. A key challenge in creating these models is the need to efficiently sample the system state during training and testing. Too few (or poorly distributed) training simulations yields inaccurate models, while too many results in wasted effort. However, the exact number required is difficult to determine ahead of time, as it is both context and model dependent. Here, we present a bootstrap strategy to estimate the number of simulations required to develop the reduced order model. The strategy uses an iterative approach in conjunction with a quasi-random sampling technique to ensure consistent coverage of the sample space. We demonstrate the approach by developing a ROM describing the coupled hydro-mechanical response of a section of an open-cut coal mine. Excellent agreement is achieved between the resulting reduced order model and the original high-fidelity simulation for significantly less computational effort.