Sand distributed eco-hydrological two-way coupling modeling system

Abstract

The present application relates to the technical field of numerical simulation of ecological hydrological process, and discloses a sand-distributed ecological hydrological two-way coupling modeling system, which comprises initialization and discretization, ecological growth calculation and source item generation, rhizosphere biochemical environment evolution and hydraulic constitutive correction and fluid transport solving modules, a finite difference grid is constructed in the system, and a root exudate release source item is calculated based on a current water state; exudate migration is simulated to reconstruct an effective dynamic viscosity of fluid and a dynamic contact angle of a sand particle surface; the soil unsaturated hydraulic conductivity and the matric potential are corrected according to the reconstructed parameters, a fluid transport control equation is solved to update the water distribution, and the state is fed back to the ecological module to form a closed loop. Through simulation of the rheological resistance and infiltration shielding effect induced by root exudate, the present application solves the problem that the traditional model ignores the dynamic evolution of the physical properties of the medium, and improves the physical authenticity of the simulation of the hydrological process in the arid region.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for eco-hydrological processes, specifically a two-way coupled modeling system for distributed eco-hydrological processes in sandy areas. Background Technology

[0002] In ecohydrological research in arid and semi-arid regions, numerical simulation is a key tool for understanding the interaction between vegetation and water resources. Existing ecohydrological models of sandy areas typically employ a one-way driving mechanism, focusing on calculating how soil moisture distribution limits or promotes vegetation growth. In the construction of these models, the porous soil medium is often simplified as an abiotic, inert system, with its hydraulic physical parameters, such as saturated hydraulic conductivity and matrix potential parameters, set as constants or functions that only change with water content, neglecting the bidirectional biophysical feedback between vegetation roots and the soil medium. This assumption of treating the physical properties of the medium as static prevents current techniques from reflecting the adaptive mechanisms by which vegetation actively alters the physical properties of the rhizosphere microenvironment by releasing substances into the rhizosphere when subjected to water stress.

[0003] In the calculation of fluid transport, conventional hydrological models are usually based on the Richards equation and assume that the pore fluid is pure water or a Newtonian fluid with constant viscosity. However, vegetation root exudates typically contain high molecular weight polymers, and their accumulation can significantly alter the rheological properties of pore fluids. Existing techniques fail to consider the nonlinear effect of root exudate concentration on the effective dynamic viscosity of the fluid and neglect the hindering effect of high-viscosity fluids on water conduction. Due to the lack of description of this non-Newtonian fluid characteristic, existing models struggle to quantitatively simulate water retention in the rhizosphere caused by increased fluid viscosity, resulting in insufficient accuracy in simulating rhizosphere water movement processes.

[0004] Furthermore, in simulating the properties of media interfaces, traditional theories typically presuppose that the sand grain surface remains hydrophilic, i.e., the contact angle is less than 90 degrees. Existing models do not include a description of the dynamic evolution of sand grain surface wettability, particularly failing to cover the physical process by which root exudates solidify into a biofilm during drying, leading to a reversal of the media surface from a hydrophilic to a hydrophobic state. Due to the lack of consideration for this mechanism, current techniques cannot calculate the repulsive potential generated by soil hydrophobicity, thus failing to accurately simulate the infiltration shielding effect and the significantly reduced water re-infiltration rate observed in dry, hydrophobic soils, resulting in biased predictions of rainfall infiltration and redistribution processes in arid regions. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a two-way coupled modeling system for distributed eco-hydrology in sandy areas. This system solves the problem that existing sandy hydrological models typically treat soil as an abiotic inert medium, neglecting the dynamic correction of soil pore fluid rheological properties and medium interface wetting performance by vegetation root exudates. As a result, they cannot accurately simulate the complex biophysical feedback mechanism between sandy ecosystems and hydrological processes under drought stress.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] A distributed eco-hydrological bidirectional coupled modeling system for sandy areas includes an initialization and discretization module, an ecological growth calculation and source term generation module, a rhizosphere biochemical environment evolution module, and a hydraulic constitutive correction and fluid transport solution module.

[0008] The initialization and discretization module is configured to construct a finite difference grid system in a three-dimensional physical space domain, and sets the initial state parameters of the system, including the effective dynamic viscosity of the fluid, the dynamic contact angle of the sand grain surface, the unsaturated hydraulic conductivity of the soil, and the matrix potential.

[0009] The ecological growth calculation and source term generation module is configured to generate an environmental water stress signal based on the current water distribution state within the finite difference grid system, and convert the environmental water stress signal into a root exudate release source term. The function of this module is to calculate the material flux input from the vegetation to the rhizosphere environment based on the current water conditions.

[0010] The rhizosphere biochemical environment evolution module is configured to receive the root exudate release source, simulate the exudate transport process, calculate the spatiotemporal distribution of biochemical substances in soil pores, and reconstruct the effective dynamic viscosity of the fluid and the dynamic contact angle of the sand grain surface based on the concentration distribution.

[0011] The hydraulic constitutive correction and fluid transport solution module is configured to correct the unsaturated hydraulic conductivity of the soil using the effective dynamic viscosity of the fluid and the matrix potential using the dynamic contact angle. Based on the corrected parameters, the fluid transport control equation is solved, the water distribution state is updated, and the updated water distribution state is fed back to the ecological growth calculation and source term generation module, forming a dynamic closed-loop coupling between the sandy ecological process and the hydrological process.

[0012] Preferably, the specific execution process of the initialization and discretization module is as follows: the continuous three-dimensional physical space domain is divided into an orthogonal finite difference grid system, and each grid cell is assigned initial volumetric water content, initial matrix potential, saturated water content, residual water content, initial saturated hydraulic conductivity, initial effective dynamic viscosity of the fluid, and initial dynamic contact angle of the sand grain surface.

[0013] Preferably, the specific calculation logic of the ecological growth calculation and source term generation module is as follows: extract the current volumetric water content from the current water distribution state, calculate the effective saturation using the current volumetric water content, map the effective saturation to the environmental water stress signal, and calculate the root exudate release source term in combination with the vegetation root length density distribution. In the calculation model setting, as the effective saturation representing the environmental water stress signal decreases, the root exudate release source term exhibits a non-linear exponential growth, thereby reflecting the physiological response of vegetation to increase exudate release under drought conditions.

