Machine Learning Methods for Thermodynamic Energy Flow Pattern Prediction at Regional Scales with High-Dimensional Climate Data

2.1 Overall framework design concept

The core logic of the dual-path physics-data hybrid framework is to achieve deep collaborative integration between thermodynamic physical laws and high-dimensional climate data, in order to overcome the adaptation limitations of traditional models under complex conditions. The framework uses the regional energy balance equation as the physical skeleton and integrates thermodynamic conservation laws into the model structure through differentiable modeling, ensuring that all prediction results strictly follow basic physical principles. At the same time, high-dimensional climate data drives the model through a spatiotemporal neural parameterization unit, which adaptively generates spatiotemporally heterogeneous key parameterization terms, accurately capturing the unique characteristics of complex underlying surfaces and non-uniform environments. The complete link of the framework is as follows: inputs include initial and boundary conditions such as surface temperature, radiation forcing, and high-dimensional climate data like wind, temperature, humidity, pressure, and surface types; the core modules consist of differentiable physical kernels and spatiotemporal neural parameterization units, which are collaboratively optimized through gradient backpropagation; outputs include high-precision energy flow spatial distribution predictions, regional thermodynamic parameterization fields, and corresponding uncertainty quantification results. The core innovation of this design lies in elevating physical laws from external fixed constraints in traditional models to discoverable and interpretable structures within the model, achieving an organic unity between physical skeleton and data-driven dynamics. The overall framework architecture is shown in Figure 1.

Figure 1. System flow architecture diagram

2.2 Design of differentiable forward physics simulator

Based on the thermodynamic energy conservation principle, the core control equation for energy balance per unit area of the surface at the regional scale can be expressed as:

${{\text{R}}_{\text{n}}}\text{=H+LE+G}$    (1)

where, ${{R}_{n}}$ is the net radiation flux at the surface, representing the total radiation energy received by the unit area of the surface per unit time; $H$ is the sensible heat flux, describing the sensible heat exchange between the surface and the atmosphere; $LE$ is the latent heat flux, corresponding to the latent heat transfer associated with water vapor phase change; $G$ is the soil heat flux, reflecting the energy conduction process between the surface and the soil. This equation forms the core constraint of the physical simulation, clarifying the basic balance relationship of regional energy flow transfer. To adapt to regional-scale simulation requirements, the finite volume method is used for spatiotemporal discretization. This method leverages the inherent conservation properties of integration to ensure that thermodynamic conservation laws are not violated during discretization. In space, the study area is divided into uniform grids, and the energy balance equation is integrated over each grid unit to obtain the integral form of the control equation at the grid scale; in time, a first-order implicit finite difference scheme is used to discretize the time derivative term, resulting in the discretized control equation:

$\frac{1}{\Delta t} \int_{\Omega_i}\left(T_{t+\Delta t}-T_t\right) \rho c~d \Omega=\int_{\Omega}\left(R_n-H-L E-G\right) d \Omega$    (2)

where, ${{\text{ }\!\!\Omega\!\!\text{ }}_{i}}$ is the control volume of the $i$-th grid unit, $\Delta t$ is the time step, $T$ is the surface temperature, and $\rho$ and $c$ are the density and specific heat capacity of the surface medium, respectively. The numerical solution steps are as follows: first, initialize the initial temperature and boundary fluxes for the grid units; then, substitute the parameterization terms generated by the neural parameterizer to solve for the sensible and latent heat fluxes; and iteratively update the energy balance state of each grid unit until the numerical solution meets the preset convergence criteria, ensuring that the discretized numerical solution maintains physical reasonableness.

