Production forecast and estimation of the recovery factor of the Los Humeros geothermal field, Mexico

Production data from the Los Humeros geothermal field in Mexico are used to develop a forecast method for operation of the field in the future. This method supports understanding the limitations of sustainably exploiting a geothermal reservoir. Using such forecasts, fluid extraction could be scheduled in order to fulfill the steam demand of the installed capacity. Moreover, the forecast of the extraction can diminish the commercial risk involved and therefore clear the way for new investments. Herein, we determined the total extracted heat from the Los Humeros geothermal system, with the aim to forecast the next phase of fluid extraction. We took the information from 29 producing wells to analyze their historic geothermal fluid production. From their statistical distribution of data, we estimated the amount of extracted fluids for a future phase. The size of the heat reserve previous to exploitation was calculated from the forecasted and the extracted heat. Thus, the recovery factor of the system was calculated as the ratio of the accumulated production of heat and the initial heat in place. The general recovery factor for 60 years of production was calculated as 37–44%. Due to the heterogeneity of the system, this study conducted individual assessments. The forecasted heat extraction for the next 30 years of production, amounts to 580 PJ assuming a constant extraction rate between 6 and 56 ton·[h]−1 depending on the particular well. The results of this study potentially complement the existing models of the geothermal exploitation and offer an estimation of the future production of the Los Humeros, which could be a good basis for decision-making and management of the site.

In Los Humeros, Mexico, a conventional geothermal system has been exploited since the early 1980s. During 30 years of continuous production, several surveys have been conducted about the response and evolution of the geothermal system as a consequence of fluid extraction (Aragon-Aguilar et al. 2017, 2019García-Gutiérrez et al. 2002;Arellano-Gómez et al. 2008). Moreover, as Los Humeros is considered to be a potential superhot geothermal system, it is one of the study sites of the Mexican-European collaborative research project GEMex (Jolie et al. 2018). The objective of GEMex was the identification of critical parameters to utilize the superhot geothermal resources. Although the use and development of a superhot geothermal system could happen in the near future, the assessment and monitoring of the conventional geothermal system of the Los Humeros caldera is needed in order to extend and preserve its use for future generations.
Previous studies about the evolution of the volcanic system exploitation show changes on temperature, pressure and thus enthalpy at reservoir depth. In general, these changes do not seem to affect the performance of the wells (Aragon-Aguilar et al. 2005; Arellano-Gómez et al. 2008). The analysis of productivity data suggests the existence of two reservoirs separated by a tuff layer (Gómez et al. 2003). The performance of the wells could be linked with both the grade of their connection with the deeper and hotter reservoir and with the physicochemical parameters of the brines (Aragon-Aguilar et al. 2005;Arellano-Gómez et al. 2008). In this context, the evaluation of the size of the reserve based on the production, which is the focus of this paper, substantially complements the characterization of the system.
The main goals of this work are first to estimate the original amount of heat in the currently exploited reservoir of Los Humeros and second to forecast the mass extraction for the next 30 years. The extrapolation time is taken from the time frame suggested by development banks to these types of energy projects. The results of the goals allow us to calculate the heat recovery factor of the system for 60 years of continuous production. The forecasted production is assumed as economically exploitable, as the necessary technology is already installed. Hence, we consider this quantity to be the heat reserve of the Los Humeros (Muffler and Cataldi 1978;Rybach 2015).

