## 1 Introduction  
As a measure of turbulence exchanges in the atmospheric surface layer, atmospheric stability can significantly affect the wind turbine wake deficit and its recovery rate. In general, turbulence exchanges between the wake and atmosphere are depressed under stable conditions. High wake deficits and slow wake recovery, thus, are usually observed in the stable stratification boundary layer [1; 2]. As wakes play critical roles in wind farm energy production and the fatigue loads of wind turbines, there is an increasing interest in studying the effects of atmospheric stability on wakes and developing wake models for non-neutral conditions. The impact of atmospheric stability on wakes is widely observed in wind tunnel measurements of small-scale models of wind turbines [3; 4; 5] and full-scale field experiments [1; 2; 6; 7; 8]. According to the wind tunnel measurements in [PERSON] et al. [3], the stronger inlet velocity gradient in the stable case leads to a slightly stronger turbulence intensity and extends the region of enhanced turbulence intensity from a distance of approximately 4 to 5.5 rotor diameters to 3 and 6 rotor diameters downwind of the turbine location. In [PERSON] et al. [4], a 15% smaller velocity deficit at the wake center, a more rapid momentum recovery due to an enhanced radial momentum transport, and 20% higher peak turbulence intensity were observed in the unstable case, as compared to the wake of the same wind turbine under neutral conditions. According to the field experiments of wakes using LiDARs or masts, [PERSON] et al. [6], [PERSON] et al. [7], and [PERSON] et al. [2] also found that the velocity deficit decreases more slowly under stable conditions and more quickly under unstable conditions. Those observations suggest that atmospheric stability should be considered for improved wake models and predictions of wind power harvesting. Besides field measurements and wind tunnel experiments of wind turbine wakes under non-neutral conditions, numerical models based on the computational fluid dynamics (CFD) have also been used to investigate the effect of atmospheric stability on the wind turbine wakes [1; 9; 10; 11; 12; 13; 14]. High-fidelity large eddy simulations (LES) with the actuator line model [9; 10] or the actuator disk model [1; 11] are commonly used to study the structure and dynamics of wind turbine wakes in varying stability cases. This is mainly because the LES approaches allow capturing the near wake structure and resolving the interactions of tip vortices with large-scale eddies of the ambient flow and wake meandering. However, the high-fidelity approaches are computationally expensive for wind energy engineering applications. There are some efforts towards developing turbulence models for wind turbine wakes under non-neutral conditions in reynolds-averaged navier-Stokes (RANS) to reduce computational costs. A widely used turbulence model for the thermal stratified boundary layer was developed by [PERSON] and [PERSON] [15], where a coefficient of the buoyant terms in the transport equation of the turbulence kinetic energy dissipation is calibrated with atmospheric stability. However, [PERSON] et al. [16] showed that this model cannot keep the flow homogeneity in a large domain under unstable conditions and thus developed a turbulence model consistent with monin-obukhov similarity theory (MOST). In [13], an additional buoyancy production based on the Richardson number, \(Ri\), was added to the turbulent kinetic energy equation to model wind turbine wakes under stable conditions. [PERSON] et al. [14] further applied this model with an additionally dissipation source term to investigate the wake behavior at different atmospheric stability conditions. Note that these turbulence models are based on MOST and could fail in modeling wind turbine wakes when MOST breaks down under stable conditions for the flux Richardson number \(R_{f}>0.25\).