An atmospheric numerical simulation method based on multiphase water material conservation constraints and application

Abstract

The present application relates to the field of atmospheric science and numerical calculation technology, and discloses a kind of atmospheric numerical simulation method and application based on the conservation constraint of multiphase water substance, and the evolution calculation of each phase water substance under the condition of discrete grid and discrete time step is oriented to the cloud microphysical parameterization of atmospheric numerical model.The method obtains the density and velocity field of each component in the grid unit, constructs the mass proportion variable normalized with wet air density;Define the growth rate per unit volume per unit time and represent the phase change source and sink with the component continuity equation;According to the mass proportion evolution equation, the coupling discrete updating rule is established, the total water substance growth rate residual error and the proportion residual error corresponding to the equivalent conservation expression are calculated, and the corrected growth rate that satisfies the conservation constraint is obtained through consistent correction to rewrite and update the mass proportion field.The method can inhibit the drift of water substance balance and reduce the phase distribution error under the condition of mixed phase cloud, and is suitable for numerical weather prediction and regional numerical simulation.

Description

Technical Field

[0001] This invention belongs to the field of atmospheric science and numerical computing technology, and relates to cloud microphysical parameterization schemes and the construction of multiphase water material conservation constraints in atmospheric numerical models, wet physical source-sink consistency correction and discrete time advancement error suppression. Specifically, it is an atmospheric numerical simulation method and application based on multiphase water material conservation constraints. Background Technology

[0002] Atmospheric numerical simulation and climate models are crucial technical tools for modern meteorological research and operational weather forecasting. Their simulation accuracy directly determines the accuracy of precipitation forecasts and the reliable description of atmospheric moisture cycles. In atmospheric numerical model systems, multiphase variables such as water vapor, cloud water, cloud ice, and precipitation particles are typically forecast on a grid, and subgrid processes such as phase transitions, collisional coalescence, deposition, and turbulent mixing are described using cloud microphysics parameterization. To adapt to dynamic frameworks, models can employ characterization methods such as specific humidity, mixing ratio, or mass concentration, and couple momentum, thermodynamic, and water substance equations in altitude or pressure coordinates to form a discrete solution system.

[0003] However, the water mass budget under three-phase coexistence conditions involves the coupling of multiple source and sink terms and energy exchange, and its conservation relationship is easily disrupted by various factors at the discrete level. First, cloud microphysics schemes often split the source and sink terms according to the process or update them item by item. The mass transfer between different phases may not be completely symmetrical numerically, leading to the accumulation of the total water mass budget residual over time. Second, precipitation deposition and ground flux introduce open boundary effects, meaning that the water mass within the control volume does not simply satisfy the approximation of a constant local proportion. At the same time, atmospheric compressibility causes changes in density and volume. If the selection of the denominator of the conservation quantities and the treatment of the transformation relationship are inconsistent, the non-closure error will be amplified. Third, processes such as spectral distribution, precipitation velocity, reevaporation, and refreezing rely heavily on empirical relationships, and the parameters are sensitive to time step and resolution, which can easily cause inconsistencies in the transfer between scales.

[0004] To mitigate numerical problems caused by imbalances, existing techniques typically employ flux discretization, mass-conserving advection schemes, source-sink term limiting, and numerical filtering, and assess wet process errors through water vapor budget diagnosis and budget closure checks. However, under conditions of strong convection, mixed-phase clouds, and persistent precipitation, frequent phase transitions and large source-sink term amplitudes can still lead to phenomena such as total water mass drift, local conservation violations, and inconsistencies between energy and water evolution, resulting in precipitation magnitude deviations and long-term integral drift.

[0005] In summary, existing atmospheric numerical simulations still have shortcomings in characterizing the conservation of multiphase water matter and in discrete budget closure, leading to problems such as decreased stability and accuracy. Therefore, how to effectively improve the numerical consistency and budget closure capability under the constraint of conservation of multiphase water matter in atmospheric numerical simulations is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0006] (a) Purpose of the invention

[0007] To address the aforementioned deficiencies and shortcomings of existing technologies, this invention aims to provide an atmospheric numerical simulation method and application based on the conservation constraints of multiphase water matter. By introducing verifiable conservation constraints and residual indices into grid discretization and time step progression, and implementing consistency coordination and adaptive correction for source and sink terms such as advection, turbulence, sedimentation, and phase transition, this method achieves multiphase water matter budget closure, reduces long-term integral drift, and improves numerical stability. It also provides recalculated diagnostic basis for wet process assessment in numerical weather prediction and climate simulation. Simultaneously, while ensuring computational efficiency, it reduces reliance on empirical limits and manual parameter tuning, and outputs updated multiphase fields and their conservation diagnostic results.

[0008] (II) Technical Solution

[0009] To achieve the objective of this invention and solve its technical problems, the present invention adopts the following technical solution:

[0010] The first objective of this invention is to provide an atmospheric numerical simulation method based on the conservation constraint of multiphase water matter, used in atmospheric numerical models to numerically calculate the evolution of multiphase atmospheric water matter, including water vapor, liquid water matter, and solid water matter, under discrete grid and discrete time step conditions. The method includes at least the following steps:

[0011] SS1. State Variable Acquisition and Construction: At the beginning of each time step, acquire the density of each component of the multiphase mixed atmosphere in each grid cell of the region to be simulated. r λ and velocity field subscript l include d , v , w , i , representing dry air, water vapor, liquid water, and solid water respectively, calculate the density of moist air. r = r d + r v , r d The density of dry air, r vGiven the water vapor density, calculate the mass percentage of water in each phase. q μ = r μ / r ,in q μ , r μ These represent the mass percentage and density of water substances in each phase, respectively, with subscripts. m include v , w , i ;

[0012] SS2. Growth Rate Definition and Source-Sink Expression: Controlling Volume in Mesh Cells α Define the mass of water substances in each phase. m μ and m μ = r μ · α And define the corresponding growth rate per unit volume per unit time. s μ and subscript m include v , w , i , t For time, while constraining the dry air growth rate s d and s d =0; Based on the component continuity equation, establish the density evolution relationship of water substances in each phase. The growth rate per unit volume per unit time s μ As mass source and sink terms caused by phase transitions and phase state transformations in the component continuity equation;

[0013] SS3. Discrete Coupling Update Rule: at time step Δ t Internally, based on the precise evolution equation of the mass percentage of water substances in each phase. Construct the corresponding discrete update quantity Δ q μ and And based on this, obtain the corresponding prediction quality percentage. and ;

[0014] SS4. Constructing Conservation Residuals: Calculating the Total Water Material Growth Rate Residual R s and R s = sv + s w + s i ,in s v , s w , s i The growth rate per unit volume per unit time for water vapor, liquid water, and solid water are respectively defined, and the proportion residual is constructed based on the equivalent conservation expression of the water mass conservation equation. R q and The equivalent conservation expression is as follows: ,in q v , q w , q i These represent the mass percentages of water vapor, liquid water, and solid water, respectively.

[0015] SS5. Consistency Correction and Write-back: Based on Constraints R s =0 and R q =0 determines the corrected growth rate of water substances in each phase. subscript m include v , w , i and use Alternative s μ Recalculate Δ q μ The corresponding corrected discrete update quantity is obtained. Then, the mass percentage of water substances in each phase after the time step was completed was obtained. and ;

[0016] SS6. Loop and Output: Repeat steps SS1~SS5 for subsequent time steps, and output the mass percentage of each phase of water within the target simulation period. q μ The spatiotemporal distribution and corresponding residuals R s , R q .

[0017] The second objective of this invention is to provide a computer program product, including computer instructions, for executing the steps of the above-described atmospheric numerical simulation method based on the constraint of conservation of multiphase water matter.

[0018] The third objective of this invention is to provide a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the above-described atmospheric numerical simulation method based on the constraint of conservation of multiphase water matter.

[0019] (III) Technical Effects

[0020] Compared with existing technologies, the atmospheric numerical simulation method and application based on the conservation constraint of multiphase water matter in this invention have the following beneficial and significant technical effects:

[0021] (1) Based on the mass ratio definition system of normalized humid air density, this invention constructs a characterization of the growth rate per unit volume per unit time around the mass of water substances in each phase and the volume of the control volume within the grid cell. The total water substance growth rate conservation relationship is used as a hard constraint. The source and sink terms corresponding to phase change and phase transformation processes such as condensation, evaporation, sublimation, freezing and melting are organized and verified in a conservation-consistent manner. In the discrete time step progression, through the calculation of conservation residuals and consistency correction, the water substance budget drift caused by the non-closure of source and sink terms and its cumulative error with the number of integration steps are significantly reduced, thereby improving the conservation, numerical stability and long-term integration reliability of cloud microphysics calculation.

[0022] (2) This invention uses the mass ratio precise evolution equation to establish the nonlinear coupling update rule of liquid and solid components on water vapor growth rate, and combines the equivalent conservation expression to construct the ratio residual, so as to realize the quantitative diagnosis and closed-loop correction of the systematic deviation introduced by the traditional linear approximation, thereby reducing the phase distribution error under mixed phase cloud conditions and improving the consistency of cloud water, cloud ice and precipitation particle evolution.

[0023] (3) The present invention outputs the mass proportion field and the conserved residual diagnostic quantity in the discrete grid and discrete time step of the numerical model, which can be used as the basis for online numerical inspection and error location of the microphysics parameterization scheme. It supports the recalculation evaluation of the conserved error under different time steps, resolutions and physical processes, reduces the dependence on empirical limiting and manual parameter tuning, and improves the efficiency of model debugging and cross-scenario portability. Attached Figure Description

[0024] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0025] Figure 1 This is a flowchart of the atmospheric numerical simulation method based on the conservation constraint of multiphase water matter according to the present invention;

