A seawater recharge type geothermal reservoir temperature attenuation prediction and dynamic evaluation system

Abstract

The application discloses a seawater recharge type geothermal reservoir temperature attenuation prediction and dynamic evaluation system, relates to the technical field of geothermal resource exploration and development and utilization, and comprises a seawater boundary module, an isotope correction module and an evaluation module. The seawater boundary module is used for collecting real-time monitoring data of a target geothermal reservoir and a near-shore sea area thereof, and identifying the instantaneous pressure conduction coefficient and the mixed proportion fluctuation range of the seawater recharge boundary under the driving of tidal dynamic load according to the real-time monitoring data. The seawater recharge type geothermal reservoir temperature attenuation prediction and dynamic evaluation system corrects the real-time value of the thermal storage property parameter from the chemical evolution mechanism level, and eliminates the process distortion of the equivalent thermal parameters of the thermal storage deviating from the real state caused by the one-way decoupling of the physical field and the chemical field.

Description

Technical Field

[0001] This invention relates to the field of geothermal resource exploration and development technology, specifically a system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs. Background Technology

[0002] In the field of geothermal resource exploration and development, accurate prediction of reservoir temperature decay and scientific assessment of reservoir dynamic stability under mining conditions have always been core technical aspects determining whether geothermal fields can achieve long-term sustainable development. For conventional deep hydrothermal geothermal systems, domestic and foreign scholars have developed relatively mature reservoir engineering evaluation methods. By establishing a coupled numerical model of heat conduction and convection, combined with long-term dynamic monitoring data, it is possible to provide relatively accurate predictions of temperature and pressure change trends in geothermal fields on a macroscopic scale. However, when the research object shifts to seawater-recharged geothermal systems that are commonly developed in coastal areas, the above-mentioned traditional prediction and assessment technology systems are facing fundamental applicability challenges. The essential characteristic that distinguishes seawater-recharged geothermal systems from inland land-based geothermal systems is that their reservoir boundaries are not closed or supplied by stable lateral runoff, but are directly or indirectly connected to open sea areas. Seawater periodically intrudes into the reservoir under the combined drive of density difference, head difference, and tidal dynamics, forming dynamic boundary conditions of alternating freshwater and saline water mixing. It is this unique hydraulic connectivity that causes the physicochemical properties of the fluids inside the thermal reservoir to exhibit spatiotemporal variability far greater than that of inland thermal systems.

[0003] However, a closer examination of existing prediction and assessment technologies reveals that, whether it's early analytical models or the widely used multi-field coupled numerical simulation platforms, almost without exception, they treat the seawater end-member as a stable recharge source with constant temperature and salinity when characterizing the recharge boundary. They also treat the chemical response of water-rock interaction as a late-stage superposition effect decoupled unidirectionally from the temperature and flow fields. This generalization of boundary conditions directly leads to the model's inability to effectively capture the nonlinear disturbances to the reservoir temperature field caused by the mixing of variable-temperature seawater during tidal cycles, and it also fails to reflect the shift and readjustment of the water-rock exchange equilibrium system caused by the gradual increase in the proportion of seawater mixing due to long-term depressurization during mining. Because the model lacks an adaptive response mechanism to dynamic changes in the recharge boundary, its medium- and long-term simulation results of the temperature decay curve often deviate from actual monitoring data, with prediction biases showing a cumulative amplification trend, severely weakening the practical guiding value of the assessment results for mining scheme formulation. At the same time, existing assessment systems generally follow an offline steady-state calculation paradigm, lacking the technical means to dynamically embed real-time monitoring data into the prediction model and perform rolling corrections. Therefore, they cannot implement graded early warning and dynamic control of reservoir temperature decay risks based on the current mining status.

[0004] In summary, the technical problem to be solved by this invention is how to construct a dynamic evaluation system that can fully consider the tidal dynamics of the seawater replenishment boundary and use the water-rock chemical response process as an embedded correction factor for temperature decay prediction. Summary of the Invention

[0005] The purpose of this invention is to provide a system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs, in order to solve the problems mentioned in the background art.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution: a system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs, including a seawater boundary module, an isotope correction module, and an evaluation module;

[0007] The seawater boundary module collects real-time monitoring data of the target thermal reservoir and its nearshore waters, identifies the fluctuation range of the instantaneous pressure transmission coefficient and mixing ratio of the seawater replenishment boundary under the driving force of tidal dynamic load, and outputs a set of dynamic boundary conditions modulated by the tidal cycle.

[0008] The isotope correction module receives isotope detection data from water samples from geothermal wells. Based on the isotope detection data and a pre-set isotope fractionation response library, it inverts the rate of change in the mixing ratio of freshwater and seawater within the geothermal reservoir and the resulting mineral dissolution-precipitation equilibrium shift. Then, based on the rate of change in the mixing ratio and the mineral dissolution-precipitation equilibrium shift, it calculates the dynamic reduction coefficient of the water-rock exchange process on the equivalent thermal conductivity and heat capacity of the geothermal reservoir. The dynamic reduction coefficient is then output as an internal parameter correction term of the temperature field prediction model.

[0009] The evaluation module is connected to the seawater boundary module and the isotope correction module respectively, receives the dynamic boundary condition set and the internal parameter correction term, performs multi-field coupled temperature decay prediction and outputs the evaluation result.

[0010] Preferably, the seawater boundary module includes a monitor and a spectrum analysis unit; the monitor is deployed in the key structural zone of the target geothermal reservoir and on the nearshore seabed to continuously collect geothermal water level, fluid temperature, fluid pressure, and real-time nearshore tide data; the spectrum analysis unit performs main frequency extraction and phase difference analysis on the long-period tidal effect response sequence and short-period tidal effect response sequence of the geothermal water level, and identifies the fluctuation range of the instantaneous pressure transmission coefficient and mixing ratio based on the analysis results.

[0011] Preferably, the spectrum analysis unit performs a discrete Fourier transform on the nearshore real-time tide level time series to extract the dominant frequency components of tidal level fluctuations and records them as the tidal fundamental frequency.

[0012] Preferably, the spectrum analysis unit performs a discrete Fourier transform on the geothermal water level time series, extracts the frequency components corresponding to the energy peaks in the frequency band near the tidal fundamental frequency, records them as the main frequency of the geothermal water level tidal response, and extracts the frequency components corresponding to the energy peaks in the integer multiples of the tidal fundamental frequency, records them as the multiples of the geothermal water level tidal response.

[0013] Preferably, the spectrum analysis unit subtracts the phase at the dominant frequency of the geothermal water level tidal response from the phase at the tidal fundamental frequency to obtain the fundamental frequency phase lag; subtracts the phase at the harmonic frequency of the geothermal water level tidal response from the phase at an integer multiple of the tidal fundamental frequency to obtain the harmonic phase lag; divides the amplitude at the dominant frequency of the geothermal water level tidal response from the amplitude at the tidal fundamental frequency to obtain the fundamental frequency amplitude attenuation ratio; and divides the amplitude at the harmonic frequency of the geothermal water level tidal response from the amplitude at an integer multiple of the tidal fundamental frequency to obtain the harmonic amplitude attenuation ratio.