To achieve collaborative training of the physical kernel and the neural parameterizer, automatic differentiation is used to convert the full numerical computation process of the physics simulator into a differentiable computational graph. By tracking the gradient propagation paths of each numerical operation, the gradients of the parameterization terms can be backpropagated along the numerical solution chain to the parameter space of the neural parameterizer. Specifically, the core steps of the finite volume method, including integral computation, flux interpolation, and iterative solving, are all encapsulated as differentiable operators to ensure that small changes in the parameterization terms are accurately transmitted and converted into update signals for the neural parameterizer. To balance numerical stability and gradient transmissibility, a step-size adaptation strategy is designed: the time step is dynamically optimized based on the residual size of the current iteration step. When the residual is large, the step size is reduced to improve stability, and when the residual converges, the step size is increased to enhance computational efficiency. At the same time, the QUICK scheme is used as the spatial interpolation scheme, as its inherent higher-order smoothing properties can reduce the gradient discretization error caused by numerical diffusion, ensuring the continuity of gradient propagation. The simulator’s inputs include surface initial temperature, top-of-atmosphere radiation forcing, near-surface meteorological variables, and other initial and boundary conditions. The outputs include energy flow fields such as sensible heat flux, latent heat flux, surface temperature, and state variables. Among them, input variables must satisfy thermodynamic dimensional consistency constraints, such as radiation flux with units of $\text{W}\cdot {{\text{m}}^{-2}}$ and temperature in Kelvin; output energy flow fields must meet the residual threshold requirements of the energy balance equation to ensure the physical validity of the output results. Figure 2 provides an illustrative diagram of the regional energy balance control volume analysis based on the finite volume method, showing the energy input and output terms (R_n, H, LE, G) within the grid unit and the time/space discretized grid structure, in conjunction with the physical meanings of Eqs. (1) and (2).

Figure 2. Regional energy balance control volume analysis diagram based on finite volume method

2.3 Design of spatiotemporal neural parameterizer

The spatiotemporal neural parameterizer adopts an encoder-decoder architecture. Its core function is to adaptively learn from high-dimensional climate data and generate parameterization fields that are fully matched with the physics simulator’s grid, achieving precise characterization of spatiotemporal heterogeneity of key parameters such as turbulence exchange coefficients and surface impedance. Figure 3 provides a detailed structure diagram of the encoder-decoder. It shows how 3D convolution layers extract spatiotemporal features, Bi-GRU captures long-term dependencies, the spatiotemporal attention module (Attention) is connected, and the embedding of physical prior constraints at the end of the network. The encoder part consists of 3 layers of 3D convolution layers and 1 layer of bidirectional gated recurrent units (Bi-GRU). The 3D convolution layers use 5 × 5 × 3 spatiotemporal convolution kernels to extract local spatiotemporal features from wind, temperature, humidity, pressure, and surface type data, while the Bi-GRU units capture long-term dependencies of climate variables along the time dimension. The decoder uses a fusion design of transpose convolution layers and spatiotemporal attention mechanisms. The transpose convolution layers upsample the feature maps to restore the spatial resolution consistent with the physical grid, while the spatiotemporal attention module calculates the correlation weights of different spatiotemporal positions in the feature map to enhance critical information. The attention weight calculation formula is:

$\alpha_{i, j, t}=\frac{exp \left({sim}\left(F_{i, j, t}, F_{i^{\prime}, j^{\prime}, t^{\prime}}\right)\right)}{\sum_{t^l, j^l, t^l} exp \left({sim}\left(F_{i, j, t}, F_{i^{\prime}, j^{\prime}, t^{\prime}}\right)\right)}$    (3)

where, ${{\alpha }_{i,j,t}}$ is the attention weight at grid point (i, j) at time $t$, ${{F}_{i,j,t}}$ is the corresponding feature vector at that position, and $\text{sim}\left( \cdot  \right)$ uses cosine similarity for computation. The network output layer directly outputs the parameterization field corresponding to each grid point of the physics simulator, through a combination of the Sigmoid activation function and linear transformation. To ensure the thermodynamic validity of the output parameters, two physical prior constraints are embedded in the network design: first, a constraint on the range of parameter values, by adding hard threshold activations to the output layer, which limits the turbulence exchange coefficient to the range of ${{10}^{-3}}$ to ${{10}^{2}}\,{{\text{m}}^{2}}\cdot {{\text{s}}^{-1}}$ and surface impedance to the range of ${{10}^{2}}$ to ${{10}^{4}}\,\text{s}\cdot {{\text{m}}^{-1}}$; second, spatial smoothness initialization, where the weights of the decoder convolution kernels are initialized to a Gaussian distribution. The spatial smoothing property of the Gaussian kernel guides the spatial continuity of the parameter field, which is consistent with the spatial gradient characteristics of thermodynamic parameters.

