Investigation of Land Atmosphere Coupling during the Extreme Rainstorm on 20 July 2021 over Central East China

1. Introduction

Land atmosphere coupling in weather and climate is important but inherent challenging and complex [1,2]. The soil moisture (SM) as one of the key factors in land atmosphere coupling is especially evidential for affecting rainfall during both model [3] and observational [4] studies. And these evidences have been concluded through various typical ways (or feedback chains) but with broad uncertainties because of coupling advantages and/or deficits [5]. Therefore, quantifying these varied land-atmosphere coupling characteristics during local typical rainstorms can be of great significance for enhancing the insights into local typical feedback chains.

As indicated by well documented global observation network, early studies have shown that rainfall has been highly coupled with SM in many “hot spot” (or typical) regions [6], where local morning SM can significantly affect afternoon convection [7]. That coupling evidence has highlighted the relation chain base of the Global Energy and Water Cycle Experiment (GEWEX), which promotes the completely interpreted relations among land surface, planetary boundary layer (PBL) (e.g., clouds or convection), and rainfall. However, under the background of shifted local surface caused by global warming [8], these rainfall-related couplings have been found to be complexly intersected with various local characteristics such as soil moisture gradient [9], and soil thickness [10]. Nevertheless, many model-relative studies have shown that typical regional couplings can be both negative and positive, indicated by the complexly intersected characteristics among soil moisture, land cover, land surface model (LSM), resolved scale of coupling model, and so on [11,12,13,14,15,16]. Generally, these complex local land atmosphere couplings (LoCo) have been characterized by immeasurable uncertainties in observations or model deficits to some extent.

In order to quantify various typical LoCo characteristics rooted in observable process-level scales, the variations of surface soil moisture, evapotranspiration, PBL states, PBL top entrainment, near-surface atmospheric thermal conditions, precipitation and/or clouds resulting from surface changes have been investigated through energy budget analysis and/or state covariations to evaluate various typical surface advantages and/or deficits, e.g., morning PBL moistening or deepening [17,18], locally and non locally physical items resulted into PBL top humidity changes [19,20], diagnostic states regarding well-mixed turbulence [21,22,23], surface and atmospheric flux contributions under surface thermal and moisture coevolution space [24], atmospheric state responses to SM [25,26,27], the persistence of SM anomalies [28], and so on. And these quite different measurements have been organized into one overall LoCo process chain, and further collected into the Coupling Metrics Toolkit (www.coupling-metrics.com), in order to foster grassroots participation [29].

LoCo metrics have been applied for regional evaluation of both models and datasets recently, regarding the coupling characteristics of their observational physics or statistics [31,32,33]. However, their performances are likely semi-coupling and require point samples of long-time series. Meanwhile, the uniformly scaled metrics of fixed space could be challenging because of insufficient PBL observations and/or coupling theory deficits (e.g., the varied scale dependence) [6,7,8,24]. While different from that, ensemble simulations within perturbed SM can achieve the overall relation chains (e.g., correlations) within regional studies [34,35], but usually non-observational and qualitative (e.g., negative and/or positive). Therefore, these various coupling indicators have been limited for broader applications by scale dependence and observational sample deficits.

In general, the local couplings indicated by observation-based metrics and ensemble statistical correlations lack of mutual confirmation mainly due to the unsolved spatiotemporal scale issues. Especially, the coupling characteristics indicated by ensemble statistical correlations lack verification of their reliability. To fill this gap, this study has evaluated the LoCo characteristics of the recent extreme rainstorm occurs on 20 July 2021 over central east China [36,37] on the resolved model grids by using various LoCo metrics and ensemble statistical relations, intending to highlight the potential perspectives of LoCo metrics on the local weather application.

2. Model and Data

This study has taken the Unified Noah Land Surface Model (Noah LSM) [38] coupled with the advanced weather research and forecasting model (WRF, Version 3.9.1) [39] to conduct the land-atmosphere coupling simulation of the typical rainstorm on 20 July 2021. The Noah LSM can run on finer temporal and spatial steps, and it has been widely used in numerical weather prediction and well-verified in local coupling studies [24,30,31,32].

As seen in Figure 1, two domains (e.g., D01 and D02) are one-way nested and centered at 113.45 ^oE,33.85 ^oN with 51 model levels topped at 50 hPa and run within a 60-second time step. The outer and inner domain resolutions are 12 and 4 Km with 100×100 and 165×159 grids respectively. During both domains, the microphysics chose the Thompson scheme [40], the short and long wave physics chose the RRTMG schemes [41], the cumulus parameterization physics chose the Kain-Fritsch scheme [42], and the surface layer physics chose the Revised MM5 scheme [43]. Especially, the Yonsei University PBL scheme (YSU) has adopted a non-local entrainment treatment method more suitable for nocturnal PBL within lower heights [44]. The lower 21 levels below 1 Km are conducted to achieve finer PBL profiles within YSU physics.

The atmospheric forcing datasets are collected from the first generation global atmospheric/land-surface reanalysis project (CRA40, 1978-2018), which has a resolution of 34 Km, 64 pressure levels, and 6 and 3 hours intervals for the atmospheric and surface layer respectively [45,46]. Also, the Noah LSM driven by the Global Land Data Assimilation System reanalysis (GLDAS) has a resolution of 0.25^o, 3-hour intervals, and 4 layers [47]. In addition, the 3 s terrain derived from the United State Geological Survey (USGS) and the 15 s land cover derived from the Moderate resolution Imaging Spectro radiometer (MODIS) product have been taken as the static underlying datasets.

The observations used in this study mainly include the automatic weather stations (AWS), nine S-band Doppler weather radars, and two soundings (Figure 1) derived from CMA Henan Meteorological Bureau. The AWS has a resolution of around 0.1^o and 1-hour interval and each radar are observed at a resolution of around 250 m, 11 elevation numbers, and 6-minute interval. The sounding sites Zhengzhou (ZZ) and Nanyang (NY) are taken as the two typical points of this domain to investigate the rainstorm-affected typical factors and non-typical factors respectively. Also, the grid precipitation product of the Land Surface Data Assimilation System (CLDAS) grid precipitation product of CMA, which has a resolution of 0.1^o×0.1^o and 1 hour [48] during the local operational application, is collected to fulfill the accumulative rainfall validation. Noted that the local hourly CLDAS rainfall products are used for calculating the longer time series accumulated precipitation (e.g., 12 hour rainfall) in order to avoid the discontinuity in AWS observations at different times. Also, the AWS datasets are supplied alongside to ensure the objectivity of instantaneous observations.

3. Experiments

