Numerical models can be utilized to understand and anticipate the future behavior of a geothermal reservoir, and hence aid in the development of efficient reservoir engineering strategies. However, as each system has a unique geological context, individual characterization is required. In this research, the Nevados de Chilla<acute accent>n Geothermal System (NChGS) in the Southern Volcanic Zone of the Andes is considered. The NChGS is controlled by the geology of the active Nevados de Chilla<acute accent>n Volcanic Complex (NChVC) including their basement units (Miocene lavas and volcaniclastic layers from Cura-Mall & iacute;n Formation and the Miocene, Santa Gertrudis granitoids) as well as the key structural control from crustal scale faults, all of which combine to influence the reservoir characteristics. The presence of faults acts to generate a high secondary permeability which favors the circulation of hydrothermal fluids. Based on previous studies in the NChGS, we designed a thermo-hydraulic model in COMSOL Multiphysics (R) combining equations of heat transfer and Darcy's law in order to determine the distribution of isotherms and surface heat flux. The boundary conditions of the model were informed by a conceptual model of depth 3 km and width of 6.6 km which considers a highly fractured granitic reservoir, a clay cap behavior of Miocene lavas and volcaniclastic units, and transitional zones between a regional zone and the reservoir. A lowangle reverse fault affecting the clay cap unit was also incorporated into the models. Results indicate convective behavior in the reservoir zone and a surface heat flux of 0.102 W/m2 with a local peak up to 0.740 W/m2 in the area affected by the low-angle reverse fault zone. The models suggest hydrothermal fluid residence times of around 9-15 thousand years are required to reach a steady-state thermal configuration, which is consistent with the deglaciation age proposed for the NChVC latitude of the complex (c. 10-15 ka). Permeability in the fractured reservoir is one of the most complex parameters to estimate and the most sensitive and hence requires further constraint. Finally, using the volumetric method and the results obtained in this research, we estimate a geothermal potential of 39 +/- 1 MWe for the NChGS.