The Los Humeros production history
The Los Humeros caldera is a liquid-dominated system. It is located in the Central part of Mexico, at the eastern part of the Trans Mexican volcanic belt (Fig. 1). The drilling activities began in 1986. After 4 years of exploration, drilling, well testing and constructing infrastructure and facilities, the power plant began operation in 1990 (Méndez et al. 2011). The power plant started with 5 MW of installed capacity and the reinjection of brine commenced in 1992 (Cendejas 1992). Today, the total capacity has increased to 95 MWe, distributed in three main production clusters: north, center and south, each one with different production units (Garcia-Gutierrez et al. 2015;Gutiérrez-Negrín 2019). The injected mass amounts to 16% of the extracted mass (Arellano-Gómez et al. 2008), and the reinjected temperature is 90 °C on average according to the available operation data provided by the Federal Commission of Electricity (CFE), the proprietor company. Nowadays, there are five injection wells operating inside the production zone. In spite of a temperature drop between production and reinjection of more than 100 K, the reinjection does not seem to affect the production of the hot fluid (Cendejas 1992; Iglesias et al. 2012).
The hot brine extraction network consists in 29 wells, with an average depth of 2100 m. Production wells use a pump to extract the brine from the reservoir. Then, the brine is conducted through a separator unit, splitting the flow into liquid and steam. Resulting brine is reinjected, while the steam flows to the turbines by different pipe arrangements (López and César 2006). Available monitoring data, collected on site since the beginning of operation, include the well head pressure, the extracted mass as a sum of brine and steam, and steam enthalpy. Although the reports state the use of pumps, they do not offer any description of them. Our objective is to, based on this data, estimate the amount of heat that has been produced from the system, the prospects for future production and the implications for the sustainability of extraction.

Production forecast and estimation of geothermal heat reserve
The production forecast using historical production data includes three main categories: decline analysis, lumped parameter and numerical simulations. Some authors have applied these procedures analyzing the productivity data for different geothermal systems. Decline analysis takes a time series of one variable and looks for trends to adjust it and forecast the behavior systems (Atkinson et al. 1978;Nathenson 1975). This procedure helps to identify signs of reservoir's pressure depletion and therefore a consecutive decline in the future production by a trend analysis. The lumped parameter is a tool that considers the reservoir as a closed tank. It takes the average properties of the reservoir in order to find a declining trend between the joint analysis of cumulative production and reservoir pressure (Westwood and Castanier 1981;Ciriaco et al. 2020). Each methodology depends on the quality of the data and the goal of the estimation. For our case, Fig. 1 Location of the Los Humeros, Mexico. Los Humeros caldera is located at the eastern part of the trans Mexican volcanic belt (Ferrari et al. 2012). Inside the TMVF are 2 more power plants under exploitation. The Cerritos Colorados power plant is under construction. The total installed capacity of electricity generation using geothermal sources in Mexico is 975.5 MWe (Gutiérrez-Negrín 2020). DEM: INEGI 2020.*Power plant under construction the available data were not sufficient to carry out these approaches. The time series of extraction rates used in this study do not show any kind of trend (Fig. 2).
The estimation using a volumetric heat in place model for the Los Humeros has already been made. In that study, the authors performed a numerical simulation including the heat flow and fluid transport equations (Bonté et al. 2020). From the geological model of the Los Humeros structure, they consider all the layers of rock with the same heat potential and with their respective physical characteristics. Their result offers an estimation of the energy content of the whole caldera. The estimation proposed here of the size of the Los Humeros reserve is limited to the productivity zone and it is based on the production data; moreover, it considers the layer of rocks along the feed zones reported on the well logs by the CFE. Gómez et al. (2003) suggested a pressure model to estimate the undisturbed distribution of pressure along the reservoir; later, Arellano-Gómez et al. (2008) calculated a decline pressure rate with values between 0.9 and 1.03 bar/year. The use of lumped parameter demands direct measurements of pressure (Brigham and Neri 1979). Unfortunately, the Los Humeros lacks an observation well and therefore direct records of pressure changes. Although the use of Arellano's approach might lead to a pressure distribution during exploitation, this estimation must be validated against direct measurements. Moreover, extraction data do not show a clear declining rate for the majority of the wells (Fig. 2), making the decline analysis difficult to perform.
In this sense, the method to forecast heat recovery presented in this study is based on the statistical evaluation of the data. In order to perform it, two quartiles of the heat production history of each well are calculated. The first quartile represents a pessimistic scenario, while the third quartile stands for an optimistic scenario of production performance. This approach considers the variability of production (as a range between those two quartiles) while ignoring irrelevant outliers that determine the absolute minimum and the maximum, the obvious characterizing the range of values. The production was forecasted assuming the Fig. 2 Extraction history of the most productive wells of the system established extraction regime for the next 30 years without any recharge by conduction or convection from below. Furthermore, the initial amount of heat with a sum of the reserve plus the accumulated production of heat is estimated. These two quantities allow the calculation of the heat recovery factor.