The flow chart of this study is shown in Figure 2. Firstly, all the datasets used for various coupling metric calculation are prepared as follows, the simulation driven by CMA datasets is compared with various local observations to ensure one basis of this study that the main land-atmosphere state characteristics have been well documented in this modeled rainstorm (See Section 5.1), base on which, the ensemble forecast within perturbated soil moisture has been conducted to achieve another basis of this study that the reliable soil ensemble forecasts. Secondly, the gridded land-PBL profiles derived from the well-modeled rainstorm and the ensemble states within soil perturbations derived from the soil ensemble forecast are collected for the calculation of LoCo metrics and ensemble statistical relations to achieve the LoCo diagnoses and the paired state relations respectively. Finally, the ongoing mentioned statistics are analyzed on their main spatial and temporal characteristics to achieve the overall quantitative coupling insights during this event.

Noted that the simulation period covers from 00:00 to 24:00 on 20 July (UTC), which includes the extreme rainfall event of ZZ city around late afternoon on 20 July. The inner domain (Figure 1; D02) is taken as the target region for land-atmosphere coupling investigation. Three dry (wet) members assembled by systematically reduced (added) soil moisture at an interval amplitude with the self magnitude of 30% for the target region are conducted to generate soil ensemble forecast [35]. And the land atmospheric coupling chains indicated by various metrics have been compared to investigate the main coupling characteristics on their applicability in local weather forecasts.

4. Methods

4.1. The local coupling metrics

Table 1 has shown various LoCo diagnoses during the CoMeT (Table 1). For various surface advantages or deficits for potential convection indicated by CHF are classified by the observation-based thresholds suggested by the former study (detail described in Section 5.2.1) [17,18]. For various deficits and advantages indicated by RHT [19,20], HCF [21,22,23], and MDT [24] are measured within hourly steps. Meanwhile, we take four 6-hour ranged periods to conduct both the TCP [25,26,27] and SMM [28], the diagnoses of which have been calculated from two adjacent 3-hour forecast during one period. Both sensible and latent heat fluxes for TCP (e.g., TCP_SH and TCP_LH) are considered, and the threshold for memory loss in SMM is set to ${\mathrm{e}}^{-1}$ (about 0.5). Noted that we take a 3-hour interval simulation of 20 July 2021 to assemble the necessary samples for calculating the LoCo diagnoses, and further analysis on both local daytime (00:00~12:00) and nighttime (12:00~24:00) are conducted.

The typical LoCo chains can be completed by using predefined variation $\mathsf{\Delta }$ within temporally and/or spatially issued diagnoses to hopefully address their closure [24]. Noted that variations $\mathsf{\Delta }$ resulting from the surface (see Table 1) have been implicitly defined as the spatial difference of an instantaneous moment (e.g., surface-upper layer height differences in CHF or HCF), or both temporal and spatial differences (e.g., MDT, RHT, TCP and SMM). Especially, the variation $\mathsf{\Delta }$ of surface has been defined in the form of time lags for TCP and SMM.

4.2. The ensemble statistical metrics

This study conducts ensemble forecasts with dry or wet perturbations (DP or WP) by utilizing the initially increased or decreased soil moisture ($∆\mathrm{S}\mathrm{M}$) of one well-documented land-atmosphere simulation (or control run; $\mathrm{c}\mathrm{t}\mathrm{l}$). Following, the differences between DP (WP) and control run are defined as variations $∆\mathrm{D}$($∆\mathrm{W}$) that resulted from $∆\mathrm{S}\mathrm{M}$. Then, the temporal correlations (${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$) among various $∆\mathrm{D}$($∆\mathrm{W}$) are analyzed to identify ensemble statistical relations. The formulas for dry and wet variations are calculated as follows.

$\overline{\mathrm{D}}=\frac{1}{\mathrm{N}}{\sum }_{\mathrm{k}=1}^{\mathrm{k}=\mathrm{N}} {\mathrm{D}}_{\mathrm{k}}$$,\overline{\mathrm{W}}=\frac{1}{\mathrm{N}}{\sum }_{\mathrm{k}=1}^{\mathrm{k}=\mathrm{N}} {\mathrm{W}}_{\mathrm{k}}$	(1)
$∆\mathrm{D}=\frac{1}{\mathrm{N}}{\sum }_{\mathrm{k}=1}^{\mathrm{k}=\mathrm{N}} \left({\mathrm{D}}_{\mathrm{k}}-\mathrm{c}\mathrm{t}\mathrm{l}\right)$$,∆\mathrm{W}=\frac{1}{\mathrm{N}}{\sum }_{\mathrm{k}=1}^{\mathrm{k}=\mathrm{N}} \left({\mathrm{W}}_{\mathrm{k}}-\mathrm{c}\mathrm{t}\mathrm{l}\right)$	(2)

In equation (1) and (2), $\mathrm{k}$ and $\mathrm{N}$ represent the ensemble members and total number ($\mathrm{N}=3$), respectively, $\overline{\mathrm{D}}$ and $\overline{\mathrm{W}}$ are the averaged state at one specific layer during the study area for the dry and wet states (${\mathrm{D}}_{\mathrm{k}}$ and ${\mathrm{W}}_{\mathrm{k}}$) respectively, and $\mathrm{c}\mathrm{t}\mathrm{l}$ represents the control state (e.g., the well-documented modeled rainstorm). Noted that the state average is usually desired to remove the unrealistic land-atmosphere response at the beginning [35].

Additionally, the relation strength based on the Spearman rank correlation [49] among the averaged states ($\overline{\mathrm{D}}$ or $\overline{\mathrm{W}}$) or variations ($∆\mathrm{D}$ or $∆\mathrm{W}$) of different layers are calculated as follows.

$\mathrm{I}={\mathrm{C}\mathrm{C}}_{\mathrm{T}}\left(\mathrm{a}\mathrm{b}\mathrm{s}\right({\mathrm{C}\mathrm{C}}_{\mathrm{T}})>0.5),$ ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}=\frac{{\sum }_{\mathrm{i}} \left({\mathrm{x}}_{\mathrm{i}}-\stackrel{\mathrm{‾}}{\mathrm{x}}\right)\left({\mathrm{y}}_{\mathrm{i}}-\stackrel{\mathrm{‾}}{\mathrm{y}}\right)}{\sqrt{{\sum }_{\mathrm{i}} {\left({\mathrm{x}}_{\mathrm{i}}-\stackrel{\mathrm{‾}}{\mathrm{x}}\right)}^2{\sum }_{\mathrm{i}} {\left({\mathrm{y}}_{\mathrm{i}}-\stackrel{\mathrm{‾}}{\mathrm{y}}\right)}^2}},\mathrm{i}\in \mathrm{T}$

(3)

In equation (3), $\mathrm{i}$ represents for the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ index of the overall time steps $\mathrm{T}$, $\mathrm{x}$ and $\mathrm{y}$ represent for two sorted state such as PBL thermal or energy variables, and the superscript horizontal line represents for temporal average. Therefore, the coupling intensity ($\mathrm{I}$) is identified when ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ value is larger than 0.5 or less than -0.5, and its reliability is indicated by the significant level (p) when p is less than a threshold (e.g., 0.05).

