The Effect of Different Implementations of the Weak Temperature Gradient Approximation in Cloud Resolving Models

The Weak Temperature Gradient (WTG) approximation is a popular method used to couple convection in limited‐area domain simulations with large‐scale dynamics. However, several different schemes have been created to implement this approximation, and these different WTG schemes show a wide range of different results in an idealized framework. Our investigation shows that different model behavior is caused by the treatment of the different baroclinic modes by the different WTG schemes. More specifically, we hypothesize that the relative strengths of the baroclinic modes plays a large role in these differences, and show that modifying these schemes such that they treat the baroclinic modes in a similar manner accounts for many of the significant differences observed.

10.1029/2023GL104350 2 of 10 equilibrium (RCE) profile produced by the model itself, similar to previous studies (e.g., Emanuel et al., 2014;Sessions et al., 2010).If this state is stable with WTG, applying the WTG framework will cause the model to converge to it, given the right initial conditions.However, if the RCE state is unstable with WTG, deviations can be expected to occur under the WTG framework, either in the form of the abovementioned bifurcation of the resulting model states in the previous paragraph, or in some form of oscillatory behavior if no stable state exists.
Over time, three main schemes that implement the WTG approximation in models have emerged: the (a) Temperature Gradient Relaxation (TGR) implementation (Raymond & Zeng, 2005), the (b) Damped Gravity Wave (DGW) implementation (Blossey et al., 2009;Kuang, 2008a), and the (c) Spectral (SPC) Weak Temperature Gradient implementation (Herman & Raymond, 2014).We provide more elaboration on these schemes in Section 2. Despite the prevalence of these schemes in modeling work for tropical climate, they often produce noticeably different results.For example, several studies (e.g., Daleu et al., 2015;Romps, 2012a;Romps, 2012b) show that the TGR implementation results in a vertical profile that is more top-heavy than the DGW implementation (Figure S1 in Supporting Information S1).
Our study builds upon previous work done to quantify these discrepancies in model results (e.g., Daleu et al., 2015) by attempting to understand why they arise in the first place.In Section 2 we will discuss these three main implementations of the WTG approximation in models, explain how we implement them in Section 3 and then show in Section 4 that these schemes give markedly different results even in idealized setups.In Section 5, we perform a vertical-mode decomposition of the WTG schemes, and discuss our results in the framework of Gross Moist Stability in Section 6.

Weak Temperature Gradient Implementations in Models
There are three major schemes enforcing the WTG approximation that are widely used in single-column and small-domain cloud resolving modes.