Methodology
This study involves the productivity analysis of 29 of 54 wells with the aim to estimate the heat reserve and the heat recovery factor of the Los Humeros productive geothermal system. Only 29 wells on Los Humeros were used for fluid extraction. The data were collected during 30 years of continuous production and provided by CFE. The cumulated extracted heat was calculated from the fluid and enthalpy reported per month per well. The production of the next 30 years was forecasted from the statistical analysis of the data as described above. The recovery factor is the ratio of the extracted quantity and the initial amount of heat. Heat loss through the pipe wall and heat flow from the reservoir to the fluids were neglected. Furthermore, the possible interference of the injection wells was not considered.

Extraction of heat
The productivity data provided by the CFE include mass extraction rate (ṁ), enthalpy (h), and well pressure. The data show two frequencies: from 1986 to 2000 it was recorded daily, but from 2000 to 2017 it changed into a monthly record. In order to have the same frequency, we calculated a monthly average of the daily data. Then, we described the data distribution to make the extraction forecast. The estimation of the total available heat at the surface began with the calculation of the heat flow ( Q ), multiplying the mass flow rate (ṁ) and the enthalpy (Δh). The mass flow rate was reported as a sum of steam (ṁ′′) and brine (ṁ′). However, the provided enthalpy describes to the vapor phase (h′′). Without pressure data to calculate the corresponding enthalpy of the liquid phase, the model only considered the steam fraction, which represents more than 95% of the total mass flow for the majority of the wells. The last member in the model to get the available heat at the surface per well (Q e ) was the time difference between data points (Δt): where h ′′ = h ′′ − h 0 , and h 0 is the value of enthalpy of water at Los Humeros average annual temperature, i.e., 21 °C (Vidal-Zepeda, 1990).
Then, we estimated the area of the reservoir affected by the fluid extraction. Since we neglect heat loss, the available heat at the surface equals the same quantity in the reservoir (Q R ); i.e., Q R = Q e . To estimate the energy stored, the volumetric heat in place consider an effective thickness ( �ξ ) times an area (A), i.e., an effective volume (Limberger et al. 2018). In our case, we are assuming that this area is the use that each well makes in the reservoir. Hence, this calculation was carried out by the combination of the volumetric heat in place and Eq. (1): where T fz is the temperature difference between the feed zone and the injected brine; �ξ is the effective thickness of the feed zone. In this model ρc p is the weighted average of the volumetric heat capacity of the rocks. This includes the isobaric heat capacity (c p ), the density (ρ) and the porosity ( φ ) for each different rock type (r) along the feed zone of the well: where w is water; ( N , n ) are the number of different rock types along the feed zone, and the weight ω is the depth percentage of the respective rock type around the corresponding well ( N n=1 ω n = 1 ). Therefore, the weighted average is: We calculated the depth percentage of rocks from the borehole reports of the CFE. The values to feed each member of ρc p are shown in Table 1, it also includes the porosity of the rocks ( φ ) measured at laboratory conditions. It has to be noted that these values are vulnerable to change due to the pressure exerted by ca. 1.8 km of rocks (which corresponds to the average depth of the wells' feed zones) (Petford 2003). Moreover, the Los Humeros caldera is within a highly complex geological site, resulting in a heterogeneous productivity zone. With this in mind, the criteria to choose the petrophysical values were based on the location of the samples reported on the source (Weydt et al. 2020). Those samples inside the productivity zone were chosen. Then, the stratigraphic classification of the rock samples was tied with the CFE boreholes reports in order to approximate the best candidate which fits with the feed zone depth. In order to measure how this reduction in porosity affects the calculation on the area, a sensitivity analysis was carried out. This analysis took three minor values that are reported in Table 1: 80%, 50%, and 20%.
As previously mentioned, the difference T fz = T fz − T inj , takes the temperature at the feed zone depth (fz) and the plant's brine reinjection (inj). Finally, the term �ξ is the effective thickness of the aquifer. In this particular case, the effective thickness was considered identical to the permeable zone per well. We took the location and thickness of the permeable zones from injection test reports of the CFE. According to the reports, the majority of the wells have more than one permeable zone, although it is not clear which ones are the actual feed zones. Lacking this information, we assume a sum of the thicknesses of the permeable zones as well as an average of T fz (Fig. 3). (3) ω n ρ wn c wn p φ n . Table 1 Physical parameters of rocks, σ is the standard deviation (Weydt et al. 2020)

Forecast of production and estimation of the original reserve of heat
The following forecast was based on the productivity rate of each well. The equation that was applied is: where the subscripts C1, C3 are first and third quartile of the production flow rate, which define the optimistic (opt) and pessimistic (pes) scenario. t ex is the extrapolation time, in this case 30 years. Q ex is the production forecast for the next period. In our case, we are considering 30 years as the duration of a production period, therefore Q ex → Q 30 , and t ex = 30 years. The theoretical used area was also calculated with Eq. (2). The initial reserve of heat ( Q 0 ) was calculated with the sum of the extracted ( Q e ) plus the forecasted heat ( Q 30 ): To validate Eq. (6), we performed a test using the known production data. We applied this methodology to the first 20 years of production to forecast the last ten. After calculating the results, the forecast was compared with the monitoring data. Theoretical heat content The theoretical heat content was calculated by combining the results obtained by Bonté et al. (2020) and the accumulated heat production ( Q 30 ). From Q 30 , the used area was calculated with Eq. (2). Then, this area was multiplied by the average heat density value reported in the map from Bonté et al. (2020). This result was considered as the theoretical heat content.

Recovery factor of the Los Humeros
The last parameter estimated was the recovery factor ( R f ). The R f is defined as the ratio between the available heat produced at the wellhead ( Q wh ) and the theoretical heat potential ( Q R ) (Gringarten 1978;Muffler and Cataldi 1978;Williams 2007;Garg 2010;Grant 2014): In this research, we did not focus on the theoretical potential of the system. Instead, R f was calculated as the ratio of the actual extracted heat ( Q e ) and the original reserve of heat ( Q 0 ): The R 60 is an estimation of how much heat can be extracted from the system under the particular circumstances of 60 years of production, 29 wells and holding a certain rate of extraction.

Results
Some of the statistics of the data are shown in Fig. 4b. The violin charts show the distribution of the data including the first and the third quartile, as well as the median of the data, being the second quartile. The steam productivity ranges from 6 to 56 ton·h −1 . The first and third quartiles of each well were used in the pessimistic and optimistic scenarios. Individual assessment of the extracted heat and the forecast is shown in Table 2. Results displayed in Table 2 show a huge production variability among the wells. The validation of the forecast model shows that 75% of the production lies within the forecasted value (shown in Fig. 5).
After 30 years of continuous production, the extracted heat from the system is calculated as approximately 340 PJ until 2017. The forecasted total heat production for the next 30 years is between 430 (pessimistic) and 580 PJ (optimistic). Therefore, the total reserve size amounts to 770-920 PJ. These are the resulting sums of the individual assessments for each well. The corresponding recovery factor in 2017 is 37-44%.
The area calculated with Eq.
(2) was plotted in the map of Fig. 6b. This allows us to appraise the extension of the area of the reservoir that each well has affected. For this estimation, no recharge was considered. A single sensitivity calculation was carried out in order to demonstrate the influence of the porosity on the results of Eq. (2).
The map in Fig. 6b includes the volumetric heat in place calculated during the GEMex project by Bonté et al. (2020). The accumulated production of heat ( Q 30 ) until 2017 is 340 PJ. The equivalent area covered by the reservoir beneath is supposed to be 4.1 km 2 , calculated with Eq. (2). The average value of the heat in place is 1750 PJ km −2 . Therefore, the resulting Q R is 7175 PJ and the R f is 4.75% (Table 3).   6). The observed value for the extraction is represented as a dot. Forecast using our model with 80% of the data for training is shown as a range whisker. The model is considered to be appropriate when the observed value is inside the range. Wells enclosed in parentheses operated for 22 years and they are not considered for validation as they do not have enough data Fig. 6 a Full extension of the Los Humeros caldera color indicating the density of the heat in place (Bonté et al. 2020). b Resulting area affected by the production of the heat reserve of the Los Humeros nor the extraction of heat in the future. This study addresses these topics and related arguments. The energy extraction forecast based on the history of production and the estimate of the theoretical spatial size of the original reserve of heat was discussed in brief. The forecast for energy was 430-580 PJ with a model accuracy of 75%. The original reserve is 770 and 920 PJ depending on a pessimistic or optimistic scenario, respectively. Assuming no recharge, the recovery factor for 60 years and 29 wells, under a production rate of 6-56 ton·h −1 , suggests that less than 45% of the system's heat can be harvested.

