Improvement of Atmospheric Objective Analysis Over Sloping Terrain and Its Impact on Shallow‐Cumulus Clouds in Large‐Eddy Simulations

Surface topography strongly impacts the regional atmospheric circulation dynamically and thermodynamically. Over the Great Plains in the United States, the gently tilted slope is an important factor that impacts clouds, convection, and regional circulations. This study enhances an atmospheric constrained variational analysis by using a terrain‐following sigma vertical coordinate and applies to the data collected at the Southern Great Plain (SGP) site of the Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) program. Sensitivity studies are performed to examine the impact of the terrain effects on the derived large‐scale atmospheric forcing fields and the simulated shallow‐cumulus clouds in large‐eddy simulation (LES) driven by the forcing. We found that the terrain impacts on the large‐scale forcing fields are mainly at lower levels and are strongly controlled by up/downslope winds. The response of the derived forcing to the slope of the terrain is monotonic, but the response of the simulated shallow‐cumulus clouds is more complex and depends on several factors. Overall, the terrain impact is small over SGP due to the small slope angle. However, the flat‐surface assumption may cause larger biases in the large‐scale forcing fields at locations with steeper terrain. The new sigma coordinate algorithm, with its consideration of surface slope, should be more suitable to derive large‐scale objective analysis over regions with steep terrain for application to force single‐column models, cloud resolving models, and LES models.


Introduction
Surface topography strongly impacts the regional atmospheric circulation via dynamically forced upslope or downslope airflow and thermally heated elevated surfaces. Dynamically, the solid lower-boundary condition forces air flowing upslope/downslope and thus changes the atmospheric vertical motion and horizontal advections. The terrain-caused ascending motion and cold/moist advections often initiate clouds and trigger convection when the atmosphere is unstable (e.g., Banta, 1990;Houze, 2014), and affect the properties of cloud systems such as the fraction of shallow-cumulus clouds (Rabin & Martin, 1996) and the transition of shallow-to-deep convection (Panosetti et al., 2016;Tian & Parker, 2003). Thermodynamically, the diurnal heating and cooling of the terrain slope is an important factor (Holton, 1967) that contributes to the nocturnal low-level jet (Fast & McCorcle, 1990;Fedorovich et al., 2017;Parish, 2017;Parish & Oolman, 2010;Shapiro et al., 2016;Shapiro & Fedorovich, 2009), which is further associated with the initiation of deep convective systems and propagation of mesoscale convective systems (Gebauer et al., 2018;Reif & Bluestein, 2017Stelten & Gallus, 2017;Zhang et al., 2019). Overall, surface topography changes the regional circulation, clouds and convection, and further impacts the radiative budget and the hydrological cycle over the topography and the surrounding regions.
Although the topography was found to be important for clouds and convection, its consideration in objective atmospheric analysis has not been investigated. For example, the current constrained variational analysis algorithm (Zhang et al., 2001;Zhang & Lin, 1997) that is used to derive the large-scale forcing from the observational data collected at the Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) research facility for modeling clouds and convection using single-column models (SCM), cloud-resolving models (CRM), and large-eddy simulations (LES) was designed for a cylindric atmospheric column using a pressure coordinate. This algorithm assumes a flat surface at the lower boundary with the surface pressure equal to the pressure at the center of the constrained variational analysis domain.
The current constrained variational analysis has been applied on many field campaigns over different types of underlying surfaces such as ARM97 (Xie et al., 2002;Xu et al., 2002;Zhang et al., 2016), ARM0003 (Xie et al., 2005;Xu et al., 2005), the Mixed-Phase Arctic Cloud Experiment (MPACE) (Xie et al., 2006), the Tropical Warm Pool-International Cloud experiment (TWPICE) (Xie et al., 2010), the Midlatitude Continental Convective Clouds Experiment (MC3E) (Xie et al., 2014), the Green Ocean Amazon Experiment (GOAmazon) (Tang et al., 2016), and so forth, as well as continuously from 1999 to present at ARM's fixed site at the Southern Great Plains (SGP) (Tang et al., 2019;Xie et al., 2004). These products are widely used in model simulation studies with SCM (see Zhang et al., 2016, for a complete review of SCM studies with the ARM forcing data) and with LES such as the DOE LES ARM Symbiotic Simulation and Observation Workflow (LASSO) (Gustafson et al., 2020) project for shallow-cumulus clouds. Considering the importance of terrain, specifically the slope of terrain, on the large-scale atmospheric environment, it is important to investigate how sloping terrain would impact the derived large-scale forcing data and therefore the simulated cloud systems. In this study, we develop a new version of constrained variational analysis algorithm using a terrain-following sigma (σ) coordinate to capture the effect of surface terrain and investigate the impact on the derived large-scale forcing fields and simulated shallow-cumulus clouds over SGP.
The paper is organized as follows: Section 2 introduces the enhancement of constrained variational analysis from pressure coordinate to sigma coordinate; section 3 shows a sensitivity study to investigate the impact of surface terrain on the derived large-scale forcing; section 4 further examines the impact of the slope-induced changes in the large-scale forcing on the shallow-cumulus clouds simulation in LES; section 5 is the summary and discussion.

