2.1. LOS HUMEROS GEOTHERMAL FIELD (LHGF)

LHGF is situated in Puebla, at the eastern end of the Mexican volcanic belt (19° 40' north latitude; 97° 25' west longitude), approximately 200 km from Mexico City, at an elevation of 3000 meters above sea level (as shown in Figure 1). LHGF marks the northern boundary of the Serdán-Oriental basin (Figure 1). This region is renowned for its bimodal, monogenetic volcanism, characterized by small, isolated cinder and scoria cones, large rhyolitic domes, and both rhyolitic and basaltic maars.

Exploratory studies in LHGF commenced in 1962, initiated by the Federal Electricity Commission (^{Mena and González, 1978}; ^{Perez, 1978}; ^{Yáñez et al., 1979}; ^{Palacios and Garcia, 1981}). The first well was drilled in 1982, and electricity production began in 1990 with the commissioning of the first 5 MW generators. Subsequently, more than 40 wells have been drilled, and by 1995, there were seven units, each with a capacity of 5 MW (^{Quijano and Torres, 1995}). Presently, the field has an installed capacity of 40 MW (^{Gutiérrez-Negrin, 2007}).

Since the early stages of exploration and production, numerous geochemical and engineering studies have been conducted on deposits in their natural state, confirming the presence of two deposits (^{Arellano et al., 1998}). The shallower deposit, located between 1600 and 1025 meters above sea level (MASL), is characterized by a liquid phase at temperatures ranging from 300°C to 330°C and a neutral pH. In contrast, the deeper deposit, situated below 850 MASL, exhibits convective behavior, with steam typically rising from the reservoir and condensing into fluid. This deposit has a low liquid saturation and operates at temperatures between 300 °C and 400 °C (^{Barragán et al., 2008}). LHGF, with the exception of well H-01, has a high enthalpy of production, and the fluid is predominantly vapour. Heat transfer dominates the fluids in the wells, which in some cases are fed from the shallowest deposit, sometimes from the deepest, and sometimes from both (^{Torres, 1995}).

Wells such as H-01, H-13, H-15, H-16, H-17, H-19, H-33, and H-24 (Figure 2), among others, experienced reduced production due to mineral precipitation and corrosion. In some instances, extended downtime was necessary for repairs (^{Cornejo, 1996}; ^{Gutiérrez and Viggiano, 1990}).

The calculation of flow within wells requires the consideration of momentum transport, which includes inertia and viscous forces between fluids and between fluids and the well walls. Flows within wells are characterized by high speeds, often turbulent behavior, significant pressure gradients, and phase changes. Therefore, the calculation of flow within wells must primarily yield realistic values for density, pressure, and temperature as functions of the position along the drilling shaft. The mathematical foundations for this computational development are outlined below.

3.1. ANALYTICAL APPROACH AND CALCULATIONS

In the present study, we employ the equations of mass, momentum, and energy conservation to characterize the flow within wells. There exist correlation-based and flow regime-based models that have proven effective in describing this type of flow, obviating the need for separate simulations of two- phase dynamics. We utilize a non-linear system of first-order differential equations, originally developed by ^{Gimon-Bastidas et al. (2018)} and integrate it with the algorithms in our study. This system simulates the kinetic dissolution and precipitation of minerals in aqueous electrolyte solutions, and it is applicable across a temperature range of 20 °C to 320 °C, pressures from 1 to 1000 bar, and ionic strengths of up to 6 mol/kg (Gimon-Bastidas et al., 2018). Their numerical tool provides accurate estimations of the following: a) the types and concentrations of minerals dissolved or precipitated over time; b) the ions and cations generated and released during the kinetic dissolution or precipitation of one or more minerals in the solution (Gimon-Bastidas, et al. 2018).

In this article, we simulate the wells using a direct modeling approach based on the work of ^{Gimon-Bastidas et al. (2018)}. We establish the conditions of pressure, temperature, enthalpy, and specific volume at the depth where amorphous silica filling occurs.

3.1.1. DETERMINATION OF THE TIME AT WHICH AN AMORPHOUS SILICA RING ADHERES TO THE INNER WALLS OF THE PIPELINES