The preprocessing of high-dimensional climate data focuses on ensuring thermodynamic consistency, and a full process handling scheme is constructed, including "spatiotemporal matching - missing value filling - physical consistency normalization." Spatiotemporal matching uses bilinear interpolation to resample climate data from different sources and resolutions to a grid scale consistent with the physics simulator, aligning the time dimension to the same time step. Missing value filling discards traditional interpolation methods and adopts a generative filling strategy based on thermodynamic equations. Using the temperature-pressure static balance relationship and the water vapor pressure-temperature saturation relationship as constraints, missing data is filled through least squares optimization, ensuring the filled data satisfies basic thermodynamic associations. The normalization process adopts a variable-specific physical normalization method. For variables with clear physical baselines such as temperature and pressure, normalization is performed based on standard atmospheric state parameters. Specifically, ${{T}_{\text{norm}}}=\frac{T-{{T}_{0}}}{{{T}_{\text{max}}}-{{T}_{0}}}$, where ${{T}_{0}}$ is the standard atmospheric sea-level temperature, and ${{T}_{\text{max}}}$ is the maximum temperature in the study region. This avoids the destruction of thermodynamic variable correlations caused by traditional normalization methods. Feature fusion uses a cross-modal attention fusion module, which dynamically allocates the feature contribution of different variables by calculating the correlation weight between each climate variable and the target parameterization term. Specifically, a variable weight vector $\omega =\left[ {{\omega }_{1}},{{\omega }_{2}},\ldots ,{{\omega }_{n}} \right]$ is learned, where $\mathop{\sum }_{k=1}^{n}{{\omega }_{k}}=1$, to achieve weighted fusion of multivariable features. At the same time, a physical feature extraction submodule is embedded to specifically extract derived physical features, such as wind shear and temperature gradient, which are directly related to turbulence exchange. Feature concatenation enhances the model’s ability to identify key thermodynamic driving factors, providing high-quality feature inputs for the precise learning of parameterization terms.

Figure 3. Detailed structure diagram of the spatiotemporal neural network encoder-decoder

2.4 Multi-objective progressive learning strategy

The multi-objective progressive learning adopts a two-phase progressive process: "supervised pretraining - physics-enhanced training." The core innovation lies in achieving an ordered optimization from "accuracy-first" to "balancing physical consistency and mechanism discovery" through phased target decomposition. The first phase is supervised pretraining, where the core loss function is the MSE between the outputs of the physics simulator and the observed energy flow values, expressed as:

$L_{M S E}=\frac{1}{N} \sum_{i=1}^N\left\|\widehat{Q}_i-Q_i^{o b s}\right\|^2$    (4)

where, $N$ is the number of observation samples, ${{\hat{Q}}_{i}}$ is the energy flow value output by the simulator, and $Q_{i}^{\text{obs}}$ is the corresponding observed value. The AdamW optimizer is selected, with an initial learning rate of ${{10}^{-4}}$ and a weight decay coefficient of ${{10}^{-5}}$. The training termination condition is when the MSE on the validation set decreases by less than ${{10}^{-5}}$ over 10 consecutive epochs. The core goal is to guide the neural parameterizer through data supervision to generate a parameter field that roughly aligns with real-world observations, laying the foundation for subsequent physics-enhanced training. The second phase is physics-enhanced training, where three physical constraint losses are introduced to build the total loss function on top of the MSE loss:

$L_{{total }}=L_{M S E}+\lambda_1 L_{p h y}+\lambda_2 L_{ {sparse }}+\lambda_3 L_{ {disc }}$    (5)

where, ${{\lambda }_{1}},{{\lambda }_{2}},{{\lambda }_{3}}$ are loss weights. The physical residual regularization loss ${{L}_{\text{phy}}}$ combines energy conservation and momentum conservation constraints. The energy conservation residual is obtained by calculating the deviation of the simulated energy flow field from the energy balance equation:

$L_{{energy }}=\frac{1}{M} \sum_{j=1}^M\left|R_{n, j}-\widehat{H}_j-\widehat{L} E_j-\widetilde{G}_j\right|$    (6)

where, $M$ is the number of grid units. The momentum conservation residual is calculated by the deviation of the simulated wind field from the geostrophic balance relationship:

$L_{ {momentum }}=\frac{1}{M} \sum_{j=1}^M\left|f \vec{u}_j-\nabla p_j\right|$    (7)

where, $f$ is the Coriolis parameter, $\vec{u}_j$ is the wind vector, and ${{p}_{j}}$ is the air pressure. Finally, ${{L}_{\text{phy}}}={{L}_{\text{energy}}}+{{L}_{\text{momentum}}}$. The sparsity penalty loss is expressed as:

$L_{{sparse}}=\frac{1}{M} \sum_{j=1}^M\left|K_{\mathrm{h}, j}\right|$    (8)

This loss uses L1 regularization to suppress redundant parameter fluctuations and encourage spatial simplicity in the parameter field, consistent with the smooth spatial gradient characteristics of parameter distribution in thermodynamic processes. The discovery loss ${{L}_{\text{disc}}}$ aims to guide the model to learn an analytically solvable parameterization form. A power function based on wind speed and stability, ${{K}_{h}}=a{{u}^{b}}R_{i}^{c}$, is selected as a simplified analytical function, where $a,b,c$ are the parameters to be fitted, $u$ is wind speed, and ${{R}_{i}}$ is the Richardson number. The loss value is the MSE between the analytical function output and the neural parameterizer output:

$L_{d i s c}=\frac{1}{M} \sum_{j=1}^M\left\|K_{\mathrm{h}, j}^{n e t}-a u_j^b R i_j^c\right\|_2^2$    (9)

To ensure the convergence of the two-stage training process, while avoiding over-suppression of the data-driven adaptive capability by physical constraints, a multidimensional stability guarantee mechanism is designed. The core innovation lies in dynamically adjusting the loss weights and fine-tuning the gradient control. The loss weights follow a progressive incremental strategy: the physical residual regularization loss weight ${{\lambda }_{1}}$ increases linearly from 0.1 at the end of pretraining to 0.5, the sparsity penalty loss weight ${{\lambda }_{2}}$ is fixed at 0.01, and the discovery loss weight ${{\lambda }_{3}}$ increases from 0.05 to 0.2. This setup ensures the model first establishes data fitting capability and gradually strengthens physical consistency and mechanism discovery objectives. Gradient control uses gradient clipping, setting a global gradient norm threshold of 1.0. Gradients exceeding this threshold are clipped to prevent gradient explosion during training, ensuring smooth parameter updates. Learning rate adjustment follows a cosine annealing decay strategy, maintaining the initial learning rate of ${{10}^{-4}}$ during pretraining, and decaying it following a cosine curve during physics-enhanced training, with a minimum value of ${{10}^{-6}}$. This guarantees rapid convergence in the early stage and stability in parameter fine-tuning later. Additionally, after each round of parameter updates, a physical validity check step is added. If the parameter field exceeds the pre-set thermodynamic validity range, the corresponding gradients are decay-corrected to avoid distortion of the physical meaning of parameterization terms, ensuring the model always converges in a physically valid direction during training.

2.5 Closed-loop integrated prediction and uncertainty quantification module