In addition, it should be noted that in order to achieve the main spatial characteristics of the gridded surface-PBL profiles, LoCo diagnoses and ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of ensemble statistics, and these three kinds of samples (number of which is equal to D02 domain gird number as 165×159) have been resampled into 100 bins and then fitted by a 2-peak Gaussian function, intending to achieve at most two main signals on the model grid resolved couplings. The Gaussian function is defined as $\mathrm{f}\left(\mathrm{s}\right)={\sum }_{\mathrm{i}=1}^{\mathrm{n}} {\mathrm{a}}_{\mathrm{i}}\mathrm{e}\left[-{\left(\frac{\mathrm{s}-{\mathrm{b}}_{\mathrm{i}}}{{\mathrm{c}}_{\mathrm{i}}}\right)}^2\right]$, where s is the sample, a is the amplitude, b is the centroid (location), c is related to the peak width, and n is the number of peaks to fit (n =2).

5. Results

5.1. Simulation overview

5.1.1. Synoptic and thermodynamics

As seen in Figure 3, for the 500 (hPa) level, a cyclonic Low-value (L) system crosses the whole central domain surrounded by three High-value (H) centers, while an intense Warm (W; around -4 ^oC) center located on the L’s west (Figure 3a). For the 700 level, the L center accompanied by an intense Cold (C; around 9 ^oC) center and strong easterly flows on its north edge is still holding but has clearly shrunk (Figure 3b). For the 850 level, the L center accompanied by strong northeast flows on its north edge is still holding but has further shrunk (Figure 3c). For the surface level, the L center accompanied by an intense C center on its northern edge is similar to the upper-level patterns (Figure 3d). Clearly, the deep low-value system holding for all levels is the primary synoptic pattern, and its range has shrunk with lowered levels, and this is especially pronounced between 500 and 700 level when the surrounding flows have clearly changed. This indicates a strong warm pool of 500 level but a strong cold invading below it. In addition, this low-value system over the relatively moist and quite warm surface ( > 25 ^oC) has lasted more than 12 hours (ending till 12:00 or later), this should be attributed to the upper systems easterly developing (e.g., the South Asia High) above the 300 hPa level (not shown).

As seen in Figure 4, at 00:00 for the ZZ site, the observed saturated atmosphere where the difference between temperature and dew point (T-T_d) is almost zero can be observed below 400 hPa, while the simulated saturated top can reach 200 hPa, and the wind rotates counterclockwise as the pressure decreases above 850 hPa for both observation and simulation (Figure 4a). At 12:00 for the ZZ site, the top of the observed saturated atmosphere reached to near 300 hPa, while the simulated saturated top decreased to 400 hPa, and the counterclockwise rotated winds between 850 and 500 hPa can be observed in both observations and simulation (Figure 4b). Clearly, for the ZZ site, below 400hPa, the regional thermodynamic conditions of this event have been well simulated, and the pressure and temperature of lifting condensation level for both simulations (sPlcl and sTlcl) and observation (oPlcl and oTlcl) are almost equitable. However, above 400 hPa for the ZZ site, the higher saturation top at 00:00 and lower saturation top at 12:00 of simulations have resulted in smaller and larger convective available potential energy (Cape) respectively, when compared to observationss. Furthermore, at 00:00 for the NY site, below about 400 hPa, both the observed and simulated temperature and dew point profiles have behaved consistently, where the counterclockwise rotated winds also can be observed, while above 400 hPa the saturated atmosphere of simulation is different from the observed dry atmosphere (Figure 4c). At 12:00 for the NY site, both the observed and simulated temperature profiles behaved consistently, and a shaking dew point profile of observation can be observed when compared to simulation, while a slightly moister atmosphere above 400 hPa in simulation can be observed when compared to observation. Also, the counterclockwise rotated winds can be observed between the surface and 500 hPa layers in both observations and simulations (Figure 4d). The low-level Plcl and Tlcl for both observation and simulation are almost equatable, while the Cape differences between simulation and observation for NY site behaviors are similar to the ZZ site.

Obviously, the mid-low layer (below 400 hPa) thermodynamic situations of both two sites have been well simulated, and the moist atmosphere holds on from 00:00 to 12:00 where the lifting condensation level (LCL) can be easily reached. The inverse wind directions indicate a large-scale cold invading, and it covers almost the whole layers at 00:00 and then shrinks to the mid-low layers at 12:00 for the ZZ site, while it has decreased from the mid-low layers at 00:00 to low layers at 12:00 for NY site. In contrast, the observed upper southern flows have been increased from 00:00 to 12:00 for both two sites, this indicates a northerly moved warm front. However, the simulation has a clear large bias on the upper-level (above 400hPa) systems. Since simulated sounding is directly initialized from the atmospheric and surface layers of CRA40 datasets, the upper layer difference can be directly attributed to the datasets’ making theory, and the upper-layer drift variation of observation and the imperfect upper-layer model physics should also not be neglectable.

5.1.2. Convection and rainfall