[0014] Preferably, the rhizosphere biochemical environment evolution module constructs governing equations based on convection, diffusion, and reaction mechanisms when simulating the exudate transport process. These governing equations consider the convection carried by water flow, the hydrodynamic diffusion driven by concentration gradients, and biochemical degradation. Specifically, an Arrhenius temperature correction function is introduced to calculate the real-time adjustment of the exudate degradation rate by ambient temperature.

[0015] Preferably, the process by which the rhizosphere biochemical environment evolution module reconstructs the effective dynamic viscosity of the fluid is as follows: A concentration-dependent nonlinear rheological model is used to map the root exudate concentration obtained from simulating the exudate transport process to the effective dynamic viscosity of the fluid. The nonlinear rheological model defines the transformation of pore fluid from a Newtonian fluid to a non-Newtonian viscous fluid as the root exudate concentration increases, thereby calculating the change in fluid flow resistance caused by the increase in solute concentration.

[0016] Preferably, the process by which the rhizosphere biochemical environment evolution module reconstructs the dynamic contact angle of the sand grain surface is as follows: Historical state variables are introduced to record the drying history of the grid cells. When the concentration of root exudates obtained from simulating the exudate transport process is higher than the critical film-forming concentration threshold and the current volumetric water content in the current water distribution state is lower than the critical drying threshold, a biofilm solidification event is determined to have occurred. This biofilm solidification event triggers the dynamic contact angle to increase with the increase of the root exudate concentration until a wettability reversal occurs from a hydrophilic to a hydrophobic state.

[0017] Preferably, the hydraulic constitutive correction and fluid transport solution module corrects the matrix potential by introducing a wettability correction factor that includes the cosine value of the dynamic contact angle to scale and adjust the original matrix potential. When the dynamic contact angle undergoes a wettability reversal from a hydrophilic to a hydrophobic state, the sign of the wettability correction factor changes, causing the corrected matrix potential to change from a negative attractive force to a positive repulsive force, forming an infiltration barrier effect that prevents water from entering the pores.

[0018] Preferably, the hydraulic constitutive correction and fluid transport solution module corrects the unsaturated hydraulic conductivity by introducing a rheological retardation factor composed of the ratio of the dynamic viscosity of pure water to the effective dynamic viscosity of the fluid. The rheological retardation factor is used to correct the relative hydraulic conductivity that only considers the pore geometry. The corrected unsaturated hydraulic conductivity is set to be negatively correlated with the effective dynamic viscosity of the fluid, thereby forming a rheological retardation effect that inhibits water conduction.

[0019] Preferably, the process of the hydraulic constitutive correction and fluid transport solution module solving the fluid transport control equations is as follows: a generalized Richards equation is constructed, and the corrected unsaturated hydraulic conductivity and the corrected matrix potential are substituted. Spatial discretization is performed using the fully implicit finite volume method, and an iterative algorithm is used to linearize and solve the discretized equations to obtain the water content distribution for the next time step. Furthermore, the hydraulic constitutive correction and fluid transport solution module integrates an adaptive time step control strategy, which automatically reduces the simulation time step when the iterative algorithm fails to converge after exceeding a preset threshold, to adapt to the nonlinear calculation requirements of the wettability inversion stage.

[0020] Through the above-mentioned technical solution, this invention integrates two physical mechanisms in the numerical model: rheological retardation induced by root exudates and infiltration shielding. It simulates the process by which sandy vegetation regulates water transport by changing the physical properties of the rhizosphere microenvironment, reflecting the physical reality of the eco-hydrological system in arid areas.

[0021] This invention provides a two-way coupled modeling system for distributed eco-hydrology in sandy areas. It has the following beneficial effects:

[0022] 1. This invention constructs a closed-loop feedback mechanism between ecological growth and hydrological processes, transforming the environmental water stress signal triggered by the current water distribution state into a source term for root exudate release, and dynamically updating soil physical parameters accordingly. This solves the defect of traditional models that treat the soil medium as a static system with constant physical properties, and can truly reflect the biophysical process by which vegetation in arid areas actively changes the rhizosphere environment through physiological activities to adapt to water stress.

[0023] 2. This invention establishes a concentration-dependent nonlinear rheological model to map the concentration of root exudates to the effective dynamic viscosity of the fluid, and uses a rheological retardation factor to correct the unsaturated hydraulic conductivity. This method introduces non-Newtonian fluid characteristics into the numerical calculation, which can quantitatively simulate the rheological retardation effect of root exudates reducing the water conduction rate by increasing the viscosity of the pore fluid, thereby improving the simulation accuracy of water retention in the rhizosphere.

[0024] 3. This invention reconstructs the dynamic contact angle of sand grains and simulates wettability flipping by introducing historical state variables to monitor the drying history and material accumulation of grid cells. By adjusting the sign of the matrix potential using a wettability correction factor, it can accurately describe the biofilm solidification and the resulting infiltration shielding effect that occur when the soil is dry and the exudate concentration is high. This corrects the prediction bias of the existing model for water infiltration rate under hydrophobic soil conditions. Attached Figure Description

[0025] Figure 1 This is a schematic diagram of the overall architecture of a distributed eco-hydrological bidirectional coupling modeling system for sandy areas according to an embodiment of the present invention.

[0026] Figure 2 This is a schematic diagram of the data processing flow of the ecological growth calculation and source term generation module according to an embodiment of the present invention;

[0027] Figure 3 This is a schematic diagram of the numerical coupling solution and discretization implementation process of a two-way coupled modeling system for distributed eco-hydrology in sandy areas according to an embodiment of the present invention.

[0028] The modules are: 101. Initialization and Discretization Module; 102. Ecological Growth Calculation and Source Term Generation Module; 103. Rhizosphere Biochemical Environment Evolution Module; and 104. Hydraulic Constitutive Correction and Fluid Transport Solution Module. Detailed Implementation

[0029] The technical solutions in 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.

[0030] See attached document Figure 1 This invention provides a distributed eco-hydrological bidirectional coupling modeling system for sandy areas. The system is configured to simulate the mechanism by which sandy vegetation, through root exudates, alters soil rheological properties and interfacial wettability, thereby inversely influencing hydrological processes. The system includes: an initialization and discretization module 101, an ecological growth calculation and source term generation module 102, a rhizosphere biochemical environment evolution module 103, and a hydraulic constitutive correction and fluid transport solution module 104. These modules interact via a data bus or memory sharing mechanism, executing cyclically according to a preset time step, collectively constructing a closed-loop, dynamically coupled simulation environment.