[0014] Preferably, the seawater boundary module constructs a phase lag set by combining the fundamental frequency phase lag and the harmonic frequency phase lag, and constructs an amplitude attenuation set by combining the fundamental frequency amplitude attenuation ratio and the harmonic frequency amplitude attenuation ratio. The phase lag set and the amplitude attenuation set are then jointly solved with pre-stored formation elastic water storage coefficient, pressure conductivity coefficient, and fracture network spatial distribution data to solve for the attenuation rate constant of tidal pressure in the fracture network, which is output as the instantaneous pressure transmission coefficient. The module also reverse-calculates the time variation range of the ratio of the increase in seawater volume entering the thermal reservoir to the original freshwater volume in the thermal reservoir per unit tidal cycle, which is output as the mixing ratio fluctuation range.

[0015] Preferably, the isotope fractionation response library pre-set inside the isotope correction module is established based on rock and mineral separation experiments in the target area and water-rock exchange kinetic experiments under different temperature and pressure conditions, and records the strontium isotope fractionation coefficient, boron isotope fractionation coefficient and corresponding mineral dissolution-precipitation equilibrium constant under different temperature, salinity and water-rock ratio conditions.

[0016] Preferably, the isotope correction module uses the strontium isotope ratio and boron isotope ratio measured in the geothermal well water sample as search inputs. It matches the nearest neighboring seawater mixing ratio identifier and mineral dissolution-precipitation equilibrium constant identifier in the isotope fractionation response library to invert the actual mixing ratio of freshwater and seawater in the current geothermal reservoir. It then calculates the rate of change of the mixing ratio over time based on the sampling time interval. Based on the rate of change, it determines the driving direction and reaction intensity of the water-rock exchange reaction. Based on the mineral dissolution-precipitation equilibrium constant identifier, it calculates the difference between the amount of dissolution and precipitation of the main mineral phases to obtain the mineral dissolution-precipitation equilibrium offset.

[0017] Preferably, the isotope correction module maps the mineral dissolution-precipitation equilibrium offset to an equivalent thermal conductivity reduction factor and a heat capacity reduction factor, which together constitute the dynamic reduction factor; the equivalent thermal conductivity reduction factor and the heat capacity reduction factor exist in the form of dimensionless product factors and are synchronously output to the evaluation module within the same time step.

[0018] Preferably, the evaluation module includes a coupled solver, a back-substitution correction unit, and a threshold determination unit. The coupled solver performs an initial multi-field coupling operation based on the initial geological model and the current mining scheme, outputting a first predicted curve of temperature field decay over time. The back-substitution correction unit substitutes the first predicted curve back to the seawater boundary module to update the tidal boundary pressure transmission coefficient for the next iteration. Simultaneously, it feeds back the temperature decay trend of the first predicted curve to the isotope correction module to update the water-rock exchange reaction progress. After being updated by the back-substitution correction unit, the threshold determination unit triggers the coupled solver to perform a second multi-field coupling operation to obtain a second predicted curve. The combined deviation rate between the second predicted curve and the first predicted curve is compared with a preset threshold. If the combined deviation rate is higher than the preset threshold, a multi-stage closed-loop iteration is triggered until the combined deviation rate of the temperature decay curves output by two adjacent iterations is lower than the preset threshold, at which point the model is determined to have converged, and the final prediction result is output.

[0019] This invention provides a system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs. It offers the following advantages:

[0020] This invention replaces the static generalization approach of treating seawater end-members as a stable supply source with constant temperature and salinity by setting up a seawater boundary module. This is achieved by dynamically outputting a set of boundary conditions that are continuously updated with the tidal phase based on real-time tidal spectrum analysis. This allows the boundary input of the multi-field coupled solver to match the time-series changes in seawater intrusion pressure and mixing ratio under actual tidal drive, incorporating nonlinear disturbance information of temperature-variable seawater mixing within the tidal cycle from the initial prediction stage. Simultaneously, the isotope correction module converts the mixing ratio change rate and mineral dissolution-precipitation equilibrium offset obtained from strontium and boron isotope inversion into dynamic reduction coefficients of equivalent thermal conductivity and heat capacity. These are used as embedded correction terms in the temperature field prediction model to solve the governing equations. This transforms the water-rock chemical response process from a post-calculation superposition analysis after temperature decay calculations to an internal constraint that iterates synchronously with the physical field. This corrects the real-time values ​​of the reservoir's physical parameters at the chemical evolution mechanism level, eliminating process distortion caused by the unidirectional decoupling of the physical and chemical fields, which leads to deviations from the true state of the reservoir's equivalent thermal parameters.

[0021] This invention utilizes the collaborative work of the back-substitution correction unit and the threshold determination unit in the evaluation module to feed back the temperature decay rate sequence and temperature value sequence extracted from the temperature decay curve output in the current round to the seawater boundary module and the isotope correction module, respectively. This updates the tidal boundary pressure conduction coefficient and water-rock exchange reaction progress parameters for the next iteration, triggering a recalculation of the multi-field coupled solver. This ensures that the prediction curves in subsequent rounds carry comprehensive information after the previous round's results have been back-corrected for boundary conditions and chemical field responses. The threshold determination unit compares the deviation rate of corresponding nodes in adjacent rounds of prediction curves with a preset threshold round by round. Only when the deviation rate is below the threshold is convergence determined and the result output; otherwise, it continues iterating until the condition is met. The final output temperature decay curve is the convergence result after multiple rounds of boundary parameter self-correction and physical property parameter self-correction. This convergence curve, together with the calculated remaining exploitable years and the matched maximum allowable extraction volume scheme, constitutes complete output information for real-time control of the mining control platform, realizing the direct conversion of prediction results into mining intensity control commands. Attached Figure Description

[0022] Figure 1 This is a flowchart illustrating the seawater-recharged geothermal reservoir temperature decay prediction and dynamic evaluation system of the present invention.

[0023] Figure 2 This is the iterative convergence judgment logic diagram of the present invention;

[0024] Figure 3 This is a diagram of the internal structure of the evaluation module of the present invention. Detailed Implementation

[0025] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0026] Example: Please refer to Figures 1 to 3 This invention provides a system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs. In this system, a seawater boundary module, an isotope correction module, and an evaluation module are sequentially connected and operate collaboratively.

[0027] The seawater boundary module deploys monitors at key structural zones of the target geothermal reservoir and along the nearshore seabed. Each monitor integrates pressure, temperature, and water level sensors to continuously collect geothermal water level, fluid temperature, fluid pressure, and real-time nearshore tidal data. The collected multiple monitoring signals are synchronously transmitted to a spectrum analysis unit within the module. The spectrum analysis unit first performs time-frequency transformation on the long-period tidal effect response sequence of the geothermal water level, extracting the dominant frequency components and their corresponding amplitude and phase information. Simultaneously, the spectrum analysis unit performs the same time-frequency processing on the short-period tidal effect response sequence of the geothermal water level to obtain the frequency response characteristics under short-period disturbances.

[0028] Based on this, the spectrum analysis unit compares and analyzes the frequency response characteristics of long and short periods with the dominant frequency of nearshore real-time tidal data to calculate the phase lag and amplitude attenuation of geothermal water level fluctuations relative to tidal loads. Based on these phase lag and amplitude attenuation values, combined with pre-stored formation compressibility and hydraulic conductivity data for the target area, the spectrum analysis unit inversely calculates the pressure transmission rate as seawater travels along the recharge channel into the reservoir, thereby identifying the instantaneous pressure transmission coefficient at the seawater recharge boundary under tidal dynamic loads. Simultaneously, based on the ratio of geothermal water level fluctuation amplitude to tidal water level amplitude, and combined with pre-input formation porosity and fracture width distribution data, the spectrum analysis unit estimates the dynamic ratio between seawater intrusion and freshwater outflow in each tidal cycle, obtaining the range of seawater mixing ratio fluctuations.

