Numerical simulation method for reactive solute transport in porous media

Abstract

The application discloses a numerical simulation calculation method for a reactive solute migration process in a porous medium and relates to the technical field of simulation calculation, and comprises the following steps: in a strong concentration front rapid advancing stage, collecting solute concentration gradient distribution data and flow velocity variation trajectory data, marking a concentration gradient sudden increase section according to a time advancing sequence, extracting a spatial section corresponding to an abrupt change range, and forming a front high sensitive area positioning result; based on the front high sensitive area positioning result, carrying out local grid refinement division on the corresponding spatial section, and simultaneously compressing a time advancing interval, and carrying out concentration evolution calculation in the refined spatial section. Through positioning the front high sensitive area and implementing local grid refinement and time compression, in combination with monotone preservation and directional consistency constraint, the application inhibits negative values and oscillation diffusion; further through dynamic boundary backtracking and neighborhood weighted replacement, the application restores the pollution front position, and improves the prediction stability and risk assessment reliability.

Description

Technical Field

[0001] This invention relates to the field of simulation calculation technology, specifically to a numerical simulation calculation method for the migration process of reactive solutes in porous media. Background Technology

[0002] Numerical simulation of reactive solute migration processes in porous media refers to constructing a mass conservation control equation that includes convection, dispersion, and reaction terms, given the porous media's porosity, permeability coefficient, velocity distribution, dispersion coefficient, and reaction kinetic parameters such as adsorption, desorption, precipitation, complexation, and biodegradation. This involves discretizing the continuous partial differential equations into a computable set of algebraic equations and then progressively solving for the evolution of solute concentration at each node over time and space within a spatial grid and time step framework. This allows for the reproduction of the entire migration and transformation process of solutes in groundwater aquifers, soil pores, or landfill linings. Its function is to quantitatively reveal the extent of pollution plume expansion, peak concentration arrival time, reaction reduction rate, and residual risk level. This provides predictable spatiotemporal distribution data for selecting remediation schemes for contaminated sites, determining the dosage of injected reagents, optimizing the location of isolation barriers, and conducting long-term environmental risk assessments. It also overcomes the limitations of field experiments, such as long cycles, limited scale, high costs, and difficulty in reproducibility, enabling controllable extrapolation and forward-looking decision support for complex coupled processes.

[0003] The existing technology has the following shortcomings: During the rapid advance of a high-concentration front, the solute concentration exhibits a steep gradient distribution along the migration direction. When the spatial grid is coarse or the node spacing is not set reasonably, the concentration changes in local areas cannot be fully analyzed, resulting in non-physical negative values ​​in the discrete calculation results at individual nodes. These negative values ​​participate in the numerical update calculations of adjacent nodes during subsequent time step iterations, gradually spreading to the surrounding areas and being amplified, forming an apparent abnormally low concentration area. This creates a false reduction phenomenon in the artificially constructed concentration field, thereby masking the actual location of the pollution front leading edge, reducing the accuracy of judging the expansion range and arrival time of the pollution plume, and affecting subsequent risk assessment and remediation decisions.

[0004] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0005] The purpose of this invention is to provide a numerical simulation method for the migration process of reactive solutes in porous media, so as to solve the problems in the background art mentioned above.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a numerical simulation calculation method for the migration process of reactive solutes in porous media, comprising the following steps: Step 1: During the rapid advance of the high-concentration front, collect solute concentration gradient distribution data and flow velocity change trajectory data. Mark the sections where the concentration gradient suddenly increases according to the time progression, extract the spatial segments corresponding to the steep change range, and form the location results of the high-sensitivity area of ​​the front. Step 2: Based on the location results of the frontal high-sensitivity area, the corresponding spatial segment is locally refined into a grid, and the time advance interval is compressed simultaneously. Concentration evolution calculations are carried out in the refined spatial segment to reduce the probability of non-physical negative values. Step 3: Around the refined spatial segment, a monotonicity maintenance control rule is introduced during the concentration update process, and the transmission range between adjacent nodes is limited according to the concentration evolution results within the refined spatial segment, so as to suppress the diffusion of numerical oscillations along the migration direction. Step 4: Based on the monotonicity maintenance control rule, apply directional consistency constraints to the flux calculation results of the migration direction, and correct the transfer direction by combining the concentration evolution results in the refined spatial segment to avoid local reverse amplification and the formation of abnormally low value segments. Step 5: Based on the directional consistency constraint results, the overall concentration distribution is dynamically reconstructed by boundary backtracking, and neighborhood weighted replacement processing is implemented for abnormally low value segments to restore the position of the pollution front and improve the overall prediction stability.

[0007] Preferably, the extraction of spatial segments corresponding to steep changes to form the location results of highly sensitive frontal zones includes the following steps: During the rapid advance phase of the high-concentration front, the solute concentration distribution of each time step in the computational domain is recorded according to the time progression sequence. The concentration difference between adjacent nodes is calculated sequentially according to the migration direction, and the directional attribute is recorded. At the same time, the velocity change trajectory data and the concentration difference data are correlated to form a time series set of concentration gradient distribution data. Based on the time series set of concentration gradient distribution data, the concentration difference of each time step is sorted according to the time progression order and compared with the gradient change threshold. Candidate nodes of gradient surge are marked, and consistent segments within a continuous time range are selected and determined as concentration gradient surge segments in combination with the direction of flow velocity change trajectory. The concentration gradient abruptly increases and extends node by node along the migration direction. Based on the concentration difference range, continuous spatial nodes are selected to form a steep change range, and the spatial node range is further defined by the direction of the velocity change trajectory. A unified numbering system was used to identify steep changes in the area, and a comprehensive data table was constructed by recording the time progression and spatial coordinate information. The spatial segment information was continuously updated during the time progression to form the location results of the high-sensitivity area of ​​the front.

