Gas-liquid interface distribution detection method for surface soil and program product

Abstract

The embodiment of the invention discloses a gas-liquid interface distribution detection method for surface soil and a program product, and the method comprises the steps: carrying out the geometric topological structure modeling of a three-dimensional structure image, obtaining a pore network topological structure, carrying out the quasi-static drainage simulation based on the pore network topological structure, and obtaining a gas-liquid interface distribution detection result. Calculating the meniscus area of the corner part and the meniscus area of the main end part, and obtaining the mapping relation between the moisture content and the interface area of the surface soil. According to the invention, quasi-static drainage simulation is carried out based on the pore network topological structure, the corner meniscus area and the main end meniscus area in the pore network topological structure can be calculated, the real gas-liquid interface distribution of the surface soil under the unsaturated condition can be accurately determined, and under the condition of knowing the moisture content of the surface soil, the real gas-liquid interface distribution of the surface soil can be accurately determined. The interface area under the unit volume can be obtained, inversion of an unobservable microstructure from an actually measured macroscopic state is achieved, and a solid structural foundation is provided for researching the related process of a gas-liquid interface. The method is widely applied to the technical field of soil detection.

Description

TECHNICAL FIELD

[0001] Embodiments of the present application relate to the technical field of soil detection, and particularly to a method and program product for detecting the gas-liquid interface distribution of surface soil. BACKGROUND

[0002] Surface soil is an important interface between the atmosphere and groundwater system. For surface soil, its pore structure is usually loose, the capillary force dominant effect is significant, and it is frequently disturbed by processes such as rainfall, evaporation, plant transpiration, etc., resulting in a highly dynamic unsaturated state. In the unsaturated state, water and gas phases exist in the pore structure of surface soil, and a complex gas-liquid contact interface is formed between the two phases. These gas-liquid interfaces not only control the water and gas migration behavior in the unsaturated zone, but also are important action sites for various environmental reactions. In recent years, it has been found that certain interfacially active pollutants, especially perfluoro and polyfluoro alkyl substances (PFAS), have a high tendency to adsorb on the gas-liquid interface, and significantly enrich on the gas-liquid interface in unsaturated soil. Compared with the water phase or the solid phase, the adsorption of such pollutants on the gas-liquid interface is more significant, resulting in that their migration ability and spatial distribution are mainly controlled by the distribution state of the gas-liquid interface.

[0003] Therefore, there is a need to detect the gas-liquid interface distribution of surface soil for the purpose of detecting soil pollutants, etc. At present, when studying the unsaturated state characteristics and water and gas distribution behavior of surface soil, experimental methods are mainly used to obtain information. Common experimental methods include soil sample collection and moisture content determination, evaporation or drainage test, and in-situ monitoring method based on pressure membrane instrument or tensiometer, which are used to evaluate the water content change and retention state at different depths or time conditions. These methods can be used to estimate parameters such as overall saturation of soil, pore water distribution and unsaturated hydraulic conductivity.

[0004] However, current experimental methods mainly focus on obtaining data at macroscopic or layered scales, and it is difficult to analyze the geometric characteristics of the gas-liquid interface inside the unsaturated soil at the microscopic scale. Especially in surface soil, due to its loose pore structure, frequent water disturbance, and dramatic dynamic evolution of water and gas interface, it is almost impossible to directly measure the position, shape and area of the gas-liquid interface under non-destructive conditions. Although the interface evolution can be indirectly inferred from the change in saturation, there is a lack of accurate expression matching the microscopic structure of soil. SUMMARY

[0005] In view of the technical problems in the current detection technology of the gas-liquid interface distribution of surface soil, such as the difficulty in accurately expressing the gas-liquid interface distribution information of surface soil matching the microscopic structure of soil, the purpose of embodiments of the present application is to provide a method and program product for detecting the gas-liquid interface distribution of surface soil.

[0006] In one aspect, embodiments of the present application include a method for detecting air-liquid interface distribution of surface soil, the method for detecting air-liquid interface distribution of surface soil comprising:

[0007] obtaining a three-dimensional structure image of the surface soil;

[0008] modeling a geometric topology of the three-dimensional structure image to obtain a pore network topology, the pore network topology comprising a plurality of pore bodies and a plurality of pore throats;

[0009] performing a quasi-static drainage simulation based on the pore network topology;

[0010] calculating a corner meniscus area and a main end meniscus area according to a result of the quasi-static drainage simulation;

[0011] obtaining a water content-interface area mapping relationship of the surface soil according to the corner meniscus area and the main end meniscus area.

[0012] Further, the obtaining a three-dimensional structure image of the surface soil comprises:

[0013] obtaining a columnar surface soil sample;

[0014] performing X-ray computed tomography on the columnar surface soil sample;

[0015] preprocessing a result of the X-ray computed tomography to obtain the three-dimensional structure image.

[0016] Further, the preprocessing a result of the X-ray computed tomography to obtain the three-dimensional structure image comprises:

[0017] removing edge region images and non-target region images from the result of the X-ray computed tomography to obtain a middle region image;

[0018] performing Gaussian filtering on the middle region image;

[0019] performing binarization on a result of the Gaussian filtering to obtain the three-dimensional structure image.

[0020] Further, the modeling a geometric topology of the three-dimensional structure image to obtain a pore network topology comprises:

[0021] performing Euclidean transformation and watershed segmentation on the three-dimensional structure image to determine a boundary range of each of the pore bodies and to identify the pore throats between each of the pore bodies and other pore bodies;

[0022] calculating geometric parameters of each of the pore bodies and each of the pore throats;

[0023] record the topological connection relationship between each of the pore bodies and each of the pore throats.

[0024] Further, the quasi-static drainage simulation based on the pore network topology structure comprises:

[0025] taking the capillary pressure as a driving variable;

[0026] setting a condition that the external capillary pressure gradually increases, simulating the gas-water two-phase distribution when the pore network topology structure forms a non-saturated state under the trap mechanism.