The core logic of the closed-loop integrated prediction is to construct a “prediction-feedback-update” autoregressive operation loop to achieve continuous scenario simulation of short-term thermodynamic energy flow evolution. In the initialization phase, the initial surface temperature, radiative forcing, and other boundary conditions, along with high-dimensional climate data, are input into the model. The neural parameterizer generates the initial parameter field, and the physical simulator outputs the energy flow field and surface state variables at time $t$. When transitioning to the simulation at time $t+1$, the energy flow field and surface temperature at time $t$ are used as new boundary conditions, and the real-time high-dimensional climate data for $t+1$ are input into the model. The neural parameterizer then adapts the parameter field based on the updated inputs, and the physical simulator performs a new round of energy balance calculations under the new parameters and boundary constraints. This closed-loop design dynamically updates input boundary conditions and parameter fields, allowing the simulation process to adaptively capture the spatiotemporal evolution patterns of energy flow. To ensure thermodynamic consistency in the autoregressive process, a dual-check procedure is performed after each simulation step: the energy balance equation residual $\mid {{R}_{n}}-H-LE-G\mid$ is computed. If the residual exceeds the preset threshold of $10\,\text{W}\cdot {{\text{m}}^{-2}}$, the output parameters of the neural parameterizer are recalculated by backtracking and adjusting; simultaneously, the value range and spatial smoothness of the parameter field are checked to ensure the physical validity of the parameterization terms does not distort during the iterative process, thus maintaining the physical validity of continuous simulations.

Uncertainty quantification is implemented using the Monte Carlo Dropout technique. The core innovation is to deeply bind uncertainty quantification with the parameterization process, enabling accurate characterization and spatial tracing of cognitive uncertainty. The specific configuration is as follows: Dropout layers are embedded in the key layers of the encoder and decoder of the spatiotemporal neural parameterizer, with a Dropout rate set to 0.15. During the prediction phase, the Dropout layers are kept active, and 50 independent forward propagations are performed for the same input sample, generating 50 independent parameterization fields and energy flow prediction results. Cognitive uncertainty is quantified through the statistical dispersion characteristics of the parameterization fields. For each grid cell, the uncertainty of the turbulence exchange coefficient is represented by the standard deviation of the 50 sampling results:

$\sigma_{K_{\mathrm{h}}}(i, j)=\sqrt{\frac{1}{s-1} \sum_{s=1}^S\left(K_{\mathrm{h}}^s(i, j)-\bar{K}_{\mathrm{h}}(i, j)\right)^2}$    (10)

where, $S=50$ is the number of samples, $K_{h}^{s}\left( i,j \right)$ is the turbulence exchange coefficient for the $s$-th sample at grid point (i, j), and $\bar{K}_h(i, j)$ is the mean of the 50 samples. Uncertainty tracing is achieved by mapping ${{\text{ }\!\!\sigma\!\!\text{ }}_{{{\text{K}}_{\text{h}}}}}\text{(i, }\!\!~\!\!\text{ j)}$ to the spatial grid of the physical simulator, generating an uncertainty spatial distribution map. High standard deviation regions are identified as "uncertainty hotspot areas." Further analysis is performed to establish the intrinsic relationship between these hotspot regions and thermodynamic features: the proportion of different surface types in the hotspot regions is statistically calculated to determine the contribution of complex surfaces to uncertainty; the wind speed gradient and Richardson number distribution in the hotspot regions are computed to reveal the positive correlation between turbulence exchange mechanisms and uncertainty. This provides data support for locating fuzzy regions in thermodynamic mechanisms and focusing subsequent research efforts.

2.6 Model output definition