[0029] After the above processing, the seawater boundary module arranges the fluctuation range of instantaneous pressure transmission coefficient and seawater mixing ratio according to time series, and outputs a dynamic boundary condition set that is continuously updated with the tidal phase. This set contains the time-varying sequence of seawater intrusion pressure and the time-varying sequence of seawater mixing ratio corresponding to the tidal phase.

[0030] The isotope correction module receives water samples periodically transported from the wellhead sampling device of the geothermal well. After preprocessing, the samples are sent to the isotope mass spectrometry unit to determine the strontium and boron isotope ratios in the water samples, forming isotope detection data of the water samples under the current mining conditions. This module has a pre-built isotope fractionation response library, established based on previous rock-mineral separation experiments and water-rock exchange kinetic experiments under different temperature and pressure conditions in the target area. This library records the fractionation coefficients of strontium and boron isotopes and the corresponding mineral dissolution-precipitation equilibrium constants under different temperatures, salinities, and water-rock ratios.

[0031] The isotope correction module uses the strontium and boron isotope ratios from the water sample isotope detection data as input. It performs a matching search in a pre-set isotope fractionation response library to deduce the actual mixing ratio between freshwater and seawater within the current hydrothermal reservoir. Simultaneously, it calculates the rate of change of the mixing ratio over time, taking into account the time interval between the current and previous sampling. A rate greater than zero indicates an increasing seawater mixing ratio, while a rate less than zero indicates a decreasing ratio. Based on this, the isotope correction module determines the driving direction and intensity of the water-rock exchange reaction at the current stage, and further calculates the difference between the dissolution and precipitation of the main mineral phases in the reservoir at the current mixing ratio, using the corresponding mineral dissolution-precipitation equilibrium constant from the pre-set isotope fractionation response library. This mineral dissolution-precipitation equilibrium offset reflects the changing trend of the reservoir's pore space and fluid channels due to the altered seawater mixing ratio.

[0032] Subsequently, the isotope correction module maps the mineral dissolution-precipitation equilibrium offset to reduction factors for the equivalent thermal conductivity and heat capacity of the thermal reservoir. Specifically, when mineral dissolution dominates, porosity increases, equivalent thermal conductivity decreases, and heat capacity also decreases; when mineral precipitation dominates, porosity decreases, equivalent thermal conductivity increases, and heat capacity also increases. The reduction factors for equivalent thermal conductivity and heat capacity together constitute the dynamic reduction factor. The isotope correction module uses this dynamic reduction factor as a correction term for the internal physical property parameters of the temperature field prediction model and outputs it to the evaluation module according to a time step synchronized with the tidal phase.

[0033] The evaluation module is internally configured with a coupled solver, a back-substitution correction unit, and a threshold determination unit. The coupled solver embeds a multi-field solver, which uses the heat transfer equation, fluid dynamics equation, and chemical reaction kinetics equation as the governing equations. It uses the stratigraphic geometry and lithological parameters stored in the initial geological model as the basic inputs, the current mining flow rate and pressure as the operating condition inputs, the dynamic boundary condition set output by the seawater boundary module as the solution boundary conditions, and the dynamic reduction coefficient output by the isotope correction module as the adjustment values ​​for physical property parameters at each calculation step during the solution process, performing the first multi-field coupled numerical calculation. After the first calculation, the coupled solver outputs the first predicted curve of the reservoir temperature field decaying with mining time. After completing the aforementioned back-substitution correction, the coupled solver re-executes the second multi-field coupled numerical calculation using the corrected instantaneous pressure conduction coefficient and the corrected water-rock exchange reaction progress parameters to obtain the second predicted curve. The second predicted curve is then sent to the threshold determination unit. The threshold determination unit stores a preset threshold, which is determined by historical monitoring data of the target area and the critical value for safe mining calibrated by numerical experiments.

[0034] The threshold determination unit calculates the absolute value of the temperature difference between the second and first prediction curves at corresponding time points, divides it by the temperature value at the same point on the first prediction curve, and obtains a deviation rate sequence. The maximum value in this deviation rate sequence is compared with a preset threshold. When the maximum value is lower than the preset threshold, the threshold determination unit determines that the model has converged, outputs the second prediction curve as the final prediction curve, and calculates the remaining recoverable years of the reservoir under the current mining scheme based on the ratio of the difference between the final prediction curve's end temperature value and the initial temperature value, combined with pre-stored reservoir volume and specific heat capacity parameters. It then uses this remaining recoverable years to back-match a pre-stored extraction volume-temperature decay correspondence table and outputs the corresponding dynamic control scheme for the maximum allowable extraction volume. When the maximum value is not lower than the preset threshold, the threshold determination unit sends the second prediction curve as a new back-substitution input back to the back-substitution correction unit, repeating the above back-substitution correction and re-calculation process until the maximum deviation rate between the temperature decay curves output from two adjacent iterations is lower than the preset threshold. The final prediction curve and the corresponding dynamic control scheme for the maximum allowable extraction volume, output after the above closed-loop iteration and threshold determination, are transmitted to the geothermal extraction control platform to adjust the actual extraction flow rate and extraction pressure.

[0035] After receiving the first prediction curve, the back-substitution correction unit executes two parallel back-substitution correction paths. In the first back-substitution path, the back-substitution correction unit extracts the temperature decay rate sequence from the first prediction curve, i.e., the ratio sequence of temperature difference between adjacent output nodes to time intervals, and sends this temperature decay rate sequence as feedback to the spectrum analysis unit of the seawater boundary module. After receiving this feedback, the spectrum analysis unit uses the temperature decay rate sequence as a constraint condition to readjust the weight allocation of phase lag and amplitude decay of the long-period tidal effect response sequence and the short-period tidal effect response sequence, so that the identification result of the instantaneous pressure conduction coefficient in the next iteration matches the current temperature decay rate. In the second back-substitution path, the back-substitution correction unit extracts the temperature value sequence corresponding to each time output node from the first prediction curve and sends this temperature value sequence as a temperature constraint vector to the isotope correction module. After receiving the temperature constraint vector, the isotope correction module uses the temperature value corresponding to each time node as the search condition to rematch the strontium isotope fractionation coefficient and boron isotope fractionation coefficient at the corresponding temperature in the preset isotope fractionation response library. Based on the rematched fractionation coefficient, it updates the mixing ratio change rate and mineral dissolution-precipitation equilibrium offset under the current iteration round, thereby refreshing the specific values ​​of the effective thermal conductivity reduction term and heat capacity reduction term in the dynamic reduction coefficient.

[0036] After completing the two back-substitution correction paths, the back-substitution correction unit sends a recalculation command to the coupled solver. The coupled solver uses the instantaneous pressure transmission coefficient corrected by the first back-substitution path and the dynamic reduction coefficient corrected by the second back-substitution path as updated input parameters, keeping other parameters in the initial geological model, mining scheme, and dynamic boundary condition set unchanged, and performs a second multi-field coupled numerical operation to obtain the second prediction curve.

[0037] After the second predicted curve is generated, the coupled solver transmits it to the threshold determination unit. Upon receiving the first and second predicted curves, the threshold determination unit compares the temperature values ​​at corresponding nodes on the two curves one by one according to the same time output nodes. The threshold determination unit internally stores a preset threshold, which is a dimensionless boundary value determined based on historical monitoring data of the target area and the critical value for safe mining calibrated by numerical experiments. The threshold determination unit calculates the absolute value of the difference between the temperature value of the second predicted curve and the temperature value of the first predicted curve at each time node, divides this absolute value by the temperature value of the same node on the first predicted curve to obtain the deviation rate of each node, and then takes the maximum value of the deviation rates of each node as the comprehensive deviation rate of the current iteration round. This comprehensive deviation rate is the deviation rate between iteration rounds.