[0008] Preferably, based on the location results of the high-sensitivity frontal zone, the corresponding spatial segment is locally refined into a finer mesh and the time advance interval is simultaneously compressed, including the following steps: Based on the spatial starting point and spatial ending point coordinates recorded in the location results of the highly sensitive area of ​​the front, spatial segments are located in the spatial node arrangement sequence of the computational domain. The original spatial units are uniformly divided along the migration direction, and the continuous spatial units at the leading edge of the pollution front are further divided to form a continuous high-density spatial node sequence. Based on the start time step number and end time step number recorded in the location results of the high-sensitivity frontal zone, the time advance interval within the time range is divided into continuous time advance sub-intervals, and a correspondence is established between the time advance sub-intervals and the refined spatial node spacing. Around the refined spatial segments and time-progression sub-intervals, the concentration transfer amount is updated node by node according to the migration direction within each time-progression sub-interval, limiting the concentration change to cover the range of the refined sub-units, and continuously saving the node concentration values ​​to advance the concentration evolution calculation.

[0009] Preferably, the original spatial unit is uniformly divided into sub-units along the migration direction, and the continuous spatial unit at the leading edge of the pollution front is further divided to form a continuous high-density spatial node sequence. The time-progression sub-interval is divided according to the start time step number and the end time step number, and is set in accordance with the refined spatial node spacing. The concentration transfer is limited to the coverage area of ​​the refined sub-unit.

[0010] Preferably, introducing a monotonic hold-up control rule during the concentration update process includes the following steps: At the start of the time-progression sub-interval, the concentration values ​​of spatial nodes in the refined spatial segment are read in the order of migration direction to construct a node concentration sequence. Non-monotonic change nodes are identified by comparing three consecutive nodes to form a list of monotonically maintained control nodes and record the corresponding concentration range. Based on the monotonicity-preserving control node list, the corresponding spatial nodes are subjected to the transfer range limitation processing between adjacent nodes within the time-progressing sub-interval. The concentration change range is limited to the concentration range between the preceding and following nodes. Spatial nodes not included in the monotonicity-preserving control node list are updated according to the predetermined concentration evolution order. Based on the concentration evolution results at the end of the time-progression sub-interval, the list of monotonic maintenance control nodes is adjusted to include non-monotonic change nodes that recur within the continuous time range and adjacent spatial nodes within the range of the transfer amount interval, forming a continuous control segment along the migration direction.

[0011] Preferably, the transfer range limitation processing uses the concentration range of the preceding and following nodes recorded in the monotonically maintained control node list as a constraint condition, and dynamically updates the monotonically maintained control node list in each time advancement sub-interval by combining the concentration evolution results in the refined spatial segment, thereby forming a continuous control segment in the migration direction.

[0012] Preferably, applying directional consistency constraints to the flux calculation results in the migration direction includes the following steps: When the concentration is updated and the monotonicity maintenance control rule is executed in the time-progression sub-interval, the flux calculation results between adjacent nodes in the refined spatial segment are extracted, a flux record table is established, and a candidate list of directional consistency constraints is formed by identifying fluxes with inconsistent directions in combination with the concentration sequence. Based on the candidate list of directional consistency constraints, the directional inconsistency flux is adjusted, the reverse flux value is corrected to zero and redistributed to adjacent nodes according to the migration direction, while maintaining the transmission range between adjacent nodes as defined by the monotonicity maintenance control rule. Based on the concentration evolution results within the refined spatial segments, the transmission direction of spatial segments with reverse flux records within a continuous time range is uniformly adjusted, and the updated flux records are used as the basis for concentration updates in the next time-progressing sub-interval.

[0013] Preferably, during the flux adjustment process for inconsistent directions, continuous spatial segments with reverse flux records are uniformly identified, and flux redistribution is completed within the range of transfer between adjacent nodes defined by the monotonicity maintenance control rule, so as to ensure that the flux direction in the refined spatial segment is consistent with the concentration evolution result.

[0014] Preferably, the dynamic boundary backtracking reconstruction of the global concentration distribution based on the directional consistency constraint results includes the following steps: When the concentration update is completed in the time-progression sub-interval, the global concentration distribution sequence is read and the flux record in the directional consistency constraint result is combined to identify abnormal low value segments, and the spatial start node number and spatial end node number are recorded to form the positioning information. Based on the location information of abnormal low value segments, the starting point and ending point of the backtracking are determined in the migration direction. A concentration transition interval is constructed between the starting point and the ending point, and the concentration values ​​of intermediate nodes are reconstructed according to the relative position ratio. The reconstruction results are written into the global concentration distribution sequence. Neighborhood weighted replacement processing is implemented around the abnormally low value segment in the migration direction. The node concentration distribution is adjusted according to the concentration values ​​of adjacent nodes and the proportion of their spatial locations to form a continuous concentration sequence for subsequent time progression.

[0015] Preferably, the location of abnormally low value segments is determined by the flux record in the directional consistency constraint results and the concentration change direction within the continuous time-progressing sub-interval. Furthermore, the concentration transition interval reconstruction and neighborhood weighted replacement processing are both continuously implemented along the migration direction to maintain the unidirectional progressive relationship of the global concentration distribution sequence.