Past
The violin charts in Fig. 4 show the historical distribution of the production data. There is high variability in the extraction range of the wells. The highest value which is an outlier, belongs to well H31. But the highest median belongs to well H12. The violin chart of the well H12 presents two peaks, which could be related to two different regimes of production, just like its neighboring wells (H6 and H39). The extraction history for the well H12 (Fig. 2), shows a lot of variation for older records than the year 2000. Although both plots were made after the frequency reduction, if this bimodality were a consequence of this data treatment, all the wells would present this distribution. Moreover, the well H39 still present this bimodality and this well did not suffer a frequency reduction. This bimodality might be related with local changes in this cluster or to their geographical position.
One example of geographical influence is the faulting. These three wells (H6, H12, H39) are in the south, they are close to each other, and, they are next to Maztaloya fault. Therefore, the permeability at this zone caused by the fault might have a seasonal effect linked with the rainfall recharge in the aquifer. Although, the bimodality could be also a consequence of different surface processes, as different pumping rates for example. If this change was the consequence of the bimodality, these three wells should to experience similar changes on their surface processes. But, Aragon-Aguilar et al. (2005) suggested a geographical classification of the wells, according to them those wells close to the faults have a similar productivity index and extraction rate. Following this idea, those wells with a bimodal distribution (H1, H8, H16, H6, H12, H9, H39) are at the left of the Maztaloya fault. Yet, this feature is not conclusive as more wells are at the left of this fault and they do not present a bimodality. Maybe the bimodal wells not only are to the left of Maztaloya, but also, their feed zones are within the same rocks.