[0038] The threshold determination unit compares the comprehensive deviation rate with a preset threshold: when the comprehensive deviation rate is less than the preset threshold, the threshold determination unit determines that the model has converged, marks the second predicted curve as the final predicted curve, and terminates the iteration process; when the comprehensive deviation rate is not less than the preset threshold, the threshold determination unit determines that the model has not yet converged, and sends the second predicted curve, along with its temperature decay rate sequence and temperature value sequence, back to the back-substitution correction unit. The back-substitution correction unit then performs the back-substitution correction operation again according to the first and second back-substitution paths described above, and subsequently triggers the coupled solver to perform the next multi-field coupled numerical operation. This process is repeated until the comprehensive deviation rate between the temperature decay curves output by two adjacent iterations is less than the preset threshold, at which point the iteration terminates. In each iteration, the back-substitution correction unit uses the predicted curve output by the latest round as the sole input for the next round of back-substitution; the historical back-substitution parameters from the previous round are not included in the correction calculation for the current round.

[0039] The threshold determination unit in the evaluation module stores preset thresholds. These preset thresholds are determined as follows: Before the system is put into operation, a retrospective analysis is conducted on existing historical geothermal extraction monitoring data for the target area, including long-term continuous records of geothermal water levels, fluid temperatures, and fluid pressures at different extraction stages, as well as nearshore tidal observation data for the corresponding time periods. From this historical monitoring data, the deviation sequence between the actual temperature decay curve of the geothermal reservoir and the theoretical decay curve simulated according to the initial boundary conditions within multiple completed extraction cycles is extracted. The distribution characteristics of this deviation sequence are statistically analyzed, and the upper limit of the deviation sequence is used as the first reference boundary.

[0040] Simultaneously, in the numerical experimental environment, multiple sets of forward modeling examples under different mining intensities were built based on the geological model of the target area. Each set of examples was set with different seawater mixing ratio growth rates and different initial reservoir temperatures. The simulation was stopped when the reservoir temperature dropped to the minimum economic temperature for geothermal utilization. The temperature decay curve deviation evolution trajectory of each set of examples before reaching the minimum economic temperature was recorded. The maximum allowable deviation value that ensures the reservoir temperature does not fall below the minimum economic temperature within the preset service life was taken as the second reference limit. The threshold determination unit determined the smaller value between the first reference limit and the second reference limit as the preset threshold and stored it in the internal register of the threshold determination unit in dimensionless numerical form.

[0041] When the threshold determination unit determines that the model has converged, that is, when the comprehensive deviation rate between the temperature decay curve output in the current iteration and the temperature decay curve output in the previous iteration is less than the aforementioned preset threshold, the threshold determination unit outputs a convergence confirmation signal to the system. This convergence confirmation signal triggers a synchronous output process: the threshold determination unit reads the temperature decay curve output in the final convergence iteration, compares the average temperature value of the geothermal reservoir corresponding to the end time node of the curve with the original geothermal reservoir temperature value stored in the initial geological model, and obtains the total temperature reduction; then, based on the temperature change rate on the temperature decay curve, it calculates the time span required for continued mining from the current end time node until the temperature drops to the minimum economic temperature for geothermal utilization, and outputs this time span as the remaining recoverable years of the geothermal reservoir.

[0042] Simultaneously, the threshold determination unit matches the remaining recoverable years of the aforementioned geothermal reservoir with a pre-stored table of extraction volume-temperature decay correspondence within the system. This table was established during the numerical experiment phase by running forward simulations on the geological model of the target area for different combinations of extraction flow rate and extraction pressure levels. Each combination of extraction flow rate and extraction pressure in the table corresponds to a simulated recoverable year when the temperature decays to the minimum economic temperature. The threshold determination unit uses the currently calculated remaining recoverable years of the geothermal reservoir as the target value, retrieves the extraction flow rate and extraction pressure combination corresponding to the recoverable years closest to the target value from the table, determines the extraction flow rate value in this combination as the maximum allowable extraction volume, and outputs it in the form of an extraction flow rate control command. This maximum allowable extraction volume is used to constrain the upper limit of the frequency converter output of the geothermal extraction pump during actual extraction operations. When the actual extraction flow rate exceeds this maximum allowable extraction volume, the extraction control platform triggers a current limiting protection action.

[0043] Following the aforementioned synchronous output process, the threshold determination unit ultimately outputs three pieces of information to the outside of the system: the final prediction curve after convergence confirmation, the remaining recoverable years of the thermal reservoir, and the dynamic adjustment scheme for the maximum allowable extraction volume. The dynamic adjustment scheme for the maximum allowable extraction volume includes a stepped extraction volume adjustment recommendation sequence corresponding to different prediction time nodes.

[0044] The seawater boundary module consists of two parts: a monitor and a spectrum analysis unit, which are connected by a data acquisition cable or a wireless data transmission link.

[0045] The monitoring system comprises multiple downhole monitoring modules deployed within the key structural zones of the target geothermal reservoir, and subsea monitoring modules deployed on the nearshore seabed. The downhole monitoring modules are installed in pre-drilled observation wells within the geothermal reservoir. Each module integrates a water level sensor, a temperature sensor, and a pressure sensor, with the sensing ends of all three probes located within the main water-producing section of the reservoir. The water level sensor detects real-time changes in groundwater level within the observation well, the temperature sensor detects real-time temperature of the reservoir fluids, and the pressure sensor detects real-time pore water pressure. The subsea monitoring modules are installed at a predetermined depth below the nearshore seabed surface and integrate a tide level sensor unit for continuously sensing real-time tide levels in the nearshore waters. Both the downhole and subsea monitoring modules synchronously acquire their respective physical quantities at a set sampling frequency. The acquired signals are converted from analog to digital and uploaded in real-time as data frames.

[0046] The spectrum analysis unit receives geothermal water level time series, fluid temperature time series, fluid pressure time series, and nearshore real-time tide level time series uploaded by the monitor. The spectrum analysis unit first performs a discrete Fourier transform on the nearshore real-time tide level time series to extract the dominant frequency component of tidal level fluctuations, and records this dominant frequency component as the tidal fundamental frequency.

[0047] Subsequently, the spectrum analysis unit performs the same time-frequency transformation on the geothermal water level time series, extracts the frequency components corresponding to the energy peaks located in the frequency band near the tidal fundamental frequency from the transformation results, and uses them as the dominant frequency of the geothermal water level tidal response; at the same time, the spectrum analysis unit further extracts the frequency components corresponding to the energy peaks located in the integer multiples of the tidal fundamental frequency from the transformation results of the geothermal water level time series, and uses them as the multiples of the geothermal water level tidal response.

[0048] Based on this, the spectrum analysis unit subtracts the phase at the main frequency of the geothermal water level tidal response from the phase at the tidal fundamental frequency to obtain the fundamental frequency phase lag of the geothermal water level relative to the tidal water level; and subtracts the phase at the harmonic frequency of the geothermal water level tidal response from the phase at an integer multiple of the tidal fundamental frequency to obtain the harmonic frequency phase lag of the geothermal water level relative to the tidal water level.