[0016] The technical effects and advantages provided by the present invention in the above technical solution are as follows: This invention precisely locates the highly sensitive region of a strong concentration front during its rapid advance and implements local mesh refinement and synchronous compression of the time progression interval within this region. This allows concentration evolution to unfold on a finer spatiotemporal scale, reducing the probability of non-physically negative values ​​from the outset. Simultaneously, by introducing monotonicity-maintaining control rules and directional consistency constraints during concentration updates, the invention ensures a continuous progression of concentration along the migration direction, preventing numerical oscillations and thus improving the stability and physical consistency of the concentration field evolution process.

[0017] After constraining directional consistency, this invention further reconstructs the global concentration distribution through dynamic boundary backtracking and applies neighborhood-weighted replacement processing to abnormally low-value segments, accurately restoring the position of the pollution front's leading edge and avoiding interference from spurious attenuation in the prediction results. Through these processes, the determination of the pollution plume's expansion range and arrival time more closely approximates the actual migration state, thereby improving the overall prediction stability and providing a more reliable spatiotemporal distribution basis for subsequent risk assessment and remediation decisions. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0019] Figure 1 This is a flowchart of the numerical simulation calculation method for the migration process of reactive solutes in porous media according to the present invention. Detailed Implementation

[0020] Exemplary embodiments will now be described more fully with reference to the accompanying drawings. However, these exemplary embodiments can be implemented in many forms and should not be construed as limited to the examples set forth herein; rather, they are provided so that the description of this disclosure will be more complete and fully convey the concept of the exemplary embodiments to those skilled in the art.

[0021] This invention provides, for example Figure 1 The numerical simulation method for the migration process of reactive solutes in porous media, as shown, includes the following steps: Step 1: During the rapid advance of the high-concentration front, collect solute concentration gradient distribution data and flow velocity change trajectory data. Mark the sections where the concentration gradient suddenly increases according to the time progression, extract the spatial segments corresponding to the steep change range, and form the location results of the high-sensitivity area of ​​the front. The spatial segments corresponding to the steep changes in elevation are extracted to form the location results of the highly sensitive area of ​​the front. The specific steps are as follows: During the rapid advance phase of the high-concentration front, the solute concentration distribution at each time step in the computational domain is recorded layer by layer according to a predetermined time progression sequence. Within each time step, all spatial nodes are arranged sequentially according to the migration direction, and the concentration values ​​between adjacent nodes are read in turn. The concentration difference between adjacent nodes is calculated and its directional attribute is recorded. After recording the concentration difference values ​​of all spatial nodes within a single time step, the velocity change trajectory data and concentration difference data corresponding to that time step are matched one-to-one. Specifically, for each spatial node, the velocity direction, velocity magnitude, and its component in the migration direction at the current time step are recorded, and the velocity component and the concentration difference between the two nodes before and after that node are stored together in the same time layer data sequence. During the time progression, the magnitude of the concentration difference change between the previous time step and the current time step is compared node by node to form a time series set of concentration gradient distribution data covering the entire computational domain. At the same time, the velocity change trajectory data within the corresponding time range are arranged in chronological order to form a complete synchronous recording sequence of concentration gradient distribution data and velocity change trajectory data.

[0022] After obtaining the time series set of concentration gradient distribution data under continuous time progression, the concentration difference data within each time step is scanned layer by layer according to the time progression order. The concentration difference of each spatial node is sorted by magnitude and compared with a pre-set gradient change threshold. When the concentration difference of a spatial node in the current time step exceeds the gradient change threshold, the node is marked as a candidate node for gradient surge, and its spatial coordinates and corresponding time step number are recorded. After marking all candidate nodes within a single time step, the candidate nodes of that time step are compared with those of the adjacent time steps before and after. Inter-overlap comparison: If candidate nodes for a sudden increase in concentration gradient appear in the same spatial segment in multiple consecutive time steps, and the direction of the velocity change trajectory within the segment remains consistent in consecutive time steps, then the spatial segment is identified as a concentration gradient surge segment. During the labeling process, the direction of the concentration difference value and the direction of the velocity change trajectory of each candidate node are checked one by one. Only when the direction of the concentration difference value and the direction of the velocity component are the same and there is no reversal of direction in consecutive time steps, is the corresponding segment confirmed as a concentration gradient surge segment, thereby ensuring the continuity and directional consistency of the labeling results in the time progression sequence.

[0023] After completing the temporal continuous marking of concentration gradient surge segments, each confirmed concentration gradient surge segment is spatially extended and extracted. Specifically, starting from the first and last spatial nodes of the marked concentration gradient surge segment, the concentration difference change is checked node by node along the migration direction. When the concentration difference of an adjacent node is still within the top 20% of all concentration differences in the current time step, that node is included in the same spatial segment. Extension stops when the concentration difference falls below the aforementioned range. Then, the same checking process is performed node by node in reverse along the migration direction. Nodes meeting the criteria are incorporated into the same spatial segment. After completing the spatial extension within a single time step, the extension results of the corresponding spatial segments in consecutive time steps are superimposed. Spatial nodes that appear repeatedly in multiple time steps are retained, while isolated nodes that appear only in a single time step are deleted, thus forming spatial segments corresponding to stable steep change ranges. During the spatial extension process, the velocity change trajectory data are checked simultaneously to ensure that the velocity component direction of each node included in the spatial segment is consistent with the overall migration direction, thereby ensuring that the extracted spatial segment completely covers the transition zone of the pollution front.