Table 3 Sensitivity of the area against changes in porosity
From the lack of a proper porosity-depth model for Los Humeros, this comparison took 3 different values of φ . r(φ 100% ) represents the porosity average for this well considering the value at lab conditions reported by Weydt et al. (2020). The entries r(φ 80% ), (φ 50% ), (φ 20% ) show the sensitivity on the radius calculation when the porosity decreases by 80, 50 and 20%, respectively. It can be seen that the radius is sensitive to these changes and the highest level of variation is when only 20% of the original value is considered H12 100 + 2 + 6 + 10 H9 100 The extraction rate of the wells is quite different. The three most productive wells have one feed zone below 800 m (a.s.l.). According to Gómez et al. (2003) below this level, a second and hotter aquifer is found. Moreover, the only feed zone of the well H9 is below this level and it is the most productive well in the field. Thus, the high productivity rates could be a consequence of the location of this feed zone, or the deepest reservoir is more productive. This observation agrees with Arellano-Gómez et al. (2008). In their work, they showed a change in the physical parameters of the fluids. They linked the decrease of the temperature of the produced fluids to a bad connection with the deeper aquifer. Nevertheless, a feed zone located below the 850 m (a.s.l) is not the only requirement to have a high productivity rate as the wells H16, H11, and H17 also have a feed zone located below the 850 m, but their production is much lower than well H7 for example.