The model adopts a dual-output design, with the core advantage being the simultaneous consideration of prediction accuracy and physical interpretability, enabling the coordinated output of “prediction results - physical mechanisms.” The first type of output is a high-precision energy flow spatial distribution prediction map, covering the spatiotemporal distribution of sensible heat flux, latent heat flux, and surface temperature. The output results undergo thermodynamic consistency checks and strictly satisfy the energy balance constraints, and can be directly used for regional energy flow evolution scenario analysis and energy flow forecasting under extreme weather conditions. The second type of output is region-specific parameterization results with clear physical meaning, along with corresponding uncertainty estimates, including the turbulence exchange coefficient field, surface impedance field, and their spatial uncertainty distribution maps. The two types of outputs are tightly thermodynamically linked: the parameterization results are the core physical foundation for energy flow prediction, and their spatiotemporal heterogeneity directly determines the energy flow transport efficiency and distribution pattern. The accuracy of energy flow predictions, in turn, verifies the rationality of the parameterization mechanism. In terms of application value, the energy flow spatial distribution predictions provide high-precision data support for climate simulation and ecosystem assessment; the parameterization results and uncertainty estimates provide direct observational calibration for improving traditional regional thermodynamic parameterization schemes. Particularly, the mechanism analysis of uncertainty hotspot areas can guide targeted field observation experiment design and promote the iterative optimization of regional thermodynamic models.

4.1 Energy flow prediction accuracy verification

The spatiotemporal scale verification and extreme scenario verification results for energy flow prediction accuracy are shown in Tables 1 and 2. The proposed model demonstrates significant advantages under different spatiotemporal scales and extreme conditions, with the core advantage stemming from the deep integration of physical consistency and data-driven adaptability.

Table 1. Comparison of energy flow prediction accuracy of different models under different spatiotemporal scales and surface types (Units: MSE: W²·m⁻⁴; MAE: W·m⁻²; R: dimensionless)

Table 2. Comparison of energy flow prediction accuracy of different models under extreme weather scenarios (Units: MSE: W²·m⁻⁴; MAE: W·m⁻²; R: dimensionless)

From Table 1, it can be seen that on the daily scale, the proposed model reduces the MSE of sensible heat flux by 38.7% compared to the traditional PGML and by 30.8% compared to CNN-LSTM for homogeneous farmland. For the heterogeneous urban-rural fringe, the MSE of sensible heat flux is reduced by 40.6% compared to traditional PGML and by 29.7% compared to CNN-LSTM. Although the weekly scale prediction accuracy slightly decreases, the advantage remains stable. Spatially, the accuracy improvement for heterogeneous surfaces is significantly higher than for homogeneous surfaces, indicating that the neural parameterizer effectively captures the heterogeneity of the underlying surface. Regarding energy conservation error, the proposed model consistently stays below 5.7 W·m⁻², far better than the 8.6~10.5 W·m⁻² of traditional PGML and 11.2~14.3 W·m⁻² of CNN-LSTM, which verifies the role of embedded physical consistency in ensuring the reliability of predictions.

The extreme scenario validation results in Table 2 show that the accuracy degradation rate of the proposed model under high temperature and strong wind conditions is only 13.5%~17.3%, far lower than the 28.4%~33.7% for traditional PGML and 33.1%~38.2% for CNN-LSTM. In the high-temperature scenario, traditional PGML fails to adapt to the strong convection conditions due to fixed turbulence exchange coefficients, leading to a sharp increase in MSE for sensible heat flux prediction. CNN-LSTM, lacking physical constraints, exhibits numerical extrapolation bias. In contrast, the proposed model dynamically adjusts the turbulence exchange coefficients using the neural parameterizer and maintains energy balance under physical conservation constraints, effectively controlling prediction bias. In the strong wind scenario, the strict adherence to momentum conservation further reduces the deviation in energy flow simulation, verifying the core role of embedded physical consistency in enhancing model generalization capability.

4.2 Innovation component effectiveness analysis

The effectiveness verification results of the innovative components based on the ablation experiments are shown in Table 3. The synergistic effect of the neural parameterizer, multi-objective learning strategy, and uncertainty quantification module is the key to ensuring model performance. The independent contributions of each component are significant and complementary.