[0024] After establishing spatial segments corresponding to continuous and stable steep changes, these segments are uniformly numbered across the entire computational domain. Their corresponding start and end time step numbers, spatial origin coordinates, and spatial end coordinates are centrally recorded, constructing a comprehensive data table containing information on the time progression sequence, concentration gradient distribution data characteristics, and velocity change trajectory direction. During subsequent time progression, each time a new time step is calculated, the same steps are followed to re-identify segments with abrupt increases in concentration gradient. The newly identified spatial segments are compared with existing numbered segments. If spatial overlap exists, they are merged into the same numbered segment, and their time range records are updated. If no overlap exists, a new numbered segment is generated and added to the comprehensive data table. By continuously updating the time and spatial range information, a location result covering the entire rapid advance phase of the high-concentration front is ultimately formed. This location result maintains continuous recording in both time and space dimensions and fully reflects the spatial segment distribution corresponding to the steep changes, providing clear spatial boundaries and temporal order for subsequent local mesh refinement and time-progression interval compression based on this high-sensitivity front location result.

[0025] Step 2: Based on the location results of the frontal high-sensitivity area, the corresponding spatial segment is locally refined into a grid, and the time advance interval is compressed simultaneously. Concentration evolution calculations are carried out in the refined spatial segment to reduce the probability of non-physical negative values. Based on the location results of the high-sensitivity frontal zone, the corresponding spatial segment is locally refined into a mesh. The specific steps are as follows: Based on the spatial start and end coordinates recorded in the location results of the frontal high-sensitivity area, the starting and ending nodes of the spatial segment are located one by one in the original spatial node arrangement sequence of the computational domain, and each original spatial unit within the segment is re-divided. Specifically, if the node spacing of the original spatial unit in the migration direction is a fixed length, each original spatial unit is evenly divided into four sub-units, with each sub-unit occupying one-quarter of the original unit length along the migration direction, while maintaining the same horizontal node arrangement order as the original order. After completing the first round of even division, a second division is performed on the three consecutive original spatial units located at the center of the frontal high-sensitivity area location results, and each sub-unit in these three segments is then evenly divided again. The system is divided into two smaller units, resulting in a local spatial node distribution at the leading edge of the pollution front with a density eight times that of the original nodes. During the splitting process, each newly added spatial node is assigned a unique number, and the refined spatial node list is reorganized according to the migration direction, so that the refined spatial segment forms a continuous and uninterrupted high-density node sequence. At the front and rear boundaries of the refined segment, the outermost sub-unit of the refined system is docked with the original spatial unit of the original unrefined area. The boundary position is split only once proportionally to avoid abrupt changes in node spacing. Through the above specific division method, the corresponding spatial segment forms a spatial node distribution structure that gradually transitions from the original spacing to one-eighth of the original spacing in the migration direction, thereby completing the local mesh refinement.

[0026] After refining the local mesh for the corresponding spatial segment, the time advance interval within that time range is synchronously compressed based on the start and end time step numbers recorded in the frontal high-sensitivity area location results. Specifically, if the original time advance interval is a fixed length, then within the time range defined by the frontal high-sensitivity area location results, each original time step is divided into three consecutive time advance sub-intervals, each sub-interval being one-third the length of the original time advance interval. During the time advance process, when the time advance enters the start time step number recorded in the frontal high-sensitivity area location results, it switches to the compressed time advance interval. The time-advancing interval sequence is compressed, and after the end time step numbering, it is restored to the original time-advancing interval. While compressing the time-advancing interval, each compressed time-advancing sub-interval is arranged in correspondence with the refined spatial node spacing to ensure that within a single time-advancing sub-interval, the concentration change in the migration direction only covers the length of one or two refined minimum sub-units. During the time-advancing process, the compressed time-advancing sub-intervals are numbered one by one, and corresponding records are established with the refined spatial node numbers in sequence to keep the time-advancing order consistent with the spatial node order, thereby forming a advancement structure in which time and space shrink synchronously within the corresponding spatial segment.

[0027] After completing the local grid refinement and time-progression interval compression, concentration evolution calculations are performed within the refined spatial segments. Specifically, at the beginning of each compressed time-progression sub-interval, the concentration values ​​of all spatial nodes within the refined spatial segment at the end of the previous time-progression sub-interval are read. The concentration transfer is updated node by node from upstream to downstream according to the migration direction, limiting the concentration change between two adjacent refined sub-units to the propagation distance corresponding to the current time-progression sub-interval. During the update process, for high-density node sequences located at the leading edge of the pollution front, the concentration difference is calculated node by node, and the transfer within a single time-progression sub-interval is distributed to multiple refined sub-units, allowing the concentration gradient to unfold progressively across multiple continuous small-scale spatial units. At the end of a time-progression sub-interval, the concentration of all nodes within the refined spatial segment is calculated. The numerical values ​​are saved as a new time-layer state, and the process continues from this state at the start of the next time-advance sub-interval. Within the time limit defined by the location results of the highly sensitive frontal region, the concentration evolution is continuously advanced according to the compressed time-advance sub-intervals. When the time advance exceeds the end time step number, the refined spatial segment structure remains unchanged, but the time-advance interval is restored to its original length. By refining the local mesh in the corresponding spatial segment, synchronously compressing the time-advance interval, and performing concentration evolution calculations sub-interval by sub-interval within the refined spatial segment, the concentration changes during the rapid advance of the strong concentration front are progressively transmitted in smaller-scale spatial nodes and shorter time-advance sub-intervals. This reduces the probability of non-physical negative values ​​and provides continuous and smooth concentration evolution conditions for the subsequent introduction of monotonicity-maintaining control rules and directional consistency constraints within the refined spatial segment.