The Temperature Gradient Relaxation Implementation
The TGR implementation directly links local buoyancy anomalies to large-scale vertical motion.Differences in buoyancy between the single-column or small-domain cloud-resolving model and the large-scale environment over a time-scale τ are balanced by the vertical advection of potential temperature w∂ z θ, such that at a height z i in the free troposphere the WTG-induced vertical velocity w wtg is given by: where z t is the height of the tropopause, θ is the model potential temperature and θ 0 is the reference large-scale potential temperature. (⋅) represents the domain-average of the variable (•).This implementation was first done by Raymond and Zeng (2005), and has been used in a number of other studies (e.g., Daleu et al., 2012;Sessions et al., 2010).In contrast to Raymond and Zeng (2005) who fixed z t = 15 km, in our runs we allowed z t to vary by setting it to be the level of the cold-point tropopause.We decided to let this level fluctuate over time for two reasons: (a) for consistency in our comparison with the setup of Blossey et al. (2009), and (b) the mean-state tropopause height in our experimental runs can change depending on the mean-state of the model when the WTG approximation is enforced-a model in a moist, highly-precipitating state will have a higher tropopause height compared to a model in a dry, non-precipitating state (Figure S1 in Supporting Information S1).To prevent unrealistically large values of w wtg , we set a lower-bound on static stability,  ( ∕ Raymond and Zeng (2005).

The Damped Gravity Wave Implementation
In contrast to the TGR implementation, the link between buoyancy and temperature anomalies to large-scale vertical motion in the DGW implementation is derived from the damping of gravity wave perturbations in the momentum equations (without Coriolis force) using a Rayleigh damping coefficient a m : where the other variables have their usual meteorological meaning.(•)′ represents the perturbation of the variable (•) from the large-scale reference profile.Assuming (a) steady state, (b) that a m is constant with height, and (c) using the ideal gas, hydrostatic balance and mass conservation laws, the momentum equations are transformed into the following governing equation for WTG-induced pressure velocity ω wtg in pressure-coordinates: where R d is the dry gas constant, T v is the virtual temperature, and k is the horizontal wavenumber of the gravity wave.As mentioned above,  (⋅) and (•)′ respectively denote the domain average of (•) and its perturbation from the large-scale reference profile.The strength of the implementation is controlled by k 2 /a m .As varying either will change model behavior in a similar manner, we keep k = 2π/λ constant, taking λ = 2,600 km and a m = 1 day −1 as in Blossey et al. (2009), and multiply k 2 /a m by a dimensionless constant α.Kuang (2008a) also derived a similar form using height coordinates instead of pressure coordinates.However, we used Equation 4 to maintain consistency with Blossey et al. (2009).Furthermore, while we used virtual temperature T v to be consistent with previous studies (e.g., Blossey et al., 2009), we have also verified by replacing T v with absolute temperature T that the virtual effect does not contribute significantly to the differences we see across the different implementations.Herman and Raymond (2014) published an updated version of the TGR implementation of Raymond and Zeng (2005).Instead of assuming that gravity waves of all vertical wavelengths are equally effective in redistributing buoyancy/temperature anomalies, the relaxation time τ j for the jth vertical mode is τ j = j • τ, where τ is the relaxation timescale of the first vertical mode.Therefore, we perform a vertical decomposition of both vertical velocity and scaled potential temperature anomaly as follows:

The Spectral Weak Temperature Gradient
where the vertical modes are of the form: 10.1029/2023GL104350 4 of 10 where similar to the TGR implementation as above, we decided to let z t fluctuate over time.The Spectral Weak Temperature Gradient implementation then assumes that strength of the vertical mode of vertical velocity as a function of the vertical mode of the scaled potential temperature anomaly is given by w j = θ j /τ j , such that the spectral WTG vertical velocity is given by We take n = 32 and neglect higher-order modes as importance decreases as the order increases.

Model Description
We used the System for Atmospheric Modeling (SAM) (Khairoutdinov & Randall, 2003) version 6.11.8.The model solves the anelastic continuity, momentum, and tracer conservation equations, with total nonprecipitating water (vapor, cloud water, cloud ice) and total precipitating water (rain, snow, graupel) included as prognostic thermodynamic variables.Simulations are run in three dimensions with doubly-periodic boundaries and a horizontal resolution at 2 km to permit clouds, with a horizontal domain of 128 km by 128 km.There are 64 vertical levels in our model, with the vertical spacing increasing from 50 m at the boundary layer to around 500 m at the tropical tropopause, to a total height of ∼27 km with a rigid upper-bound.Damping is applied to the upper third of the model domain to reduce reflection of gravity waves.A simple Smagorinsky-type scheme is used for the effect of subgrid-scale motion.
In all our experiments, the SST is fixed at 300 K, spatially uniform and time-invariant.We run two versions of the model: (a) the default version of SAM with the RRTM radiative scheme (Mlawer et al., 1997), and (b) the idealized radiative scheme of Pauluis and Garner (2006) that uses a fixed radiative-cooling rate of −1.5 K day −1 in the troposphere and Newtonian relaxation when the temperature is less than 205 K with a relaxation timescale of 5 days.

Obtaining the Large-Scale Reference Profiles for WTG Simulations
As mentioned in Section 1, all simulations involving the WTG approximation require coupling of the model to a large-scale profile of the relevant buoyancy-variable (for e.g., in the DGW implementation (Equation 5) this would be virtual temperature T v ).These reference profiles were obtained by spinning a 10-member ensemble to RCE over 2,000 days, taking the last 500 days for statistics, with separate profiles constructed for full-radiation and idealized-radiation simulations.We then take the average of the vertical profiles of temperature and specific humidity of these ensemble members to construct the large-scale reference profiles.
When each model run is initialized, SAM reads in a sounding file containing vertical heights, pressure levels, and the profiles of potential temperature and specific humidity in order to construct the initial state of the atmosphere.
If the profile is close to RCE that is in balance with the time-invariant SST, then the state of the equilibrated atmosphere after 1,000 days should be close to the initial profile.We reinitialize the model with the equilibrated sounding profiles of temperature and specific humidity from our 10-member ensemble run and repeat this cycle until the root-mean-squared difference between the initial and final ensemble-mean temperature profiles was <0.01 K.

Implementing the Different Schemes Into SAM
Once the models have been spun-up to RCE, we take the average temperature and humidity vertical profiles of the 10-member ensemble as the large-scale reference profiles.We then enforce the WTG approximation over a range of τ or α (depending on the scheme used) values, and run a 5-member ensemble over a period of 250 days for each of the configurations, taking statistics every hour over the last 100 days.For each member in the ensembles, perturbations were made to the initial state of the model, resulting in a mix of wet and dry final states.In order to make it easier to obtain both wet-and dry-states of the multiple-equilibria regime, we perturbed the large-scale reference profile uniformly in the vertical by −0.05 K for another 5-member ensemble, and +0.05 K for a final 5-member ensemble respectively. 10.1029/2023GL104350 5 of 10 In order to showcase the difference between RCE and WTG states, we implement a smooth transition from a pseudo-RCE state (α(t = 0) = τ(t = 0) = ∞) to a WTG state (α = α 0 or τ = τ 0 ), where α 0 and τ 0 are the final strength of the WTG approximation at t = t wtg .In all our experiments, we take t wtg = 25 days, which means that in our experimental runs the WTG implementations will reach maximum strength at 25 days from model startup.

Divergence in Model Behavior With Different WTG Schemes Under an Idealized Model Framework
Applying the WTG approximation to small-domain models with interactive radiative schemes results in multiple-equilibria (see Figure 2i), with permanent wet and dry model states both being possible outcomes regardless of the WTG scheme.Results from the different WTG schemes are qualitatively similar to each other and to the results of Emanuel et al. (2014) using the MITgcm in single-column mode, but have significant quantitative differences.As the strength of the WTG adjustment increases, the model eventually enters an oscillatory regime where the model rapidly alternates between wet and dry states (see whiskers in Figure 2, and daily-averaged time-series plots in Figure S2 in Supporting Information S1).However, we note that the magnitude of these oscillations is very small in TGR simulations compared to when the DGW and SPC implementations are used.
Model behavior varies even more markedly between the different WTG schemes (Fig. 2ii) when the idealized-radiation framework described in Section 3 is used.In the DGW framework, while the multiple-equilibrium regime is greatly reduced compared to realistic-radiation simulations, it is still significant and leads into an oscillatory regime (see the timeseries of daily-averaged precipitation in Figure S3 in Supporting Information S1), and the results found by Sessions et al. (2016).However in the SPC framework, the bifurcation between the wet-and dry-states of the multiple-equilibrium regime is reduced until it is almost indistinguishable from the RCE-mean (though the presence of yellow and blue dots in Figure 2cii indicates that it is not entirely gone).A significant oscillatory regime still exists when the strength of the implementation is large (τ < 10 hr).
In the TGR framework the oscillatory regime does not even become significant until τ approaches values that are not physical (e.g., τ < 0.5 hr).
We see that these differences in model behavior upon the implementation of different WTG schemes is larger in a simple model framework with idealized radiation (Figure 2).The implementation of full interactive radiation serves to mask the differences in model behavior by amplifying the multiple-equilibria regime, similar to how Since the contrast between WTG schemes is best shown in model frameworks with idealized radiation, the model results in the sections below are therefore limited to experimental setups with idealized radiation.Nonetheless, because the model results from the DGW and SPC implementations are qualitatively more similar to each other than between the DGW and TGR implementations across different radiation schemes, we believe that our discussions in Sections 5 and 6 would still be applicable to model frameworks with fully-interactive radiation.

Revisiting the Different WTG Schemes Using a Vertical Mode Decomposition
We now seek to understand the differences between these implementations.Similar to Kuang (2008b); Herman and Raymond (2014), we decompose both the left-and right-hand-side components of Equation 4 into linear combinations of the vertical modes G j (see Equation 6): Noting that the equations in the DGW implementation solve not for ω′, but for ∂ zz ω′, we see that ω j and T j are related to each other as follows: 2 , and since fluctuations in c depend only on z t , which can be assumed to be constant compared to the range of α explored, we can assume that c is constant as well.
A similar analysis of the TGR implementation gives: Lastly, analysis of the SPC implementation gives (see Section 2.3): A comparison of Equations 10, 12 and 13 show that the higher-order modes in vertical velocity associated with the respective higher-order vertical modes of local buoyancy-temperature anomalies are different in the different WTG schemes.For a given buoyancy-temperature perturbation, the resulting higher-order modes in vertical velocity decrease in strength in order of (a) DGW, (b) SPC and (c) TGR respectively.Therefore, the vertical structure of vertical velocity will be different across the different WTG schemes, where profiles from the TGR implementation are likely to have stronger higher-order modes compared to the profiles from the DGW or SPC implementations, and this has been well-documented (Romps, 2012b;Daleu et al., 2015, see also Figure S1 in Supporting Information S1).

Bringing the Different WTG Schemes Together Using the Gross Moist Stability Framework
Previous studies have shown that the basic dynamics of convectively coupled tropical waves can largely be captured by models which contain the first two baroclinic modes of the vertical structure of the tropical atmosphere (e.g., Haertel & Kiladis, 2004;Khouider & Majda, 2006;Kuang, 2008b;Majda & Shefter, 2001;Mapes, 2000).Using the first two modes and ignoring all higher-order terms, we analyze our vertical mode decomposition of the various WTG implementations in the context of the GMS framework.Following Raymond et al. (2009); Kuang (2011); Inoue andBack (2015, 2017) we define: where angle brackets represent vertical averages.This is the ratio of the lateral export of moist static energy h to the vertical export of dry static energy s.W 1 , and W 2 are the first and second modes of vertical velocity.Taking idealized vertical profiles of the dry and moist static energies shown in Figure 3, we see that Equation 14 can be reduced to: Thus, any change to the GMS is ultimately dominated by the relative strengths of the first two baroclinic modes.However, as we have discussed previously, the response of higher-order baroclinic modes to a given buoyancy perturbation is different across the WTG implementations.For example, because the SPC and TGR implementations result in stronger second baroclinic modes, and thus stronger second-order modes of vertical velocity, it would favor higher GMS magnitudes than the DGW implementation and thus larger magnitudes of export (or import) of moist static energy.This is in line with the characterization of GMS as a quantity that describes the (de)stabilization mechanisms of convective disturbances in the atmosphere (e.g., Inoue & Back, 2015, 2017;Raymond et al., 2009).We believe that the ratio w r = w 2 /w 1 therefore constrains how rapidly these convective disturbances are magnified/reduced.
As an example, we consider a moist environment with stronger-than-RCE deep convection.Such a moist and strongly-convecting environment will often have temperature profiles that are warmer in the upper troposphere and cooler in the lower troposphere.Assuming Rayleigh friction, a limited area warm anomaly aloft and an underlying cool anomaly will result in an adjustment circulation (in a non-rotating environment) that has ascent at upper levels and descent at lower levels, that is, a circulation with positive second mode structure shown in Figure 3a.As elaborated by Raymond et al. (2009); Inoue andBack (2015, 2017) and many other studies, this stratiform profile of convection tends to export GMS and return the domain-mean back to RCE.The greater the value w r , the stronger this tendency.As the TGR implementation's greater emphasis on higher-order baroclinic modes naturally results in higher values of w r , we see that in the idealized-radiation framework there is no visible bifurcation or multiple-equilibria (Figure 2aii) when the TGR implementation is used.In contrast, higher-order baroclinic modes are weak in the DGW implementation, which results in a multiple-equilibria regime and a noticeable bifurcation in the resulting wet and dry states (Figure 2cii).
We therefore hypothesize that the discrepancies in model behavior when different WTG schemes are used can be attributed to the differences in treatment of the baroclinic modes between these schemes.If we modify the TGR and SPC implementations such that the response strength of higher baroclinic modes is reduced, the multiple-equilibria regime may appear.To test this hypothesis, we modified the DGW and TGR implementations such that only the response of the first two baroclinic modes impact the system (note that in such a case, the form of the TGR and SPC implementations would be the same), and calculated the WTG-induced vertical velocities for the DGW and TGR implementations respectively to be: where c 1 and c 2 vary vertical velocity associated with the first and second baroclinic modes to the first and second vertical modes of the temperature perturbation.
We vary different configurations of c 1 and c 2 as follows: Similar to Section 3.3, to obtain both wet-and dry-states of the multiple-equilibria regime, we perturbed the large-scale reference profiles, but this time by ±0.1 K.We used the idealized radiation scheme of Pauluis and Garner (2006), and plot the results for α = 10 and τ = 25 hr in Figure 4.As postulated above, the presence and strength of multiple-equilibria is indeed tied to the ratio of c r = c 2 /c 1 , with smaller values of c r resulting in stronger bifurcation into the wet and dry equilibrium states.When c 1 = 0, there is no bifurcation between wet and dry equilibrium states, nor any oscillatory behavior, even at much lower values of τ.
We also note the discrepancy when c 2 = 0.5, which is when the idealized TGR implementation is equivalent to the SPC implementation (if n = 2 in Equation 6 of the SPC implementation, see Section 2.3).In Figure 2 we see that the SPC implementation's multiple-equilibria regime is weaker than in Figure 4 for an equivalent τ.This is presumably due to the effect of higher-order baroclinic modes beyond the second-order.We are able to verify this by running a modified version of the SPC implementation where Equation 14 is modified to: and our results (Figure S4 in Supporting Information S1) show that the multiple-equilibria regime is now visible.
Lastly, in this paper we are focused on understanding the behaviors of different implementations of the WTG approximation under the assumption of Rayleigh damping.In observations, a warm anomaly aloft with an underlying cool anomaly are sometimes found to correspond to a more bottom-heavy convective mass flux (Fuchs-Stone et al., 2020;Raymond & Fuchs-Stone, 2021;Raymond et al., 2014Raymond et al., , 2015)), opposite to what is shown in Figure 3a.Whether this discrepancy is due to the assumption of Rayleigh friction or the assumption of the WTG balance warrants additional study, but this is beyond the scope of our paper.