The effectiveness of the neural parameterizer is significant: when removed, the MSE of sensible heat flux increases by 75.7%, the MSE of latent heat flux increases by 69.3%, and the energy conservation error doubles, with the parameter reasonability coefficient dropping from 0.87 to 0.65. This indicates that fixed empirical parameters cannot adapt to the spatiotemporal heterogeneity of complex surfaces, whereas the neural parameterizer, driven by high-dimensional climate data, can precisely generate parameter fields that match regional features, providing core support for high-precision simulations. In the multi-objective learning strategy, removing physical residual regularization increases the energy conservation error from 4.8 to 8.5 W·m⁻², and the parameter reasonability coefficient drops by 0.11, verifying its key role in strengthening physical consistency. Removing the discovery loss decreases the parameter reasonability coefficient by 0.15, indicating that this loss term effectively guides the model to learn a simplified and physically intuitive parameterization form. The effectiveness of the uncertainty quantification module is verified by the confidence interval coverage, with the model achieving a coverage rate of nearly 94%, far exceeding the approximately 82% coverage rate of the model without the neural parameterizer. This demonstrates its ability to accurately depict cognitive uncertainty. The CRPS metric shows that the uncertainty quantification accuracy of the model is over 35% higher than that of the model without the neural parameterizer, and the uncertainty hotspot areas are primarily concentrated in regions with complex underlying surfaces such as urban-rural fringe zones and mountain-plain transitions, as well as areas with sparse observational data, which aligns highly with regional thermodynamic characteristics, achieving precise association between uncertainty and mechanisms.

Table 3. Impact of each component on model performance in ablation experiments (Units: MSE: W²·m⁻⁴; MAE: W·m⁻²; R: dimensionless; energy conservation error: W·m⁻²; Confidence interval coverage: %; CRPS: W·m⁻²)

4.3 Discovery results of regional thermodynamic parameterization mechanism

The discovery results of the regional thermodynamic parameterization mechanism are the core innovation of this study, including the spatial distribution characteristics of parameterization terms, the derivable parameterization relationships, and the core thermodynamic driving mechanisms. The related results are shown in Table 4 and Table 5. The corresponding Figure 4 displays the regional turbulence exchange coefficient field ${{K}_{h}}$, surface impedance field ${{r}_{s}}$, and the corresponding uncertainty spatial distribution map generated by the model inversion. It highlights the parameter gradient changes in heterogeneous regions such as the urban-rural fringe.

Table 4 shows that the spatial distribution characteristics of parameterization terms are significant and consistent with thermodynamic principles: The turbulence exchange coefficient ${{K}_{h}}$ follows the distribution pattern of waterbody > bare land > forest > farmland > urban built-up area, which is negatively correlated with surface roughness—water bodies have smooth surfaces and low momentum exchange resistance, yielding the largest ${{K}_{h}}$; urban areas with dense buildings have high roughness, resulting in the smallest ${{K}_{h}}$. The surface impedance ${{r}_{s}}$, in contrast, shows an opposite trend, with urban built-up areas having the largest ${{r}_{s}}$ due to low vegetation cover and weak evaporation, while water bodies have the smallest ${{r}_{s}}$. Correlation analysis shows that ${{K}_{h}}$ is significantly positively correlated with both wind speed and temperature gradient, with coefficients higher than 0.85 for bare land and water, indicating that turbulence exchange in these regions is primarily driven by wind speed. For farmland and forest, the correlation with temperature gradient is higher, reflecting the important contribution of thermal convection to turbulence exchange. These distribution characteristics are consistent with the physical understanding of regional thermodynamic processes.

Figure 4. Regional spatial distribution heatmap

Table 4. Statistical results of turbulence exchange coefficient and surface impedance for different surface types (Units: ${{K}_{h}}$: m²·s⁻¹; ${{r}_{s}}$: s·m⁻¹)

Table 5. Comparison of fit performance between newly discovered parameterization relationships and existing empirical formulas

Figure 5. Fit comparison of newly discovered parameterization relationship with existing empirical formulas