Constrained Variational Analysis in Sigma Coordinate
The general idea of constrained variational analysis is to adjust the atmospheric state variables (u, v, s, q) to satisfy the following conservations of mass, water vapor, heat, and momentum within an atmospheric column: where s = C p T + gz is dry static energy; E is surface evaporation; P is surface precipitation; L v is latent heat of vaporization; LWP is cloud liquid water path; R TOA and R SFC are net downward radiations at the top of atmosphere (TOA) and surface, respectively; SH is surface sensible heat flux; f is the Coriolis parameter; k → is the unit vector in vertical direction; ϕ is geopotential; and τ s → is surface wind stress. Other variables (g; V h ! ; p; t; q) are as commonly used in meteorology. Ice processes are ignored for simplicity. The terms on the right-hand sides of the equations are constraint variables usually from observations measured at the surface or the TOA of the atmospheric column, while those on the left-hand sides of the equations are vertical integrals from surface to TOA within the column.
The above equations are written in pressure coordinate, with the surface assumed at constant pressure. In sigma coordinate, vertical coordinate sigma (σ) is defined as where π ≡ p SFC − p TOA and p TOA is set as 100 hPa in this study. For simplicity, we are not considering the momentum constraints in this study. The mass, moisture, and energy constraint equations in sigma coordinate are written as Comparing to the constraint equations in pressure coordinate (Equations 1-3), the equations in sigma coordinate have similar format but with an additional variable π and the change of vertical integration from dp to dσ. Therefore, the sigma coordinate constrained variational analysis (σCVA) is implemented based on the original code structure and numerical algorithm of the pressure coordinate constrained variational analysis (pCVA). In the constrained variational analysis, the divergence terms ð∇ · V h ! ; ∇ · ðq V h ! Þ and ∇ · ðs V h ! ÞÞ are calculated as line integrals along the analysis domain using Green's theorem. Therefore, only the terrain along the boundary of the domain (i.e., the slope) is considered in the revised constrained variational analysis algorithm ( Figure 1). As the large-scale forcing fields represent the average of the analysis domain, subgrid terrain variations within the domain cannot be resolved and are not the subject of this study.
After performing the constrained variational analysis in sigma coordinate, the large-scale forcing fields are converted back into pressure coordinate for the purpose of comparing with original large-scale forcing data and driving SCM/CRM/LES: Surface heights at all boundary points are the surface heights in reality 2x The difference between surface height and the domain-mean height doubles −2x The difference between surface height and the domain-mean height doubles and reverses In the next section we will investigate the terrain effects on the large-scale forcing using a set of sensitivity studies.