[0049] Meanwhile, the spectrum analysis unit calculates the ratio of the amplitude at the main frequency of the geothermal water level tidal response to the amplitude at the tidal fundamental frequency to obtain the fundamental frequency amplitude attenuation ratio; it also calculates the ratio of the amplitude at the harmonic frequency of the geothermal water level tidal response to the amplitude at integer multiples of the tidal fundamental frequency to obtain the harmonic frequency amplitude attenuation ratio.

[0050] The aforementioned fundamental frequency phase lag and harmonic phase lag together constitute the phase lag set, and the aforementioned fundamental frequency amplitude attenuation ratio and harmonic amplitude attenuation ratio together constitute the amplitude attenuation set. The spectrum analysis unit performs joint calculations on the phase lag set and amplitude attenuation set with the pre-stored formation elastic water storage coefficient, pressure conductivity coefficient, and fracture network spatial distribution data within the system.

[0051] Specifically, the spectrum analysis unit calculates the propagation time required for pressure waves to travel from the nearshore seabed to the location of the thermal reservoir observation well based on the correspondence between the phase lag value and the tidal period. Then, it calculates the actual propagation speed of the pressure waves in the fractured medium by combining the propagation path length. Based on the deviation between the actual propagation speed and the theoretical elastic wave speed, and combined with the pre-stored fracture width statistics, it calculates the attenuation rate constant of tidal pressure in the fracture network and outputs this attenuation rate constant as the instantaneous pressure transmission coefficient.

[0052] Simultaneously, based on the amplitude attenuation ratio and the spatial distribution data of formation porosity, the spectrum analysis unit inversely calculates the ratio between the increase in seawater volume entering the geothermal reservoir through the recharge channel per unit tidal cycle and the original freshwater volume within the reservoir. The range of this ratio's variation over time is output as the mixing ratio fluctuation range. The spectrum analysis unit discretizes and arranges the instantaneous pressure transmission coefficient and mixing ratio fluctuation range according to the time step corresponding to the tidal fundamental frequency, using integer multiples of the tidal fundamental frequency as time nodes, forming a set of serialized boundary condition data that is continuously updated with the tidal phase—this is the dynamic boundary condition set.

[0053] The isotope correction module is equipped with an isotope fractionation response library. This library was pre-established for geothermal reservoir rock samples and regional seawater samples in the target area before the system was put into use.

[0054] The process of establishing the isotope fractionation response library is as follows: Core samples at different depths were collected from typical geothermal exploration wells in the target area, seawater samples at different tidal stages were collected from the nearshore waters of the area, and geothermal water samples were collected from the outflow wells of the target geothermal field. The collected core samples were classified according to their mineral composition, and single mineral components such as plagioclase, potassium feldspar, quartz, calcite, and clay minerals were separated.

[0055] For each of the aforementioned single mineral components, multiple temperature gradient conditions and multiple pressure gradient conditions were set up. Each temperature gradient and each pressure gradient was combined to form multiple sets of temperature and pressure experimental conditions. Under each set of temperature and pressure experimental conditions, the single mineral component was mixed with seawater samples and geothermal water samples at different solid-liquid mass ratios, and water-rock exchange kinetic experiments were conducted in a closed reaction vessel.

[0056] During the experiment, the reaction solution was extracted from the reaction vessel at set time intervals, and the strontium and boron isotope ratios in the reaction solution were measured. Simultaneously, the strontium and boron concentrations in the reaction solution were also measured. Based on the measurement results, for each mineral type, temperature condition, pressure condition, solid-liquid mass ratio, and reaction time, the strontium isotope fractionation coefficient, boron isotope fractionation coefficient, and the corresponding mineral dissolution rate constant and mineral precipitation rate constant were recorded. All the above records were organized and stored according to a multi-dimensional index structure of mineral type-temperature-pressure-solid-liquid ratio-reaction time, forming an isotope fractionation response library.

[0057] During actual system operation, the isotope correction module receives water samples from the wellhead sampling device of the geothermal extraction well. The wellhead sampling device automatically collects water samples from the outlet of the geothermal extraction well at set time intervals. After on-site filtration and acidification, the water samples are sent to the isotope mass spectrometry detection unit configured within the module. The isotope mass spectrometry detection unit performs thermal ionization mass spectrometry analysis and multi-receiver inductively coupled plasma mass spectrometry analysis on the water samples, determining the strontium isotope ratio and boron isotope ratio, respectively. The strontium isotope ratio is recorded as the ratio of strontium-87 to strontium-86, and the boron isotope ratio is recorded as the ratio of boron-11 to boron-10.

[0058] After the measurement is completed, the isotope correction module uses the strontium and boron isotope ratios of the current water sample as search inputs to find the experimental record closest to the current measurement value in the isotope fractionation response library using a nearest neighbor matching algorithm. During the matching process, the isotope correction module first performs a coarse screening based on the major order of magnitude of the strontium isotope ratio to obtain a candidate record set, and then performs a fine screening based on the boron isotope ratio within the same candidate record set to obtain the optimal matching record. From this optimal matching record, the seawater mixing ratio identifier and the mineral dissolution-precipitation equilibrium constant identifier under the corresponding conditions are extracted to complete the inversion analysis of the current water sample's isotope data. This inversion analysis result serves as the basic input data for the isotope correction module to subsequently calculate the mixing ratio change rate and the mineral dissolution-precipitation equilibrium offset.

[0059] The evaluation module has a coupled solver with an embedded solver.

[0060] The solver uses the heat transfer equation, fluid dynamics equation, and chemical reaction kinetics equation as the governing equations. The heat transfer equation describes the conduction, convection, and exchange of heat between the solid-phase framework and the liquid-phase fluid within the geothermal reservoir. It is a heat energy conservation equation with a time derivative term, where the heat conduction term is expressed using Fourier's law, the heat convection term is expressed as the product of the fluid's specific heat capacity and Darcy's velocity, and the heat exchange between the solid and liquid phases is introduced into the source term on the right-hand side of the equation in the form of the volumetric heat transfer coefficient. The fluid dynamics equation describes the flow of geothermal water within the geothermal reservoir's fracture network. It is a combination of mass and momentum conservation equations with time derivative terms. The mass conservation equation is expressed using the continuity equation for compressible fluids, and the momentum conservation equation is expressed using a generalized form of Darcy's law. The permeability tensor in the equations is dynamically updated according to changes in the geothermal reservoir porosity. Chemical reaction kinetic equations are used to describe the rates and equilibrium relationships of dissolution and precipitation reactions between geothermal fluids and geothermal reservoir rocks and minerals. The equations use the volume fraction of each mineral phase as state variables, the reaction rate constant and ion activity product as control parameters, and the reaction progress as output variables.

[0061] In the heat transfer equations described above, the rate of change of temperature over time is jointly determined by the heat conduction, convection, and solid-liquid heat exchange terms. The heat conduction term is related to the spatial second derivative of temperature, the convection term is related to the fluid velocity vector and temperature gradient, and the solid-liquid heat exchange term is related to the temperature difference between the two phases and the heat transfer area. In the fluid dynamics equations, the rate of change of fluid pressure over time is jointly determined by the change in fluid density and flow divergence, while the fluid velocity is jointly controlled by the pressure gradient and permeability tensor. In the chemical reaction kinetics equations, the rate of change of the volume fraction of each mineral phase is jointly controlled by the ratio of the reaction rate constant, the ion activity product, and the equilibrium constant. The greater the deviation of this ratio from the equilibrium state, the faster the reaction rate.