The observed composite reflectivity (CR) and 3-hour accumulative rainfall (${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) are compared with simulation to validate the modeled performance of convection (defined as when CR>10 dBZ) and rainfall. The diurnally observed convection is relatively stronger (> 26 dBZ) when compared to the simulation (> 24 dBZ), and this difference is especially pronounced at the beginning (e.g. from 03:00 to 06:00) (Figure 5a). And the dominant convection of observation is also stronger (around 28 dBZ), but with a lower frequency of around 2000, which is almost half of the simulation (Figure 5b). The ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ varied between 0.3 and 0.5 on 20 July. And the most pronounced ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ at 12:00 (when the extreme rainfall ends) can reach nearly 0.6. Meanwhile, RMSE has shown the opposite behavior of ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$. Especially, a significant decrease (from 17 to 13 dBZ) at the beginning and a significant increase (from 13 to 16 dBZ) at the end can be found (Figure 5c). The significant high ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ and low RMSE (when ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$>0.3 and RMSE<14 dBZ) indicate the main convection can been well captured from local noon to early night (06:00~15:00), while before and following that, the simulated convection is likely missed. Additionally, the daytime ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ has shown significantly negative values around the ongoing mentioned main rainfall area, with positive values in the peripheral area which can be observed (Figure 5d). Clearly, the beginning simulation deviation should due to the model warm-up deficits resulting from insufficient datasets, and this also contributes partly to the daytime ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ around the rainstorm area (e.g., the significant increase of CR in simulation), while the ending simulation deviation should mostly due to model physics deficits.

The domain averaged 3-hour rainfall (${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) of simulation has behaved a significant increase from local noon to early night, which is consistent with observations but mostly with lower values (Figure 6a). Nevertheless, significant rainfall biases can be found at the beginning and ending periods. Meanwhile, the frequency of ${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$ is lager in simulation, which is especially pronounced for light rainfall (e.g., ${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$< 20mm), though the frequency distributions of both simulation and observation behaved consistently (when ${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$ is larger than 10mm) (Figure 6b). Additionally, the daytime rainfall amount has a ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ of about 0.62, and the shapes of heavy rainfall areas are especially similar (Figure 6c and d). Also, the rainfall extremes (>200 mm) are both located at the ZZ city. This indicates the main rainfall spatial characteristics and location of extremes have been well simulated. However, the northeastern edge of this developing rainfall in simulation has been clearly missed when compared to observation. This could be related to the latter poorly simulated convection that developed in the north.

Overall, though this modeled convection and rainfall are relatively weaker than observations, both the main convection and rainfall characteristics have been well captured from local noon to early night. Recall the fact that during simulation, the long-lasting daytime low-value systems of mid-low synoptic, and the mid-low layer significant moistening and low layer strong cold invading of thermodynamics have shown quality promised consistency with observations. Therefore, the ongoing simulation advantages can be likely attributed to the fair quality of both datasets and model laws. In contrast, the simulated convection and/or rainfall deficits during the beginning and ending periods that are possibly related to the model spins-up/physical-dynamic law deficits or datasets deficits can be too difficult to conclude because of insufficient observations. Nevertheless, the largely biased northern flows of the upper layer above 400hPa can be directly related to the biased northerly developed convection and rainfall. And this can be likely attributed to the dataset’s deficits because of the sounding drifts at the stratosphere.

5.1.3. Underlying surface characteristics

The domain averaged ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ simulation has shown a consistent diurnal variation with observation, but is warmer than the latter, which is especially pronounced around local noon time (Figure 7a). However, the frequency of the temporally averaged observed ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ has been mostly distributed above the simulation (Figure 7b). Clearly, the warmer diurnal variation and colder frequency spatial distribution indicate the overall simulated ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ characteristics have fallen into two different phases. But the observation is below the simulation at around 302K with a frequency of around 5000, which is almost twice the former. Additionally, the ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ has a significant decrease at the beginning, then a relatively flat variation of around 0.2 from local morning to late afternoon, and a significant decrease at the ending period (Figure 7c). Meanwhile, RMSE has shown an opposite variation compared to ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$, and the significant RMSE around 8 K occurs at local noontime. Furthermore, the daytime ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ that exceeds 0.5 with confidence (whose significance level is lower than 0.05) can be observed in most areas, while significant negative relations can be observed over the Southern Taihang mountain (Figure 7d). Clearly, the significantly negatively biased daytime surface temperature over mountain areas has mainly contributed to the phase differential surface characteristics. However, since well-documented convection and rainfall can be achieved over the phase differential biased thermal surfaces, the mountains thermal effects likely contributed little to the daytime upper atmosphere.

In general, the long-lasting low-value system with the upper warm flow and lower strong cold invading, the mid-low layer thermodynamic situations, the convection and rainfall spatiotemporal characteristics, and the diurnal surface thermal characteristics are consistent with the available observations. However, the stratospheric (higher than 400 hPa) thermodynamics that related to the northern developed rainfall and convection, and the temperature during mountain areas have been found biased. Except for the mountain areas, the main characteristics during the low atmosphere and the surface can be well documented in this modeled event, therefore, this is taken as the base of the surface-atmosphere state for further investigation.

5.2. The local coupling evaluation

5.2.1. CHF

The thresholds of CHF for typical characteristic discrimination are defined as follows. For atmospherically controlled advantage (ACA), TCP varies between 0 and 300 (J/Kg) and HI_low varies between 0 and 5 (℃). For wet soil advantage (WSA), TCP varies between 0 and 230 (J/Kg) and HI_low varies between 5 and 10 (℃). For dry soil advantage (DSA), TCP varies between 180 and 300 (J/Kg) and HI_low varies between 10 and 15 (℃). For trans region (TR), TCP varies between 100 and 180 (J/Kg) and HI_low varies between 10 and 15 (℃). For stable regions when too dry for rainfall (SRD), HI_low is larger than 15 (℃), and for regions when too stable atmosphere for rainfall (SRC), CTP is smaller than 0 (J/Kg).

At 00:00 on 20 July (local morning time), most grid points have fallen into the ACA and WSA regions during the CHF, while fewer grid points have fallen into the TR, DSA, and SRC regions (Figure 8a). Especially, one SRC region lies in the northwest of HP which locates near the center of the upper low-value system (Figure 3a and b), the adjacent TR and DSA areas have located on the northwest corner of the D02 domain, two WSA areas have located in the northeast corner and southeast of D02 domain respectively, and most of D02 domain has been occupied by ACA (Figure 8b).

Overall, the CTP values of SRC are mostly negative, and HI_low values of SR are quite small, this indicates the atmosphere is likely too stable to rain near the storm center. Meanwhile, the moderate CTP and high HI_low values of the TR area indicate a neutral boundary layer and dry surface, and additional CTP is needed to meet with convection on the northwest of the D02 domain. Clearly, except SRC and TR, other regions are favorable for potential convection induced by atmospheric control (ACA), wet soil (WSA), and dry soil (DSA).

5.2.2. RHT

Noted that the relative humidity tendency (RHT) describes the RH tendency of PBL top that resulted from surface evaporation (SE), boundary layer growth (BLG), boundary layer warming (BLW), and dry air entrainment (DE), also the non-evaporation (NE) has been included for broader factor investigation. They are point wised indicators for local cloud development [20]. Therefore, the main characteristics of RHT have been analyzed by using temporally averaged diagnoses.

The PBL height of the ZZ site has varied consistently with RH at the PBL top, and a significant decrease and increase during the beginning and ending periods respectively can be observed (Figure 9a). Especially, the PBL height can decrease to near 0 at local noontime while RH at the PBL top is mostly larger than 0.8. However, the PBL height of the NY site has behaved oppositely with the RH at the PBL top, and significant daytime increase and decrease in PBL height and RH respectively can be observed while this behaves inversely at night(Figure 9b). Especially, growing PBL height and decreasing RH can reach around 1000m and 0.65 respectively during local noontime while significant moistening and decreasing PBL occurs in the local late afternoon. Additionally, the significant RHT tendency of the ZZ site at the local early noon time with positive BLW, EF, and BLG can be observed (Figure 9c). While besides the early noon, the positive RHT in the afternoon can be observed for the NY site (Figure 9d). Especially EF increases at early night for both ZZ and NY sites, while the intermittent RHT, BLG, BLW, and absent DE for both ZZ and NY sites have indicated that the PBL clouds developing resulting from the local surface are highly nonlinear.

RHT is mostly distributed around 4.0 with a significant frequency of around 3000 (Figure 10a and b). And SE is mostly distributed around 0.85 with a significant frequency of around 1500. Then NE is mostly distributed around 1.9 with a significant frequency of around 600. While DE is narrowly distributed around 0 with a significant frequency of around 1500. Especially, BLW is broadly distributed around 1.9 with a significant frequency of around 500. Nevertheless, BLG is narrowly distributed around 0.02 with a significant frequency of around 800. Clearly, the SE and BLW have contributed mostly to RHT, while NE has shown relative significance during this event. Additionally, the spatial distribution of NE has behaved consistently with BLW, and strong SE surrounding the main rainfall area, weaker BLG, and weakest (mostly negative) DE can be identified. It is obvious that except for the main rainfall area, EF and BLW have played the dominant roles in the local PBL cloud development, while the consistent BLW and NE distribution also indicates the significant environmental factors except for the SE that directly account for surface factors.

5.2.3. HCF

The diurnal ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ for ZZ has shown a sharp increase at the beginning accompanied by increasing ${\mathsf{\theta }}_{\mathrm{B}\mathrm{M}}$, weakened rainfall and developing QCloud (Figure 11a). Especially, ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ has met with ${\mathsf{\theta }}_{\mathrm{B}\mathrm{M}}$ at local noontime. Followed by the heavily developed rainfall and QCloud accompanied by sharply decreased ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ and ${\mathsf{\theta }}_{\mathrm{B}\mathrm{M}}$ at local afternoon. Nevertheless, ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ and rainfall have been slowly decreased during local night time, while ${\mathsf{\theta }}_{\mathrm{B}\mathrm{M}}$ sharply increases at early night, but decreases to meet with ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ at late night when rainfall and QCloud develop. ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ for the NY site has shown a consistent diurnal variation but with 3 hour lagged as the ZZ site, while ${\mathsf{\theta }}_{\mathrm{B}\mathrm{M}}$ varies consistently with ZZ during local daytime, but with 3 hours advanced during local nighttime (Figure 11b). Furthermore, afternoon clouds and late afternoon rainfall can be observed. Clearly, the difference between ${\mathsf{\theta }}_{\mathrm{B}\mathrm{M}}$ and ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ (or ${\mathsf{\theta }}_{\mathrm{d}\mathrm{e}\mathrm{f}}$), and QClouds developing during daytime has indicated the following rainfall increase or not for both ZZ and NY sites. Additionally, for the domain averaged diurnal variation, minimal differences between ${\mathsf{\theta }}_{\mathrm{B}\mathrm{M}}$ and ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ occur around local noon time, while ${\mathsf{\theta }}_{\mathrm{a}\mathrm{d}\mathrm{v}}$ has shown a bimodal phase, which is pronounced at local early noon and local late afternoon respectively (Figure 11c). Especially, ${\mathsf{\theta }}_{\mathrm{B}\mathrm{M}}$, ${\mathsf{\theta }}_{\mathrm{a}\mathrm{d}\mathrm{v}}$, and ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ have met with each other at local noon time. ${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$ has shown a pronounced decrease during the daytime, while ${\mathrm{h}}_{{\mathrm{E}}_{\mathrm{a}\mathrm{d}\mathrm{v}}}$ has shown an increase at the beginning and a decrease at local afternoon during the daytime (Figure 11d). Especially, $\mathrm{M}\mathrm{E}\mathrm{D}$ and ${\mathrm{E}}_{\mathrm{a}\mathrm{d}\mathrm{v}}$ are quite pronounced in local afternoon.

The $h_{\mathrm{p}\mathrm{b}\mathrm{l}}$ distribution is pronounced at around 1000 m with a frequency of around 1100 (Figure 12a). And ${\mathrm{E}}_{\mathrm{a}\mathrm{d}\mathrm{v}}$ has shown an inverted L-shaped distribution which is pronounced around 90 with a top frequency. While ${\mathsf{\theta }}_{\mathrm{d}\mathrm{e}\mathrm{f}}$ distribution is narrowly centered around 274K with a frequency of around 2800. Clearly, the L-shaped frequency distribution indicates the chances of ${\mathrm{E}}_{\mathrm{a}\mathrm{d}\mathrm{v}}$ are significantly high or low, and the most frequent ${\mathsf{\theta }}_{\mathrm{d}\mathrm{e}\mathrm{f}}$ is quite small (e.g., about 1 K larger than the zero degrees Celsius). Additionally, except in the northwest region where the upper low-value system lasts, its surrounding $h_{\mathrm{p}\mathrm{b}\mathrm{l}}$ is mostly around 1000m or larger (Figure 12b). The high ${\mathrm{E}}_{\mathrm{a}\mathrm{d}\mathrm{v}}$ area lies on the east while the low is on the west, and three small missing regions during the high ${\mathrm{E}}_{\mathrm{a}\mathrm{d}\mathrm{v}}$ area. Nevertheless, except for the marginal areas, low ${\mathsf{\theta }}_{\mathrm{d}\mathrm{e}\mathrm{f}}$ can be observed in most areas.

Overall, point wised lower ${\mathsf{\theta }}_{\mathrm{d}\mathrm{e}\mathrm{f}}$ (or ${\mathsf{\theta }}_{\mathrm{B}\mathrm{M}}-{\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$) with developing clouds could trigger the convectional rainfall at local daytime (e.g., ZZ and NY sites). Especially, that the domain averaged ${\mathsf{\theta }}_{\mathrm{B}\mathrm{M}}$,${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$, and ${\mathsf{\theta }}_{\mathrm{a}\mathrm{d}\mathrm{v}}$ all intersect at local noon time indicates the quite pronounced advantage and potential on triggering local convection for the whole domain. Meanwhile, the domain averaged $h_{\mathrm{p}\mathrm{b}\mathrm{l}}$ and $h_{{\mathrm{E}}_{\mathrm{a}\mathrm{d}\mathrm{v}}}$ intersects at the local late afternoon when significant ${\mathrm{E}}_{\mathrm{a}\mathrm{d}\mathrm{v}}$ and $\mathrm{M}\mathrm{E}\mathrm{D}$ occurs. This indicates the a significant local energy transformation occurs during the whole PBL when the main rainstorm ends. Furthermore, the 1000m boundary condensation layer height and lower temperature deficit, and significant energy transformation have dominated the local PBL cloud developing result from surface during this event. Also, the western lower $h_{\mathrm{p}\mathrm{b}\mathrm{l}}$, ${\mathrm{E}}_{\mathrm{a}\mathrm{d}\mathrm{v}}$ and ${\mathsf{\theta }}_{\mathrm{d}\mathrm{e}\mathrm{f}}$ but higher on the east indicate the clear regional differences of coupling advantage or deficit during the local daytime.

5.2.4. MDT

During local daytime, ${\mathrm{H}}_{\mathrm{s}\mathrm{f}\mathrm{c}}$, ${\mathrm{L}\mathrm{H}}_{\mathrm{s}\mathrm{f}\mathrm{c}}$, and ${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ for ZZ and NY sites are 3.3 W/m², 63.4 W/m² and 0.05, and 56.06 W/m², 252.1 W/m² and 0.22 respectively, meanwhile, ${\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$, ${\mathrm{L}\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$, and ${\mathsf{\beta }}_{\mathrm{e}\mathrm{n}\mathrm{t}}$ for the ZZ and NY sites are 47.84 W/m² 53.21 W/m² and 0.9, and 181.73 W/m², 37.7 W/m² and 4.8 respectively (Figure 13a and b). The surface flux slopes are smaller than the entrainment flux slopes for both ZZ and NY sites, and the flux slopes for the ZZ site are smaller than the NY site, while the entrainment flux slope for NY site is quite pronounced. However, during local nighttime, ${\mathrm{H}}_{\mathrm{s}\mathrm{f}\mathrm{c}}$, ${\mathrm{L}\mathrm{H}}_{\mathrm{s}\mathrm{f}\mathrm{c}}$, and ${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ for ZZ and NY sites are -12.3 W/m², 19.4 W/m² and -0.63, and -7.12 W/m², 4.78 W/m² and -1.49 respectively, meanwhile, ${\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$, ${\mathrm{L}\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$, and ${\mathsf{\beta }}_{\mathrm{e}\mathrm{n}\mathrm{t}}$ for ZZ and NY sites are 3.71 W/m², -76.76 W/m² and -0.05, and -23.11 W/m², -10.1 W/m² and 2.29 respectively. The surface latent flux has contributed more to PBL than heat flux for both ZZ and NY sites, while lighter entrainment heat flux for ZZ contribution can be identified. Clearly, the inverted flux contributions between day and night can be identified, and the pronounced entrainment contribution for NY indicates an intense PBL energy transformation. Additionally, regarding QClouds and entrainment relations, the daily, daytime, and nighttime fluxes are centered at around 0 W/m², 100 W/m², and -50 W/m² respectively, while nighttime fluxes have shown relatively narrower distribution than the daytime and daily fluxes (Figure 13c and d), while ${\mathrm{L}\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$ distribution is broader than ${\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$. Clearly, nighttime negative flux contribution indicates the energy reduction resulted from the phase change of clouds. Also, two local factors of the relation between PBL height and daily EF can be identified, e.g., the heavy (>0.6) and light (<0.3) daily EF resulted from the wet and intermediate soil respectively (Figure 13e).

The daytime averaged $\mathrm{S}\mathrm{M}$ has a narrow distribution centered around 0.4 m3/m3 with a frequency of around 2200 (Figure 14a). ${\mathrm{H}}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ and ${\mathrm{L}\mathrm{H}}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ have been most frequently observed around 60 W/m² and 280 W/m² respectively. Meanwhile, ${\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$ and ${\mathrm{L}\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$ have been most frequently observed around 90 W/m² and 100 W/m² respectively. And the latent flux distributions are clearly broader than heat fluxes. Furthermore, except in the northwest and northeast corners, quite moist soil ($\mathrm{S}\mathrm{M}$>0.3m3/m3) can be observed in most areas during the daytime (Figure 14b). However, that both lower surface heat and latent fluxes can be observed along with the main rainstorm area has indicated that besides soil moisture the atmospheric forcing has greatly shaped the surface fluxes distribution during daytime. Meanwhile, low ${\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$ and relatively higher ${\mathrm{L}\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$ have dominated the PBL flux contributions resulting from local.

Overall, the quite different surface latent or heat flux contribution to PBL for ZZ (typically affected by the rainstorm) and NY (non-typical) sites indicate the magnitude of flux slope has been suppressed by the strong atmospheric forcing during local daytime. The strong ${\mathsf{\beta }}_{\mathrm{e}\mathrm{n}\mathrm{t}}$ with weak ${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ has indicated the quite effective energy transformation of PBL, and it’s more pronounced under weak atmospheric forcing (e.g., NY site). However, except for latent flux, most nighttime flux contributions are negative. And in relation to PBL clouds, the daytime positive and nighttime negative entrainment fluxes and that broader latent entrainment flux distribution indicate the dominant effects of daytime PBL latent energy on local clouds. Also, the quite pronounced daily EF across almost all the PBL height extremes has identified the dominant wet evaporation during PBL. In addition, besides soil moisture, the atmospheric forcing has greatly shaped the surface fluxes distribution during the daytime, which is further characterized by low ${\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$ and relatively higher ${\mathrm{L}\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$.

5.2.5. TCP

TCP_LH has shown distinguished bimodal distributions during the time periods of 00~06, 06~12, and 12~18 when both positive and negative centers with significant frequency can be observed, while the distribution during 18~24 is negatively centered (Figure 15a). Except for the small scattered blocks on the west domain, low negative can TCP_LH be observed during 00~06. Following, except the whole middle east region, significant positive TCP_LH can be observed during 06~12. Then, the positive TCP_LH shrunk to the corners on the south during 12~18. Finally, little scattered positive TCP_LH can be observed on the domain margin during 18~24. Clearly, the surface coupling indicated by the positive TCP_LH is pronounced from the local afternoon to early night. Furthermore, the TCP_SH frequency distributions for different time periods have behaved consistently with TCP_LH, but narrower (Figure 15b). Also, similar spatial distributions of TCP_SH but with weaker values compared to TCP_LH for different time periods can be observed. In general, the surface coupling indicated by TCP_SH and TCP_LH are both significant from the local afternoon to early night, and they are especially pronounced in the surrounding areas of the whole middle east domain in the local afternoon.

Recall the fact that the greatly developed low-level systems during the middle east domain (Figure 3), where the daytime surface coupling has been clearly suppressed. Overall, surface coupling surrounding the middle east domain occurs in the local afternoon are significant during this event.

5.2.6. SMM

SMM for different periods has shown quite different narrow distributions with sharp frequency, whose center has increased with time increase (Figure 16). Furthermore, SMM has shown a significant northwest-southeast gradient distribution during 00~06, and this is opposite to the initial SM gradient. Following, scattered high SMM blocks lay surrounding the whole middle east region during 06~12, this is consistent with the TCP of the same period. However, a significant increase of overall SMM for the 12~18 and 18~24 periods can be found. Especially, decreased SMM in the northern region has resulted in a clear south-north gradient distribution during 18~24. This is obviously related to the northerly developed rainfall late at night.

Overall, the overall SMM has increased along with time, daytime SMM is weaker than nighttime. And the weak SMM of daytime should attribute to the moist advantage in soil evaporation and/or moist deficits in soil infiltration. Generally, the consistency of SMM with initial soil gradient at the beginning and developed rainfall at the end have indicated that both the surface soil and upper rainfall have shaped its spatial distribution.

5.3. The ensemble statistical relations

The ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ between various variations of surface and PBL states including surface flux slope (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$), surface temperature and water vapor (${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{q}}_{\mathrm{s}\mathrm{f}\mathrm{c}}$), 30m temperature and water vapor (${\mathsf{\theta }}_{30\mathrm{m}}$,${\mathrm{q}}_{30\mathrm{m}}$), 500m moist static energy and PBL height(${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$), temperature and water vapor at the height averaged between 0 and 1.5Km (${\mathsf{\theta }}_{\mathrm{h}}$,${\mathrm{q}\mathrm{v}}_{\mathrm{h}}$), and 3-hour rainfall (${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) resulted from initial soil moisture deviations are conducted to investigate the soil atmosphere relations chain during this event. Noted that the correlations with high confidence where the paired states had the same or opposite physics, or measured at the same layer or space are excluded to identify the significance in stepwise correlated state pairs (e.g., adjacent states in LoCo chain). Also, both the domain averaged (representing the overall domain estimation) and point-wised relations are evaluated. And the stepwise relation chains with significant correlation (e.g., 0.5) and confidence are selected and analyzed.

5.3.1. Spatially averaged relations

For dry perturbation (DP), significant negative correlations with high confidence can be found in the distinguished state pairs as (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$), (${\mathsf{\theta }}_{30\mathrm{m}}$,${\mathrm{q}\mathrm{v}}_{\mathrm{s}\mathrm{f}\mathrm{c}}$), (${\mathsf{\theta }}_{30\mathrm{m}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$), (${\mathsf{\theta }}_{30\mathrm{m}}$,${\mathrm{q}\mathrm{v}}_{\mathrm{h}}$), (${\mathrm{q}}_{30\mathrm{m}}$,${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$), (${\mathrm{q}}_{30\mathrm{m}}$,${\mathsf{\theta }}_{\mathrm{h}}$), (${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{q}\mathrm{v}}_{\mathrm{h}}$), (${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$) and (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathsf{\theta }}_{\mathrm{h}}$), however, significant positive correlations with high confidence can be found in the distinguished state pairs as (${\mathrm{q}}_{30\mathrm{m}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$), (${\mathrm{q}\mathrm{v}}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$) and (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{q}\mathrm{v}}_{\mathrm{h}}$) (Figure 17a). The paired states of different layers or space with opposite physics (e.g. humidity and temperature) are further excluded, then ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ that exceed the strength of 0.5 for (${\mathsf{\theta }}_{30\mathrm{m}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$), (${\mathsf{\theta }}_{30\mathrm{m}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$), (${\mathrm{q}}_{30\mathrm{m}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$) and (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{q}\mathrm{v}}_{\mathrm{h}}$) during LoCo chains have been identified (Figure 17b). For wet perturbation (WP), significant negative correlations with high confidence can be found in the distinguished state pairs as (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$), (${\mathrm{q}}_{30\mathrm{m}}$,${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$), (${\mathrm{q}}_{30\mathrm{m}}$,${\mathsf{\theta }}_{\mathrm{h}}$), and (${\mathrm{q}\mathrm{v}}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathsf{\theta }}_{\mathrm{h}}$), however, significant positive correlations with high confidence can be found in the distinguished state pair as (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathsf{\theta }}_{\mathrm{h}}$) (Figure 17c). Furthermore, ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ that exceeds the strength of 0.5 for (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathsf{\theta }}_{30\mathrm{m}}$), (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{q}}_{30\mathrm{m}}$), (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{q}\mathrm{v}}_{\mathrm{h}}$), and (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) during LoCo chains have been identified (Figure 17d).

Clearly, ${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ and ${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$ of DP have fewer correlations (or are less coupled) than those of WP, while other states of DP are more informative than those of WP in contrast. ${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$ has played a more significant role than ${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$ during the stepwise LoCo chains for both DP and WP. And ${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$ relations are only significantly on surface flux slope for DP. This indicates the significantly different response mechanisms of the atmosphere between DP (e.g., higher PBL with steeper surface flux slope) and WP (e.g., lower PBL with smoother surface flux slope).

5.3.2. Point wised relations

${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$) has shown a distribution centered around -0.4 with a frequency around 600, and ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$) has shown a distribution with a sharp unimodal centered around 0.2, then ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) has shown a distribution centered around 0.3 with a frequency around 500, while ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) has shown a distribution centered at -0.2 with the frequency around 550 (Figure 18). Clearly, (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$) and (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) have shown the dominant negative correlation, while (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$) and (${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) have behaved oppositely. Furthermore, significant convinced negative ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ (or relation intense I of equation 3.) of (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$) has mainly distributed on the east domain along the front of the weather system (Figure 3). Meanwhile, significant convinced positive ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$) within many small strip-like aggregations have been mainly distributed in the moving path of the synoptic system, which is likely related to the upper and middle layers of the dominant flows (see Section 5.1). Moreover, significant convinced positive ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$of (${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) within many strip aggregations have been mainly distributed on the west domain surrounding the main rainstorm area (null). In addition, significant convinced negative ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) have been mainly distributed in a limited distance surrounding the rainstorm area, while convinced positive ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ within many small strip aggregations have been distributed on the west away from the rainstorm area.

For DP, ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$) are mainly negative and mostly pronounced in the east region, which is the front end of the precipitation system movement direction, while ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$) mainly lies in the regions within the path of the rainstorm southwest inflows. ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) are mainly positive and mostly pronounced on the west of the rainfall area, while ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) with both negative significance near the rainfall area and positive significance far away from the rainfall area can be observed.

${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$) has shown a distribution centered at -0.3, and ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$) has shown both a broader distribution centered around 0.1 and a sharp unimodal distribution centered around 0.4, then ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) has shown both a broader distribution centered around 0.4 and a sharp unimodal distribution centered around -0.2, while ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) has shown a distribution centered at -0.1 (Figure 19). Clearly, the dominant ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ for WP is similar to DP but much larger, while the two sharp unimodal distributions have indicated the relative significance of their centered values for both (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$) and (${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$). Moreover, the ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ spatial distribution of the four pairs for WP is quite similar to DP, but the relation intensity of WP is likely narrower and smaller than that of DP.

The dominant ${\mathrm{C}\mathrm{C}}_{\mathrm{T}}$ of (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$), (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$), (${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$), and (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$) for DP and WP are -0.4, 0.2, 0.3, and -0.2, and -0.3, 0.1, 0.3, and -0.2 respectively. And this indicates the point-wised relation chains are consistent for both DP and WP, and also the relation intense of DP is relatively larger than WP. Moreover, the significant relation intensity spatial differences indicate this relation chain has been highly shaped by the atmospheric forcing (e.g., environmental flows). This is specially pronounced for the chains characterized by the relation intensity among surface flux, PBL height, and rainfall (e.g., east region (${\mathsf{\beta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$,${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$) and west region (${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$,${\mathrm{r}}_{3\mathrm{h}\mathrm{r}}$)).

6. Conclusions

This study has simulated the typical rainstorm on 20 July 2021 over central east China by using the CRA40 and GLDAS datasets, and the Noah land surface model coupled with the advanced weather research and forecasting model (WRF, version 3.9.1). Based on this, the gridded planetary boundary layer profiles and ensemble states within soil perturbations are collected to investigate the main land-atmosphere coupling characteristics during the modeled rainstorm by using various local coupling metrics and ensemble statistical relations. The results are concluded as follows.

• The long-lasting low-value system with the upper warm flow and lower strong cold invading, the mid-low layer thermodynamic situations, the convection and rainfall spatiotemporal characteristics, and the diurnal surface thermal characteristics are consistent with the available observations. However, the stratospheric (higher than 400 hPa) thermodynamics that related to the northern developed rainfall and convection, and the ${\mathsf{\theta }}_{\mathrm{s}\mathrm{f}\mathrm{c}}$ during mountain areas have been found biased. Except for the mountain areas, the main characteristics during the low atmosphere and the surface can be well documented in this modeled event.
• In CHF, the PBL near the west of the storm center is likely too stable to rain (SRC), and thee PBL on the northwest needs additional CTP to trigger convection (TR) while other regions have shown different advantaged (e.g., ACA, WSA, and DSA) and are favor of afternoon convection. In RHT, great contributions to RHT from the surface evaporation (SE), PBL warming (BLW), and non-evaporative factors (NE) have indicated their dominant roles in the local PBL clouds developing before noontime, during which SE around 0.8 and NE around 2 are especially significant. In HCF, the noontime lower buoyant mixing temperature deficits (e.g, around 274K) with developing clouds could trigger the convection except in the SRC region, while the significant energy transform of PBL occurs when the main rainstorm ends and these have dominated the daytime PBL cloud developing but with regional differences. In MDT, both the daytime PBL and surface latent energy contributions around 100 and 280 W/m² respectively have dominated the surface relations to PBL clouds, nevertheless, soil moisture and atmospheric forcing have greatly shaped the daytime surface fluxes distribution characterized as low entrainment heat flux (${\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$) but high entrainment latent flux (${\mathrm{L}\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$). In TCP, surface coupling surrounding the middle east domain occurs in the local afternoon are significant during this event. In SMM, it has increased along with time, and the comparable distributions of both initial SM and developed rainfall at the end have indicated that both the surface soil and upper rainfall have shaped the spatial distribution of SMM.
• Moist static energy (${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$) is more significant than PBL height (${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$) during the stepwise relation chains for both DP and WP. Deeper PBL with steeper surface flux slope in DP and shallower PBL with smoother surface flux slope in WP are significantly different. However, the point-wised relation chains interfaced by ${\mathrm{E}}_{\mathrm{m}\mathrm{s}\mathrm{e}}$ or ${\mathrm{h}}_{\mathrm{p}\mathrm{b}\mathrm{l}}$ are consistent for both DP and WP, while the relation intensity of DP is larger than WP. Nevertheless, the point-wised relation chains have been highly shaped by atmospheric forcing (e.g. environmental flows). This is especially pronounced for the chains characterized by relation intensity among surface flux, PBL height, and rainfall.

Overall, the WRF model initialized with the CRA40 and GLDAS datasets can well simulate various characteristics of synoptic, convection, rainfall, and surface during this event, especially for the mid-low layers. The local coupling during this event can be quantitatively identified by the ongoing mentioned metrics as follows. Wet soil contributes greatly to the local surface moisture evaporation of daytime, which then contributes to various PBL advantages for shallow convection development (e.g., increase the early noon PBL warming and/or enhance the noontime buoyant mixing) through local flux allocation. However, the weak latent flux contribution but more significant than the weaker heat flux contribution when compared to former studies [30,31] have indicated an overall suppression of the upper environmental factors (e.g., southwestern low-level inflows of the rainstorm). In contrast, the strong rainstorm forcing has significantly shaped the spatial distribution of soil moisture memory, terrestrial coupling parameters, and ensemble statistical relations during this study.

7. Discussion

Wet Soil Advantage (WSA) on the east and Atmospherically Controlled Advantage (ACA) of other regions have dominated the convection triggering potential, this indicates the dominant surface overall deficits and moistening advantage in the morning [17,18]. Both strong moisture evaporation and atmospheric interference (e.g., great daytime SE and NE in RHT) before noon have enhanced the noontime buoyant mixing (e.g., lower daytime ${\mathsf{\theta }}_{\mathrm{d}\mathrm{e}\mathrm{f}}$ in HCF), which indicates the favorable shallow clouds development surrounding the SRC region [19,20]. However, great entrainment latent fluxes (e.g., large daytime ${\mathrm{L}\mathrm{H}}_{\mathrm{e}\mathrm{n}\mathrm{t}}$ in MDT) that are mostly pronounced over the main rainstorm areas (Figure 14b) and likely controlled by the upper flows have taken over the cloud developing in the early afternoon [21,22,23,24]. The greatly developed upper systems have resulted in the huge energy transformation of the whole east domain during the late afternoon (e.g., increase ${\mathrm{E}}_{\mathrm{a}\mathrm{d}\mathrm{v}}$ in HCF when the extreme rainfall ends) [21,22,23]. Moreover, the significant TCP between noontime and afternoon has also emphasized the significance of soil-surface coupling during the rainstorm surrounding areas. However, SMM increases with rainfall develops and it is likely matched with TCP in the late afternoon period. In addition, the relation intensity for DP is greater than that for WP should be due to the initial wet soil surface, and the consistent point-wised relation chains interfaced by moist static energy or PBL height for both WP and DP indicate the overall dominant role of strong atmospheric forcing.

For real-world applications, the daytime HCF diagnoses have shown to be the reliable indicators for incoming rainfall or convection of a specific site (e.g., ${\mathsf{\theta }}_{\mathrm{d}\mathrm{e}\mathrm{f}}$ for ZZ and NY sites), and point scale weather service could possibly benefit from this; also, the point wised CHF and other metric diagnoses within various advantages can be directly referenced to for local convection initiation decisions, while the ensemble statistic relations and nighttime LoCo diagnose should be further studied. In addition, the radar [50,51] and satellite [52,53,54,55,56,57,58,59,60] studies within PBL layers but not completely observed could benefit from the LoCo metrics for the possibly comprehensive data or method evaluation.

Overall, the abundant evaporation of wet surfaces is heavily suppressed during the main rainstorm periods. Though the LoCo metrics qualify various local factors by using point wised simulation profiles under resolved model grid (4 Km) during this study, the actual land-atmosphere interaction in nature occurs on much smaller scales, e.g., the land surface feature sizes and the atmospheric turbulence scales, and the nonlinear systems over a single surface may respond quite differently among different regions, also, the spatiotemporal characteristics of paired states in both DP and WP have exhibited significantly regional differences (chaotic distributions), which should be possibly induced by external factors (e.g., human activities) [37,50]. And more investigation on local coupling among different regions should be carried out in the future to enhance the LoCo insights on typically varied underlying characteristics.