[0031] The initialization and discretization module 101 is configured to define the three-dimensional physical space domain and initial state of the sand system to be simulated. The initialization and discretization module 101 divides the continuous three-dimensional physical space domain into an orthogonal finite difference grid system, with each grid cell identified by a unique spatial coordinate index. The initialization and discretization module 101 assigns initial hydrophysical parameters to each grid cell, including initial volumetric water content, initial matrix potential, saturated water content, residual water content, and initial saturated hydraulic conductivity of the sand medium. Furthermore, the initialization and discretization module 101 sets the fundamental physical constants of the fluid, including the dynamic viscosity of pure water and the intrinsic contact angle of the sand grain surface. The initialization and discretization module 101 is also responsible for setting the time control parameters of the numerical simulation, including the initial time step, maximum allowable time step, minimum time step, and total simulation duration, and initializes all system state variables to the initial state.

[0032] The ecological growth calculation and source term generation module 102 is configured to calculate the vegetation growth dynamics and root exudate release flux at the beginning of each simulation time step, based on the current water distribution state of the system. The ecological growth calculation and source term generation module 102 is connected to the hydraulic constitutive correction and fluid transport solution module 104, receiving soil moisture content data from the previous time step or iteration step, and calculating the vegetation transpiration rate and root water uptake rate accordingly. The ecological growth calculation and source term generation module 102 updates the root length density distribution in three-dimensional space according to a preset vegetation growth model. Furthermore, the ecological growth calculation and source term generation module 102 is further configured to calculate root exudate source terms. Root exudate source terms characterize the mass of mucus released from roots into soil pores per unit time and unit volume. The calculation process of root exudate source terms considers the nonlinear coupling relationship between root length density and water stress, i.e., by monitoring whether the water content within the grid is below a specific threshold to determine whether drought stress is triggered, and then dynamically adjusting the release intensity of exudates.

[0033] The rhizosphere biochemical environment evolution module 103 is configured to receive root exudate source terms output by the ecological growth calculation and source term generation module 102, simulating the material transport process of exudates and the resulting changes in the physical properties of the medium in porous media. The rhizosphere biochemical environment evolution module 103 first calculates the root exudate concentration field at each grid node based on convection, diffusion, and reaction mechanisms. This process considers convection carried by water flow, diffusion driven by concentration gradients, and biodegradation. Subsequently, based on the calculated concentration field and the wet / dry history of the grid, the rhizosphere biochemical environment evolution module 103 reconstructs the rheological parameters of the fluid and the surface physical parameters of the medium in real time. Regarding rheological parameters, the rhizosphere biochemical environment evolution module 103 treats pore water as a non-Newtonian fluid and dynamically updates the effective dynamic viscosity of the fluid according to the exudate concentration. Regarding surface physical parameters, the rhizosphere biochemical environment evolution module 103 introduces historical state variables to record the drying history of the grid. When a high concentration of secretions is detected and the water content is lower than the critical drying threshold, a biofilm solidification event is determined to have occurred, and then the dynamic contact angle is calculated to realize the wettability reversal logic from hydrophilic to hydrophobic state.

[0034] The hydraulic constitutive correction and fluid transport solution module 104 is configured to receive the effective dynamic viscosity and dynamic contact angle output by the rhizosphere biochemical environment evolution module 103, perform dual corrections on the fundamental constitutive relations controlling soil moisture movement, and solve the corrected governing equations to obtain a new flow field distribution. The hydraulic constitutive correction and fluid transport solution module 104 performs two core correction operations: First, it corrects the unsaturated hydraulic conductivity so that it depends not only on saturation but also is inversely proportional to the effective dynamic viscosity of the fluid, thus reflecting the flow obstruction effect caused by viscous thickening; second, it corrects the matrix potential by adjusting the capillary force using the dynamic contact angle. When the contact angle reverses, it simulates the transformation of matrix suction into hydrophobic potential, thus reflecting the infiltration barrier effect caused by biofilm. The hydraulic constitutive correction and fluid transport solution module 104 substitutes the corrected hydraulic conductivity and matrix potential into the unsaturated fluid transport equation (generalized Richards equation) and uses a numerical iterative algorithm to solve for the water content distribution and pressure head distribution at the next time step. After the solution is completed, the hydraulic constitutive correction and fluid transport solution module 104 feeds back the updated water state data to the ecological growth calculation and source term generation module 102 as the input condition for the next time step, thereby completing the two-way coupled cyclic simulation of the sandy land eco-hydrological process.

[0035] See attached document Figure 1 and attached Figure 2 The ecological growth calculation and source term generation module 102, based on the principles of plant physiological ecology, converts environmental water stress signals into quantitative biochemical release signals, providing material source term inputs for subsequent soil rheology and wettability evolution simulations.

[0036] The ecological growth calculation and source term generation module 102 first obtains the current simulation time. The moisture status of each grid cell is calculated. The ecological growth calculation and source term generation module 102 reads the moisture status of each spatial grid cell. The volumetric water content data within the plant is used to calculate the effective saturation, which serves as a standardized indicator for measuring the degree of drought in the environment in which plant roots reside. The calculation formula is:

[0037] ;

[0038] in, It represents the effective saturation, which is a dimensionless value between 0 and 1, reflecting the relative position of the current soil moisture content between the residual moisture content and the saturated moisture content; This indicates the volumetric water content of the current grid cell; Residual moisture content refers to the maximum amount of water that remains adsorbed on the surface of soil particles after natural evaporation or drainage and cannot be discharged. Saturated water content refers to the volumetric water content when the soil pores are completely filled with water.

[0039] While determining the water status, the ecological growth calculation and source term generation module 102 performs dynamic updates on the spatial distribution of the root system. Considering the non-uniform vertical distribution of roots in psammophytes, the ecological growth calculation and source term generation module 102 uses an exponential decay model to determine the root length density at each depth layer. :

[0040] ;

[0041] in, This represents the root length density at the calculated node, which is the total root length per unit volume of soil. This represents the maximum root length density at the ground surface, and this parameter is updated according to the vegetation growth curve over time during the simulation. Represents the natural exponential function; This represents the attenuation coefficient of the vertical distribution of the root system, used to characterize the rate at which the root system decreases with increasing depth; This indicates the vertical depth of the calculation node from the Earth's surface.