Future
The model used to forecast the production of the data is very simple. It is based on a constant future extraction rate set by the two quartiles. Despite its simplicity, the validation of the model shows that 75% of the forecasted heat production values are within the quartile range (Fig. 5). Nevertheless, when comparing this result with the forecast calculated with the decline analysis and lumped parameter method, the approach of this study looks inaccurate. Both procedures were applied in the geothermal field Cerro Prieto, having an accurate prediction of the future values and a good match between the production and the trend line proposed (Hector and Campbell 1990;Westwood and Castanier 1981). However, there are some remarks to be made about the presented model herein.
Since the presented model is based on the mass flow rate, the forecast shows the quantity of heat available at the surface if the extraction rate remains constant. For example, the well H7 has an accumulated heat production of 34.4 PJ and the forecast for the optimistic scenario is 34.3 PJ; both numbers after 30 years of continuous extraction. This result might offer a clue about the extraction rate needed to produce a similar amount of heat. The assumption of a constant rate of extraction is a bit naïve, especially after the violin charts show a huge heterogeneity in the data, but it offers a rough perspective about where the limits of the production could take place in the future. Still, this forecast is not entirely meaningless. It helps to estimate the reserve of the heat of the Los Humeros.
The reserve magnitude calculated for the Los Humeros equals the size of the system's technical potential as defined by Rybach (2015), i.e., it is an exploitable amount of heat under the actual legal frame and using the installed technology (Rybach 2015). The cumulative heat production per well tells us how much heat has been already taken (heat loss is neglected therein). By the premise of a constant production rate, the model forecasts how much heat can be produced. The sum of these quantities offers an approximate insight of the original heat in place assuming a "closed tank", i.e., without recharge. Nevertheless, the limits of the production are not considered at all. Since the mass extraction has a close relationship with the system pressure, this assumption is far from reality. More extraction could reduce reservoir pressure. In this sense, a better approximation requires precise pressure monitoring.