[0028] Step 3: Around the refined spatial segment, a monotonicity maintenance control rule is introduced during the concentration update process, and the transmission range between adjacent nodes is limited according to the concentration evolution results within the refined spatial segment, so as to suppress the diffusion of numerical oscillations along the migration direction. A monotonic hold-through control rule is introduced during the concentration update process. The specific steps are as follows: At the beginning of each time-advancement sub-interval, the concentration values ​​of all spatial nodes within the refined spatial segment at the end of the previous time-advancement are read sequentially according to the migration direction. These concentration values ​​are then arranged from upstream to downstream according to the node number, forming a complete node concentration sequence. Within this node concentration sequence, three consecutive nodes are grouped together, and the concentration values ​​of the preceding node, the intermediate node, and the following node are extracted sequentially. The magnitude relationship between the concentration values ​​of the intermediate node and the preceding and following nodes is compared group by group. When the concentration value of the intermediate node is higher than the larger value between the preceding and following nodes, or lower than the smaller value between the preceding and following nodes, the intermediate node is marked as a non-monotonic node, and its specific position number in the refined spatial segment is recorded. After completing the scanning of all three-node combinations, all marked non-monotonic nodes are summarized to form a list of monotonicity-maintaining control nodes. At the same time, the concentration value range between the preceding and following nodes corresponding to each non-monotonic node is recorded as reference data for subsequently limiting the transmission range between adjacent nodes, thereby clarifying the specific target of the monotonicity-maintaining control rule within the refined spatial segment.

[0029] After obtaining the list of monotonic hold control nodes, when updating the concentration node by node within the current time advance sub-interval, a range limitation process for the transfer amount between adjacent nodes is implemented for each node in the list of monotonic hold control nodes. Specifically, when updating the concentration transfer amount between a monotonic hold control node and its predecessor node, the concentration values ​​of the monotonic hold control node and its predecessor node at the start of the current time advance sub-interval are first read. The direction of the concentration difference between the two is calculated, and the maximum allowable transfer amount is limited to not exceeding the full amplitude of the concentration difference, while the minimum allowable transfer amount is limited to zero. When the transfer amount in the original update result exceeds the above allowable range, the transfer amount is corrected to the range. The boundary values ​​are calculated, and the corrected transfer amount is rewritten into the concentration update record between the node and the previous node. Subsequently, the same processing procedure is performed on the concentration transfer amount between the monotonic control node and its next node to ensure that the concentration change amplitude of the monotonic control node is always within the range formed by the concentration values ​​of its preceding and following nodes within the same time advancement sub-interval. When performing the above-mentioned transfer amount range limitation processing between adjacent nodes, for other nodes not included in the monotonic control node list, they are still updated node by node according to the predetermined concentration evolution order within the refined spatial segment, without additional limitation on their transfer amount. This ensures the continuity of the overall concentration evolution while implementing targeted control at locations where numerical oscillations may occur.

[0030] After updating the concentration of all nodes within the current time-progression sub-interval, the updated concentration values ​​of all nodes within the refined spatial segment are arranged again according to the migration direction, and the concentration difference between adjacent nodes is recorded node by node. The node concentration sequence at the end of the current time-progression sub-interval is compared node by node with the node concentration sequence at the end of the previous time-progression sub-interval. If a node is found to be a non-monotonic node in two consecutive time-progression sub-intervals, at the beginning of the next time-progression sub-interval, this node and its three consecutive nodes before and after it are included in the list of monotonic hold-up control nodes, and the interval limitation processing of the transfer amount between adjacent nodes is performed synchronously within the continuous segment of these three nodes. When a non-monotonic node is found to be a non-monotonic node, the node is considered a non-monotonic node. When a changing node shows a forward movement trend in the migration direction, the two adjacent refined nodes in front of that node are included in the monotonicity-preserving control node list to form a control segment covering five consecutive nodes, thereby forming a continuous transmission range limit zone in the migration direction. By dynamically expanding or shrinking the monotonicity-preserving control node list according to the concentration evolution results in the refined spatial segment after each time advancement sub-interval, the monotonicity-preserving control rule is always adjusted around the region where numerical oscillations may occur, forming a continuous and progressive transmission range limit structure in the refined spatial segment, thereby suppressing the diffusion of numerical oscillations along the migration direction and maintaining the orderly change of concentration along the migration direction during the rapid advancement phase of the strong concentration front.

[0031] Step 4: Based on the monotonicity maintenance control rule, apply directional consistency constraints to the flux calculation results of the migration direction, and correct the transfer direction by combining the concentration evolution results in the refined spatial segment to avoid local reverse amplification and the formation of abnormally low value segments. Apply directional consistency constraints to the flux calculation results in the migration direction. The specific steps are as follows: After each time-advance sub-interval completes concentration updates and executes monotonic maintenance control rules, the flux calculation results between all adjacent nodes within the refined spatial segment are extracted one by one, and an ordered flux record table is established from upstream to downstream according to the migration direction. In this flux record table, each flux data is clearly marked with its starting node number, target node number, and transmission direction identifier. Flux from upstream nodes to downstream nodes is recorded as forward transmission, and flux from downstream nodes to upstream nodes is recorded as reverse transmission. Subsequently, the concentration of each node within the refined spatial segment at the end of the current time-advance sub-interval is read. The flux values ​​are calculated and arranged in order of migration direction to form a concentration sequence. For each flux data in the flux record table, the concentration values ​​of the corresponding two nodes are compared one by one. When the flux direction is consistent with the concentration gradient direction, that is, the concentration is transferred from the high-value node to the low-value node, the original flux result is retained. When the flux direction is opposite to the concentration gradient direction, that is, the concentration is transferred from the low-value node to the high-value node, the flux record is marked as a flux with inconsistent direction, and the specific node number and flux value are marked in the flux record table to form a candidate list of directional consistency constraints. In this way, all flux locations with directional conflicts are identified in the refined spatial segment.