[0026] Figure 2 The figure shows a comparison curve of the evolution of the conserved residuals with time steps in the warm cloud gas-liquid two-phase example. In the figure: (a) is the instantaneous growth rate residual | R s |Evolution curve over time step; (b) is the instantaneous mass percentage residual| R q |Evolution curve over time step; (c) is the cumulative growth rate residual Σ R s The evolution curve over time steps; (d) represents the cumulative mass percentage residual Σ R q Evolution curve over time step;

[0027] Figure 3 The figure shows a comparison curve of the evolution of the conserved residuals with time steps in the mixed-phase three-phase example. In the figure: (a) is the instantaneous growth rate residual | R s |Evolution curve over time step; (b) is the instantaneous mass percentage residual| R q |Evolution curve over time step; (c) is the cumulative growth rate residual Σ R s The evolution curve over time steps; (d) represents the cumulative mass percentage residual Σ R q Evolution curve over time steps. Detailed Implementation

[0028] This invention aims to provide an atmospheric numerical simulation method and application based on the constraint of multiphase water material conservation. To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions in the embodiments of this invention will be described in more detail below with reference to the accompanying drawings. The described embodiments are some, but not all, embodiments of this invention, and are exemplary and intended to explain this invention, and should not be construed as limiting this invention. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0029] Example 1: Atmospheric Numerical Simulation Method

[0030] As a specific example, such as Figure 1As shown in the embodiments of the present invention, the atmospheric numerical simulation method based on the conservation constraint of multiphase water matter is used to numerically calculate the evolution of multiphase atmospheric water matter, including water vapor, liquid water matter, and solid water matter, under discrete grid and discrete time step conditions in the cloud microphysics parameterization scheme of atmospheric numerical models. Its core lies in constructing a dual diagnostic constraint using growth rate conservation and equivalent conservation expressions, and suppressing the accumulation of conservation bias with integration through correction write-back. Its implementation mainly includes the following steps:

[0031] SS1. State acquisition and construction:

[0032] At the beginning of each time step, the densities of each component of the multiphase mixed atmosphere in each grid cell of the region to be simulated are obtained. r λ and velocity field subscript l include d , v , w , i , representing dry air, water vapor, liquid water, and solid water respectively, calculate the density of moist air. r = r d + r v , r d The density of dry air, r v Given the water vapor density, calculate the mass percentage of water in each phase. q μ = r μ / r ,in q μ , r μ These represent the mass percentage and density of water substances in each phase, respectively, with subscripts. m include v , w , i .

[0033] Preferably, in step SS1, after obtaining the density of each component... r λ And calculate the density of moist air. r Then, the atmospheric density of multiphase mixtures is defined. r s and r s = r d + r v + r w + r i ,in r d , r v , r w , r i The densities of dry air, water vapor, liquid water, and solid water are respectively defined as the total water mass density. r t and r t = r v + r w + r i Based on this, the total water mass percentage was constructed. q t and q t = r t / r At the same time, construct the dry air quality ratio q d and q d = r d / r Based on the defined relationship of quality proportion, establish a consistency constraint on quality proportion, so that... q d + q v =1、 q t = q v + q w + q i ,in q d , q v , q w , q i The mass percentages of dry air, water vapor, liquid water, and solid water were respectively determined, and the mass percentage of the multiphase mixed atmosphere relative to moist air was further constructed. q s and q s = q d + q v + q w + qi and q s =1+ q w + q i and will q s Additional constraint variables used for subsequent conservation diagnostics, error analysis, or numerical checks, at each time step. q s The deviation is calculated, and when the deviation exceeds a preset tolerance, [the following is taken]: r d , r v , r w , r i Or obtained from its calculation q v , q w , q i Perform a standardization process to ensure that each mass percentage variable remains consistent under the same denominator system and to reduce nonlinear coupling errors caused by inconsistencies in the denominator system.

[0034] It should be noted that step SS1 not only completes the reading of state variables and the construction of proportions, but also serves to form a recalcible normalization benchmark and consistency check caliber. Specifically, it uses the density of moist air as an example. r = r d + r v Non-multiphase mixed atmospheric density r s As a percentage of quality q μ (including) q d , q t , q s The unified denominator is the core design decision for achieving multiphase coupling conservation in this method. This choice effectively avoids the disturbance of the denominator system by changes in condensate density, thus... q μ Maintaining the inherent consistency of the nonlinear coupling structure during phase transitions allows for rapid identification of inconsistencies in proportions caused by input field interpolation, boundary updates, or numerical rounding at each time step. This lays the foundation for denominator stability in the rigorous derivation of the coupling subtraction term in step SS3. Simultaneously... q s =1+ q w + qi The additional constraints provide a quantitative criterion for online verification of denominator consistency at each time step, significantly improving the numerical reliability of the multiphase density field initialization. Furthermore, the recording of deviations can serve as prior information for subsequent residual source attribution, thereby enhancing the interpretability of conservation diagnostics and strengthening the robustness of numerical implementation.

[0035] SS2. Definition of growth rate and source / sink expression:

[0036] In the control volume of the grid cell α Define the mass of water substances in each phase. m μ and m μ = r μ · α And define the corresponding growth rate per unit volume per unit time. s μ and subscript m include v , w , i , t For time, while constraining the dry air growth rate s d and s d =0; Based on the component continuity equation, establish the density evolution relationship of water substances in each phase. The growth rate per unit volume per unit time s μ As a mass source and sink term caused by phase transitions and phase state transformations in the component continuity equation.

[0037] Preferably, in step SS2, the growth rate of each phase of water substance is... s μ The phase transition and phase state transition mass fluxes output by the cloud microphysics parameterization scheme are constituted and organized according to the mass conservation transfer pairs. Water vapor to liquid and liquid to water vapor, as well as water vapor to solid and solid to water vapor, are represented by source and sink terms with opposite signs. Simultaneously, the freezing and melting processes between liquid and solid are written with opposite signs. s w and s i To ensure that within the same grid cell, the following conditions are met. s v + s w + s i =0, ensuring that phase transitions only change the phase components without introducing additional mass sources.

[0038] Furthermore, in step SS2, in addition to establishing the component continuity equations for each phase of water, a total continuity equation is also established based on the total mass conservation of the multiphase mixed atmosphere, thus ensuring the density of the multiphase mixed atmosphere... r s Satisfy the velocity field The conservation relationship of change, and based on this, The calculation results are checked for consistency. When the divergence terms of the overall continuity equation and the component continuity equations are inconsistent, the same divergence term is used to update the continuity equations of each component synchronously.

[0039] It should be noted that the key point of step SS2 is to convert the microphysical process output into a growth rate source-sink expression that satisfies the conservation structure. This is achieved by converting the output of the cloud microphysical parameterization scheme into a growth rate... s μ The control volume organization ensures that source and sink terms maintain dimensional consistency under coordinate transformation and adaptive resolution conditions, thus achieving strict decoupling from divergence transport terms and enhancing the method's portability to different mesh frameworks. By transferring the organization and cancellation constraints, the risk of source-sink non-closure introduced when different physical processes are superimposed within discrete walks can be reduced. Simultaneously, using the overall continuity equation for divergence term consistency verification avoids spurious quality changes in the component equations due to inconsistent divergence apertures, thereby providing a stable source-sink basis for subsequent proportion coupling updates and residual diagnosis. Furthermore, dry air growth rate... s d The constraint of 0 not only directly guarantees the conservation of dry air mass, but also ensures the velocity divergence through the derivation of the continuity equation in which it participates. The physical self-consistency of the expression eliminates potential sources of divergence error for the accurate derivation of the mass ratio coupling evolution equation in step SS3.

[0040] SS3. Discrete Coupling Update Rules:

[0041] At time step Δ t Internally, based on the precise evolution equation of the mass percentage of water substances in each phase. Construct the corresponding discrete update quantity Δ q μ and And based on this, obtain the corresponding prediction quality percentage. and .

[0042] Preferably, in step SS3, the precise evolution equation of the mass ratio is expanded into three phase update rules in a coupling manner with the water vapor growth rate, so that the water vapor mass ratio satisfies The mass ratio of liquid water meets the requirements. The solid water mass ratio meets the requirements. And at time step Δ tEach discrete increment is constructed to reflect the nonlinear coupling constraint between the liquid phase and the solid phase on the growth rate of the water vapor phase.

[0043] It should be noted that step SS3 introduces a coupled subtraction term into the exact evolution equation of the proportion. s v q μ / r This is not an empirical correction term, but rather a rigorous derivation obtained jointly from the denominator design in step SS1 and the component continuity equation in step SS2. Its physical meaning lies in quantifying the feedback effect of the time-varying variability of the moist air reference density with respect to the water vapor phase change on the evolution of the proportions of each phase. The introduction of this coupling term allows the discrete updates of the liquid and solid phases to be transmitted via the water vapor growth rate. s v This implicit constraint fundamentally overcomes the systematic nonlinear conservation error caused by neglecting the time-varying denominator in traditional phase-by-phase independent update schemes. Furthermore, the introduction of this coupling term allows the update of non-gaseous components to automatically sense the normalization effect of the water vapor component, thereby avoiding the inversion of Δ... q w With Δ q v The systemic bias caused by simply setting it as the opposite quantity is eliminated; at the same time, the predicted proportion, as an intermediate state, provides a unified input for subsequent residual calculation and correction write-back, which facilitates the decoupling of the prediction stage and the correction stage at the implementation level, and improves the maintainability and recalculation of the numerical solution process.

[0044] SS4. Conservative Residual Construction:

[0045] Calculate the residual of total water material growth rate R s and R s = s v + s w + s i ,in s v , s w , s i The growth rate per unit volume per unit time for water vapor, liquid water, and solid water are respectively defined, and the proportion residual is constructed based on the equivalent conservation expression of the water mass conservation equation. R q and The equivalent conservation expression is as follows: ,in q v , q w , qi These represent the mass percentages of water vapor, liquid water, and solid water, respectively.