Influence of the wells in the reservoir
To estimate the reservoir expansion projected onto the topographic surface, a geometrical concept was applied using Eq. (2). This expression calculates an equivalent heat area, i.e., the area needed to have equality between the cumulative heat production ( Q 30 ) and the stored heat in place ( Q 0 ). This area is presented as a circle, which in turn is the 2D surface projection of a 3D subsurface cylinder (Fig. 3). The circles traced per well show the resulting area calculated with Eq. (2) assuming a radial and steady flow (Fig. 6). Equation (2) describes the growing area as a function of extracted heat ( Q e ). These circles will grow until production is stopped. Thus, the circle's area might show the limit of the maximum level of extracted heat for the accumulated time of production.
When two cylinders' outer boundaries meet, the heat production of the implicated wells could be affected. If two cylinders meet, their shapes merge, decreasing the heat flow rate in the overlapping zone to the wells. To compensate for this reduction in the flow, the other parts of both cylinders will have to grow faster.
The results of the model show that the wells H15, H30, and H31 might be good candidates to represent this effect. Figure 4a shows the feed zone of these three wells. Since it is at the same depth these wells could share the same aquifer. The deformation of the cylinder is a much more complicated case since the feed zones are not at the same depth nor the same thickness. The cylindric deformation could include a stretching or shortening of the feed zone. This effect is beyond the presented consideration.
As previously stated, recharge is not considered. However, when the recharge of fluid is included, this effect will be different. The encounter of cylinders might derivate into a decrease of the water level in the effective thickness, i.e., a drop of pressure. In response to this pressure drop, a higher rate of the aquifer fluid might enter into the cylinder. Thus, the deformation of the shape could be slower or it could go to one well. For example, if the pumping rate increases in well Hy (Fig. 7), this change on fluid demand could take fluid supply from the well Hx, diminishing the whole performance of Hx. In any case, this encounter effect results in an extra demand of fluid.
Here is another consideration. Although the injection wells do not seem to affect the temperature of fluid (Cendejas 1992;Iglesias et al. 2012), to satisfy this demand the cold front of injection could move faster towards these wells. If this new fluid is colder, then the temperature might be decreased affecting the performance of the implicated wells. Although no decrease in heat production has been recorded yet, it could happen within the next 30 years of production. Fig. 7 Sketch of the circles' deformation when two wells share the same aquifer. The size of the arrows is proportional to the heat flow into the well. In the central sketch, these arrows are larger. The right sketch represents the circles' deformation as a consequence of a higher extraction rate. The name tag makes a reference to any couple of wells, sharing the same aquifer and whose circles meet The CFE reports do not mention from which feed zone the fluid is coming into the well and there are no pressure measurements regarding the flow. Therefore, it was assumed that the geothermal brine is streaming into the well in all the feed zones. This model is extremely simple, but it might offer a very rough visualization of the possibility of wells interruption. Therefore, there is no evidence of this interruption measurable at the surface since the productivity of all the wells has not decreased. Thus, if there is any interruption in extraction between the wells, nothing at the surface could indicate the effectiveness of the recharge of hot fluids.
Finally, the porosity changes with depth have an implication on the results of Eq.
(2). Taking the three most productive wells as an example, it can be seen that the area of each well increases up to 13% when the porosity is reduced by 80% of its measured value at lab conditions. The inverse of this relation is set by Eq. (2), where:

Technical potential extension
The size of the circles represents the used area, it also could be seen as the affected surface area of the system. The affected surface area is the portion of the reservoir that actually contains the hot fluid that was extracted by the wells. So, the area of Eq. (6) is the Los Humeros equivalent reserve extension [ A(Q 0 ) ]. The difference between A(Q 0 ) and A(Q 30 ) is the potential area to be claimed with the installed equipment. In other words, it is the extension of the technical potential area of the Los Humeros.
The theoretical potential is 20 times bigger than the technical one. Bonté et al. (2020) consider a much larger reservoir thickness (ca. 5 km). In their work, each layer of rock has the same potential to be a reservoir of heat. In this study, it is not the case. It is just considered the effective thickness of the aquifer with less than 1 km (H12). In any case, this comparison just reflects the general idea about the technical potential of a geothermal resource, from which the entire potential cannot be extracted (Nathenson and Muffler 1975;Muffler and Cataldi 1978;Rybach 2015).

Technical recovery factor
The recovery factor based on 60 years of continuous production gives an insight into the maximum producible heat under these particular conditions (29 wells, extraction rate etc.). This amount of heat is expressed in the technical potential of the system. In Los Humeros, all the technology needed to harvest the heat is already installed. Thus, this reserve is the technical portion of the system heat. Ergo, the R 60 delimitates the maximum level of heat that can be harvested from the technical potential stored.
The heat recovery factor ( R f ) is closer to the relation of what is theoretically contained in the reserve vs. what can be extracted from it (Gringarten 1978;Muffler 1979;Williams 2007;Garg and Combs 2015;Grant 2014). In this sense, R f can help to delimitate the technical potential from the theoretical potential of the Los Humeros. It could define the maximum amount of heat attainable from the theoretical potential reserve.
The individual results show that the older wells have close up R 60 . Although each well is unique, this parameter can be found within a range imposed by their technical boundaries. Therefore, it gives an idea of the production limits. The R f , calculated with the accumulated production and the theoretical heat in place, is similar to the value reported in the literature. The R 60 shows how much energy can be extracted under a given time frame. Equation (6) implies that the size of the reserve as demonstrated in this study is defined by the extraction rate and the production time. In this sense, when this production phase ends in 30 years, another recovery factor could be defined (in that case R 90 ).