[0032] After obtaining the candidate list of directional consistency constraints, each directional inconsistency flux in the list is processed one by one. Specifically, for a reverse flux record pointing from a downstream node to an upstream node, the concentration values ​​of the downstream and upstream nodes at the beginning and end of the current time-progression sub-interval are first read. The concentration change trends of the two nodes within the current time-progression sub-interval are analyzed. When the concentration of the downstream node decreases within the current time-progression sub-interval while the reverse flux shows transmission from that node to the upstream, the reverse flux value is directly adjusted to zero, and the original reverse flux value is proportionally distributed between that node and its downstream nodes. In the forward flux records between adjacent nodes, when the concentration of an upstream node increases within the current time advance sub-interval and the reverse flux shows that it is being transmitted from downstream to that node, the reverse flux value is also adjusted to zero, and the original value is redistributed between adjacent node pairs in the same direction of migration. When redistributing flux values, the transmission range between adjacent nodes defined in the monotonicity control rules is always kept within the allowable transmission range. Through the above step-by-step processing, all flux calculation results in the refined spatial segment maintain a unified direction in the migration direction, eliminating the influence of reverse flux on concentration updates.

[0033] After completing the directional consistency constraint processing and updating the flux record table, the transfer direction is corrected as a whole based on the concentration evolution results within the refined spatial segment. Specifically, after the current time-progression sub-interval ends, the concentration values ​​of all nodes within the refined spatial segment are rearranged according to the migration direction. A local segment consisting of three consecutive nodes is scanned. If the concentration value of the intermediate node is found to be lower than that of the nodes before and after it for two consecutive time-progression sub-intervals, and if the intermediate node had a reverse flux record before the directional consistency constraint processing, this three-node segment is marked as an abnormally low-value segment. For this abnormally low-value segment, all flux records within the two most recent time-progression sub-intervals are read, and all flux directions related to this segment are corrected. The flow direction is adjusted to a single direction from upstream to downstream, and the original reverse flux values ​​are redistributed to adjacent nodes according to the current concentration difference ratio. After the direction adjustment and value redistribution are completed, the updated flux records are rewritten into the flux record table, and the updated flux results are used as the basis for concentration updates in the next time interval. By applying directional consistency constraints to the flux calculation results in the migration direction under the premise of monotonicity maintenance control rules, and combining the concentration evolution results in the refined spatial segment, the flow direction is corrected segment by segment. This ensures that the concentration continues to advance along the migration direction during the rapid advance of the strong concentration front, avoiding the formation of abnormally low value segments by local reverse amplification, thereby maintaining a stable and coherent concentration evolution state in the refined spatial segment.

[0034] Step 5: Based on the directional consistency constraint results, the global concentration distribution is dynamically reconstructed by boundary backtracking, and neighborhood weighted replacement processing is implemented for abnormally low value segments to restore the position of the pollution front leading edge and improve the overall prediction stability. Based on the directional consistency constraint results, the global concentration distribution is dynamically reconstructed through boundary backtracking. The specific steps are as follows: After the current time-progression sub-interval ends, the concentration values ​​of all nodes in the entire computational domain are read in spatial node number order to form a complete global concentration distribution sequence. This global concentration distribution sequence is then compared node-by-node with the global concentration distribution sequence saved at the end of the previous time-progression sub-interval, recording the direction and magnitude of concentration change for each node within two consecutive time-progression sub-intervals. Using the migration direction as the main axis, the system scans node by node from upstream to downstream. When a node in a continuous segment shows a decrease in concentration value in two consecutive time-progression sub-intervals, and the direction of this decrease is inconsistent with the expected concentration decrease trend under the migration direction, this continuous segment is initially marked as a candidate range for an abnormally low-value segment. After marking, the flux record table in the direction consistency constraint results is retrieved, and the flux direction records of each node within the candidate range in the most recent two time-progression sub-intervals are searched one by one. When it is found that there is reverse flux from downstream to upstream participating in the update history within this segment, the segment is finally confirmed as an abnormally low-value segment, and its spatial start node number, spatial end node number, and corresponding time-progression sub-interval number are recorded, providing clear spatial and temporal location information for subsequent processing.

[0035] After identifying the abnormally low-value segment, a dynamic boundary backtracking reconstruction of the overall concentration distribution is performed. Specifically, using the spatial starting node of the abnormally low-value segment as the boundary, a backtracking search is conducted upstream along the migration direction until a node is found where the concentration change direction within the last two time intervals is consistent with the migration direction and no continuous decrease has been observed. This node is then designated as the backtracking starting point. Similarly, using the spatial ending node of the abnormally low-value segment as the boundary, a backtracking search is conducted downstream along the migration direction until a node is found where the concentration change direction is consistent with the migration direction and no continuous decrease has been observed. This node is then designated as the backtracking ending point. A new concentration distribution is then constructed between the backtracking starting point and the backtracking ending point. In the transition interval, all nodes within this interval are rearranged according to the migration direction. Using the concentration values ​​at the starting and ending points of the backtracking as the two end control values, the concentration values ​​of intermediate nodes are reconstructed by decreasing or increasing node by node along the migration direction, ensuring that the concentration change within this interval remains unidirectional. During the reconstruction process, the concentration value of each intermediate node is allocated according to its relative position to the starting and ending points of the backtracking, ensuring that the reconstructed concentration sequence is spatially continuous and without abrupt changes. After the concentration reconstruction within this interval is completed, the updated concentration values ​​are written into the global concentration distribution sequence, replacing the concentration records in the original abnormal low-value segment, thereby achieving dynamic boundary backtracking reconstruction.