As amorphous silica constitutes the predominant hydrothermal mineral phase at LHGF, silica is employed to model the incrustation processes. The precipitation of amorphous silica is most frequently encountered in pipe seals (^{Ocampo-Díaz et al., 2005}; ^{Miranda-Herrera, 2012}), and the literature extensively documents that the sealing of pipes with amorphous silica (^{Cornejo, 1996}; ^{Gutiérrez and Viggiano, 1990}) is the primary issue facing the geothermal industry (see Figure 3). Over time, the accumulation of amorphous silica in pipelines and fluid handling equipment proves detrimental to geothermal power production. The observations made in LHGF are included in Appendix A, providing ample observational data for adjusting and validating the proposed model. Studies conducted at the Hellisheidi geothermal power plant in SW-Iceland on solution composition, saturation indices, and silica precipitations have indicated silica precipitation rates exceeding 0.35-0.75 g m^-2 day^-1 (^{Van den Heuvel et al., 2018}). The determination of the time required for amorphous silica scale formation relies on the following assumptions: 1) The mass flow within the pipe section remains constant; 2) The temperature, pressure, and enthalpy within this section remain consistent; 3) The pipe maintains a uniform diameter; 4) In LHGF, the predominant phase undergoing precipitation is amorphous silica; 5) Scale formation occurs regularly and uniformly; 6) Kinetically precipitated amorphous silica originates from a film whose volume represents a percentage of the water flowing through the pipe; 7) The brine is in equilibrium before entering the pipe section, but after entering, the fluid's temperature decreases, leading to precipitation.

Figure 3 (a) Sealed pipe of the H-43 well at a depth of 2500 m, where the precipitation of amorphous silica can be observed (c and d), and amorphous silica with a hematite patina (b). Note the banded texture of progressive growth toward the center of the pipe. Boxes e, f, and g correspond to the areas where SEM-EDS analyses were performed, shown in microphotographs. SEM-EDS analyses reveal the compositions of the components: amorphous silica in e) and f); hematite in g) from the red patina (b). 

3.1.2. DETERMINATION OF TOTAL MOLES OF AMORPHOUS SILICA DEPOSITED PER KILOGRAM OF WATER

Brines found in geothermal wells contain a large amount of dissolved ionized salts, primarily chlorides and sulfates, which are highly corrosive (^{Valdez et al., 2009}). The chemical analysis of the brine extracted from the Humeros wells is presented in Table 1 (^{Arellano et al., 2001}). The pH of the various brines varies, ranging from acidic (pH 5) to basic (pH 8.9), depending on the location of the well and changes over time. In the case of well H-16 brine, with a pH of 5.1, reports of corrosion exist despite having samples with a neutral pH (^{Gutiérrez and Viggiano, 1990}). These solutions were used to determine the moles of kinetically deposited amorphous silica per kilogram of water Nt0. The characteristics of the rock in contact with the brine are presented in Table 2. The rock has a porosity of 15%.

It can take thousands of years for the interactions between brine and rock to reach a state of equilibrium with the surrounding minerals, such as potassium chloride, sodium chloride, etc., or saturation with certain salts like silica, bicarbonate, sulfide, and sulfate (^{Mercado et al., 1989}). ^{Pérez et al., (2012)} utilized the temperatures listed in Table 3 and the thermodynamic equilibrium code SPCALC to determine the state of equilibrium between the brine (Table 1) and the rock (Table 2). After calculating the concentrations of species at equilibrium and the amount of thermodynamically precipitated amorphous silica, the kinetic precipitation of the amorphous silica phase was measured at a temperature 5 °C lower than the equilibrium temperature. Following the determination of kinetically precipitated moles at each of the evaluated temperatures, correlations corresponding to each well were established (Table 3).

3.1.3. KINETIC PARAMETER CALCULATION Nt0

To determine the moles of kinetically precipitated amorphous silica Nt0 over a specified time, the algorithm utilizes the following input variables: 1) The concentration of all aqueous species obtained through speciation at equilibrium, as determined by a thermodynamic speciation code (^{Pérez et al., 2012}); 2) The mass of amorphous silica that has thermodynamically precipitated, also derived from the same thermodynamic speciation; 3) The assumed specific surface area of amorphous silica, set at 1000 cm²/g; 4) Temperature; 5) The activation energy and pre-exponential factor for the dissolution or precipitation of amorphous silica, as per ^{Palandri and Kharaka (2004)}.

Appendix B provides the equation defining the reaction rate for this system. The pressures used in both thermodynamic speciation and kinetic precipitation correspond to vapor pressure. Results from kinetic precipitation are presented in Figure 4 through 8 (5, 6, 7).

Figure 4 Kinetic precipitation of amorphous silica at a temperature of 290^o C for the different wells studied. 

Figure 5 Kinetic precipitation of amorphous silica at a temperature of 260^o C for well H-16. 

Figure 6 Kinetic precipitation of amorphous silica at a temperature of 280^o C for well H-16. 