Sustainability of the system
The presented approach requires a stationary concept of operation. Since neither heat nor fluid recharge was considered, the production cannot be claimed as sustainable in the actual sense of the word. Nevertheless, the geothermal system probably experiences recharge. This hint is shown by the wells H15, H30, and H31. If this recharge were in equilibrium with the extraction, the Los Humeros would be a fully renewable system, as described by Stefansson (2000). In contrast, O'Sullivan et al. (2010) argue that full renewability is rarely reached. Yet the production of wells H15, H30, and H31 is prominent, this could be related to other reasons. For example, the original reserve of energy is bigger than the actual calculation (i.e., a larger tank). In this case, a colder body of water has not been reached by the cylinder extensions.
In this work, we assumed 30 years as the duration time per production phase. The forecast of energy offers an assessment of the productivity of the system over this particular or other defined time spans. Based on these results, the system can sustain another period of production under comprehensive management of the site (reasonable extraction rates, enough injection etc.). Based on the volcanic nature of the Los Humeros location, the availability of enough heat to ensure the longevity of the exploitation is very likely, thus the main challenge is keeping the balance between extraction of fluids and the natural recharge of hot fluids, including injection. Although the injection of colder water does not seem to affect the whole production, precise pressure monitoring is needed, even the use of observation wells to enhance this activity is highly recommended. Axelsson et al. (2004) presented a maximum level of production, which defines the limit of sustainable production. According to them, this limit is reached when the water level is stable under its actual exploitation regime and, below this limit the production can be sustained for 100 years or more. Besides, Rybach (2003) argues a long productivity period is the main characteristic of geothermal sustainability. Both papers emphasize the need for permanent monitoring of the aquifer.
In this sense, at this stage it is difficult to define the sustainable limit of production of the Los Humeros. In this study estimation, the future production is achieved under a fixed operation limit. Furthermore, the visualization of the use of the reservoir by the wells is a useful tool to enhance the management of the site and ensures the future production.

Outlook
A full insight of the sustainability of the heat extraction includes an economic assessment plus an evaluation of the environmental impact. Based on the fluid extraction and the energy production, a profitability assessment can be carried out which could help with the management of the system. Finally, a life cycle assessment could close an integrated evaluation of the Los Humeros exploitation.

Conclusions
We presented a forecasting method of the feasible heat extraction of the Los Humeros and the estimation of the heat recovery factor. Our findings define a limit of the future heat production levels (assuming a constant extraction rate). Based on this limit, the definition of the original reserve of the system delimitates the extension of the technical potential of the Los Humeros.
These results offer: • A comprehensive overview of the future production limits. These values show the quantity needed in order to extract similar amounts of heat and maintain the power production. • The operation limits can also be combined with pressure monitoring to keep the extraction rates within the lowest production depletion. • Stakeholders can take this information to clarify the limitations of the production and therefore meet the risk associated with new investment in the system. • A data-based methodology, which is a good basis for the applying of machine learning methodologies.
Finally, the technical recovery factor states the maximum level of heat production attained for 30-year-old wells under a defined extraction rate. Although these results are not enough to claim something about the sustainability of the system, they open important topics to put the focus on. For example, the overlapping of circles of the resulting used area. The exploitation of the same aquifer by near wells could be an aggressive strategy, which might increase the pressure drop and compromising the heat production.