[0062] During operation, the solver constructs the three governing equations into a coupled system. The temperature field calculation results of the heat transfer equation affect the chemical reaction kinetic equation through the temperature-dependent reaction rate constant. The output of the chemical reaction kinetic equation is fed back to the permeability parameter in the fluid dynamics equation by changing the porosity through the amount of mineral dissolution and precipitation. The output of the fluid dynamics equation affects the distribution of the convection term in the heat transfer equation through the direction and velocity of fluid flow. A two-way coupling relationship is formed between the three equations, and the output variable of any one equation is used as part of the input parameters of the other two equations in the calculation.

[0063] When the system actually performs prediction calculations, the solver first reads the initial geological model stored inside the system and obtains the spatial geometric grid information of the thermal reservoir, the lithological classification identifier of each grid cell, the initial temperature field distribution, the initial pressure field distribution, and the initial porosity distribution from the model.

[0064] Subsequently, the solver applies the set of dynamic boundary conditions output by the seawater boundary module to the outer boundary of the solution domain and the inner boundary where the supply channel is located in sequence according to the time step. At each time step, the boundary assignment is updated with the seawater intrusion pressure value and seawater mixing ratio value corresponding to the current time.

[0065] Simultaneously, the solver applies the dynamic reduction coefficients output by the isotope correction module to the equivalent thermal conductivity calculation unit and the heat capacity calculation unit at each time step: In the equivalent thermal conductivity calculation unit, the equivalent thermal conductivity reduction term in the dynamic reduction coefficient is multiplied by the initial equivalent thermal conductivity in the form of a product factor to obtain the actual equivalent thermal conductivity of each grid cell at the current time step; in the heat capacity calculation unit, the heat capacity reduction term in the dynamic reduction coefficient is multiplied by the initial volumetric heat capacity in the form of a product factor to obtain the actual volumetric heat capacity of each grid cell at the current time step.

[0066] The actual equivalent thermal conductivity and actual volumetric heat capacity mentioned above replace the physical property parameter values ​​stored inside the solver immediately after each time step update, and are used as the reference values ​​for the calculation of the next time step.

[0067] Following the aforementioned coupling mechanism, the solver sequentially performs temperature field, pressure-velocity field, and reaction progress field solutions within each time step. These three solutions iterate internally within the same time step according to their respective convergence criteria until the residuals within that time step meet the convergence requirements, before proceeding to the next time step. After all time steps are completed, the solver summarizes the spatial distribution data of the temperature field corresponding to each time step, extracts the overall temperature variation sequence of the thermal reservoir over time calculated using a spatial averaging method, and outputs a temperature decay prediction curve.

[0068] The seawater boundary module is equipped with a boundary condition arrangement unit. After the spectrum analysis unit completes the identification of the instantaneous pressure transmission coefficient and the estimation of the mixing ratio fluctuation range for each tidal cycle, the arrangement unit serializes the above two parameters according to the tidal phase to generate a dynamic boundary condition set.

[0069] The tidal phase is defined based on the tidal fundamental frequency of the nearshore real-time tide level time series. The boundary condition arrangement unit first obtains the start and end times corresponding to the tidal fundamental frequency within the current tidal cycle from the spectrum analysis unit. The time span between the start and end times is divided into multiple equally spaced time nodes according to the set discretization step size. Each time node corresponds to a tidal phase value. The tidal phase value starts from the zero phase at the start time and increases along the time axis to the complete periodic phase at the end time.

[0070] For each of the above time points, the boundary condition orchestration unit performs the following serialization process:

[0071] For the time-varying sequence of seawater intrusion pressure, the boundary condition arrangement unit reads the instantaneous value of the conduction coefficient corresponding to the time node from the instantaneous pressure conduction coefficient output by the spectrum analysis unit. Simultaneously, it extracts the tide height value at the same time node from the nearshore real-time tide level time series, converts this tide height value into an equivalent seawater column pressure value, and then multiplies it by the instantaneous conduction coefficient value to obtain the seawater intrusion pressure value at that time node. This process is repeated for all time nodes within the current tidal cycle, arranging each time node and its corresponding seawater intrusion pressure value in chronological order to form a time-varying sequence of seawater intrusion pressure that continuously changes with time within the current tidal cycle.

[0072] For the time-varying sequence of seawater mixing ratio, the boundary condition orchestration unit reads the instantaneous estimated value of the mixing ratio corresponding to the time node from the mixing ratio fluctuation range output by the spectrum analysis unit. Simultaneously, it reads the effective connected pore volume corresponding to the current tidal phase at that time node from the pre-stored geothermal reservoir pore volume distribution data. The instantaneous estimated value of the mixing ratio is correlated and converted with the effective connected pore volume to obtain the instantaneous proportion of the increase in seawater volume entering the geothermal reservoir through the replenishment channel to the total fluid volume within the geothermal reservoir at that time node. This process is repeated for all time nodes within the current tidal cycle, arranging each time node and its corresponding instantaneous proportion value in chronological order to form a time-varying sequence of the seawater mixing ratio that continuously changes over time within the current tidal cycle.

[0073] The boundary condition arrangement unit aligns and arranges the aforementioned time-varying sequences of seawater intrusion pressure and seawater mixing ratio according to the same tidal phase time nodes. This ensures that the corresponding seawater intrusion pressure value and seawater mixing ratio value are simultaneously present at the same time node, thus forming a complete pair of boundary condition sequences. This set of dynamic boundary conditions for the current tidal cycle is then output to the evaluation module. When the next tidal cycle arrives, the boundary condition arrangement unit repeats the above process, generating a new set of dynamic boundary conditions to replace the set of the previous cycle, achieving real-time rolling updates of boundary conditions with the tidal phase.

[0074] The isotope correction module is equipped with a dynamic reduction factor calculation unit. This calculation unit receives the mixing ratio change rate and mineral dissolution-precipitation equilibrium offset obtained from the isotope mass spectrometry detection unit as input data, and performs a time-step calculation process for the dynamic reduction factor.

[0075] The rate of change of the mixing ratio indicates the direction and speed of the increase or decrease in the proportion of seawater mixed in the reservoir under the current mining conditions. The dynamic reduction coefficient calculation unit uses the positive or negative sign of this rate of change of the mixing ratio as the basis for determining the driving direction of the water-rock exchange reaction: when the rate of change of the mixing ratio is positive, it indicates that the proportion of seawater mixed in is in the rising stage, the salinity and ionic composition of the fluid in the reservoir are continuously changing, breaking the original water-rock chemical equilibrium state, and driving the mineral dissolution reaction in the positive direction. At this time, the driving force points in the direction of mineral dissolution. When the rate of change of the mixing ratio is negative, it indicates that the proportion of seawater mixed in is in the falling stage, the fluid salinity tends to decrease, and the supersaturated state drives the mineral precipitation reaction to dominate. At this time, the driving force points in the direction of mineral precipitation. At the same time, the absolute value of the rate of change of the mixing ratio is used to determine the intensity level of the driving force. The larger the absolute value, the greater the change in the mixing ratio per unit time, the deeper the water-rock exchange reaction deviates from the equilibrium state, and the stronger the driving force of the reaction.

[0076] The mineral dissolution-precipitation equilibrium offset records the difference between the total amount of dissolution and precipitation of each major mineral phase in the thermal reservoir at the current mixing ratio. After reading this offset, the dynamic reduction factor calculation unit extracts the dissolution and precipitation values ​​for each mineral type. For silicate minerals, represented by plagioclase and potassium feldspar, the dissolution process releases silicate ions and metal cations, while consuming hydrogen ions in the fluid, leading to changes in fluid acidity and increased pore space. For carbonate minerals, represented by calcite, the dissolution process releases calcium ions and bicarbonate ions, resulting in increased porosity. For secondary minerals, represented by clay minerals, the precipitation process fills pore throats, leading to decreased porosity. The dynamic reduction factor calculation unit weights and summarizes the dissolution and precipitation amounts of each mineral type according to their respective volume fractions to obtain the direction and magnitude of change in the overall porosity of the thermal reservoir. This direction of change is consistent with the sign of the mineral dissolution-precipitation equilibrium offset, and the magnitude of change is positively correlated with the absolute value of the offset.