Figure 7 Kinetic precipitation of amorphous silica at a temperature of 300^o C for well H-16. 

Figure 8 Kinetic precipitation of amorphous silica at different temperatures over a time of 2s and 4s for well H-16. 

The kinetics of amorphous silica precipitation for the various wells under study exhibit consistent behavior with increasing temperature. As temperature rises, the quantity of kinetically precipitated amorphous silica also increases. This trend is evident in Figure 4, where the ordinate axis is an order of magnitude greater than that in Figure 5.

For the H-16 well, an investigation into the effect of pH on the moles of kinetically precipitated amorphous silica was conducted. Differences in the moles of amorphous silica precipitated in solutions with varying pH levels increase at temperatures between 280 °C and 300 °C and decrease at lower (260 °O280 °C) and higher (> 300 °C) temperatures, as depicted in Figure 4, 5, 6, 7, 8. However, the sealing time for the H-16 well was determined based on the solution with a pH of 5.1, as this solution best represents the behavior of the well during the study period (^{Arellano et al., 2001}; ^{Cornejo, 1996}; ^{Gutierrez and Viggiano, 1990}).

The algorithm summarizing the methodology used to determine the value is depicted in Figure 9. For each of the wells, a correlation was established between the moles of kinetically precipitated amorphous silica and the variables of time and temperature at the sealing point (Table 3).

Figure 9 Algorithm used in the determination of kinetically precipitated moles of amorphous silica per kg of water. Aqueous speciation was performed for each of the chemical fluids of the wells studied in a temperature range of 260^o C to 300^o C. 

Another parameter required for input into the Nt0 calculation algorithm (Figure 9) is the residence time (t_res) of the fluid in a specific section, calculated as:

Where, V_W is the volume of water entering the pipe section under study in the first iteration and Q_W is the flow of water.

3.5. DETERMINING THE TIME AT WHICH THE AMORPHOUS SILICA RING FORMS

The two main assumptions regarding the shape of the scale and the origin of the kinetically precipitated amorphous silica are as follows: 1) The scale forms as a uniform cylindrical ring and adheres to the inner wall of the pipe (see Figure 3); 2)The amorphous silica that precipitates kinetically originates from a water film, the volume of which corresponds to a percentage of the total water present in the cylindrical section at the residence time of the water body; 3) This percentage of volume is known as the ‘k’ parameter and is an adjustable parameter.

In determining the time when the amorphous silica ring forms, the following input variables are utilized: a) Pressure conditions (p), temperature (T), enthalpy (H), and mass flow rate (W) in the section of the pipe under investigation. The first three conditions are obtained from well simulations, while the mass flow rate is reported by ^{Cornejo (1996)}, ^{Quijano and Torres (1995)}, ^{Torres (1995)}, and the Institute of Electrical Research; b) The percentage of pipe clogging (%Obt), which is either derived from literature data (^{Gutiérrez and Viggiano, 1990}) or determined based on increased pressure in the well; c) The percentage of the total volume that constitutes the water film (k). This parameter is adjustable; d) The internal radius of the pipe (R); e) Moles of kinetically deposited amorphous silica per kilogram of water Nt0; f) Water flow rate (Q_w), calculated using Equation (2); g) Liquid density (ρ_l), calculated using the function of temperature. The mass fraction of fluid vapor depends primarily on temperature, pressure, enthalpy, and is determined using ^{IAPWS tables (2007)}.

The radius obtained after ring formation (Equation 3) remains constant throughout the calculation of the ring formation time or shutter time:

Every time the cylindrical section is filled with liquid an iteration occurs. (i), at each iteration (i) The resulting radius is calculated (r_i) after the kinetic precipitation of amorphous silica using Equation 4. The algorithm for calculating training time of the amorphous silica ring, stops when:

The variation of the internal radius of the pipe (Δri-Equation 5) It is calculated from equalizing the volumes of the cylindrical ring formed by precipitated amorphous silica, calculated from the geometry of the ring (Equation 6) and the molar volume of amorphous silica (Equation 7):

Where, dx is the height of the cylindrical section, which is calculated as dx = L/n_x. L is the total length of the pipe and n_x is the number of subdivisions.

Moles of amorphous silica deposited (Nt,i) at each iteration (Equation 8):

Where Nt0 are the moles of amorphous silica per kilogram of water deposited kinetically and m_w,i is the body of water in the cylindrical section (Equation 9):

Where Vw,i is the volume of water that fills the cylin der at each iteration and is directly proportion al to the resulting radius (r_i) squared (Equation 10):

