Geothermal reservoir simulation often considers the coupled thermo-hydro-mechanical physics, so the computational cost is remarkably expensive, which brings challenges for rapid reservoir optimization for geothermal management. In this work, we developed a parsimonious thermal decline model with only 3 parameters, namely 𝐻𝑦𝑝𝑒𝑟𝑅𝑒𝐿𝑈 model. It can accurately predict the produced fluid temperature behavior in geothermal recovery, which captures both the early thermal breakthrough and the later decline behavior. Further, a forward surrogate model based on deep neural network is developed to map the reservoir parameters to the 𝐻𝑦𝑝𝑒𝑟𝑅𝑒𝐿𝑈 model parameters and the ultimate total net energy. The forward model is integrated with a multi-objective optimizer (MOO) based on Non-dominated Sorting-based Genetic Algorithm II (NSGA-II), which considers reservoir uncertainties of rock properties and subjects to nonlinear engineering constraints for robust reservoir optimization. The 𝐻𝑦𝑝𝑒𝑟𝑅𝑒𝐿𝑈 model is validated through processes including enhanced geothermal recovery (EGS) and geothermal recovery from hot sedimentary aquifers (HSA) without fracturing. The mean relative error of the 𝐻𝑦𝑝𝑒𝑟𝑅𝑒𝐿𝑈 model is less than 1%. We also examined the deep neural network to predict 4 parameters including the total energy and 3 𝐻𝑦𝑝𝑒𝑟𝑅𝑒𝐿𝑈 model parameters in EGS, with decent 𝑅2 scores 0.998, 0.998, 1.000 and 0.946, respectively. The MOO converges well to achieve the optimum total energy, and solutions with different (low, median, high) risk levels are consistent with the results based on reservoir simulation. The decision variables including injection temperature and rate, extraction well pressure and well distance are provided based on the MOO framework. The number of forward model evaluations during optimization is 20000, and the average CPU time of MOO based on the forward surrogate model is 28.32 s, while the optimization based simulation is estimated to be around 600 min. Therefore, the newly proposed workflow is highly scalable and ready for field or regional scale geothermal optimization.