[0027] Further, the calculation of the corner meniscus area comprises:

[0028] For any pore body, when the saturation of the pore body is greater than or equal to the drainage critical water saturation, the corner meniscus area A is calculated according to the formula

[0029]

[0030] obtaining the cross-sectional area A occupied by the liquid phase in the pore body w , wherein r wd is the liquid equivalent radius of the pore body at the initial time of the quasi-static drainage simulation, and areCoe is the area coefficient of the pore cross section of the pore body;

[0031] according to the formula

[0032] A nw = A cross -A w

[0033] obtaining the cross-sectional area occupied by the gas phase in the pore body, wherein A cross is the cross-sectional area of the pore body;

[0034] based on the volume conservation, according to the formula

[0035]

[0036] calculating the liquid column length l w , wherein l is the total length of the pore body, s w is the saturation of the pore body,

[0037] according to the formula

[0038]

[0039] calculating the corner perimeter P per unit length nw , wherein β i is the corner half angle of the pore body, and θ is the contact angle,

[0040] according to the formula

[0041] A AMs = P nw · l nw + A cross

[0042] calculating the corner meniscus area A AMs of the pore body.

[0043] Further, the corner meniscus area is calculated by:

[0044] For any pore body, when the saturation of the pore body is less than the drainage critical water saturation, the effective radius r

[0045]

[0046] obtaining the water film radius r w in the pore body, wherein A w is the cross-sectional area occupied by the liquid phase in the pore body, A cross is the cross-sectional area of the pore body, s w is the saturation of the pore body, and areCoe is the area coefficient of the pore cross-section of the pore body;

[0047] calculating the corner meniscus area A of the pore body according to the formula:

[0048]

[0049] AMs ; wherein l is the total length of the pore body, β i is the corner half-angle of the pore body, and θ is the contact angle.

[0050] Further, the main end meniscus area is calculated by:

[0051] For any pore body-pore throat connection pair, when the pore body and the pore throat in the pore body-pore throat connection pair are in different phases, the effective radius r

[0052]

[0053] of the pore body-pore throat connection pair is calculated according to the formula: eff cross ; wherein A c is the cross-sectional area of the pore body, and P is the cross-sectional perimeter of the pore body.

[0054] calculating the capillary pressure at both ends of the pore body in the pore body-pore throat connection pair and the capillary pressure at both ends of the pore throat

[0055] according to the formula:

[0056] .where σ is the surface tension, θ is the contact angle, and a is the half angle;

[0057] According to the formula

[0058]

[0059] The pressure difference ΔP between the interfaces is calculated c ;

[0060] According to the formula

[0061]

[0062] The radius of curvature r of the main end meniscus type interface is calculated MTM ;

[0063] According to the formula

[0064]

[0065] A flag = min(A body , A throat ) · f overlap

[0066] The connection surface area A of the pore body and the pore throat in the pore body-pore throat connection pair is calculated flag ; wherein f overlap is a geometric overlap factor estimated according to the angular opening angle, A body is the cross-sectional area of the pore body, A throat is the cross-sectional area of the pore throat, a body is the average value of the half angle of the pore body, and a throat is the average value of the half angle of the pore throat;

[0067] According to the formula

[0068]

[0069] The curvature correction coefficient ρ is calculated, wherein is the hydraulic radius of the pore body in the pore body-pore throat connection pair, is the hydraulic radius of the pore throat in the pore body-pore throat connection pair;

[0070] According to the formula

[0071] A MTM = A flag · ρ

[0072] The main end meniscus area A of the pore body-pore throat connection pair is calculated MTM .

[0073] Further, the obtaining the water content-depth mapping relationship of the surface soil according to the corner meniscus area and the main end meniscus area comprises:

[0074] The total effective volume is obtained by accumulating the volumes of all the pore bodies and all the pore throats in the pore network topology structure.

[0075] The total corner meniscus area is obtained by accumulating all the corner meniscus areas corresponding to the pore network topology structure.

[0076] The total main end meniscus area is obtained by accumulating all the main end meniscus areas corresponding to the pore network topology structure; and the total main end meniscus area is obtained according to the formula

[0077]

[0078] The total interface area a is obtained by calculation according to the formula Wherein, The total corner meniscus area is The total main end meniscus area is

[0079] The total interface area a is obtained according to the formula

[0080]

[0081] The interface area a per unit volume of the surface soil is obtained by calculation according to the formula nw ; wherein V total The total effective volume is The porosity of the surface soil is

[0082] The depth-water content curve of the surface soil is obtained.

[0083] The water content of the surface soil is determined according to the depth-water content curve.

[0084] The mapping relationship is established according to the water content and the interface area per unit volume.

[0085] In another aspect, the embodiments of the present application also include a computer program product comprising a computer program, which, when executed by a processor, implements the method for detecting the gas-liquid interface distribution of the surface soil in the embodiments.

[0086] The embodiment of the present application has the beneficial effect that the method for detecting the gas-liquid interface distribution of the surface soil in the embodiment can calculate the corner meniscus area and the main end meniscus area in the pore network topology structure by carrying out quasi-static drainage simulation based on the pore network topology structure, so as to accurately determine the real gas-liquid interface distribution of the surface soil under the unsaturated condition, which is significantly better than the traditional method based on the empirical formula, finally obtains the moisture content-interface area mapping relationship of the surface soil, so that the interface area per unit volume of the surface soil can be obtained in the case that the moisture content of the surface soil is known, thereby realizing the inversion of the unobservable microstructure from the measured macrostate, and providing a solid structural foundation for the research on the gas-liquid interface related process. BRIEF DESCRIPTION OF DRAWINGS

[0087] Figure 1 The step schematic diagram of the method for detecting the gas-liquid interface distribution of the surface soil in the embodiment;

[0088] Figure 2 The schematic diagram of the corner meniscus in the embodiment;

[0089] Figure 3 The schematic diagram of the main end meniscus in the embodiment;

[0090] Figure 4 The schematic diagram of the depth-moisture content curve in the embodiment. DETAILED DESCRIPTION