[0077] After determining the change in porosity, the dynamic reduction coefficient calculation unit converts the change in porosity into the direction and magnitude of permeability change based on the functional relationship between porosity and permeability. Specifically, the conversion is as follows: permeability increases when porosity increases, and permeability decreases when porosity decreases. The magnitude of permeability change and the magnitude of porosity change follow a predefined nonlinear mapping relationship. These changes in porosity and permeability together constitute the dynamic changes in the physical properties of the thermal reservoir.

[0078] The dynamic reduction factor calculation unit then converts the dynamic change of this physical property into a reduction factor for thermal parameters. Specifically, the conversion is as follows: the change in porosity is mapped to the rate of change of equivalent thermal conductivity; the change in permeability is mapped to the rate of change of thermal convection contribution; and the combined effect of the changes in porosity and permeability is mapped to the rate of change of volumetric heat capacity. These rates of change—equivalent thermal conductivity, thermal convection contribution, and volumetric heat capacity—are expressed as dimensionless product factors, and these three factors together constitute the dynamic reduction factor output. In this dynamic reduction factor, the rate of change factor of equivalent thermal conductivity is used to correct the heat conduction term in the temperature field solution process; the rate of change factor of thermal convection contribution is used to correct the convection term of heat carried by fluid flow; and the rate of change factor of volumetric heat capacity is used to correct the thermal energy storage term in the energy conservation equation. The three factors are updated synchronously within the same time step and act synchronously on the corresponding calculation modules of the solver.

[0079] The back-substitution correction unit and threshold determination unit in the evaluation module collaboratively execute multi-order closed-loop iterative operations according to the iteration round sequence. During each iteration, the back-substitution correction unit uses the temperature decay curve output in the current iteration round as the sole back-substitution input, extracts the temperature decay rate sequence and temperature value sequence carried in the curve, and sends these sequences to the seawater boundary module and isotope correction module respectively. After updating and correcting the boundary conditions and physical property parameters for the current round, it triggers the coupled solver to execute the next multi-field coupled numerical operation to obtain the temperature decay curve for the next iteration round. After the temperature decay curve of each new round is output, the threshold determination unit immediately calculates the comprehensive deviation rate between the new round curve and the previous round curve and compares this comprehensive deviation rate with a preset threshold. When the overall deviation rate reaches a state less than the preset threshold, the threshold determination unit stops outputting the iteration command, and the system records the temperature decay curve of the current round as the final prediction result; when the overall deviation rate does not reach a state less than the preset threshold, the threshold determination unit sends a continue iteration command to the back-substitution correction unit, and the back-substitution correction unit then uses the temperature decay curve of the current round as the only back-substitution input for the next round of iteration, and repeats the above back-substitution correction, multi-field coupling solution and threshold comparison process.

[0080] In the aforementioned multi-stage closed-loop iteration process, each iteration replaces the back-substitution input data used in the previous round with the latest output temperature decay curve. The temperature decay curves and back-substitution parameters corresponding to earlier rounds are no longer involved in the correction calculations of subsequent rounds. Specifically, when the iteration reaches the Nth round, the back-substitution correction unit only reads the temperature decay curve output in the Nth round as the basis for the back-substitution correction in the N+1th round. The temperature decay curves output in the N-1th round and earlier rounds are marked as historical data within the back-substitution correction unit and are no longer used. The temperature decay rate sequence and temperature value sequence they carry no longer have any impact on the identification of instantaneous pressure conduction coefficient and the calculation of dynamic reduction coefficient in subsequent rounds.

[0081] At the start of each iteration, the coupled solver uses the basic parameters stored in the initial geological model as the starting point for calculation, and adds the boundary condition corrections and physical property parameter corrections that have been updated up to the current iteration, rather than directly adding increments based on the calculation results of the previous iteration. This ensures that the solution process in each iteration is based on the complete parameter state rather than relative changes.

[0082] The above iterative process terminates after the threshold determination unit determines that the comprehensive deviation rate is less than the preset threshold. After that, the system will no longer perform new multi-field coupling operations, and the temperature decay curve corresponding to the final convergence round is output as the final prediction result of the entire closed-loop iterative process.

[0083] After determining that the model has converged, the evaluation module outputs a set of evaluation results to external devices through its internal result output interface. This set of evaluation results includes the temperature decay prediction curve and the corresponding maximum allowable mining volume scheme.

[0084] The temperature decay prediction curve is recorded in the form of time series data, including the average reservoir temperature value corresponding to each output node on the prediction time axis, calculated from the current calculation start time. The curve is provided simultaneously in two forms: a data table and a graphical curve. The data table lists the prediction time and the corresponding average reservoir temperature value row by row. The graphical curve plots the continuously changing temperature decay trajectory with the prediction time on the horizontal axis and the average reservoir temperature on the vertical axis. The temperature decay prediction curve includes a horizontal reference line corresponding to the minimum economic temperature for geothermal utilization. The intersection of this reference line and the temperature decay trajectory indicates the critical time point when the reservoir temperature drops to the point where further exploitation is no longer possible under the predicted conditions.

[0085] The maximum allowable extraction volume scheme is output in conjunction with the aforementioned temperature decay prediction curve. This scheme is recorded as time-series data, containing the recommended upper limit of extraction flow rate for each output node on the same time axis as the temperature decay prediction curve. The generation rule for the maximum allowable extraction volume scheme is as follows: In the early stages of the temperature decay prediction curve, the recommended upper limit of extraction flow rate is set to the flow rate value of the initial extraction scheme; as the prediction time progresses, when the slope of the temperature decay prediction curve changes, the recommended upper limit of extraction flow rate is adjusted downwards in a stepwise manner in the opposite direction to the change in slope. That is, when the temperature decay rate accelerates, the upper limit of extraction flow rate is lowered; when the temperature decay rate slows down, the upper limit of extraction flow rate is maintained or moderately adjusted back. This stepwise adjustment is based on multiple consecutive time intervals, with each time interval corresponding to a constant recommended upper limit of extraction flow rate. The recommended upper limit of extraction flow rate between adjacent time intervals is adjusted stepwise according to the above rule, forming a dynamically changing stepwise extraction volume adjustment sequence that progresses with the prediction time.

[0086] The temperature decay prediction curve and the maximum allowable extraction volume scheme are used together to formulate the extraction-replenishment balance control strategy for seawater-recharged geothermal resources. The strategy is formulated as follows: the predicted temperature value and its trend direction at the current moment are extracted from the temperature decay prediction curve; the recommended upper limit of extraction flow rate at the current moment is extracted from the maximum allowable extraction volume scheme; the extracted predicted temperature value and its trend direction are used as the basis for determining the reservoir status; the extracted upper limit of extraction flow rate is used as the target constraint value for the operating parameters of the extraction equipment; the two are combined to form a control command set, which is then sent to the geothermal extraction central control platform.