[0042] Based on the above calculations, the effective saturation is obtained. With root length density The ecological growth calculation and source term generation module 102 further calculates the source term of root exudate release. This step quantifies the rate at which plant roots release mucilage into the soil matrix under different moisture conditions. This computational model not only considers the basic contribution of root biomass but also specifically incorporates the excitation effect of water stress, simulating the physiological mechanism by which plants alter the rhizosphere microenvironment by increasing exudate release under drought conditions. Root exudate release source term. The calculation formula is as follows:

[0043] ;

[0044] in, This represents the source term of root exudate release, defined as the increase in the mass of exudate per unit volume of soil per unit time; The basal secretion rate coefficient per unit root length represents the background secretion intensity under conditions of sufficient water (no stress). This represents the root length density at the calculated node, which is the total root length per unit volume of soil. This represents the base ratio under no-stress conditions; This represents the maximum stress gain coefficient, indicating the maximum increase in secretion rate relative to the baseline rate under extreme drought conditions; Represents the natural exponential function; This represents a moisture-sensitive factor that controls the secretion rate as a function of effective saturation. Reduce the sensitivity of the response. The higher the value, the more severe the plant's response to mild drought; It represents the effective saturation, which is a dimensionless value between 0 and 1, reflecting the relative position of the current soil moisture content between the residual moisture content and the saturated moisture content.

[0045] In this embodiment, the specific values ​​of the model parameters depend on the specific vegetation type (such as Artemisia arenaria, Caragana korshinskii, etc.) and soil texture. Specifically, the parameters... , as well as The data was obtained through root box control experiments. In these experiments, different soil moisture gradients were set, and the carbon secretion of the roots per unit time was collected and measured. The resulting data was then obtained by fitting the exponential stress response equation using a nonlinear least squares method. The root vertical attenuation coefficient was also determined. The coefficient is determined through statistical regression analysis of field root excavation data, and it is dynamically adjusted as the vegetation grows in different stages to reflect the growth strategy of roots penetrating deeper into the soil.

[0046] Using this formula, the ecological growth calculation and source term generation module 102 generates a three-dimensional source term field that changes in real time with spatial location and moisture conditions. The three-dimensional source term field data is then transmitted to the rhizosphere biochemical environment evolution module 103 as the source term input in the material transport equation, thereby driving the evolution of the physical properties of the underground rhizosphere microenvironment.

[0047] See attached document Figure 1 The rhizosphere biochemical environment evolution module 103 is configured to receive root exudate release source term data output by the ecological growth calculation and source term generation module 102, and simulate the dynamic transport process of the exudate as a solute in an unsaturated porous medium. This process quantifies the accumulation, diffusion, and biodegradation behavior of root exudates in the soil space, thereby outputting a root exudate concentration field that changes over time. .

[0048] Module 103, which describes the evolution of the rhizosphere biochemical environment, constructs a solute transport governing equation based on convection, diffusion, and reaction mechanisms. This equation considers not only water-driven convective transport but also hydrodynamic diffusion transport driven by concentration gradients, as well as source-sink terms induced by biochemical processes. In a three-dimensional spatial grid system, the governing equations describing changes in root exudate concentration are as follows:

[0049] ;

[0050] in, It represents the rate of change of the mass of root exudates per unit volume of soil over time, reflecting the local accumulation or consumption of substances. This indicates the volumetric water content of the current grid cell; The concentration of root exudates in soil pore water is defined as the mass of exudates per unit volume of pore water. Indicates the simulated time; This represents the Hamiltonian operator (gradient operator). The hydrodynamic dispersion coefficient tensor is a second-order tensor that encompasses both molecular diffusion coefficient and mechanical dispersion coefficient, and is used to characterize the anisotropic diffusion capacity of solutes in sandy media. Represents the spatial gradient of concentration; This represents the hydrodynamic dispersion term, which describes the diffusion flux of solutes under the combined effects of flow rate differences and molecular diffusion in micropores. This represents the Darcy velocity vector, which is the amount of water passing through a unit cross-section per unit time, and its direction is consistent with the direction of the hydraulic gradient. This represents the convective transport term, which describes the flux of solute as it migrates along the average pore water flow. This represents the source term of root exudate release, defined as the increase in the mass of exudate per unit volume of soil per unit time; This represents the first-order biodegradation rate constant, which depends on soil temperature and microbial activity and reflects the persistence of exudates in the soil environment. This indicates a biodegradation sink, describing the process by which root exudates are decomposed and consumed by soil microorganisms.

[0051] It is worth noting that, considering the dramatic diurnal temperature variations characteristic of sandy environments, the biodegradation rate constant... It is not a fixed constant, but is updated in real time using the Arrhenius temperature correction function.

[0052] By numerically solving the above partial differential equations, the rhizosphere biochemical environment evolution module 103 updates the concentration distribution throughout the entire simulation domain at each time step. The calculated concentration field data not only reflects the local enrichment effect of exudates around the rhizosphere, but also their leaching and migration process as water infiltrates into deeper soil layers, providing a material basis for subsequent calculations of fluid viscosity changes and determination of interfacial wettability inversion.

[0053] See attached document Figure 1 After completing the transport calculation of root exudate concentration, the rhizosphere biochemical environment evolution module 103 is further configured to reconstruct the rheological physical parameters of pore fluids based on real-time concentration data at each grid node. This step aims to quantify the physical process by which root exudates dissolve in soil water, causing the pore fluid to transform from an ideal Newtonian fluid to a non-Newtonian viscous fluid.

[0054] The rhizosphere biochemical environment evolution module 103 reads the root exudate concentration field data calculated in the previous step. Since root exudates are mainly composed of high-molecular-weight polysaccharides and glycoproteins, their accumulation in aqueous solutions significantly increases the internal friction of the fluid. Therefore, the rhizosphere biochemical environment evolution module 103 abandons the assumption in traditional hydrological models that the pore water viscosity is constant at the level of pure water, and adopts a concentration-dependent nonlinear rheological model to dynamically calculate the effective dynamic viscosity of each grid cell.

[0055] Effective dynamic viscosity The calculation follows the following rheological constitutive equation:

[0056] ;

[0057] in, This represents the effective dynamic viscosity of fluids with secretions in porous media, and this parameter directly determines the flow resistance of the fluid in porous media. This represents the dynamic viscosity of pure water at the current soil temperature, serving as a benchmark reference value for rheological calculations. This represents the baseline state when no secretions are present. The rheological thickening coefficient is a parameter that depends on the molecular weight and molecular structure of the secretion and is used to characterize the ability of a specific plant secretion to change the viscosity of a fluid. The concentration of root exudates in soil pore water is defined as the mass of exudates per unit volume of pore water. The concentration power law exponent is used to describe the non-linear, sharp increase in viscosity as concentration increases.

[0058] Through the calculation using this formula, the rhizosphere biochemical environment evolution module 103 will convert the scalar root exudate concentration field Viscosity field mapped to tensor or scalar form Viscosity field data accurately reflects the significant differences in fluid properties between the rhizosphere microregion (high concentration zone) and the non-rhizosphere region (low concentration zone). In the rhizosphere region, due to the presence of high-concentration exudates, the calculated... The value is significantly higher than This localized high viscosity characteristic will manifest as a strong impediment to water movement in subsequent hydraulic calculations, thus physically simulating the behavior of plants locking in rhizosphere water through secretion of mucus. The calculated and updated effective dynamic viscosity data is stored in the system state matrix for subsequent hydraulic constitutive correction and fluid transport solution module 104 to correct the hydraulic conductivity parameters.

[0059] See attached document Figure 1 While simultaneously updating the rheological properties, the rhizosphere biochemical environment evolution module 103 is further configured to monitor the wet-dry alternation history of the soil pore media and determine whether a biofilm altering the wetting properties has formed on the sand grain surface, thereby calculating the dynamic contact angle. This step simulates the physicochemical mechanism by which root exudates transform from a hydrophilic gel state to a hydrophobic glass state after undergoing a dehydration and drying process.

[0060] The rhizosphere biochemical environment evolution module 103 first introduces a historical state variable for recording the physical state history of grid cells. Historical state variables are used to identify whether the sand surface at a specific spatial location has undergone a drying event that induces biofilm formation. At the initial moment of the simulation, all grid cells... The default value is set to 0, representing that the sand surface is in its original clean state. At each time step, the rhizosphere biochemical environment evolution module 103 simultaneously reads the volumetric water content of the current grid. and root exudate concentration And execute the following state determination logic: if the current grid's volumetric water content is... Below the preset critical drying threshold And at the same time, the concentration of root exudates within this grid Higher than the preset critical film-forming concentration threshold If the biofilm solidification event has occurred in the grid, then the historical state variables will be determined. Updated to 1. Once Set to 1, this state will remain in the subsequent simulation steps until a specific long-term soaking and rewetting event occurs (optional reset condition) or the simulation ends, reflecting the hysteresis effect after biofilm formation.

[0061] Based on the updated historical state variables Based on the current concentration field data, the rhizosphere biochemical environment evolution module 103 calculates the dynamic contact angle of each grid cell. The dynamic contact angle is a key parameter for measuring the wetting performance of the fluid-solid interface. The rhizosphere biochemical environment evolution module 103 uses the following piecewise function model to describe the dynamic evolution of the contact angle:

[0062] ;

[0063] in, Indicates the dynamic contact angle, in radians or degrees; The intrinsic contact angle of a sand grain corresponds to the surface of clean sand that is not contaminated by secretions, and is usually less than 90 degrees. This represents the limiting contact angle after the plant is completely covered by a biofilm of dried secretions. Based on the hydrophobic properties of the secretions of psammophytes, this value is set to be greater than 90 degrees. This represents a historical state variable that records the physical state history of a grid cell. The concentration of root exudates in soil pore water is defined as the mass of exudates per unit volume of pore water. This represents the half-saturation constant, used to control the contact angle as it approaches a limiting value with increasing concentration. The rate.

[0064] This calculation shows that a high concentration of secretions alone is not sufficient to immediately change the wettability of the medium (if...). The dynamic contact angle remains unchanged. It must undergo a film-forming drying process. Only when the ratio is flipped to 1 will the hydrophobic modification potential of the secretion be activated. After activation, the dynamic contact angle... Size and current concentration of residual root exudates They are positively correlated. When When the calculated angle exceeds 90 degrees, it physically signifies a reversal of the soil medium from hydrophilic to hydrophobic. The calculated dynamic contact angle data... The output is sent to the hydraulic constitutive correction and fluid transport solution module 104 to correct the relationship curve between matrix potential and saturation, thereby introducing the capillary force direction reversal and infiltration shielding effect caused by biofilm into the model.

[0065] See attached document Figure 1The hydraulic constitutive correction and fluid transport solution module 104 is configured to receive dynamic contact angle data output by the rhizosphere biochemical environment evolution module 103, and use this data to correct the matrix potential model describing the soil moisture and energy state. This process aims to map the wettability changes of the micro-interface onto the macro-water characteristic curve, thereby reproducing the soil water repellency and hydraulic hysteresis loop distortion phenomena caused by biofilm in the numerical model.

[0066] The hydraulic constitutive correction and fluid transport solution module 104 executes correction logic based on capillary pressure theory. In porous media physics, matrix potential (or capillary suction) is inversely proportional to the pore radius and directly proportional to the cosine of the contact angle at the fluid-solid interface. The hydraulic constitutive correction and fluid transport solution module 104 first calculates the original matrix potential unaffected by biofilms based on the current effective saturation using the standard Van Genuchten model. The original matrix potential represents the energy state required to maintain the current moisture content in a clean, hydrophilic sandy medium.

[0067] Subsequently, the hydraulic constitutive correction and fluid transport solution module 104 introduces a wettability correction factor to scale and adjust the original matrix potential, thus obtaining the corrected matrix potential. The mathematical expression of this correction process is as follows:

[0068] ;

[0069] in, This represents the corrected matrix potential, expressed in units of length, and in unsaturated soil physics, it typically corresponds to pressure head. This represents the inverse parameter of the intake suction in the Van Genuchten model, which is related to the average pore size of the soil. It represents the effective saturation, which is a dimensionless value between 0 and 1, reflecting the relative position of the current soil moisture content between the residual moisture content and the saturated moisture content; Indicates the shape parameters of the pore distribution; Denotes the porosity distribution index parameter, and satisfies The constraint relationship; Indicates the dynamic contact angle, in radians or degrees; Indicates the intrinsic contact angle of the sand grain surface; This represents the cosine trigonometric function.