[0091] Term explanation:

[0092] 1. Gas-liquid interface: the interface region formed by the contact of the gas phase and the liquid phase in the pore structure in the porous medium, which is divided into two types of corner meniscus and main end meniscus;

[0093] 2. Corner meniscus: refers to the local gas-liquid interface with high curvature formed at the corner region of the pore body or pore throat wall surface and controlled by capillary force;

[0094] 3. Main end meniscus: refers to the main interface formed by the liquid at the end of the pore throat or pore body, which is the main meniscus dividing the gas-liquid two phases and controls the entry and exit process of the liquid;

[0095] 4. Trapping mechanism: refers to the phenomenon that part of the wetting phase cannot be drained because it is surrounded by the non-wetting phase during the drainage process, so that the wetting phase is trapped.

[0096] For the topsoil, its pore structure is usually loose, and the capillary force dominates the effect significantly. The processes of rainfall, evaporation, and plant transpiration frequently disturb the topsoil, leading to its highly dynamic unsaturated state. In this state, water does not completely disappear during the drainage process, but remains in the form of residual water in the corners, slits, and surfaces of the pore structure, and exists in the form of thin film (main end meniscus) or corner water film (corner meniscus). These discontinuous water bodies and the gas phase form a large number of widely distributed gas-liquid interfaces. Therefore, if the residual water state (i.e. residual saturation) can be used as the input variable, combined with the three-dimensional pore structure information of the soil, a set of physically mechanism-driven model can be established to accurately deduce the geometric shape and area of the gas-liquid interface under a specific state, and then the enrichment position and content of the pollutants on the interface can be further calculated, so as to realize the quantitative prediction of the spatial distribution and environmental risk of the pollutants in the topsoil.

[0097] Although experimental methods can measure the overall water content of the soil, it is difficult to distinguish the spatial distribution and geometric characteristics of the gas-liquid interface at the microscale. Especially in the unsaturated state, the soil water is distributed in the form of residual water film or corner water body in the complex pore, and the area of the gas-liquid interface formed by it and its connection state are currently difficult to be directly quantified by experimental methods.

[0098] With the development of high-resolution imaging technology, X-ray computed tomography has become an important tool for studying soil pore structure, which can realize three-dimensional reconstruction and obtain the geometric shape, size distribution, and connectivity of real pores. If the method based on the pore network model is introduced, the drying and wetting processes under the capillary pressure can be reproduced in the virtual environment, and the distribution and area of the gas-liquid interface can be deduced by constructing the position of the main end meniscus (MTM) and the corner meniscus (AMs).

[0099] If only the identification of the pore scale interface is focused on, and the gas-liquid interface formed in the pore throat region is not considered, the obtained results will have problems such as insufficient expression of the gas-liquid interface. For example, in the pore network model of the topsoil, although the volume of a single pore throat is small, its number is much larger than that of the pore body. Ignoring the pore throat part will ignore its interface contribution, resulting in a significant underestimate of the overall gas-liquid interface area.

[0100] Moreover, if the trapping mechanism in the unsaturated drainage process is not fully embodied, but the water phase is assumed to be freely drained under the driving of capillary pressure, since this drainage mode is only applicable to deep layers or regular media with high pore connectivity, and the surface soil has complex pore structure and irregular invasion path of the gas phase, the water phase is easily surrounded by the non-wetting phase to form a trapping area, so this drainage mode is not applicable, which may lead to underestimation of the interface distribution under the residual saturation state. Therefore, when simulating the gas-liquid distribution of the surface soil, the trapping mechanism cannot be ignored, and if the ideal drainage assumption is still used, the gas-liquid interface area will be systematically overestimated, especially at low saturation.

[0101] To solve the above problems, a technical path for reconstructing the pore structure of the surface soil based on the two-dimensional CT image and combining the quasi-static drainage simulation to calculate the gas-liquid interface can be used. In the simulation, the trapping mechanism is introduced to accurately describe the behavior that the water body is trapped by the gas phase and cannot be drained during the drainage process, and the gas-liquid interface contribution inside the pore throat is explicitly calculated, thereby improving the accuracy of the interface area estimation. At the same time, the correspondence between the water content of the surface soil and the residual water content in the quasi-static drainage simulation is established. Since the water content of the surface soil is close to the drainage final state relative to the saturation zone, it can be regarded as the residual water state in the simulation. Based on this consistency, by inputting the measured water content of the surface soil under the actual working condition, the gas-liquid interface distribution and area under the drainage final state can be deduced, and the quantitative mapping from the macroscopic water content state to the microscopic interface structure is realized.

[0102] Based on the above principle, an embodiment of the present application provides a gas-liquid interface distribution detection method for surface soil. Referring to Figure 1 , the gas-liquid interface distribution detection method for surface soil includes the following steps:

[0103] S1. Obtain a three-dimensional structure image of the surface soil;

[0104] S2. Model the geometric topological structure of the three-dimensional structure image to obtain a pore network topological structure;

[0105] S3. Perform quasi-static drainage simulation based on the pore network topological structure;

[0106] S4. Calculate the corner meniscus area and the main end meniscus area according to the results of the quasi-static drainage simulation;

[0107] S5. Obtain a water content-interface area mapping relationship of the surface soil according to the corner meniscus area and the main end meniscus area.

[0108] In this embodiment, when performing step S1, that is, obtaining the three-dimensional structure image of the surface soil, the following steps can be performed:

[0109] S101. Obtain a columnar topsoil sample;

[0110] S102. Perform X-ray computed tomography on the columnar topsoil sample;

[0111] S103. Preprocess the results of the X-ray computed tomography to obtain a three-dimensional structure image.