[0046] Preferably, in step SS4, in addition to calculating the residual of total water matter growth rate... R s With proportion residual R q In addition, the percentage of total water mass is calculated. q t At time step Δ t The change Δ within q t and Δ q t and R s Combined for conservation diagnosis, among which R s Δ is used to characterize the conservation of water mass transfer per unit volume per unit time. q t Used to characterize the percentage change with moist air density as a reference; when Δ q t Not zero and R s When the value is zero, it is determined that there is a change in the proportion caused by density change or transport process at this time step, rather than a mass non-conservation caused by phase transition.

[0047] It should be noted that in step SS4, the residual of the total water mass growth rate is constructed in parallel. R s With proportion residual R q These can respectively indicate the source-sink closure error and the normalized coupling error. In engineering implementation, when... R s abnormal R q Under normal conditions, priority should be given to investigating microphysical source and sink structures; otherwise, priority should be given to investigating... r Update and adjust the calculation method for percentages; further introduce Δ q t It can be used to distinguish between proportionate changes caused by transport and mass changes caused by phase transitions, thereby improving the localization efficiency and interpretability of residual diagnosis. Furthermore, proportionate residuals... R q The construction form originates from the equivalent algebraic transformation of the total water mass conservation condition within the mass proportion framework, where Δ q v / (1- q v ) and (Δ q w +Δ q i ) / (1+q w + q i The weighted structure quantifies the contribution of the proportion changes of the water vapor phase and the non-gaseous water phase, respectively. The sum of the two is zero, which indicates that the conservation of total water matter is satisfied in the humid air reference denominator system. The dual residual system elevates the conservation judgment at the growth rate level to a refined conservation diagnosis at the mass proportion level, significantly enhancing the sensitivity to the identification of implicit nonlinear conservation errors and providing a reliable quantitative input for the accurate correction of step SS5.

[0048] SS5. Consistency Correction and Write-back:

[0049] Based on constraints R s =0 and R q =0 determines the corrected growth rate of water substances in each phase. subscript m include v , w , i and use Alternative s μ Recalculate Δ q μ The corresponding corrected discrete update quantity is obtained. Then, the mass percentage of water substances in each phase after the time step was completed was obtained. and .

[0050] As a preferred option, steps SS3 and SS5 are set for warm cloud gas-liquid two-phase conditions. q i =0, and at time step Δ t The update constraint for the gas-liquid two-phase system is constructed based on the nonlinear coupling relationship within the mass proportion definition framework, such that the increment Δ of the liquid water mass proportion is... q w The increase in the mass ratio of water vapor Δ q v It satisfies the correspondence including the denominator correction term, and the denominator correction term starts from the current time step. q v and q w Determined, the replacement will be Δ q w With Δ q v We simply approximate it as a linear approximation with opposite signs.

[0051] In addition, steps SS3 and SS5 introduce the amount of non-gaseous water under the three-phase coexistence condition of gas, liquid, and solid. qwi and q wi = q w + q i And based on the time step starting q t or q wi As a condition for determining the integration constant, q wi and q v The constraints between them are constructed within the time step, so that the corrected growth rate satisfies R s =0 and R q At the same time, the changes in the non-gaseous water mass and the water vapor ratio satisfy the coupled evolution law constrained by the initial water mass ratio condition.

[0052] Preferably, in step SS5, the growth rate is corrected. Obtained through iterative uniformity solution, with the updated solution used during the iteration process. R s and R q As a convergence criterion, the corrected growth rate is substituted into step SS3 to recalculate Δ in each iteration. q μ And update until R s , R q The iteration ends after the convergence condition is met, thereby reducing the risk of overcorrection caused by a single correction and improving the numerical stability of time step updates.

[0053] It should be noted that in step SS5, with R s =0 and R q =0 represents a double-constraint iterative solution for the corrected growth rate, on the one hand through... R s Constraints ensure source-sink closure, and on the other hand, through R qThe constraint-suppressed ratio normalization of coupling bias is mathematically equivalent to fixed-point iteration of the coupled nonlinear equation system, rather than traditional linear proportional allocation correction. Its essential advantage lies in simultaneously addressing errors at both the closed growth rate conservation and ratio conservation levels. This design enables the correction process to adaptively handle the differences in coupling structures under different phase coexistence conditions, such as gas-liquid two-phase and gas-liquid-solid three-phase scenarios, avoiding the residual nonlinear errors introduced by traditional single-constraint correction schemes that neglect the ratio coupling relationship. Thus, it balances numerical accuracy and methodological robustness in engineering practice. Furthermore, introducing two-phase constraints and non-gaseous resultant variables under warm cloud and mixed-phase cloud conditions respectively can reduce error abrupt changes during phase combination switching. Combined with iterative uniform solution, this improves numerical stability and recalculation under complex conditions.

[0054] SS6. Loops and Output:

[0055] Repeat steps SS1 to SS5 for subsequent time steps, and output the mass percentage of each phase of water within the target simulation period. q μ The spatiotemporal distribution and corresponding residuals R s , R q Preferably, in step SS6, in addition to outputting the mass percentage of water substances in each phase, q μ Spatiotemporal distribution and residuals R s , R q In addition, it also R s and R q Cumulative statistics are performed during the simulation period to form a conserved drift diagnostic value, and this conserved drift diagnostic value is then compared with... q μ The spatiotemporal distribution is written to the storage medium according to the grid index and time index, so that it can be read in subsequent time steps and used for the conservation evaluation and numerical debugging of cloud microphysical parameterization scheme.

[0056] It should be noted that in step SS6, the cumulative statistics of the conserved drift diagnostic quantity are based on the time step as the smallest unit. R s and R q Stepwise integration can distinguish the different contribution modes of short-term high-frequency phase transition processes and long-term slowly varying accumulated errors to conservation, providing quantitative spatiotemporally resolved diagnostic basis for parameter sensitivity analysis and numerical optimization of cloud microphysics parameterization schemes. Furthermore, through... R s and R qBy performing cumulative statistics, it is possible to identify whether the conservation deviation accumulates with each step and the spatial distribution characteristics of the accumulation rate, thereby providing a quantitative basis for setting the time step, combining parameterized processes, and handling boundary conditions.

[0057] Example 2: Verification Example

[0058] To verify the feasibility and technical effectiveness of the method of the present invention in engineering numerical models, based on the above steps SS1~SS6, this embodiment further selects two representative cloud microphysics computing scenarios for comparative verification, namely, the warm cloud gas-liquid two-phase calculation example (calculated for 300 time steps, time step size Δ t =1s) and a three-phase calculation example of mixed phase gas-liquid-solid (calculated over 500 time steps, time step size Δ t =1s). The vertical axis for both types of examples is taken on a logarithmic scale with base 10, and the convergence tolerance reference line is set to 10. -10 The comparison scheme uses a traditional conservation approximation scheme, while the scheme proposed in this paper uses the dual residual correction scheme. Both types of examples output the mass proportion field and residual diagnostic quantity, and use instantaneous residual and cumulative residual as the core evaluation indicators for statistical comparison.

[0059] Example 1: Verification of Warm Cloud Gas-Liquid Two-Phase Process: A warm cloud process containing only water vapor and liquid water is selected, and the mass percentage of solid water is set. q i =0, and at each time step n Obtain the mesh cells according to step SS1. Oh j Control volume α j Density of dry air, water vapor and liquid water and velocity field Calculate the density of moist air Construct a mass percentage variable normalized to moist air density. , In the microphysical process, only condensation and evaporation-type phase transition sources and sinks are retained, which are then converted into the growth rate per unit volume per unit time. And according to the transfer to the tissue, the total water matter growth rate residual is reduced. The conservation closure condition is satisfied. Time progression uses a uniform time step Δ. t The comparison scheme uses a traditional linear approximation with an incremental constraint Δ. q w ≈-Δ q v The present invention uses the proportion of step SS3 to precisely evolve discrete update amount. And obtain the predicted proportion In step SS4, construct the proportion residual. Then, in step SS5, the corrected growth rate is solved. And write back and update .

[0060] Figure 2 The figure shows a comparison curve of the evolution of the conserved residuals over time in the warm cloud gas-liquid two-phase calculation example. As can be seen from the figure, in the comparison scheme... R s |and| R q | Maintained at 10 throughout the entire calculation process -5 ~10 -3 The magnitude of the cumulative residuals for both increases linearly with time steps, reaching 10 by the 500th step. -1 The magnitude is far beyond the convergence tolerance; the solution of this invention suppresses the instantaneous residual to the convergence tolerance (10) throughout the entire process. -10 Below 10, the cumulative residual stabilizes at 10. -8 Below the order of magnitude, it is reduced by about 7 to 8 orders of magnitude compared to the comparison scheme, and no residual rebound or divergence occurs throughout the process, effectively eliminating the problem of accumulation of conservation error in the gas-liquid two-phase warm cloud scenario.

[0061] Example 2, Verification of Mixed-Phase Gas-Liquid-Solid Three-Phase Process: A mixed-phase cloud process with coexisting gas, liquid, and solid phases is selected. Mesh elements are acquired synchronously in step SS1. Oh j Density and velocity fields of the four internal components Calculate the density of moist air And construct the quality ratio Where μ = v, w, i. Microphysical sources and sinks include processes such as condensation and evaporation, sublimation, and freeze-thaw cycles. The growth rate is calculated using SS2. And according to the transfer to the organization, the residual of total water matter growth rate within the same grid cell is reduced. The mass fraction of non-gaseous water tends towards zero. To strengthen the engineering constraints of three-phase coupling, a non-gaseous water mass fraction variable is introduced into SS5. and with time step initial conditions or As the basis for determining the integration constant, a coupling constraint under three-phase conditions is constructed so that the corrected update simultaneously satisfies... With proportion residual The dual constraints. The comparative scheme still uses traditional approximate constraints and directly advances the proportion field. The present invention uses a predicted proportion. With correction increment Write-back forms the final update .

[0062] Figure 3 The figure shows a comparison curve of the evolution of the conserved residuals over time steps in a mixed-phase three-phase example. From... Figure 3As can be seen, the coupling of phase change source terms is more complex under the condition of three-phase coexistence, compared with the scheme in | R s |and| R q The envelope exhibits a bimodal oscillation pattern, influenced by both liquid-phase precipitation and solid-phase nucleation, with the peak value reaching 10. -4 The magnitude is such that the accumulated residuals exceed 10 at the end of the calculation. -1 The magnitude of the error severely compromises numerical reliability; the solution in this invention stabilizes and suppresses the instantaneous residual to 10. -11 ~10 -10 The magnitude of the double-peaked oscillation structure was completely eliminated, and the cumulative residual remained at 10 throughout the entire process. -8 Below the order of magnitude, it is reduced by about 7 to 8 orders of magnitude compared to the comparison scheme, and both the growth rate and quality ratio constraints are met.

[0063] The verification results of the two types of examples show that, regardless of whether it is a gas-liquid two-phase warm cloud scenario or a gas-liquid-solid three-phase mixed-phase scenario, the proposed solution can stably control the instantaneous residual and cumulative residual below the convergence tolerance, without generating a conservation deviation that accumulates over time. This fully demonstrates that the proposed dual residual correction scheme has good robustness and accuracy assurance capabilities in cloud microphysics calculation scenarios with different phase combinations, and can be effectively applied to the improvement of microphysics schemes in engineering numerical weather prediction models.

[0064] The objectives of this invention have been fully and effectively achieved through the above embodiments. Those skilled in the art will understand that this invention includes, but is not limited to, the contents described in the accompanying drawings and the specific embodiments described above. Although the invention has been described with reference to what is currently considered the most practical and preferred embodiments, it should be understood that the invention is not limited to the disclosed embodiments, and any modifications that do not depart from the functional and structural principles of the invention will be included within the scope of the claims.

Claims

1. A method of atmospheric numerical simulation based on the conservation of multiphase water substances, characterized in that, It should include at least the following steps: SS1. At the beginning of each time step, obtain the density of each component of the multiphase mixed atmosphere in each grid cell of the region to be simulated. ρ λ and velocity field subscript λ include d , v , w , i , representing dry air, water vapor, liquid water, and solid water respectively, calculate the density of moist air. ρ = ρ d + ρ v , ρ d The density of dry air, ρ v Given the water vapor density, calculate the mass percentage of water in each phase. q μ = ρ μ / ρ ,in q μ , ρ μ These represent the mass percentage and density of water substances in each phase, respectively, with subscripts. μ include v , w , i ; SS2. Control volume in mesh cells α Define the mass of water substances in each phase. m μ and m μ = ρ μ · α And define the corresponding growth rate per unit volume per unit time. s μ and subscript μ include v , w , i , t For time, while constraining the dry air growth rate s d and s d =0; Based on the component continuity equation, establish the density evolution relationship of water substances in each phase. The growth rate per unit volume per unit time s μ As mass source and sink terms caused by phase transitions and phase state transformations in the component continuity equation; SS3. At time step Δ t Internally, based on the precise evolution equation of the mass percentage of water substances in each phase. Construct the corresponding discrete update quantity Δ q μ and And based on this, obtain the corresponding prediction quality percentage. and ; SS4. Calculate the residual of total water mass growth rate. R s and R s = s v + s w + s i ,in s v , s w , s i The growth rate per unit volume per unit time for water vapor, liquid water, and solid water are respectively defined, and the proportion residual is constructed based on the equivalent conservation expression of the water mass conservation equation. R q and The equivalent conservation expression is as follows: ,in q v , q w , q i These represent the mass percentages of water vapor, liquid water, and solid water, respectively. SS5. Based on constraints R s =0 and R q =0 determines the corrected growth rate of water substances in each phase. subscript μ include v , w , i and use Alternative s μ Recalculate Δ q μ The corresponding corrected discrete update quantity is obtained. Then, the mass percentage of water substances in each phase after the time step was completed was obtained. and ; SS6. Repeat steps SS1 to SS5 for subsequent time steps, and output the mass percentage of each phase of water within the target simulation period. q μ The spatiotemporal distribution and corresponding residuals R s , R q .

2. The method according to claim 1, characterized in that, In step SS1, the density of each component is obtained. ρ λ And calculate the density of moist air. ρ Then, the atmospheric density of multiphase mixtures is defined. ρ s and ρ s = ρ d + ρ v + ρ w + ρ i ,in ρ d , ρ v , ρ w , ρ i The densities of dry air, water vapor, liquid water, and solid water are respectively defined as the total water mass density. ρ t and ρ t = ρ v + ρ w + ρ i Based on this, the total water mass percentage was constructed. q t and q t = ρ t / ρ At the same time, construct the dry air quality ratio q d and q d = ρ d / ρ Based on the defined relationship of quality proportion, establish a consistency constraint on quality proportion, so that... q d + q v =1、 q t = q v + q w + q i ,in q d , q v , q w , q i The mass percentages of dry air, water vapor, liquid water, and solid water were respectively determined, and the mass percentage of the multiphase mixed atmosphere relative to moist air was further constructed. q s and q s = q d + q v + q w + q i and q s =1+ q w + q i and will q s Additional constraint variables used for subsequent conservation diagnostics, error analysis, or numerical checks, at each time step. q s The deviation is calculated, and when the deviation exceeds a preset tolerance, [the following is taken]: ρ d , ρ v , ρ w , ρ i Or obtained from its calculation q v , q w , q i Perform consistency processing.

3. The method according to claim 1, characterized in that, In step SS2, the growth rate of water substances in each phase. s μ The phase transition and phase state transition mass fluxes output by the cloud microphysics parameterization scheme are constituted and organized according to the mass conservation transfer pairs. Water vapor to liquid and liquid to water vapor, as well as water vapor to solid and solid to water vapor, are represented by source and sink terms with opposite signs. Simultaneously, the freezing and melting processes between liquid and solid are written with opposite signs. s w and s i To ensure that within the same grid cell, the following conditions are met. s v + s w + s i =0, ensuring that phase transitions only change the phase components without introducing additional mass sources.

4. The method according to claim 1 or 3, characterized in that, In step SS2, in addition to establishing the component continuity equations for each phase of water, a total continuity equation is also established based on the total mass conservation of the multiphase mixed atmosphere, thus ensuring the density of the multiphase mixed atmosphere. ρ s Satisfy the velocity field The conservation relationship of change, and based on this, The calculation results are checked for consistency. When the divergence terms of the overall continuity equation and the component continuity equations are inconsistent, the same divergence term is used to update the continuity equations of each component synchronously.

5. The method according to claim 2, characterized in that, In step SS3, the mass fraction evolution equation is expanded into three phase update rules in a water vapor growth rate coupling manner, so that the water vapor mass fraction satisfies , the liquid water mass fraction satisfies , and the solid water mass fraction satisfies , and the corresponding discrete increments are respectively constructed within a time step Δ t to reflect the nonlinear coupling constraints of the liquid phase and the solid phase on the water vapor phase growth rate.

6. The method according to claim 5, characterized in that, In step SS4, in addition to calculating the residual of total water mass growth rate... R s With proportion residual R q In addition, the percentage of total water mass is calculated. q t At time step Δ t The change Δ within q t and Δ q t and R s Combined for conservation diagnosis, among which R s Δ is used to characterize the conservation of water mass transfer per unit volume per unit time. q t Used to characterize the percentage change with moist air density as a reference; when Δ q t Not zero and R s When the value is zero, it is determined that there is a change in the proportion caused by density change or transport process at this time step, rather than a mass non-conservation caused by phase transition.

7. The method according to claim 6, characterized in that, Steps SS3 and SS5 are set for the gas-liquid two-phase conditions of the warm cloud. q i =0, and at time step Δ t The update constraint for the gas-liquid two-phase system is constructed based on the nonlinear coupling relationship within the mass proportion definition framework, such that the increment Δ of the liquid water mass proportion is... q w The increase in the proportion of water vapor mass Δ q v It satisfies the correspondence including the denominator correction term, and the denominator correction term starts from the current time step. q v and q w Determined, the replacement will be Δ q w With Δ q v We simply approximate it as a linear approximation with opposite signs.

8. The method according to claim 7, characterized in that, Steps SS3 and SS5 introduce the amount of non-gaseous water mass under the condition of gas-liquid-solid three-phase coexistence. q wi and q wi = q w + q i And based on the time step starting q t or q wi As a condition for determining the integration constant, q wi and q v The constraints between them are constructed within the time step, so that the corrected growth rate satisfies R s =0 and R q At the same time, the changes in the non-gaseous water mass and the water vapor ratio satisfy the coupled evolution law constrained by the initial water mass ratio condition.

9. The method according to claim 8, characterized in that, In step SS5, the growth rate is corrected. Obtained through iterative uniformity solution, with the updated solution used during the iteration process. R s and R q As a convergence criterion; In each iteration, the corrected growth rate is substituted into step SS3 to recalculate Δ. q μ And update until R s , R q The iteration ends when the convergence condition is met.

10. The method according to claim 9, characterized in that, In step SS6, in addition to outputting the mass percentage of water substances in each phase, q μ Spatiotemporal distribution and residuals R s , R q In addition, it also R s and R q Cumulative statistics are performed during the simulation period to form a conserved drift diagnostic value, and this conserved drift diagnostic value is then compared with... q μ The spatiotemporal distribution is written to the storage medium according to the grid index and time index, so that it can be read in subsequent time steps and used for the conservation evaluation and numerical debugging of cloud microphysical parameterization scheme.

11. A computer program product comprising computer instructions, characterized in that, The computer instructions are used to execute the steps of the atmospheric numerical simulation method based on the conservation constraint of multiphase water matter as described in any one of claims 1 to 10.

12. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the atmospheric numerical simulation method based on the conservation constraint of multiphase water matter as described in any one of claims 1 to 10.

Images