[0087] The geothermal extraction control platform adjusts the operating frequency of the extraction pumps based on the recommended upper limit of the extraction flow rate in the control command set, ensuring that the actual extraction flow rate does not exceed this upper limit. Simultaneously, it determines the time interval for the next adjustment of the extraction flow rate based on the trend of the predicted temperature value: the time interval is shortened when the predicted temperature value shows a downward trend and extended when the predicted temperature value shows a stable trend. Thus, the temperature decay prediction curve provides information on the future evolution of the geothermal reservoir state, and the maximum allowable extraction volume scheme provides the control parameters for the corresponding state. Together, they achieve dynamic control of the balance between the extraction intensity and reservoir replenishment of seawater-recharged geothermal resources.

[0088] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs, characterized in that, Includes a seawater boundary module, an isotope correction module, and an evaluation module; The seawater boundary module collects real-time monitoring data of the target thermal reservoir and its nearshore waters, identifies the fluctuation range of the instantaneous pressure transmission coefficient and mixing ratio of the seawater replenishment boundary under the driving force of tidal dynamic load, and outputs a set of dynamic boundary conditions modulated by the tidal cycle. The isotope correction module receives isotope detection data from water samples from geothermal wells. Based on the isotope detection data and a pre-set isotope fractionation response library, it inverts the rate of change in the mixing ratio of freshwater and seawater within the geothermal reservoir and the resulting mineral dissolution-precipitation equilibrium shift. Then, based on the rate of change in the mixing ratio and the mineral dissolution-precipitation equilibrium shift, it calculates the dynamic reduction coefficient of the water-rock exchange process on the equivalent thermal conductivity and heat capacity of the geothermal reservoir. The dynamic reduction coefficient is then output as an internal parameter correction term of the temperature field prediction model. The evaluation module is connected to the seawater boundary module and the isotope correction module respectively, receives the dynamic boundary condition set and the internal parameter correction term, performs multi-field coupled temperature decay prediction and outputs the evaluation result.

2. The seawater-recharged geothermal reservoir temperature decay prediction and dynamic evaluation system according to claim 1, characterized in that, The seawater boundary module includes a monitor and a spectrum analysis unit; The monitors are deployed in the key structural zone of the target geothermal reservoir and on the nearshore seabed to continuously collect geothermal water level, fluid temperature, fluid pressure and real-time nearshore tide data. The spectrum analysis unit extracts the dominant frequency and performs phase difference analysis on the long-period tidal effect response sequence and short-period tidal effect response sequence of the geothermal water level, and identifies the fluctuation range of the instantaneous pressure transmission coefficient and the mixing ratio based on the analysis results.

3. The system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs according to claim 2, characterized in that, The spectrum analysis unit performs a discrete Fourier transform on the real-time tidal time series near the shore, extracts the dominant frequency components of tidal level fluctuations, and records them as the tidal fundamental frequency.

4. The system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs according to claim 3, characterized in that, The spectrum analysis unit performs a discrete Fourier transform on the geothermal water level time series, extracts the frequency components corresponding to the energy peaks in the frequency band near the tidal fundamental frequency, and records them as the main frequency of the geothermal water level tidal response. It also extracts the frequency components corresponding to the energy peaks in the integer multiples of the tidal fundamental frequency and records them as the multiples of the geothermal water level tidal response.

5. The seawater-recharged geothermal reservoir temperature decay prediction and dynamic evaluation system according to claim 4, characterized in that, The spectrum analysis unit subtracts the phase at the main frequency of the geothermal water level tidal response from the phase at the tidal fundamental frequency to obtain the fundamental frequency phase lag, and subtracts the phase at the harmonic frequency of the geothermal water level tidal response from the phase at an integer multiple of the tidal fundamental frequency to obtain the harmonic frequency phase lag. Divide the amplitude at the dominant frequency of the geothermal water level tidal response by the amplitude at the fundamental tidal frequency to obtain the fundamental frequency amplitude attenuation ratio. Divide the amplitude at the harmonic frequency of the geothermal water level tidal response by the amplitude at an integer multiple of the fundamental tidal frequency to obtain the harmonic frequency amplitude attenuation ratio.

6. The system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs according to claim 5, characterized in that, The seawater boundary module constructs a phase lag set from the fundamental frequency phase lag and the harmonic frequency phase lag, and constructs an amplitude attenuation set from the fundamental frequency amplitude attenuation ratio and the harmonic frequency amplitude attenuation ratio. The phase lag set and the amplitude attenuation set are then jointly solved with pre-stored formation elastic water storage coefficient, pressure conductivity coefficient, and fracture network spatial distribution data to solve for the attenuation rate constant of tidal pressure in the fracture network. This constant is output as the instantaneous pressure transmission coefficient. Furthermore, the module inversely calculates the time variation range of the ratio of the increase in seawater volume entering the thermal reservoir to the original freshwater volume within the thermal reservoir per unit tidal cycle, which is output as the mixing ratio fluctuation range.

7. The system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs according to claim 1, characterized in that, The isotope correction module has a pre-set isotope fractionation response library, which is established based on rock and mineral separation experiments in the target area and water-rock exchange kinetic experiments under different temperature and pressure conditions. It records the strontium isotope fractionation coefficient, boron isotope fractionation coefficient and corresponding mineral dissolution-precipitation equilibrium constant under different temperature, salinity and water-rock ratio conditions.

8. The seawater-recharged geothermal reservoir temperature decay prediction and dynamic evaluation system according to claim 7, characterized in that, The isotope correction module uses the strontium isotope ratio and boron isotope ratio measured in the geothermal well water sample as the search input. It matches the corresponding seawater mixing ratio identifier and mineral dissolution-precipitation equilibrium constant identifier in the isotope fractionation response library to the nearest neighbor. It then inverts the actual mixing ratio of freshwater and seawater in the current geothermal reservoir and calculates the rate of change of the mixing ratio over time in combination with the sampling time interval. Then, based on the rate of change, the driving direction and reaction intensity of the water-rock exchange reaction are determined. Based on the mineral dissolution-precipitation equilibrium constant value, the difference between the amount of dissolution and precipitation of the main mineral phases is calculated to obtain the mineral dissolution-precipitation equilibrium offset.

9. The system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs according to claim 8, characterized in that, The isotope correction module maps the mineral dissolution-precipitation equilibrium offset to an equivalent thermal conductivity reduction factor and a heat capacity reduction factor, which together constitute the dynamic reduction factor. The equivalent thermal conductivity reduction factor and the heat capacity reduction factor exist in the form of dimensionless product factors and are synchronously output to the evaluation module within the same time step.

10. The system for predicting and dynamically evaluating the temperature decay of seawater-recharged geothermal reservoirs according to claim 1, characterized in that, The evaluation module includes a coupled solver, a back-substitution correction unit, and a threshold determination unit; The coupled solver performs the first multi-field coupling operation based on the initial geological model and the current mining scheme, and outputs the first predicted curve of temperature field decay over time. The back-substitution correction unit substitutes the first prediction curve back to the seawater boundary module to update the tidal boundary pressure transmission coefficient for the next iteration. Simultaneously, the temperature decay trend of the first predicted curve is fed back to the isotope correction module to update the progress of the water-rock exchange reaction. After being updated by the back-substitution correction unit, the threshold determination unit triggers the coupled solver to perform a second multi-field coupling operation to obtain a second prediction curve, and compares the comprehensive deviation rate between the second prediction curve and the first prediction curve with a preset threshold. If the overall deviation rate is higher than the preset threshold, multi-stage closed-loop iteration is triggered until the overall deviation rate of the temperature decay curves output by two adjacent iterations is lower than the preset threshold. When the model is deemed to have converged, the final prediction result is output.

Images