Conclusions
Implementing different WTG schemes results in different model behavior, and these differences become more prominent under a simplified framework with idealized radiation.A multiple-equilibria regime appears when the DGW implementation is used, with persistent wet and dry states.When the WTG approximation is enhanced more strongly, the model transitions into a regime that oscillates between these wet and dry states.However, when the TGR and SPC schemes are implemented the multiple-equilibria regime either weakens or vanishes, and the oscillatory behavior only appears in the TGR scheme when the relaxation occurs over unrealistically short timescales (τ ∼ 0.1 hr).
We have shown that these discrepancies can be attributed to their different treatments of higher-order baroclinic modes.Specifically, WTG schemes with stronger higher-order baroclinic modes reduce the likelihood of the multiple-equilibria and oscillatory regimes appearing.We can understand these differences in the GMS framework, specifically in reference to how Inoue and Back (2017) characterized GMS as a measure of feedback effects to convection.By approximating GMS as the ratio of export of moist static energy to that of dry static energy (Equation 15, see also Raymond et al. (2009); Kuang (2011); Inoue and Back (2015)), we see that the choice of WTG scheme will play a significant role in the GMS of the system, particularly because the response of vertical velocity to buoyancy perturbations of the different baroclinic modes are treated differently.
As we first touched upon in our introduction, while some work has gone into quantifying the discrepancies in model results when different implementations are used (e.g., Daleu et al., 2015;Romps, 2012a;Romps, 2012b), less thought has been given to understanding why different implementations give rise to different results in the first place.We hope that this set of idealized model experiments begins to close the gap between quantifying and understanding the differences in model results when different WTG schemes are used.