Sensitivity Study on the Forcing
We choose 1 month of continuous forcing data (doi: 10.5439/1273323) in July 2015 at SGP as a test case to examine the sensitivity to the surface slope. Over the SGP region, the slope angle is about 0.15°or 1/400 (Holton, 1967;Shapiro & Fedorovich, 2009) from west down to the east. In the experiment, four sets of terrains are designed as shown in Table 1. The 1x terrain is the actual terrain over SGP region, as shown in Figure 2a. The surface heights in 0x, 2x and −2x terrains are calculated by adding the domain-mean height and the height deviations from the domain-mean multiplied by 0, 2, and −2, respectively. Surface pressure is then calculated based on hydrostatic equation from surface height and surface temperature. Other surface variables such as surface temperature, moisture, and winds are interpolated from the 3-D atmospheric fields into the level of the new terrain. The 3-D fields used in this study are obtained from the NOAA Rapid Refresh (RAP) analysis (https://rapidrefresh.noaa.gov/). For each terrain setup we perform two runs of constrained variational analysis: the "pCVA" run using the original algorithm released as version 2 of continuous forcing (Tang et al., 2019;Xie et al., 2004) and the "σCVA" run described in section 2.

Over Flat Surface
Figures 3a and 3b show divergence field derived from pCVA and σCVA, respectively, over the flat (0x) surface. Figure 3c shows the difference between the two. Similar plots are also shown for pressure vertical velocity (omega), horizontal advection of dry static energy (s), and horizontal advection of moisture (q) in Figures 4, 5, and 6, respectively. There are differences seen in divergence and horizontal advections shown mainly at the upper levels between 100 and 300 hPa (except for advection of q). Although the lower boundary conditions for pCVA and σCVA are the same for the flat surface design, we do not expect the two experiments to be exactly the same. The noticeable differences can come from the smoothing and adjustment processes used in the constrained variational analysis which is related to the atmospheric pressure. In sigma coordinate, the pressure at a given sigma level is not a constant value; this changes the weights of smoothing and adjusting. Therefore, the final large-scale forcing fields are slightly different in the sigma coordinate system than in the pressure coordinate system. This is confirmed in a sensitivity test that does not apply smoothing and uses constant adjustment weights. In this test, these noticeable differences in the upper troposphere disappear (not shown). Random errors during the iterations may also contribute to the adjustment differences. The differences are more prominent when the absolute values of the actual fields are large, since it usually requires large adjustments and thus is easier to change. Overall, the differences of large-scale forcing data over flat surface are relatively small and have little impact on the forcing patterns (Figures 3-6, panels a and b).

Sensitivity to Different Slopes
Figures 3d, 3g, 3j show the differences between divergence derived over different terrains and the divergence derived over 0x terrain using pCVA; those derived using σCVA are shown in Figures 3e, 3h, 3k. Similar plots For σCVA, although we only change the slope of the surface terrain, the entire profiles of the derived large-scale forcing are adjusted by the variational analysis. The vertical profiles of the atmospheric fields are from the RAP analysis based on the real surface terrain, which are the same for different terrain designs. However, when the surface slope is changed in the sensitivity study, the column-integrated budgets are also changed due to the different atmospheric thickness along the domain boundary. To close the budgets in the constrained variational analysis, the entire vertical profiles have to be adjusted based on the weighted least squares. The changes of adjustments for different surface terrain slopes in σCVA are related to the tilting of sigma surface and pressure surface, which decreases with height. At lower levels, the changes of adjustments show temporal oscillation diurnally, which is likely related to the diurnal oscillation of up/downslope wind (Figure 2c). More discussion on the relationship of the terrain-induced change of large-scale forcing and up/downslope wind will be given in section 3.4.

Difference Between Pressure and Sigma Coordinates Over Terrain
Figures 3f, 3i, and 3l show the differences of divergence from pCVA and σCVA over different terrains. Similar plots for omega, horizontal advection of s, and horizontal advection of q are shown in Figures 4, 5, and 6, respectively. These differences are the sum of the differences due to the change of vertical coordinates discussed in section 3.1, and the differences introduced by surface terrain discussed in section 3.2. On average, the terrain effect dominants the differences seen between pCVA and σCVA, especially at the lower levels. Over SGP (1x), the overall difference between the large-scale forcing from pCVA and σCVA is 10-20%. It does not impact the overall structure of the large-scale forcing and only has a very minor impact on the

10.1029/2020JD032492
Journal of Geophysical Research: Atmospheres simulated clouds, as shown later. This indicates that pCVA and σCVA provide similar large-scale forcing fields over regions with mild slope such as at SGP.

Vertical Structure of Derived Forcing Fields
The time-mean vertical profiles of these forcing fields from pCVA and σCVA are plotted in Figure 7. The profiles from pCVA over different terrains are almost identical and close to the σCVA in 0x terrain, except for the lowest level above the surface. The small differences at the upper levels are due to weight differences and random errors discussed in section 3.1, and the differences at the lowest level are due to the change of surface variables discussed in section 3.2 σCVA over different terrains differ prominently at the lower levels below 700 hPa. Divergence and omega differences are primarily seen below 800 hPa, where the prevailing wind is from the south or the east, which is upslope (Figure 2c). Over a sloping terrain, the pressure intervals at the upslope points are smaller than at the downslope points between the two iso-sigma levels. To ensure mass balance, stronger divergence and upward motion are forced, considering that the air cannot penetrate the lower boundary. The horizontal wind has a diurnal oscillation as shown in Figure 2c, which may be related to the oscillation of the forcing differences seen in Figures 3-6. The horizontal advections of s and q are adjusted mainly between 700 to 850 hPa to maintain the column-integrated energy and moisture budgets. The adjustments are related to the tilting sigma surface at the boundary layer top, where the horizontal gradient of dry static energy and moisture along the sigma surface are relatively large because of large pressure difference. Overall, for the large-scale forcing data derived using σCVA, the terrain effect increases with increasing terrain slope, and the effect of the opposite terrain is reversed.

Impact of Sigma Coordinate-Based Forcing on Shallow-Cumulus in LES
The terrain effect on the large-scale forcing discussed in the last section may further impact the simulation of clouds and convection when the forcing is used to drive SCM/CRM/LES. In this section we will use LES to investigate how the simulated clouds could be impacted by the inclusion or neglect of the terrain effect.   (Endo et al., 2015). The model configuration follows LASSO LES for 2016 case simulations (Gustafson et al., 2018). Physics parameterizations include the Thompson microphysics (Thompson et al., 2008), Rapid Radiation Transfer Model for Global Climate Models (RRTMG) longwave and shortwave radiation (Iacono et al., 2008), and the 1.5-order turbulent kinetic energy subgrid-scale turbulence (Deardorff, 1980).

Difference of Forcing and Clouds Between Pressure and Sigma Coordinates
We first examine the large-scale fields from the forcing data sets for the three selected cases. Figure 8 shows the large-scale relative humidity (RH) and the forcing fields from σCVA over the actual terrain, and Figure 9 shows the differences of these fields between σCVA and pCVA. The spatial patterns of these fields from pCVA are similar to those from σCVA (not shown), so the ranges of colorbars in Figure 9 are smaller than in Figure 8. All the three test cases show similar relative humidity fields but different large-scale circulation structures. All the relative humidity fields have a peak at 800-900 hPa and show time variations associated with daytime evolution of the boundary layer. The 20,160,518 case has strong westerly wind at the upper levels and easterly wind at the lower levels (Figure 8a), corresponding to upper-level downslope wind and lower-level upslope wind over SGP. As discussed in section 3.4, to maintain the mass balance in σCVA over terrain, the omega field shows weaker downward motion at the lower levels and stronger downward motion  Figure 9d is about 2 to 4 mb/hr and the actual omega in Figure 8d is in the order of 10 mb/hr. The 20,160,611 and 20,160,719 cases have much weaker upper-level winds and moderate lower-level southerly winds. Since the slope contour line is primarily northeast-southwest (Figure 2), when the horizontal wind turns from southwesterly (parallel to the slope contour) to southerly (upslope) at around 1800 UTC (Figures 8b and 8c), the wind is more upslope and the σCVA produces weaker lower-level downward motion to ensure mass balance (Figures 9e and 9f). The terrain-induced lower-level upward motion anomaly is limited below 800 hPa (~1.5 km above ground level [AGL]). As seen later, this is below the cloud base and does not directly impact the LES-simulated clouds. Figure 10 shows the time-height cross section of cloud fraction of the three cases from the ARM Active Remote Sensing of Cloud Locations (ARSCL, doi: 10.5439/1350630) product and from the LES simulations using large-scale forcing from pCVA and σCVA, as well as the differences between the two simulations. Overall, LES well captures the timing, location, and evolution of the observed shallow-cumulus clouds, with differences of cloud fraction possibly due to the mismatch of sampling scales (ARSCL data are fraction of clouds detected by point measurements in 4-s frequency while LES data represent fraction of clouds in the model domain in 30-s frequency, both averaged in 10 min). The cloud fraction differences from pCVA

10.1029/2020JD032492
Journal of Geophysical Research: Atmospheres and σCVA simulations are larger for the 20,160,518 case than the other two cases. This is partly because the 20,160,518 case has the larger terrain-induced large-scale forcing differences than the other two cases (Figures 9d-9f), and partly because the cloud base for the 20,160,518 case is between 0.4 and 1.2 km in the morning, which is within the levels of large terrain-induced upward anomaly or reduced subsidence. The reduced subsidence decreases the suppression of boundary layer deepening and moistening and thus favors the formation of shallow-cumulus clouds. For the 20,160,611 and 20,160,719 cases, the cloud bases are above the levels of large terrain-induced large-scale forcing as shown before; therefore, the simulated clouds are less impacted. Overall, for all the three cases the simulated cloud structures and evolution are not impacted by the use of large-scale forcing from pCVA and σCVA. This indicates that the flat-surface assumption in the pCVA is good enough to produce reasonable large-scale forcing data for the shallow-cumulus cases at the SGP site. Figure 11 shows the time-averaged (1200-0300 UTC) vertical profiles of large-scale forcing fields derived from pCVA and σCVA for 0x, 1x, 2x, and −2x terrains for the three cases. The major fields impacted by the surface terrain is horizontal divergence and omega. Consistent with the July 2015 case discussed in section 3, the large-scale forcing from pCVA is insensitive to the terrain change, and that from σCVA is monotonically responding to the terrain change (forcing difference increases with increasing terrain slope and reverses with opposite terrain). Among the three cases, the 20,160,518 case is more sensitive to the terrain, due to its stronger up/downslope wind (Figure 8a). The 20,160,611 and 20,160,719 cases are much less Figure 10. (From top to bottom): Cloud fractions from ARSCL observations, LESs using pCVA-based forcing, LESs using σCVA-based forcing, and the differences between the σCVAand pCVA-based simulations for the real terrain (1x) setup. The three columns from left to right are for three cases of 20,160,518, 20,160,611, and 20,160,719, respectively. Local time is UTC-6.

10.1029/2020JD032492
Journal of Geophysical Research: Atmospheres Figure 11. Time-average vertical profiles of (from up to bottom) divergence, omega, horizontal advection of s, and horizontal advection of q for (from left to right) 20, 160,518, 20,160,611, and 20,160,719. Solid lines are simulations using σCVA-based forcing, while dashed lines are simulations using pCVA-based forcing. Different colors represent different terrain configurations.

10.1029/2020JD032492
Journal of Geophysical Research: Atmospheres sensitive to the terrain, especially at the levels above 3 km (~700 hPa). The terrain effect on the large-scale forcing is highly related to the wind blowing up/downslope, as discussed in section 3.4. For more information, the time evolution of the large-scale forcing fields with different terrain configurations for the three cases is shown in supporting information Figures S1-S12. Figure 12 shows the time-averaged vertical profiles of liquid water content (LWC) and cloud fraction, as well as time series of 0-3 km integrated cloud liquid water path (LWP) and cloud fraction for the 20,160,518 case. The simulated cloud properties are consistent with the large-scale forcing changes due to the different terrains. The major differences in LWC and cloud fraction occur between 0.4 to 2 km, with the strongest upward motion anomaly (2x terrain) having the largest LWC and cloud fraction, and the strongest downward motion anomaly (−2x terrain) having the smallest LWC and cloud fraction. The large LWC and cloud fraction of −2x terrain simulation near 2 km (~760 hPa) may be related to the large vertical gradient in temperature advection at this level, which increases the stratification of the atmosphere and prevents the vertical development of the clouds. In the time series, the cloud fraction has large differences before 1900 UTC (1300 local time). During this period, the simulated cloud is low (below 1 km) and thin (a couple hundred meters deep); therefore, it is very sensitive to the forcing change. In the afternoon, the clouds rise to the levels with smaller terrain-induced forcing, so that the sensitivity to the terrain-induce forcing reduces. This indicates the complexity and nonlinearity of the cloud response to the large-scale forcing. Figures 13 and 14 show same cloud variables for the 20,160,611 and 20,160,719 cases, respectively, but the time series are based on integration from 0 to 5 km. The cloud fractions are almost identical for the forcing with different terrain configurations, while cloud LWC and LWP are slightly different and not monotonically responding to the surface terrain slopes. The much smaller sensitivity to terrain in these two cases comparing to the larger sensitivity in 20,160,518 case may come from two reasons: (1) much weaker up/downslope wind causes smaller terrain-induced changes in the large-scale forcing and (2) clouds form at the higher levels where there are small forcing differences; therefore, the larger changes in the forcing

10.1029/2020JD032492
Journal of Geophysical Research: Atmospheres below the clouds do not directly impact the cloud properties. The time-height evolution of simulated cloud fraction is given in Figures S13-S15.

Summary and Discussion
Surface terrain is an important boundary condition that alters the atmospheric circulation and impacts the initiation and evolution of clouds and convection. However, the surface is assumed as flat in the constrained variational analysis method that is currently used in ARM to derive the large-scale forcing fields for driving SCMs/CRMs/LESs. To take into account the terrain effect in the constrained variational analysis algorithm, we replaced the current pressure coordinate system with a terrain-following sigma coordinate system so that the slope of the underlying surface is considered. We also investigated the impact of the terrain on the derived large-scale forcing through a sensitivity study by artificially changing the slope of the terrain. We have shown that 1. the terrain impact on the forcing is strongly controlled by the up/downslope wind; 2. the terrain impact is stronger at lower levels, where the tilting of sigma surface and pressure surface is large; 3. the terrain impact is stronger for steeper slope; and 4. when the slope reverses, the terrain-induced forcing change also reverses.
This study has shown that the impact of terrain slope can be well captured by the proposed sigma coordinate constrained variational analysis. However, the new method is still unable to capture the subgrid terrain variation within the analysis domain. This issue could be partially addressed by carefully designing the analysis domain to avoid possible impacts from subgrid terrain variability.
The forcing data derived from the original and the new methods are tested using LES model. Although the impact of terrain on the large-scale forcing is monotonic, the simulated clouds show more complex responses. In the three shallow-cumulus cloud cases we examined, the 20,160,518 case shows more linear response while the 20,160,611 and 20,160,719 cases show more nonlinear responses. We have found that the strength of up/downslope winds and the altitude of the cloud layer impact the sensitivity of shallow-cumulus clouds simulation. The shallow-cumulus cloud cases tested in this study are primarily locally driven with weak large-scale forcing. For rainy cases, precipitation is the strongest constraint for the derived large-scale forcing in the constrained variational analysis. The terrain effect on the large-scale forcing and the corresponding model simulations may be smaller. This is subject to future study.
It is worthwhile noting that, over a flat surface, the forcing fields from the two methods only have small differences. These differences could be considered as uncertainties in the forcing fields introduced by the smoothing and adjustment processes implemented in the constrained variational analysis, and it is difficult to make the judgment which version should be more trusted because we do not have observed forcing fields to verify them. Since the overall differences between these two forcing data sets are quite small, it will not cause a big issue in using either methods to drive cloud models over a region where the surface is flat. With the mild slope at the SGP region, the terrain-induced forcing change is~10% to 20% of the actual forcing magnitude, which only has a small impact on the overall forcing pattern. As shown in the LES tests, the difference in the simulated cloud with the two forcing data sets is much smaller than that between model and observations. The relatively small difference suggests that the current large-scale forcing data sets from the original constrained variational analysis are reliable for SGP. However, when a research area has a steeper slope, we anticipate the flat-surface assumption may cause larger biases on the derived large-scale forcing fields, which may further impact the simulation of clouds and convection. The recently completed ARM field campaign, Cloud, Aerosol, and Complex Terrain Interactions (CACTI), has been conducted in an area with a much steeper slope with a slope angle of about 1/40, which is 10 times of the slope over SGP. We plan to derive large-scale forcing fields using both the original and new algorithms for CACTI and investigate the performances of the two types of forcing data on deep convective systems over steep terrain in the future.