[0112] In steps S101-S103, the three-dimensional structure image of the topsoil is obtained using X-ray computed tomography (CT) technology. Specifically, step S101 is performed to preprocess the columnar topsoil sample to ensure its geometric structure is complete and stable during the scanning process. Then, step S102 is performed to scan the columnar topsoil sample layer by layer by a CT device to obtain grayscale image data with high spatial resolution. The key parameters (such as voltage, current, resolution, and exposure time) during the scanning process are optimized according to the soil particle size and scanning accuracy requirements to ensure that the image quality is sufficient to distinguish the spatial distribution of pores and solids.

[0113] The original image data obtained in step S102 usually contains noise, boundary artifacts, and uneven grayscale distribution, which can be further processed.

[0114] When step S103 is performed, first, the original image obtained in step S102 is processed to remove edge effects and non-target regions, and only the central region image with complete structure and representative is retained; then a Gaussian filter algorithm is applied to the image to smooth the image and suppress noise, while enhancing the contrast of soil particle edges; finally, the Gaussian filtered image is binarized by a grayscale threshold segmentation method, and each pixel (voxel) in the image is classified as solid or pore, effectively separating the pore space and soil skeleton, and finally obtaining the three-dimensional structure image of the columnar topsoil sample.

[0115] In this embodiment, when step S2 is performed, that is, the three-dimensional structure image is geometrically modeled to obtain the pore network topology, the following steps can be performed:

[0116] S201. Perform Euclidean transformation and watershed segmentation on the three-dimensional structure image to determine the boundary range of each pore body and identify the pore throat between each pore body and other pore bodies;

[0117] S202. Calculate the geometric parameters of each pore body and each pore throat;

[0118] S203. Record the topological connection relationship between each pore body and each pore throat.

[0119] In step S201, on the basis of the binarized three-dimensional structure image of the surface soil, extraction of the pore network and geometric topological structure modeling are carried out. Through image analysis and network abstraction methods, the complex three-dimensional pore space is simplified into a discrete network structure composed of pore bodies and pore throats. Using network extraction tools, the spatial morphology and connectivity of the pores are identified and quantified. Specifically, the three-dimensional image is subjected to Euclidean distance transformation, and image processing algorithms such as watershed segmentation are used to determine the boundary range of each pore body and identify the connection path between the pore body and other pore bodies. The topological structure of the pore network constructed in this way retains the main geometric characteristics of the original pore structure of the surface soil and greatly simplifies the calculation complexity, making it suitable for physical process simulation such as multiphase flow and gas-liquid interface evolution.

[0120] After the extraction of the pore network topological structure in step S201 is completed, steps S202 and S203 are performed to further calculate the geometric parameters of each pore body and pore throat, including cross-sectional area, perimeter, volume, length, equivalent radius, triangular cross-sectional angle half-angle, etc., while recording the topological connection relationship between the pore body and the pore throat. These parameters provide the necessary structural basis for subsequent capillary-driven processes, meniscus calculation, and gas-liquid interface modeling.

[0121] In this embodiment, when step S3, that is, the step of performing quasi-static drainage simulation based on the pore network topological structure, is performed, the following steps can be performed:

[0122] S301. Capillary pressure is used as the driving variable;

[0123] S302. The condition of gradually increasing external capillary pressure is set, and the gas-liquid two-phase distribution when the pore network topological structure forms a non-saturated state under the trapping mechanism is simulated.

[0124] In steps S301-S302, on the basis of the pore network topological structure extracted in step S2, quasi-static drainage simulation of the gradual invasion of the gas phase into the pore space is carried out. This simulation process uses capillary pressure as the driving variable to simulate the gradual replacement of the water phase in the original saturated state by the gas phase under the condition of gradually increasing external capillary pressure, thereby forming a non-saturated state of coexistence of gas and water phases. The simulation assumes that the flow process is quasi-static, that is, each step of gas phase invasion is carried out under the premise that the system reaches local equilibrium, without considering transient inertia or dynamic effects. As the capillary pressure gradually increases, more and more pore units are sequentially invaded and occupied by the gas phase, and the simulation results can reflect the preferential invasion path and distribution process of the gas phase in the pore network.

[0125] Step S302 introduces the trapping mechanism to more realistically reflect the distribution characteristics of gas and water phases under unsaturated conditions. The trapping mechanism refers to: if a certain connected water phase cluster meets the following characteristics: all pore bodies and pore throats inside it are in a state of not being invaded by the gas phase, and are completely isolated from the outside entrance and exit, and the edge of the water phase cluster is no longer connected to other drainable areas, it is considered that the water phase cluster is surrounded by the gas phase and cannot continue to drain, i.e. the trapping phenomenon occurs.

[0126] Step S302 introduces the trapping mechanism to more realistically reflect the distribution characteristics of gas and water phases under unsaturated conditions. The trapping mechanism refers to: if a certain connected water phase cluster meets the following characteristics: all pore bodies and pore throats inside it are in a state of not being invaded by the gas phase, and are completely isolated from the outside entrance and exit, and the edge of the water phase cluster is no longer connected to other drainable areas, it is considered that the water phase cluster is surrounded by the gas phase and cannot continue to drain, i.e. the trapping phenomenon occurs.

[0127] By performing steps S301-S302, the quasi-static drainage simulation is performed based on the pore network topology, which not only reveals the basic law of the distribution of gas and liquid in unsaturated soil with the change of capillary pressure, but also provides physical basis and state data for subsequent interface type identification (such as corner wetting or main end gas-liquid interface) and interface area calculation.

[0128] In step S4, the steps of calculating the corner meniscus area (AMs) and calculating the main end meniscus area (MTM) are included.

[0129] 1. Step of calculating the corner meniscus area (AMs)

[0130] The corner meniscus area (AMs) of each pore body and pore throat in the pore network topology can be calculated according to the water saturation and geometric characteristics, i.e. the gas-liquid interface area retained in the solid corner area. This step physically simulates the gas-liquid interface formed by the liquid phase retained on the surface of the pore body and pore throat under unsaturated conditions, and the area size is controlled by factors such as local saturation, geometric parameters, material properties (such as contact angle and surface tension) and the like.