Table 5 shows that the newly discovered parameterization relationship has significant advantages: its fit quality ${{R}^{2}}$ reaches 0.89, improving by 23.6% compared to the traditional Businger-Dyer formula and by 14.1% compared to the improved Zhang formula. The RMSE is only 0.078 m²·s⁻¹, more than 50% lower than the traditional formulas. Traditional formulas are mostly suitable for single homogeneous surfaces, while the new relationship, constrained by discovery loss, integrates the thermodynamic characteristics of multiple surface types and can adapt to the turbulence exchange laws of different regions. Combined with energy flow prediction results and parameter distribution characteristics, the core driving mechanisms of regional energy flow evolution can be clearly identified: In the plain farmland region, the turbulence exchange coefficient has a high correlation with the temperature gradient, and energy flow transmission is dominated by thermal convection, with latent heat flux accounting for more than 60%. In the urban built-up area, the large surface impedance and weak turbulence exchange lead to a significant increase in sensible heat flux, which accounts for 55%~65%, forming an "urban heat island" related energy flow distribution pattern. In the mountain-plain transition zone, the surface heterogeneity causes dramatic fluctuations in the turbulence exchange coefficient, making it a hotspot area for uncertainty in energy flow prediction. In this region, energy flow transmission is jointly dominated by terrain dynamics and thermal effects. Corresponding to Table 5, Figure 5 displays the comparison of fit performance between the newly discovered parameterization relationship, the traditional Businger-Dyer formula, and observational data. The X-axis represents the physical driving variables, and the Y-axis represents the turbulence exchange coefficient, showing the data point distribution and the fitting curves of different formulas.

This study proposed a dual-path physics-data hybrid framework that achieves deep synergistic integration of high-dimensional climate data and thermodynamic physical laws. The core innovation lies in constructing a collaborative architecture of differentiable physical kernels and spatiotemporal neural parameterizers, which upgrades physical laws from traditional external fixed constraints to an internally discoverable and interpretable structure within the model. This effectively solved the challenges of unknown thermodynamic parameterizations and spatial heterogeneity adaptation under complex underlying surface conditions. The integration of multi-objective progressive learning strategies and closed-loop uncertainty quantification modules further strengthened the model's physical consistency and mechanism discovery capabilities, forming a complete technological system of "physical skeleton-data flesh-mechanism discovery."

Experimental results show that this model outperforms traditional PGML and CNN-LSTM in both spatiotemporal scale and extreme scenario validations. The MSE of sensible heat flux in heterogeneous urban-rural transition zones at the daily scale is reduced by 40.6% compared to traditional PGML. The advantage remains stable at the weekly scale, and the accuracy degradation rate under high temperature and strong wind conditions is far lower than that of the comparison models. Ablation experiments show that after removing the neural parameterizer, the MSE of sensible heat flux increases by 75.7%. The physical residual regularization ensures energy conservation, and the collaborative components enhance model performance. The discovery of regional thermodynamic mechanisms shows that the turbulence exchange coefficient follows the order: water > bare land > forest > farmland > urban built-up area. The newly discovered parameterization relationship has a goodness of fit (R²) of 0.89, significantly improving compared to traditional formulas, and clarifies the regional differences in energy flow drivers, such as thermal convection dominance in plain farmland and sensible heat dominance in urban built-up areas.

The academic contribution of this study lies in providing a new method for regional thermodynamic energy flow pattern prediction, advancing the paradigm shift in machine learning from alternative prediction to mechanism discovery and model enhancement in the field of thermodynamics. The practical value is reflected in how the discovered regional parameterization mechanism provides direct data support and theoretical basis for improving traditional climate model parameterization schemes. This is of great significance for enhancing climate simulation reliability, optimizing extreme weather warnings, and ecological environment assessments. Future research can further expand the model's scale adaptability, improve parameterization systems under complex terrain conditions, and promote its application in thermodynamic system simulations across larger areas.