Where εg is the volumetric fraction of gas and k is the percentage of the volume of liquid that forms a film of water from which amorphous silica precipitates. In each iteration the volume of water is added (Vw,i) obtained to the tot volume of water passing through the cylindrical section in such a way that ri≤robt, as presented below:

Where N is the total number of iterations, Knowing that,

The time elapsed until the given shutter radius is reached (tobt) It clears itself of the equality presented in the equation (12). The algorithm used to determine the shutter time is summarized in Figure 10.

Figure 10 Determination of the time at which a certain sealing percentage is reached. 

4.1. MODEL ADJUSTMENT PREDICTING AMORPHOUS SILICA SCALE FORMATION TIME

The parameter ‘k’ was adjusted using results obtained from simulations of six wells: H-01, H-13, H-15, H-16, H-17, and H-19. The H-01 well commenced production in 1982 (^{Cornejo 1996}) and remained productive until 1994. In 1995, maintenance was performed on the pipe, and a seal was discovered between depths of 1295 m and 1450 m in its productive zone (^{Cornejo, 1996}). Notably, this is the only well where the sealing zone aligns with initial conditions where the fluid was entirely in the liquid phase (see Figure 11(c)). In contrast, Wells H-13, H-15, H-16, H-17, H-17D, H-19, H-20, H-24, H-32, and H-33 had seals at shallower depths, occurring before instantaneous vaporization of the fluid and resulting in significant temperature and pressure losses within the formation or grid (refer to Figure 11). The Los Humeros reservoir conditions, which have governed the productive behavior since its inception, are as follows: fluids within the region of subcooled liquid near the saturation curve, low permeability, and low porosity of the system. A decrease in pressure leads to a phase change within the rock formation. Los Humeros is characterized by a low-porosity system that allows for efficient heat transfer from the rock to the fluid. Consequently, there are enthalpy variations between the wellhead and the well bottom on the order of 1200 kJ/kg. Additionally, the low permeability results in high pressure gradients between the well bottom and the wellhead (^{Quijano and Torres, 1995}).

Thermodynamic simulations provide a close approximation to the reported values of wellhead pressure and enthalpy (^{Torres, 1995}; ^{Cornejo, 1996}), with errors of less than 15%. Notably, Wells H-01 and H-15 (see Figure 12(b)), which draw from a single reservoir (^{Torres, 1995}), exhibit differences in behavior. H-01, characterized by low temperatures and minimal pressure loss, experiences limited heat transfer between the fluid and the surrounding medium. Consequently, the enthalpy difference between the wellhead and the bottom is very similar, with a difference of merely 40 kJ/kg lower at the wellhead compared to the bottom. This well exhibits the highest fraction of liquid among those studied in the Humeros region.

Figure 12 Gas fraction profiles of the wells in the LHGF. (a) Wells fed from two reservoirs. (b) Wells fed by a single reservoir. 

On the other hand, H-15, being the deepest well and closest to the central collapse, has higher bottom temperatures, significant contact with low-porosity rock after instantaneous vaporization, and a considerable temperature difference with the surrounding rock, favoring additional heat transfer. This results in a significantly higher fraction of steam and, consequently, a greater enthalpy difference between the wellhead and the bottom (^{Quijano and Torres, 1995}), on the order of 1100 kJ/kg.

The enthalpy behavior and vapor fraction of Wells H-13, H-16, H-17, and H-19 (see Figure 12(a)) align with the hypothesis that these wells draw from two reservoirs (^{Arellano et al., 1998}; ^{Torres, 1995}). Vapor fraction plots reveal that the deepest reservoir, characterized by higher temperatures, experiences early boiling and high permeability, leading to a substantial pressure loss on the order of 90 bar (^{Quijano and Torres, 1995}). This reservoir does not exhibit vapor gain beyond the instantaneous vaporization point, with vapor fractions remaining below 0.05 (see Figure 12(a) and Figure 13(c)). However, when the well is also fed by the shallower deposit, vapor gain occurs at heights near the wellhead. This behavior is not typical of dry wells with biphasic fluids since their formation (^{Arellano et al., 2011}), as observed in H-20 and H-24, which draw from a single reservoir.

Figure 13 Profiles of the gas mass fraction of the wells in the LHGF. (a) Well H-33, fed by two reservoirs, used for model adjustment. (b) Well H-24, fed by a reservoir, used for model adjustment. (c) Wells H-17D and H-32, fed by two reservoirs, used to estimate operating time. (d) Well H-20, fed by a reservoir, used to estimate operating time. 