[0131] In this embodiment, the corner meniscus area (AMs) of each pore body and pore throat in the pore network topology can be calculated based on the theoretical model of the triangular pore structure as shown in FIG. 8. Figure 2 The theoretical model considers the following core factors: current water saturation, drainage critical water saturation, residual water column radius, pore volume and cross-sectional area, corner angle and contact angle.

[0132] For any one pore body in the pore network topology, when calculating the corner meniscus area (AMs) of the pore body, one of the two calculation logics can be selected according to the relationship between the current water saturation of the pore body and the drainage critical saturation.

[0133] (1) Water column stable existing region (saturation ≥ drainage critical water saturation)

[0134] In this case, it is considered that a stable meniscus water column exists in the pore of the pore body, and the radius of the water column is the known residual water column radius. At this time, when the step of calculating the corner meniscus area in step S4 is performed, the following steps S401A-S402A can be specifically performed:

[0135] S401A. For any one pore body, when the saturation of the pore body is greater than or equal to the drainage critical water saturation, the corner meniscus area AMs is calculated according to the formula

[0136]

[0137] Obtain the cross-sectional area A occupied by the liquid phase in the pore body w , wherein r wd is the liquid equivalent radius of the pore body at the initial time of the quasi-static drainage simulation (i.e., the critical radius at the time of entering the corner drainage stage), and areCoe is the area coefficient of the pore cross section of the pore body;

[0138] S402A. According to the formula

[0139] A nw = A cross -A w

[0140] Obtain the cross-sectional area occupied by the gas phase in the pore body, wherein A cross is the cross-sectional area of the pore body;

[0141] Based on the volume conservation, according to the formula

[0142]

[0143] Calculate the liquid column length l w , wherein l is the total length of the pore body, s w is the saturation of the pore body,

[0144] At this time, the gas-liquid interface is mainly formed at the junction of the liquid column and the gas triangular corner, and according to the formula

[0145]

[0146] Calculate the corner perimeter P per unit length nw , wherein β i is the corner half angle of the pore body, and θ is the contact angle,

[0147] Since the corner meniscus area (AMs) includes the area of the sidewall interface and the area of the water column cross section, according to the formula

[0148] A AMs = P nw · l nw + A cross

[0149] The corner meniscus area A AMs of the pore body is calculated.

[0150] (2) The water column is broken but the corner film remains (water saturation < drainage critical water saturation)

[0151] When the water saturation is lower than the critical value, the water column in the pore of the pore body can no longer exist stably, but the corner water film can still remain. In this case, the new film radius can be estimated according to the current water amount. At this time, when the step of calculating the corner meniscus area in step S4 is performed, the following steps S401B-S402B can be performed in particular:

[0152] S401B. For any pore body, when the saturation of the pore body is less than the drainage critical water saturation, the corner meniscus area A

[0153]

[0154] The film radius r w in the pore body is obtained, where A w is the cross-sectional area occupied by the liquid phase in the pore body, A cross is the cross-sectional area of the pore body, s w is the saturation of the pore body, and areCoe is the area coefficient of the pore cross section of the pore body;

[0155] At this time, the entire pore is a gas phase region, so the interface area is the perimeter of the corner film multiplied by the total length, and step S402B is performed;

[0156] S402B. According to the formula

[0157]

[0158] The corner meniscus area A AMs of the pore body is calculated; where l is the total length of the pore body, β i is the corner half angle of the i-th corner meniscus in the pore body, and θ is the contact angle.

[0159] By performing steps S401A-S402A or steps S401B-S402B, the corner meniscus area (AMs) of the pore body in different cases can be calculated respectively.

[0160] 2. The step of calculating the main end meniscus area (MTM)

[0161] When one pore body and one adjacent pore throat in the pore network topology are in non-equivalent wetting states and the corner wetting phase is not connected, an interface will form a triangular or curved gas-liquid interface at the connection between the two, as shown in Figure 3 The geometry and area of the interface are determined by the curvature difference on both sides, the geometric characteristics of the connection surface, and the material properties.

[0162] Based on the model shown in Figure 3 When performing the step of calculating the main end meniscus area in step S4, the following steps can be performed:

[0163] S403. For any pore body-pore throat connection pair, when the pore body and the pore throat in the pore body-pore throat connection pair are in different phases, the effective radius r of the pore body-pore throat connection pair is calculated according to the formula

[0164]

[0165] The effective radius r of the pore body-pore throat connection pair is calculated according to the formula eff ; Wherein A cross is the cross-sectional area of the pore body, and P is the cross-sectional perimeter of the pore body;

[0166] Wherein, the pore body-pore throat connection pair includes a pore body and a pore throat connected to the pore body; the pore body and the pore throat are in different phases, including the case that the pore body is occupied by the wetting phase and the pore throat is occupied by the non-wetting phase, or the pore body is occupied by the non-wetting phase and the pore throat is occupied by the wetting phase, and when the pore body and the pore throat are in different phases, a gas-liquid interface of the MTM type will be formed between the pore body and the pore throat;

[0167] S404. According to the Young-Laplace equation, the capillary pressure at both ends of the pore body in the pore body-pore throat connection pair is calculated according to the formula

[0168]

[0169] The capillary pressure at both ends of the pore throat in the pore body-pore throat connection pair is calculated according to the formula Where σ is the surface tension, θ is the contact angle, and α is the half angle;

[0170] S405. The pressure difference ΔP between the interfaces is calculated according to the formula

[0171]

[0172] The pressure difference ΔP between the interfaces is calculated according to the formula c ;

[0173] S406. The pressure difference ΔP between the interfaces is calculated according to the formula

[0174] ​

[0175] The curvature radius r of the main end meniscus (MTM) type interface is calculated MTM ;

[0176] S407. The connection surface area A of the pore body and the pore throat in the pore body-pore throat connection pair is calculated according to the formula

[0177]

[0178] A flag = min(A body , A throat )·f overlap

