Summary A reduced‐order‐modeling (ROM) framework is developed and applied to simulate coupled flow/quasistatic‐geomechanics problems. The reduced‐order model is constructed using POD‐TPWL, in which proper orthogonal decomposition (POD), which enables representation of the solution unknowns in a low‐dimensional subspace, is combined with trajectory‐piecewise linearization (TPWL), where solutions with new sets of well controls are represented by means of linearization around previously simulated (training) solutions. The overdetermined system of equations is projected into the low‐dimensional subspace using a least‐squares Petrov‐Galerkin (LSPG) procedure, which has been shown to maintain numerical stability in POD‐TPWL models. The states and derivative matrices required by POD‐TPWL, generated by an extended version of the Stanford University Automatic Differentiation General Purpose Research Simulator (AD‐GPRS), are provided in an offline (preprocessing or training) step. Offline computational requirements correspond to the equivalent of five to eight full‐order simulations, depending on the number of training runs used. Run‐time (online) speedups of O(100) or more are typically achieved for new POD‐TPWL test‐case simulations. The POD‐TPWL model is tested extensively for a 2D coupled problem involving oil/water flow and geomechanics. It is shown that POD‐TPWL provides predictions of reasonable accuracy, relative to full‐order simulations, for well‐rate quantities, global pressure and saturation fields, global maximum‐ and minimum‐principal‐stress fields, and the Mohr‐Coulomb rock‐failure criterion, for the cases considered. A systematic study of POD‐TPWL error is conducted using various training procedures for different levels of perturbation between test and training cases. The use of randomness in the well‐bottomhole‐pressure (BHP) profiles used in training is shown to be beneficial in terms of POD‐TPWL solution accuracy. The procedure is also successfully applied to a prototype 3D example case.