The results obtained by simulating wells drawing from one or two reservoirs indicate that the algorithms developed in this study can also replicate the behavior of these wells. This expands the applicability of the mathematical model of flow and heat due to the diversity of well types in the LHGF. Another phenomenon observed aligns with ^{Ocampo-Díaz et al., (2005)} descriptions. It explains that when liquid vaporization occurs, the concen tration of SiO₂ increases, leading to oversaturation and precipitation of amorphous silica. This phenomenon is observed in all the studied wells in the LHGF (see Figure 11).

Figure 12(b) illustrates that the H-01 well draws from a shallow reservoir with low temperatures (ranging from 290 °C to 315 °C in the Mastayola sector (^{Torres, 1995})). It exhibits behavior characteristic of wells with a gas fraction (X_g) below 0.15 and minimal enthalpy variation between the formation and the wellhead. Despite drawing from a single reservoir, H-15 experiences substantial vapor gain due to fluid evaporation. In Figure 11(a), the wells initially contain a significant amount of liquid in the first meters of travel, which tends to increase at heights coinciding with the feed from the shallower reservoir (between 1400 MASL and 700 MASL, according to ^{Torres 1995}). After a few meters, heat transfer promotes steam gain within the pipe, resulting in production enthalpies higher than background values, reaching up to 1200 kJ/kg.

In Figure 13(a), it becomes evident that the H-33 well draws from two deposits. Figure 13(b) demonstrates that the H-24 well is fed by a single reservoir but exhibits substantial vapor gain, particularly after the instantaneous vaporization point. Wells H-24 and H-33 were employed to validate the predictive model. Figure 13(c) illustrates that both H-17D and H-32 wells exhibit behavior akin to wells drawing from two reservoirs. Finally, Figure 13(d) shows the behavior of the H-20 well, which draws from a single reservoir. These three wells were used to estimate amorphous silica precipitation behavior under conditions similar to those used to adjust the ‘k’ parameter."

The percentage of clogging of the wells used in the parameter adjustment varies between 35% and 85%, while the time needed for the amorphous silica ring to form ranges from 2 to 12 years. For each of the six wells used in the adjustment of the model, a different value was obtained for k (Table 5). Once these values were obtained, it was determined that the variables that have the greatest influence on the parameter are: 1) the volumetric fraction of liquid εl and 2) the product Nt0⋅ρl. From these two variables, a correlation was constructed (Equation 12) that predicts the value of the parameter, knowing the volumetric fraction of the liquid, the density of the liquid, and the moles of amorphous silica that precipitate kinetically per kilogramme of water. Table 6 presents the results of the adjustment, and it can be seen that the time needed for the scale to form is predicted with errors less than 3%.

The time in which a certain percentage of seal is reached increases as the percentage of seal increases (Figure 14). In some wells, scale forms more quickly than in others; for example, wells H-15, H-16, H-17, and H-19 seem to reach a certain percentage of seal faster than wells H-01 and H-13. This has to do with the characteristics of each well. The results of this predictive model are validated by evaluating its performance by simulating wells H-24 and H-33, for which operating time data are available.

Figure 14 Variation in the percentage of plugging over time for the wells used in the adjustment of parameter k. 

4.2. VALIDATION OF THE MODEL THAT PREDICTS THE FORMATION TIME OF AMORPHOUS SILICA SCALE

^{Cornejo (1996)} provides information on the initiation and operational periods of wells H-24 and H-33. Both wells underwent direct modeling simulations, and the errors in predicting the time when scale formation occurs for H-24 and H-33 wells (as presented in Table 8) are below 3%. This demonstrates that the predictive model is accurate and effectively replicates the time at which well sealing occurs under the following conditions: 1) Temperatures within the range of 270°C to 305°C; 2) Enthalpies ranging from 1187 to 1400 kJ/kg; 3) Volumetric fractions of liquid εl >0.5 between 0.5 and 0.99; As seen in previous cases, with the passage of time, the percentage of sealing within the pipes increases (as illustrated in Figure 15). Additionally, we selected three wells that have not experienced closure due to sealing issues, but share temperature conditions similar to those of the wells under investigation. The results are summarized in Table 7, where the estimated sealing times range from 36 to 68 years under these conditions. This estimate serves as an additional validation of the model, as, at the time of this study, none of the wells had reached the sealing time predicted by the model (refer to Figure 16). H- 17D has been in operation since 1993 following repairs (^{Cornejo 1996}), H-20 has been operational since 1988, and H-32 since 1989 (^{Torres 1995}).

Figure 15 Changes in the percentages of sealing over time for the wells used in the validation of the correlation used to determine the parameter k. 

Figure 16 Wells used to verify the model through sealing time projections.