In this study, we present an efficient and flexible adjoint-based framework for history matching and forecasting geothermal energy extraction at a large scale. In this framework, we applied the Principal Component Analysis to reduce the parameter space for representing the complex geological model. The adjoint method is implemented for gradient calculation to speed up the history-matching iteration process. Operator-based linearization (OBL) used in this framework makes the calculation of the physical state and its derivatives very efficient and facilitates the matrix assembly in the adjoint method. This study primarily focuses on history matching based on combined observation of well production and in-situ electromagnetic measurements to predict the temperature front. However, different types of misfit terms can be added to the objective function based on practical considerations. For example, our history-matching case studies include model misfit terms applied for regularization purposes. The measurement data is extracted from the true model, and realistic measurement errors are considered. Also, in this work, we propose an optimal weighting strategy for the terms of the objective function to balance their sensitivity with respect to the model control variables. The high efficiency of the framework is demonstrated for the geothermal doublet model implemented at the heterogeneous reservoir with multiple realizations. The framework allows for generating posterior Randomized Maximum Likelihood (RML) estimates of the entire ensemble of the realizations with a reasonable computational cost. Results show that the framework can achieve reliable history-matching results based on the doublets production data and the reservoir electromagnetic measurement.