[0036] After completing the dynamic boundary backtracking reconstruction, the nodes within the abnormally low value segment are further processed using a neighborhood-weighted replacement method. Specifically, within the abnormally low value segment, each node is processed sequentially from upstream to downstream. For each node, the concentration values ​​of its preceding and following nodes at the end of the current time-progression sub-interval are read, and the weighted average of the concentration values ​​of the two adjacent nodes is calculated. Nodes closer to the backtracking start point are assigned a greater weight to the concentration value of their preceding node, and nodes closer to the backtracking end point are assigned a greater weight to the concentration value of their following node. Nodes in the middle are assigned weights proportional to their distance from both ends. When the current concentration value of a node is lower than its weighted average, that node is... The node concentration values ​​are replaced with the calculated weighted average; when the current concentration value of a node is higher than its weighted average, the original value remains unchanged; after completing the neighborhood weighted replacement processing for all nodes, the concentration distribution of the entire region is arranged sequentially according to the migration direction to ensure that there is no local reverse depression pattern in the abnormally low value section; the concentration distribution of the entire region is dynamically reconstructed by backtracking based on the directional consistency constraint results, and the neighborhood weighted replacement processing is implemented node by node in the abnormally low value section, so that the position of the pollution front front is restored to a continuous advancing state in space, thereby improving the overall prediction stability and maintaining the stable evolution of the concentration field in the subsequent time advancement process.

[0037] Beneficial effect 1: This invention precisely locates the highly sensitive region of a strong concentration front during its rapid advance and implements local mesh refinement and synchronous compression of the time progression interval within this region. This allows concentration evolution to unfold on a finer spatiotemporal scale, reducing the probability of non-physically negative values ​​from the outset. Simultaneously, by introducing monotonicity-maintaining control rules and directional consistency constraints during concentration updates, the invention ensures a continuous progression of concentration along the migration direction, preventing numerical oscillations and thus improving the stability and physical consistency of the concentration field evolution process.

[0038] Benefit 2: After constraining directional consistency, this invention further reconstructs the global concentration distribution through dynamic boundary backtracking and applies neighborhood-weighted replacement processing to abnormally low-value segments, accurately restoring the position of the pollution front's leading edge and avoiding interference from spurious attenuation in the prediction results. Through these processes, the determination of the pollution plume's expansion range and arrival time more closely approximates the actual migration state, thereby improving the overall prediction stability and providing a more reliable spatiotemporal distribution basis for subsequent risk assessment and remediation decisions.

[0039] The foregoing has only described certain exemplary embodiments of the present invention by way of illustration. Undoubtedly, those skilled in the art can modify the described embodiments in various ways without departing from the spirit and scope of the present invention. Therefore, the foregoing drawings and descriptions are illustrative in nature and should not be construed as limiting the scope of protection of the claims of the present invention.

Claims

1. A numerical simulation method for the migration process of reactive solutes in porous media, characterized in that, Includes the following steps: Step 1: During the rapid advance of the high-concentration front, collect solute concentration gradient distribution data and flow velocity change trajectory data. Mark the sections where the concentration gradient suddenly increases according to the time progression, extract the spatial segments corresponding to the steep change range, and form the location results of the high-sensitivity area of ​​the front. Step 2: Based on the location results of the frontal high-sensitivity area, the corresponding spatial segment is locally refined into a grid, and the time advance interval is compressed simultaneously. Concentration evolution calculations are then carried out within the refined spatial segment. Step 3: Around the refined spatial segment, a monotonicity maintenance control rule is introduced during the concentration update process, and the transmission range between adjacent nodes is limited according to the concentration evolution results within the refined spatial segment, so as to suppress the diffusion of numerical oscillations along the migration direction. Step 4: Based on the monotonicity maintenance control rule, apply directional consistency constraints to the flux calculation results of the migration direction, and correct the transfer direction by combining the concentration evolution results in the refined spatial segment. Step 5: Based on the directional consistency constraint results, the overall concentration distribution is dynamically reconstructed by boundary backtracking, and neighborhood weighted replacement processing is implemented for abnormally low value segments to restore the position of the pollution front leading edge.

2. The numerical simulation calculation method for the migration process of reactive solutes in porous media according to claim 1, characterized in that, Extracting spatial segments corresponding to steep changes to form the location results of highly sensitive frontal zones includes the following steps: During the rapid advance phase of the high-concentration front, the solute concentration distribution of each time step in the computational domain is recorded according to the time progression sequence. The concentration difference between adjacent nodes is calculated sequentially according to the migration direction, and the directional attribute is recorded. At the same time, the velocity change trajectory data and the concentration difference data are correlated to form a time series set of concentration gradient distribution data. Based on the time series set of concentration gradient distribution data, the concentration difference of each time step is sorted according to the time progression order and compared with the gradient change threshold. Candidate nodes of gradient surge are marked, and consistent segments within a continuous time range are selected and determined as concentration gradient surge segments in combination with the direction of flow velocity change trajectory. The concentration gradient abruptly increases and extends node by node along the migration direction. Based on the concentration difference range, continuous spatial nodes are selected to form a steep change range, and the spatial node range is further defined by the direction of the velocity change trajectory. A unified numbering system was used to identify steep changes in the area, and a comprehensive data table was constructed by recording the time progression and spatial coordinate information. The spatial segment information was continuously updated during the time progression to form the location results of the high-sensitivity area of ​​the front.