Figure 1 .
Figure1.When (a) a large-domain simulation is run to RCE, the induced large-scale circulation causes self-aggregation of convection, resulting in the formation of (c) a dry, weakly/no-precipitating regimes with vertical subsidence and (d) moist, strongly precipitating regimes with vertical ascent.In (b) small-domain RCE runs, self-aggregation does not naturally occur, but previous studies have shown that implementations of the WTG approximation that parameterize the large-scale tropical circulation allow us to attain either of these two regimes.Here, w wtg and r denote the domain-mean large-scale vertical velocity and relative humidity respectively.

Figure 3 .
Figure 3.We plot an idealized profile of the (a) first two baroclinic modes of WTG-induced vertical velocity, (b) vertical profiles of (1) dry and moist static energy and (2) their vertical derivatives, and lastly (c) the product of the vertical derivatives of the static energies with the (1) first and (2) second vertical modes of vertical velocity.We see that the lateral export of moist and dry static energies are dominated by the second and first baroclinic modes respectively.

Figure 4 .
Figure 4. We show here how the strength of the bifurcation varies with the ratio of c r = c 2 /c 1 for the (a) DGW and (b) TGR implementations in experimental setups with idealized radiation.As c r decreases, the bifurcation between the wet-and dry-states of the multiple-equilibria regime increases in magnitude.