[0070] Through this corrected formula, the hydraulic constitutive correction and fluid transport solution module 104 achieves dynamic reconstruction of the soil moisture characteristic curve. When the system is in a normal hydrophilic state (i.e., The correction factor is close to 1. When the value remains negative, it indicates normal capillary suction, allowing water to be spontaneously drawn into the pores. However, when a biofilm-induced wetting flip event occurs (i.e., ... ), dynamic contact angle cosine value The value becomes negative, resulting in the corrected matrix potential. The value abruptly changes from negative to positive. Physically, a positive matrix potential means that the surface of soil particles exerts a repulsive force on water molecules; the pores no longer generate suction, and instead require external positive pressure (i.e., water pressure) to force water into the pores. This sign reversal of potential energy caused by the change in contact angle manifests on the moisture characteristic curve as a sharp separation in the shape of the desiccation curve and the hygroscopic curve, i.e., a severe distortion of the hysteresis loop. This mechanism forms a hydraulic barrier in the model, accurately simulating the physical effects of the rhizosphere biofilm preventing water re-infiltration and cutting off capillary upwelling after drying.

[0071] See attached document Figure 1 The hydraulic constitutive correction and fluid transport solution module 104, while correcting the matrix potential, is further configured to correct the unsaturated hydraulic conductivity of the soil medium based on changes in fluid rheological properties. This step aims to map the fluid viscosity changes calculated by the rhizosphere biochemical environment evolution module 103 onto the hydraulic conductivity, thereby reproducing the water transport hindrance caused by the viscosity-enhancing effect of root exudates in the numerical model.

[0072] The hydraulic constitutive correction and fluid transport solution module 104 executes correction logic based on the seepage theory of porous media. The hydraulic conductivity of soil, in its physical essence, depends not only on the pore geometry of the medium (such as porosity, pore size distribution, and connectivity) but also on the physical properties of the fluid itself (primarily fluidity, i.e., the reciprocal of viscosity). The hydraulic constitutive correction and fluid transport solution module 104 first uses the Mualem-Van Genuchten model to calculate the relative hydraulic conductivity considering only the influence of geometry, and then introduces a viscosity correction factor to construct the corrected unsaturated hydraulic conductivity. Calculation formula:

[0073] ;

[0074] in, This represents the corrected unsaturated hydraulic conductivity, expressed in length-to-time ratio. This parameter characterizes the soil medium's ability to transport water at the current water content and fluid viscosity. It represents the initial saturated hydraulic conductivity of the soil medium, corresponding to the rate at which pure water is conducted through clean sand in a saturated state; It represents the effective saturation, which is a dimensionless value between 0 and 1, reflecting the relative position of the current soil moisture content between the residual moisture content and the saturated moisture content; Indicates the shape parameters of the pore distribution; This represents the dynamic viscosity of pure water at the current soil temperature, serving as a benchmark reference value for rheological calculations. This represents the effective dynamic viscosity of fluids with secretions in porous media, and this parameter directly determines the flow resistance of the fluid in porous media. This characterizes the decrease in hydraulic conductivity caused by the reduction in the number of water-filling pores and the increase in the tortuosity of the water flow path as the soil moisture content decreases; this belongs to the geometric resistance part. This constitutes the rheological retardation factor, which characterizes the direct modulation effect of fluid viscosity changes on flow velocity.

[0075] Through this corrected formula, the hydraulic constitutive correction and fluid transport solution module 104 achieves real-time coupling of rheological and hydraulic properties. When the root exudate concentration is low, near The rheological resistance factor is close to 1, and the hydraulic conductivity is mainly controlled by water content. When high concentrations of exudates accumulate in the rhizosphere, leading to... When the value increases significantly, the rheological retardation factor becomes a decimal much smaller than 1, leading to The decrease is orders of magnitude. Physically, this means that even with high soil moisture content, the flowability of pore fluids is significantly suppressed because they become a high-viscosity gel-like substance. This mechanism accurately simulates in the model the biophysical process by which plants use mucilage to form a hydraulically viscous band around the rhizosphere to reduce water loss or alter water infiltration pathways. (Corrected) The data were then compared with the corrected matrix potential. They are substituted into the fluid transport control equations to calculate the flow field distribution at the next moment.

[0076] See attached document Figure 1 The hydraulic constitutive correction and fluid transport solution module 104 calculates the corrected unsaturated hydraulic conductivity obtained above. With the corrected matrix potential By substituting the fluid motion equations into the porous media, the two physical effects of rheological hindrance and infiltration shielding were quantified at the numerical solution level. These two effects together constitute the water regulation mechanism of the rhizosphere microenvironment of sandy vegetation, which is manifested in the following ways.

[0077] The hydraulic constitutive correction and fluid transport solution module 104 first constructs the modified Darcy's law equation describing the flux of soil moisture movement. This equation clarifies the contribution relationship of each physical quantity to the flow velocity:

[0078] ;

[0079] in, This represents the soil moisture flux vector, where the magnitude represents the amount of water passing through a unit area per unit time, and the direction represents the direction of water flow. This represents the corrected unsaturated hydraulic conductivity, expressed in length-to-time ratio. This parameter characterizes the soil medium's ability to transport water at the current water content and fluid viscosity. This represents the Hamiltonian operator (gradient operator). This represents the corrected matrix potential, expressed in units of length, and in unsaturated soil physics, it typically corresponds to pressure head. This indicates the vertical depth of the calculation node from the Earth's surface.

[0080] Within this equation framework, the rheological retardation effect is achieved through the modified unsaturated hydraulic conductivity. The viscosity term is reflected in the data. With effective dynamic viscosity They are inversely proportional. When there is a high concentration of root exudates in the rhizosphere, The value increased sharply, leading to The numerical value decreases significantly due to nonlinearity. In the flux equation, because... As a multiplication coefficient, its decrease directly affects the water flux vector. The modulus decreases proportionally. Physically, this indicates that even with a large hydraulic gradient (i.e., (Very large), due to the high viscosity of the fluid itself, the actual migration rate of water in the pores of the medium is still forcibly limited to an extremely low level. In the model, this mechanism not only simulates the physical interception of water seepage into deeper layers by rhizosphere exudates, but also simulates its inhibition of lateral water diffusion, thus effectively locking water around the roots.

[0081] The infiltration shielding effect is achieved through the matrix potential. The contact angle correction term is reflected in this. Under normal hydrophilic conditions, Negative value The direction points towards the wetting front, driving the spontaneous absorption of moisture into the dry area. However, when the rhizosphere biochemical environment evolution module 103 determines that a biofilm-induced wetting reversal event has occurred, the dynamic contact angle... exceeding 90 degrees, making The attractive force abruptly changes from a negative value to a positive value of repulsive potential. In the flux equation, this sign reversal alters the direction of the potential energy gradient. Specifically, the positive value... A potential energy barrier forms at the soil surface or the wet-dry interface. For an inward positive flux (i.e., infiltration) to occur, the external water flow must provide sufficient positive pressure head to counteract this potential energy barrier (i.e., overcome it). (The repulsive force generated). If the water pressure generated by rainfall or irrigation is insufficient to overcome this repulsive potential, water will be forced to accumulate or run off the surface and will not be able to enter the soil pores. This mechanism accurately reproduces the resistance of the dry rhizosphere biofilm to the rewetting process in the model, explaining why sandy soils experience temporary water infiltration difficulties after drought.

[0082] The hydraulic constitutive correction and fluid transport solution module 104 solves the above flux relationship in real time at each spatial grid and each time step, dynamically mapping microscopic rheological changes and interface flips into changes in macroscopic flow field distribution, thereby realizing direct physical feedback of ecological processes to hydrological processes.

[0083] See attached document Figure 3 After completing the physical process description of source term generation, material transport and constitutive relation correction, the system performs numerical discretization and coupled solution through the hydraulic constitutive correction and fluid transport solution module 104 to obtain accurate system state evolution trajectory in discrete time and space domains.

[0084] The hydraulic constitutive correction and fluid transport solution module 104 is configured to solve the generalized Richards equations, which integrate all the correction parameters calculated by the aforementioned modules. To accommodate the strong nonlinearity and hysteresis effects caused by biofilms, the hydraulic constitutive correction and fluid transport solution module 104 adopts the following modified governing equation form:

[0085] ;

[0086] in, It represents the rate of change of soil volumetric water content over time, that is, the change in water storage per unit time. This indicates the volumetric water content of the current grid cell; Indicates the simulated time; This represents the Hamiltonian operator (gradient operator). This represents the corrected unsaturated hydraulic conductivity, expressed in length-to-time ratio. This parameter characterizes the soil medium's ability to transport water at the current water content and fluid viscosity. This represents the corrected matrix potential, expressed in units of length, and in unsaturated soil physics, it typically corresponds to pressure head. This indicates the vertical depth of the calculation node from the Earth's surface; This represents the amount of water absorbed by the roots per unit time and per unit volume of soil by the plant roots, which is then lost through transpiration.

[0087] To solve the above equations, the hydraulic constitutive correction and fluid transport solution module 104 needs to incorporate the boundary conditions of the simulation domain. In this embodiment, the top boundary of the simulation domain is set as a flux boundary, the flux value of which is determined by the atmospheric rainfall intensity and potential evapotranspiration capacity, and it is allowed to automatically switch to a head boundary based on the surface water accumulation; the bottom boundary of the simulation domain is set as a free drainage boundary (unit hydraulic gradient boundary) or a constant head boundary (corresponding to the groundwater level); the lateral boundary of the simulation domain is set as a zero flux boundary to simulate independent vertical soil column elements.

[0088] To solve this partial differential equation on a computer, the hydraulic constitutive correction and fluid transport solution module 104 discretizes the three-dimensional spatial domain using the finite difference method or the finite volume method, and discretizes the time domain using a fully implicit backward difference scheme. For any node in the spatial grid system... In time step The discretized equations at time points are constructed in the form of mass conservation:

[0089] ;

[0090] in, and These represent the current time step. and the previous time step The water content of the node volume; Indicates the current time step; Indicates the volume of the control volume; This represents the summation of fluxes across all adjacent interfaces of the control volume; This represents the average corrected hydraulic conductivity at the interface. Indicates the area of ​​the control volume interface; and These represent the modified matrix potentials of the adjacent nodes and the current node at the current time step, respectively. and This indicates the corresponding gravitational potential height; Indicates the distance between the centers of the nodes; This represents the root water absorption rate at the current time step.

[0091] Given the equation and All are variables to be determined The highly nonlinear function (or pressure head) is solved linearly using the Picard iterative algorithm or the Newton-Raphson iterative algorithm in module 104, which is used for hydraulic constitutive correction and fluid transport solution. At each time step... Within the module 104, the hydraulic constitutive correction and fluid transport solution constructs a sparse coefficient matrix and performs iterative iterations. The convergence criterion for the iterations is set as follows: the norm of the head difference or the norm of the mass residual between two adjacent iterations is less than a preset tolerance error. :

[0092] ;

[0093] in, Indicates the iteration step index; Norm operations are used to represent matrices or vectors; Indicates the convergence tolerance; and These represent the current time step. Next, the Second and third The full-field corrected matrix potential vector obtained from the next iteration.

[0094] Furthermore, to balance computational efficiency and numerical stability, the hydraulic constitutive correction and fluid transport solution module 104 also integrates an adaptive time step control strategy. During the simulation, if the Picard iteration count exceeds a preset threshold (e.g., 10 iterations) and convergence is not achieved, or if the water content changes drastically within the current time step, the system will automatically reduce the time step. (For example, halving the time step); conversely, if the iteration converges rapidly, the system will attempt to increase the time step. This strategy effectively avoids numerical oscillations or computational divergence during the drastic nonlinear phases of rapid advance of the rainfall infiltration front or wetting reversal.

[0095] In the overall system-level coupling implementation, the hydraulic constitutive correction and fluid transport solution module 104 and the rhizosphere biochemical environment evolution module 103 are sequentially coupled using an operator splitting strategy. At each time step... In the initial stage, the system first calls the rhizosphere biochemical environment evolution module 103 to solve the solute transport equation and update the current root exudate concentration field. Subsequently, the effective viscosity of the fluid was calculated based on the updated concentration field. and dynamic contact angle Finally, using the updated and Corrected hydraulic parameters and The discretization and solution scheme ensures the stability of the numerical calculations and accurately captures the transient changes in soil hydraulic properties caused by ecophysiological processes.

Claims

1. A distributed eco-hydrological two-way coupled modeling system for sandy land, characterized in that, include: The initialization and discretization module is configured to construct a finite difference grid system in a three-dimensional physical space domain, and sets the initial state parameters of the system, including the effective dynamic viscosity of the fluid, the dynamic contact angle of the sand grain surface, the unsaturated hydraulic conductivity of the soil, and the matrix potential. The ecological growth calculation and source term generation module is configured to generate an environmental water stress signal based on the current water distribution state within the finite difference grid system, and convert the environmental water stress signal into a root exudate release source term; The rhizosphere biochemical environment evolution module is configured to receive the root exudate release source, simulate the exudate transport process, and reconstruct the effective dynamic viscosity of the fluid and the dynamic contact angle of the sand grain surface. The hydraulic constitutive correction and fluid transport solution module is configured to correct the unsaturated hydraulic conductivity of the soil using the effective dynamic viscosity of the fluid and the matrix potential using the dynamic contact angle. Based on the corrected parameters, the fluid transport control equation is solved, the water distribution state is updated, and the updated water distribution state is fed back to the ecological growth calculation and source term generation module, forming a dynamic closed-loop coupling between the sandy ecological process and the hydrological process.

2. The distributed eco-hydrological bidirectional coupled modeling system for sandy land according to claim 1, characterized in that, The initialization and discretization module is specifically configured as follows: The continuous three-dimensional physical space domain is divided into an orthogonal finite difference grid system, and each grid cell is assigned initial volumetric water content, initial matrix potential, saturated water content, residual water content, initial saturated hydraulic conductivity, initial effective fluid dynamic viscosity, and initial dynamic contact angle of sand grain surface.

3. The distributed eco-hydrological bidirectional coupled modeling system for sandy land according to claim 1, characterized in that, The specific configuration of the ecological growth calculation and source term generation module is as follows: Extract the current volumetric water content from the current water distribution state, calculate the effective saturation using the current volumetric water content, map the effective saturation to the environmental water stress signal, and calculate the root exudate release source term by combining the vegetation root length density distribution. The calculation model for the root exudate release source term is set to show a nonlinear exponential growth as the effective saturation representing the environmental water stress signal decreases.

4. The distributed eco-hydrological bidirectional coupled modeling system for sandy land according to claim 1, characterized in that, The rhizosphere biochemical environment evolution module is configured to simulate the secretion transport process as follows: The governing equations are constructed based on convection, diffusion, and reaction mechanisms. These equations take into account the convection carried by the water flow, the hydrodynamic diffusion driven by the concentration gradient, and the biochemical degradation. The biochemical degradation process is described by introducing an Arrhenius temperature correction function to reflect the real-time regulation of the degradation rate of secretions by ambient temperature.

5. The distributed eco-hydrological bidirectional coupled modeling system for sandy land according to claim 1, characterized in that, The process configuration for the rhizosphere biochemical environment evolution module to reconstruct the effective dynamic viscosity of the fluid is as follows: A concentration-dependent nonlinear rheological model was used to map the root exudate concentration obtained by simulating the exudate transport process to the effective dynamic viscosity of the fluid. The nonlinear rheological model defines the transformation of pore fluid from a Newtonian fluid to a non-Newtonian viscous fluid as the concentration of root exudates increases.

6. The distributed eco-hydrological bidirectional coupled modeling system for sandy land according to claim 2, characterized in that, The process of reconstructing the dynamic contact angle of the sand grain surface by the rhizosphere biochemical environment evolution module is configured as follows: A historical state variable is introduced to record the drying history of the grid cell. When the root exudate concentration obtained by simulating the exudate transport process is higher than the critical film-forming concentration threshold and the current volumetric water content in the current water distribution state is lower than the critical drying threshold, a biofilm solidification event is determined to have occurred. The biofilm solidification event triggers the dynamic contact angle to increase with the increase of the root exudate concentration until a wettability reversal occurs from a hydrophilic to a hydrophobic state.

7. The distributed eco-hydrological bidirectional coupled modeling system for sandy land according to claim 1, characterized in that, The process of correcting the matrix potential in the hydraulic constitutive correction and fluid transport solution module is configured as follows: The original matrix potential is scaled and adjusted by introducing a wettability correction factor that includes the dynamic contact angle cosine value; When the dynamic contact angle undergoes a wettability reversal from a hydrophilic to a hydrophobic state, the sign of the wettability correction factor changes, causing the corrected matrix potential to change from a negative attractive force to a positive repulsive force, thus forming an infiltration shielding effect.

8. The distributed eco-hydrological bidirectional coupled modeling system for sandy land according to claim 1, characterized in that, The process of correcting the unsaturated hydraulic conductivity in the hydraulic constitutive correction and fluid transport solution module is configured as follows: A rheological retardation factor is introduced, which is composed of the ratio of the dynamic viscosity of pure water to the effective dynamic viscosity of the fluid. The rheological retardation factor is used to correct the relative hydroconductivity that only considers the pore geometry. The corrected unsaturated hydroconductivity is set to be negatively correlated with the effective dynamic viscosity of the fluid, thus forming a rheological retardation effect.

9. The distributed eco-hydrological bidirectional coupled modeling system for sandy land according to claim 1, characterized in that, The process of solving the fluid transport control equations by the hydraulic constitutive correction and fluid transport solution module is configured as follows: The generalized Richards equation is constructed and the modified unsaturated hydraulic conductivity and the modified matrix potential are substituted into it. The spatial discretization is performed using the fully implicit finite volume method, and the discretized equation system is linearized and solved using an iterative algorithm to obtain the water content distribution at the next time step.

10. The distributed eco-hydrological bidirectional coupled modeling system for sandy land according to claim 9, characterized in that, The hydraulic constitutive correction and fluid transport solution module also integrates an adaptive time step control strategy, which is configured as follows: When the number of iterations of the iterative algorithm exceeds a preset threshold and fails to converge, the simulation time step is automatically reduced to adapt to the nonlinear calculation requirements of the wettability flipping stage.

Images