[0179] The connection surface area A of the pore body and the pore throat in the pore body-pore throat connection pair is calculated flag ; wherein f overlap is a geometric overlap factor estimated according to the angular opening angle, A body is the cross-sectional area of the pore body, A throat is the cross-sectional area of the pore throat, a body is the average of the pore body angular half-angle, a throat is the average of the pore throat angular half-angle; the connection surface area A flag calculated in step S407 represents Figure 3 the triangular interface area of the blue part in FIG. 8, that is, assuming that the pore body-pore throat connection surface is a triangular butt joint;

[0180] S408. Considering the possibility of area amplification of the curved surface interface compared with the flat surface interface, a curvature correction coefficient is introduced according to the formula

[0181]

[0182] The curvature correction coefficient p is calculated, wherein R is the hydraulic radius of the pore body in the pore body-pore throat connection pair, and R is the hydraulic radius of the pore throat in the pore body-pore throat connection pair.

[0183] S409. The main end meniscus area A MTM of the pore body-pore throat connection pair is calculated according to the formula

[0184] A MTM = A flag · p

[0185] The main end meniscus area A MTM of the pore body-pore throat connection pair is calculated according to the formula

[0186] In this embodiment, the calculation steps of the corner meniscus area (AMs) in steps S401A-S402A and steps S401B-S402B can be performed for all pore bodies in the pore network topology, so as to obtain all the corner meniscus areas; the calculation steps of the main end meniscus area (MTM) in steps S403-S409 can be performed for all pore body-pore throat connection pairs in the pore network topology, so as to obtain all the main end meniscus areas.

[0187] In this embodiment, when performing step S5, that is, obtaining the water content-interface area mapping relationship of the surface soil according to the corner meniscus area and the main end meniscus area, the following steps can be performed:

[0188] S501. Accumulate the volumes of all pore bodies and all pore throats in the pore network topology to obtain a total effective volume;

[0189] S502. Accumulate all the corner meniscus areas corresponding to the pore network topology to obtain a total corner meniscus area;

[0190] S503. Accumulate all the main end meniscus areas corresponding to the pore network topology to obtain a total main end meniscus area;

[0191] S504. Calculate according to the formula

[0192]

[0193] to obtain a total interface area wherein, is the total corner meniscus area, is the total main end meniscus area;

[0194] S505. Calculate according to the formula

[0195]

[0196] to obtain the interface area a per unit volume of the surface soil nw ; wherein V total is the total effective volume, is the porosity of the surface soil;

[0197] S506. Obtain a depth-water content curve of the surface soil;

[0198] S507. Determine the water content of the surface soil according to the depth-water content curve.

[0199] S508. Establish a mapping relationship according to the water content and the interface area per unit volume.

[0200] In step S501, the total effective volume V is obtained by accumulating the volumes of all pore bodies and all pore throats in the pore network topology. total Total effective volume V total It is the sum of the volumes of non-boundary pore bodies and pore throats, which is equivalent to the volume after excluding the volumes of the inlet and outlet pores from the total volume of the pore network topology structure, that is, only the volume of the intermediate effective unit is retained, which can be used for normalization in subsequent steps.

[0201] In step S502, the total corner meniscus area is obtained by accumulating the areas of all corner meniscus corresponding to the pore network topology.

[0202] In step S503, the total main end meniscus area is obtained by accumulating the areas of all main end meniscus corresponding to the pore network topology.

[0203] In step S504, the total corner meniscus area is calculated Total main tip meniscus area The total interface area is obtained by summing

[0204] In step S505, macro-normalization is performed according to the formula

[0205]

[0206] Calculate and obtain the interfacial area a per unit volume of surface soil nw ;in, is the porosity of the surface soil;

[0207] In step S506, you can get Figure 4 The depth-moisture content curve of the surface soil is shown. The depth-moisture content curve represents the corresponding relationship between soil depth and soil moisture content under normal conditions (e.g., in a natural environment). Therefore, in step S507, the location at a specific depth in the surface soil is determined, and the moisture content of the soil at that location can be determined based on the depth-moisture content curve.

[0208] In step S508, the moisture content determined in step S507 is compared with the interface area per unit volume of the surface soil obtained by executing steps S501-S505. nw By using the mapping relationship between moisture content and interfacial area per unit volume, the interfacial area per unit volume of the surface soil can be determined when the moisture content of the surface soil is known, thereby inferring the gas-liquid interface distribution of the surface soil from the moisture state of the surface soil.

[0209] The principle of performing steps S501-S508 in this embodiment is that, in the drainage simulation, the water phase driven out caused by the gradual increase of capillary pressure has physical consistency with the trend that the water content of the surface unsaturated soil gradually decreases with the increase of height under actual natural conditions (as shown in FIG. 6). Figure 4 Especially under the joint action of natural gravity and capillary force, the water content state of the surface soil is often in the terminal stage of the drainage simulation process, that is, close to the stable state corresponding to the "residual water" in the simulation. Therefore, under the premise that the water content of the surface soil is known, the corresponding gas-liquid interface area per unit volume can be obtained in the water content-unit volume interface area mapping curve obtained by the quasi-static drainage simulation. In this way, knowing the water content of the surface soil, the unobservable microstructure (gas-liquid interface area) can be inversely calculated from the interface area per unit volume obtained by performing steps S501-S508, thereby realizing the inversion of the unobservable microstructure (gas-liquid interface area) from the measured macrostate (water content), and providing a solid structural foundation for the study of gas-liquid interface related processes (such as the enrichment and migration of pollutants on the interface).

[0210] In summary, the gas-liquid interface distribution detection method for surface soil in this embodiment has the following effects:

[0211] 1. By extracting the pore network topology structure containing pore bodies and pore throats, and accurately constructing the geometric model of the corner meniscus (AM) and the main end meniscus (MTM), the real gas-liquid interface distribution under unsaturated conditions can be completely tracked, which is significantly better than the traditional method based on empirical formula or assumed surface;