3. The numerical simulation calculation method for the migration process of reactive solutes in porous media according to claim 2, characterized in that, Based on the location results of the high-sensitivity frontal zone, the corresponding spatial segment is locally meshed and the time advance interval is simultaneously compressed, including the following steps: Based on the spatial starting point and spatial ending point coordinates recorded in the location results of the highly sensitive area of ​​the front, spatial segments are located in the spatial node arrangement sequence of the computational domain. The original spatial units are uniformly divided along the migration direction, and the continuous spatial units at the leading edge of the pollution front are further divided to form a continuous high-density spatial node sequence. Based on the start time step number and end time step number recorded in the location results of the high-sensitivity frontal zone, the time advance interval within the time range is divided into continuous time advance sub-intervals, and a correspondence is established between the time advance sub-intervals and the refined spatial node spacing. Around the refined spatial segments and time-progression sub-intervals, the concentration transfer amount is updated node by node according to the migration direction within each time-progression sub-interval, limiting the concentration change to cover the range of the refined sub-units, and continuously saving the node concentration values ​​to advance the concentration evolution calculation.

4. The numerical simulation calculation method for the migration process of reactive solutes in porous media according to claim 3, characterized in that, The original spatial unit is uniformly divided into sub-units along the migration direction. The continuous spatial unit at the leading edge of the pollution front is further divided to form a continuous high-density spatial node sequence. The time-progressing sub-intervals are divided according to the starting time step number and the ending time step number, and are set to correspond to the spacing of the refined spatial nodes. The concentration transfer is limited to the coverage area of ​​the refined sub-units.

5. The numerical simulation calculation method for the migration process of reactive solutes in porous media according to claim 3, characterized in that, Introducing a monotonic hold-up control rule during concentration updates involves the following steps: At the start of the time-progression sub-interval, the concentration values ​​of spatial nodes in the refined spatial segment are read in the order of migration direction to construct a node concentration sequence. Non-monotonic change nodes are identified by comparing three consecutive nodes to form a list of monotonically maintained control nodes and record the corresponding concentration range. Based on the monotonicity-preserving control node list, the corresponding spatial nodes are subjected to the transfer range limitation processing between adjacent nodes within the time-progressing sub-interval. The concentration change range is limited to the concentration range between the preceding and following nodes. Spatial nodes not included in the monotonicity-preserving control node list are updated according to the predetermined concentration evolution order. Based on the concentration evolution results at the end of the time-progression sub-interval, the list of monotonic maintenance control nodes is adjusted to include non-monotonic change nodes that recur within the continuous time range and adjacent spatial nodes within the range of the transfer amount interval, forming a continuous control segment along the migration direction.

6. The numerical simulation calculation method for the migration process of reactive solutes in porous media according to claim 5, characterized in that, The transfer range limitation process uses the concentration range of the preceding and following nodes recorded in the monotonically maintained control node list as a constraint condition, and dynamically updates the monotonically maintained control node list in each time advancement sub-interval by combining the concentration evolution results in the refined spatial segment, thus forming a continuous control segment in the migration direction.

7. The numerical simulation calculation method for the migration process of reactive solutes in porous media according to claim 5, characterized in that, Applying directional consistency constraints to the flux calculation results in the migration direction includes the following steps: When the concentration is updated and the monotonicity maintenance control rule is executed in the time-progression sub-interval, the flux calculation results between adjacent nodes in the refined spatial segment are extracted, a flux record table is established, and a candidate list of directional consistency constraints is formed by identifying fluxes with inconsistent directions in combination with the concentration sequence. Based on the candidate list of directional consistency constraints, the directional inconsistency flux is adjusted, the reverse flux value is corrected to zero and redistributed to adjacent nodes according to the migration direction, while maintaining the transmission range between adjacent nodes as defined by the monotonicity maintenance control rule. Based on the concentration evolution results within the refined spatial segments, the transmission direction of spatial segments with reverse flux records within a continuous time range is uniformly adjusted, and the updated flux records are used as the basis for concentration updates in the next time-progressing sub-interval.

8. The numerical simulation calculation method for the migration process of reactive solutes in porous media according to claim 7, characterized in that, During the process of adjusting flux in case of direction inconsistency, continuous spatial segments with reverse flux records are uniformly identified, and flux redistribution is completed within the range of transfer between adjacent nodes defined by the monotonicity maintenance control rule.

9. The numerical simulation calculation method for the migration process of reactive solutes in porous media according to claim 7, characterized in that, The dynamic boundary backtracking reconstruction of the global concentration distribution based on the directional consistency constraint results includes the following steps: When the concentration update is completed in the time-progression sub-interval, the global concentration distribution sequence is read and the flux record in the directional consistency constraint result is combined to identify abnormal low value segments, and the spatial start node number and spatial end node number are recorded to form the positioning information. Based on the location information of abnormal low value segments, the starting point and ending point of the backtracking are determined in the migration direction. A concentration transition interval is constructed between the starting point and the ending point, and the concentration values ​​of intermediate nodes are reconstructed according to the relative position ratio. The reconstruction results are written into the global concentration distribution sequence. Neighborhood weighted replacement processing is implemented around the abnormally low value segment in the migration direction. The node concentration distribution is adjusted according to the concentration values ​​of adjacent nodes and the proportion of their spatial locations to form a continuous concentration sequence for subsequent time progression.

10. The numerical simulation calculation method for the migration process of reactive solutes in porous media according to claim 9, characterized in that, The location of abnormally low value segments is determined by the flux records in the directional consistency constraint results and the concentration change direction within the continuous time-progressing sub-interval. Furthermore, the concentration transition interval reconstruction and neighborhood weighted replacement processing are both continuously implemented along the migration direction.

Images