[0212] 2. The trapping mechanism is introduced in the drainage simulation process, which more truly realizes the existence state of "closed water" in unsaturated seepage, and is suitable for surface soil with complex pore structure and irregular gas phase invasion path. This method realizes effective judgment of the trapped water body by identifying the closed wetting phase area not connected with the inlet and outlet, thereby improving the physical reality of the gas-liquid distribution state in the drainage simulation, improving the accuracy of the interface area prediction, and significantly improving the accuracy of the interface area estimation under low saturation;

[0213] 3. Based on the simulation of the whole quasi-static drainage process, a quantitative mapping relationship between water content and interface area per unit volume is established, which realizes the function of inversely calculating the microstructure of the unobservable interface from the macroscopically measurable water content, and has important practical applicability;

[0214] 4. By combining the pore network topology structure and the measured water content of the surface soil, the gas-liquid interface distribution can be quantitatively inverted on the microscale, filling the gap that the internal gas-liquid interface of the surface soil cannot be directly observed in experiments.

[0215] The method for detecting the gas-liquid interface distribution of the surface soil in the embodiment can be implemented by writing a computer program for executing the method, writing the computer program into a computer device or a computer program product such as a storage medium, and executing the method when the computer program is read out and run, thereby achieving the same technical effects as the method for detecting the gas-liquid interface distribution of the surface soil in the embodiment.

[0216] It should be noted that, unless otherwise specified, when a certain feature is referred to as being "fixed", "connected" to another feature, it can be directly fixed, connected to the other feature, or indirectly fixed, connected to the other feature. In addition, the up, down, left, right and the like used in the present disclosure are only relative to the relative positional relationship of the components of the present disclosure in the drawings. The singular forms "a", "an" and "the" used in the present disclosure are also intended to include the plural forms, unless the context clearly indicates otherwise. In addition, unless otherwise defined, all technical and scientific terms used in the embodiments have the same meanings as generally understood by those skilled in the art. The terms used in the embodiments of the present disclosure are only for the purpose of describing the specific embodiments, and are not intended to limit the embodiments of the present disclosure. The term "and / or" used in the embodiments includes any combination of one or more related listed items.

[0217] It should be understood that although the terms first, second, third, etc. can be used in the present disclosure to describe various elements, these elements should not be limited to these terms. These terms are only used to distinguish one type of element from another. For example, without departing from the scope of the present disclosure, a first element can also be referred to as a second element, and similarly, a second element can also be referred to as a first element. The use of any and all examples or exemplary language (e.g., "for example", "as such", etc.) provided in the embodiments of the present disclosure is only intended to better illustrate the embodiments of the present disclosure, and unless otherwise required, does not impose any limitation on the scope of the embodiments of the present disclosure.

[0218] It should be recognized that the embodiments of the present disclosure can be implemented or embodied by computer hardware, a combination of hardware and software, or by computer instructions stored in a non-transitory computer readable storage medium. The method can be implemented in a computer program configured with a non-transitory computer readable storage medium, in which the storage medium thus configured causes the computer to operate in a specific and predefined manner according to the methods described in the specific embodiments and the drawings. Each program can be implemented in a high-level procedural or object-oriented programming language to communicate with a computer system. However, if necessary, the program can be implemented in assembly or machine language. In any case, the language can be a compiled or interpreted language. In addition, the program can be run on a programmed application-specific integrated circuit for this purpose.

[0219] Further, the operations of the processes described in this embodiment can be performed in any suitable order unless otherwise indicated herein or otherwise clearly contradicted by context. The processes described in this embodiment (or variations and / or combinations thereof) can be implemented under the control of one or more computer systems configured with executable instructions (e.g., computer programs, one or more computer programs, or one or more applications), by hardware, or combinations thereof. The computer programs include a plurality of instructions that are executable by one or more processors.

[0220] Further, the methods can be implemented in any suitable type of computing platform operatively coupled to any suitable type of computing platform, including but not limited to a personal computer, a mini-computer, a mainframe, a workstation, a network or distributed computing environment, a stand-alone or integrated computer platform, or in communication with a charged particle tool or other imaging device, and the like. Aspects of embodiments of the present invention can be implemented in machine readable code stored on a non-transitory storage medium or device, whether removable or integrated into a computing platform, such as a hard disk, an optical read and / or write storage medium, RAM, ROM, and the like, such that it can be read by a programmable computer to configure and operate the computer to perform the processes described herein when the storage medium or device is read by the computer. Further, the machine readable code, or portions thereof, can be transmitted over wired or wireless networks. The present invention encompasses these and other different types of non-transitory computer readable storage media when the instructions or programs incorporating the above steps are implemented in conjunction with a microprocessor or other data processor. Embodiments of the present invention also encompass the computer itself when programmed in accordance with the methods and techniques of embodiments of the present invention.

[0221] The computer programs are capable of applying to input data to perform the functions of embodiments of the present invention, thereby transforming the input data to generate output data that is stored to non-volatile memory. The output information can also be applied to one or more output devices, such as a display. In preferred embodiments of embodiments of the present invention, the transformed data represents a physical and tangible object, including a particular visual depiction of a physical and tangible object produced on a display.

[0222] The above merely preferred embodiments of the present invention, and the present invention is not limited to only the above-described embodiments, as long as the same means achieve the technical effects of embodiments of the present invention, any modification, equivalent replacement, improvement, etc. made within the spirit and principle of embodiments of the present invention shall be included in the scope of protection of embodiments of the present invention. The technical solutions and / or embodiments of the present invention can have various modifications and changes within the scope of protection.

Claims

1. A method for detecting the gas-liquid interface distribution of surface soil, characterized in that: The gas-liquid interface distribution detection method for surface soil includes: Obtain three-dimensional structural images of the surface soil; Performing geometric topological modeling on the three-dimensional structural image to obtain a pore network topological structure; the pore network topological structure includes a plurality of pore bodies and a plurality of pore throats; Performing quasi-static drainage simulation based on the pore network topology; Calculating the corner meniscus area and the main end meniscus area based on the results of the quasi-static drainage simulation; A mapping relationship between the moisture content and the interface area of ​​the surface soil is obtained according to the corner curved liquid surface area and the main end curved liquid surface area.

2. The method for detecting gas-liquid interface distribution in surface soil according to claim 1, characterized in that: The method of obtaining a three-dimensional structural image of the surface soil comprises: Obtaining cylindrical surface soil samples; performing X-ray computed tomography scanning on the columnar surface soil sample; The X-ray computed tomography scan result is preprocessed to obtain the three-dimensional structural image.

3. The method for detecting gas-liquid interface distribution in surface soil according to claim 2, characterized in that: The preprocessing of the X-ray computed tomography scan result to obtain the three-dimensional structural image includes: removing edge region images and non-target region images from the X-ray computed tomography scan results to obtain a central region image; Performing Gaussian filtering on the central area image; The result of the Gaussian filtering is binarized to obtain the three-dimensional structure image.

4. The method for detecting gas-liquid interface distribution in surface soil according to claim 1, wherein: The performing geometric topological structure modeling on the three-dimensional structure image to obtain the pore network topological structure includes: performing Euclidean transformation and watershed segmentation on the three-dimensional structural image to determine the boundary range of each pore body and identify the pore throat between each pore body and other pore bodies; Calculating geometric parameters of each pore body and each pore throat; The topological connection relationship between each pore body and each pore throat is recorded.

5. The method for detecting gas-liquid interface distribution in surface soil according to claim 1, characterized in that: The quasi-static drainage simulation based on the pore network topology structure includes: Capillary pressure is used as the driving variable; The condition of gradually increasing external capillary pressure is set to simulate the gas-water two-phase distribution when the pore network topology forms an unsaturated state under the closure mechanism.

6. The method for detecting gas-liquid interface distribution in surface soil according to claim 1, characterized in that: The calculating of the corner meniscus area comprises: For any pore body, when the saturation of the pore body is greater than or equal to the critical water saturation of drainage, according to the formula Get the cross-sectional area A occupied by the liquid phase in the pore body w , where r wd is the liquid equivalent radius of the pore body at the initial stage of the quasi-static drainage simulation, areCoe is the area coefficient of the pore cross section of the pore body; According to the formula A nw =A cross -A w Get the cross-sectional area occupied by the gas phase in the pore body, where A cross is the cross-sectional area of ​​the pore body; Based on the law of volume conservation, according to the formula Calculate the length of the liquid column l w , where l is the total length of the pore body, s w is the saturation of the pore body, According to the formula Calculate the corneal circumference P per unit length nw , where β i is the angular half angle of the pore body, θ is the contact angle, According to the formula A AMs =P nw ·l nw +A cross Calculate the area A of the corner meniscus of the pore body AMs .

7. The method for detecting gas-liquid interface distribution in surface soil according to claim 1, characterized in that: The calculating of the corner meniscus area comprises: For any pore body, when the saturation of the pore body is less than the critical water saturation of drainage, according to the formula Get the water film radius r in the pore w , where A w is the cross-sectional area occupied by the liquid phase in the pore, A cross is the cross-sectional area of ​​the pore, s w is the saturation of the pore body, areCoe is the area coefficient of the pore cross section of the pore body; According to the formula Calculate the area A of the corner meniscus of the pore body AMs ; Wherein, l is the total length of the pore body, β i is the angular half angle of the pore body, and θ is the contact angle.

8. The method for detecting gas-liquid interface distribution in surface soil according to claim 1, characterized in that: The calculating of the main end meniscus area comprises: For any pore-throat connection pair, when the phase states of the pore body and pore throat in the pore-throat connection pair are different, according to the formula Calculate the effective radius r of the pore body-pore throat connection pair eff Among them, A cross is the cross-sectional area of ​​the pore body, and P is the cross-sectional perimeter of the pore body; According to the formula Calculate the capillary pressure at both ends of the pore body in the pore body-pore throat connection and the capillary pressure at both ends of the pore throat Where σ is the surface tension, θ is the contact angle, and α is the half-angle; According to the formula Calculate the pressure difference ΔP between interfaces c ; According to the formula Calculate the radius of curvature r of the main end meniscus type interface MTM ; According to the formula IN flag =min(A body ,IN throat )·favorite overlap Calculate the connection surface area A of the pore body and pore throat in the pore body-pore throat connection pair flag ; Among them, f overlap is the geometric overlap factor estimated from the corner opening angle, A body is the cross-sectional area of ​​the pore, A throat is the cross-sectional area of ​​the pore throat, α body is the average half-angle of the hole body angle, α throat is the average value of the half-angle of the pore throat angle; According to the formula Calculate the curvature correction coefficient ρ, where is the hydraulic radius of the pore body in the pore-throat connection pair, is the hydraulic radius of the pore throat in the pore-throat connection; According to the formula A MTM =A flag ·r Calculate the main end meniscus area A of the pore body-pore throat connection pair MTM .

9. The method for detecting gas-liquid interface distribution in surface soil according to claim 1, characterized in that: The step of obtaining a mapping relationship between the moisture content and the interface area of ​​the surface soil according to the angled curved liquid surface area and the main end curved liquid surface area includes: Accumulating the volumes of all pore bodies and all pore throats in the pore network topology to obtain a total effective volume; Accumulating all corner meniscus areas corresponding to the pore network topology to obtain a total corner meniscus area; The total main end meniscus area is obtained by accumulating the areas of all main end meniscus corresponding to the pore network topology; according to the formula Calculate the total interfacial area in, is the total corner meniscus area, is the total main end meniscus area; According to the formula Calculate and obtain the interface area a per unit volume of the surface soil nw ; Among them, V total is the total effective volume, is the porosity of the surface soil; Obtaining a depth-moisture content curve of the surface soil; determining the moisture content of the surface soil according to the depth-moisture content curve; A mapping relationship is established according to the moisture content and the interface area per unit volume.

10. A computer program product comprising a computer program, characterized in that When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 9 are implemented.

Images