Evolution of forced shear flows in polytropic atmospheres: A comparison of forcing methods and energetics

Shear flows are ubiquitous in astrophysical objects including planetary and stellar interiors, where their dynamics can have significant impact on thermo-chemical processes. Investigating the complex dynamics of shear flows requires numerical calculations that provide a long time evolution of the system. To achieve a sufficiently long lifetime in a local numerical model the system has to be forced externally. However, at present, there exist several different forcing methods to sustain large-scale shear flows in local models. In this paper we examine and compare various methods used in the literature in order to resolve their respective applicability and limitations. These techniques are compared during the exponential growth phase of a shear flow instability, such as the Kelvin-Helmholtz (KH) instability, and some are examined during the subsequent non-linear evolution. A linear stability analysis provides reference for the growth rate of the most unstable modes in the system and a detailed analysis of the energetics provides a comprehensive understanding of the energy exchange during the system's evolution. Finally, we discuss the pros and cons of each forcing method and their relation with natural mechanisms generating shear flows.


INTRODUCTION
The relative difficulty of observing most astrophysical shear regions, such as those in the Sun, in detail makes it imperative to use analytical and numerical techniques to shed light on the motions present there. Global-scale numerical calculations of stellar interior dynamics is one approach to investigate the mechanisms maintaining differential rotation (Brun & Toomre 2002;Miesch et al. 2008). However, by using a global approach, such models can not resolve a large range of length-scales in its entirety and have to rely on artificially large transport coefficients or subgrid-scale models. Therefore, numerical investigations using a local approach, where only a small fraction of the object is simulated, can help to provide a more detailed description of the region of interest.
Previous studies of astrophysical flows looked at an assortment of different velocity profiles, depending on the problem and choice of boundary conditions. In the case of Keplerian motion, which is investigated in the context of accretion discs, usually a linear velocity profile is assumed (e.g. Brandenburg et al. 1995;Hawley, Balbus & Winters 1999;Dubrulle et al. 2005;Silvers 2008). This type of shear profile does not require an external force to balance viscous E-mail:Veronika. Witzke.1@city.ac.uk (VW); Lara.Silvers.1@city.ac.uk (LJS); Favier@irphe.univ-mrs.fr (BF) dissipation and instead they incorporate the velocity via a shearingbox approach (Goldreich & Lynden-Bell 1965;Narayan, Goldreich & Goodman 1987), where the velocity is instantaneously present. In contrast, some investigations of stellar shear flows have used polynomial functions, as for example in Tobias & Hughes (2004) and Cline, Brummell & Cattaneo (2003a) while other investigations have utilized trigonometric functions to model the velocity field (see, for example, Hughes & Proctor 2013;Cline, Brummell & Cattaneo 2003b). Such velocity profiles have a non-vanishing gradient at the boundaries. In order to minimise the effect of the boundaries on the shear layer a hyperbolic tangent profile can be used (see for example Brüggen & Hillebrandt 2001;Hughes & Tobias 2001;Vasil & Brummell 2008). Hyperbolic tangent profiles are commonly used in classical studies of Kelvin-Helmholtz instability and turbulence. However, most local numerical studies of the turbulence transition in shear flows take the approach of an unforced flow (Caulfield & Peltier 2000;Smyth & Moum 2000;Smyth & Winters 2003), which results in a finite lifetime of an initially unstable background state due to its inevitable viscous decay. However, astrophysical shear flows can be either transient features or be sustained over very long time-scales, where the physical mechanism maintaining the shear flow is usually unknown. Incorporating a forcing into the numerical model has therefore two roles first to sustain the initial state, for which a linear stability analysis can be carried out, and second to model the unknown physical processes responsible for the resulting flow in an astrophysical system. A variety of classical studies of shear driven turbulence exploit a method where a decoupled background shear flow is present (for example used by Holt, Koseff & Ferziger 1992;Jacobitz, Sarkar & van Atta 1997;Barker et al. 2012). This method requires a change of variables to incorporate a mean shear profile and does not allow for a back-reaction of the actual flow on the forcing. Whereas investigations of astrophysical shear flows exploit different methods to provide a sustained flow. For example in Miesch (2003), Vasil & Brummell (2008), and  a method to balance viscous dissipation by introducing an external force proportional to the viscous term is utilized. Another option that has been selected is the relaxation method (e.g. Prat & Lignières 2013), which incorporates an external force proportional to the difference between the actual velocity profile and the target shear profile. Furthermore, studies focusing on magnetohydrodynamical instabilities used forced shear flows in order to study magnetic buoyancy and its relevance to the formation of sunspots (Tobias & Hughes 2004;Brummell, Cline & Cattaneo 2002). To understand the role of a shear flow for magnetic field generation, several investigations on the interaction between a shear flow and initial weak structured magnetic fields have been conducted (Vasil & Brummell 2008;Miesch, Gilman & Dikpati 2007;Heifetz et al. 2015). In all cases however, the influence of the forcing method used on the evolution of the system has not been studied. In this paper, a comparative analysis of different forcing methods is performed to understand how different forcings affect the nonlinear evolution of the KH instability. The governing equations are given in Sec. 2 along with the formulation of the forcing methods and numerical methods used. Our comparison of the different forcing methods is presented in Sec 3, where the exponential growth regime is compared in Sec. 3.1 and the non-linear phase is presented in Sec. 3.2.

THREE-DIMENSIONAL MODEL
Although a linear stability analysis is a powerful tool to investigate possible instabilities of shear flows in stellar environment as previously studied in Witzke et al. (2015), understanding the complex dynamics requires three-dimensional non-linear calculations. In order set up a system that initially remains comparable to a linear stability problem, which has a non-evolving background state, external forcing is needed. The force aims to maintain the initial conditions corresponding to the equilibrium state as long as nonlinear effects are negligible. Our purpose is to investigate whether different forcing methods provide a temporally evolving system, which shows the predicted linear evolution, and in what respect the non-linear evolution depend on the method used.

Governing equations, boundary conditions, and background state
We consider a three-dimensional domain of depth d, bounded by two horizontal planes located at z = 0 and z = 1, and periodic in both horizontal directions. The fluid is assumed to be an ideal monatomic gas with the adiabatic index γ = c p /c v = 5/3 and constant dynamic viscosity µ, constant thermal conductivity κ, constant heat capacities c p at constant pressure, and c v at constant volume. The set of dimensionless differential equations we consider is: where ρ is the density, u the velocity field, T the temperature, θ denotes the uniform temperature gradient across the layer, and p is the pressure. In the dimensionless equations above, all lengths are given in units of the domain depth d. The temperature and density are recast in units of T t and ρ t , the temperature and density at the top of the layer, and we take the sound-crossing time, which is given byt = d/[(c p − c v )T t ] 1/2 , as the reference time. There are two dimensionless numbers in the set of equations above: the Prandtl number, σ = µc p /κ, which is the ratio of viscosity to thermal conductivity and the thermal dissipation parameter C k = κτ/(ρ t c p d).
The strain rate tensor in equation (3) has the form A force term F in equation (2) aims to model external forces resulting from large-scale global effects (such as Reynold stresses associated with thermal convection in global-scale calculations for example) not included in our local approach. This force can sustain a shear flow when needed and is set to zero otherwise. The different forcings that we will consider are described in Section 2.2. For the basic state a polytropic relation between pressure and density is taken. Due to the Schwarzschild criterion the fluid is stable against convection if the inequality m > 1/(γ − 1) = 1.5 holds. In this paper the polytropic index m is always chosen such that the atmosphere is stably stratified. The boundary conditions at the top and the bottom of the domain are impermeable and stress-free velocity and fixed temperature: T = 1 at z = 0 and T = 1 + θ at z = 1.
The dimensionless initial temperature and density profiles are of the form: This basic state corresponds to an equilibrium state if the fluid is at rest. In order to obtain a shear driven turbulent regime it is necessary to start with an unstable velocity profile. A hyperbolic tangent profile is assumed in order to model a localised shear layer in the middle or our domain that will minimize the effects of the boundaries. Therefore, we assume that an external force sustains the following initial background velocity profile with a shear amplitude U 0 and a scaling factor 1/L u that controls the width of the shear profile. The boundary conditions introduced in equation (5) restrict the shear profile to values of L u which will result in a low enough value of the z-derivative at the boundaries. Using this velocity profile the above basic state can still be regarded as an equilibrium state if viscosity is neglected. Some shear profiles for different widths are illustrated in Fig. 1. Our calculations are initialized by adding a small random temperature perturbation to the equilibrium state including the additional shear flow given by equation (9). In order to evolve the system in time equations (1) -(3) are solved by using a hybrid finite-difference/pseudo-spectral code (see for example Matthews, Proctor & Weiss 1995;Silvers, Bushby & Proctor 2009;Bushby 2012, 2013, andreferences therein). For the linear stability calculations, the eigenvalue-problem formulated in Witzke et al. (2015) is numerically solved on a one-dimensional grid in the zdirection that is discretised uniformly, and this method is adapted from . When using direct numerical calculations to solve the system numerically over considerable iterations an issue concerning the momentum conversation might appear, which is specific to the choice of stress-free boundary conditions (see discussion in Jones et al. 2011). Despite the fact that equation (2) together with the boundary conditions in equation (5) conserves the momentum a cumulative effect of truncation errors at each time step might lead to an unphysical change in momentum when integrating over a vast number of time steps. We have checked that both momentum and mass is indeed conserved for all the calculations presented in this paper.

Forcing methods
We will compare three different methods to sustain an initial shear flow where we distinguish between methods that are static, i.e. the force term does not change throughout the calculation, and dynamic forcing methods that might vary depending on the current flow. Furthermore, methods with a force applied to a localised region have a local force whereas for global forcing methods the force term applies on the whole domain. The main goal is to find a forc-ing method that does not significantly alter the characteristics of the linear phase of the evolution but allows to reach a quasi-steady non-linear state.

Viscous method
In order to balance the viscous dissipation associated with the initial shear flow profile given by equation (9) the force is added to the RHS of equation (2). By applying this viscous forcing the initial state given by equations (7) -(9) is in equilibrium provided that viscous heating is neglected. This method has been broadly applied in forced shear flows to model the dynamics of the solar tachocline (e.g. Miesch 2003;). Note that this method only balances for the viscous diffusion of momentum associated with the target profile and does not depend on the actual non-linear solution. In that sense, the forcing can be considered to be local and static.

Perturbation method
Our second method, the perturbation method, was previously used by Holt et al. (1992); Jacobitz et al. (1997); Barker et al. (2012). For this method a slightly different set of differential equations is solved. A decomposition of the velocity, u = U 0 +ũ into a background shear flow, U 0 , and the deviation from the background profileũ enables the maintenance of a shear flow that is independent of the unstable perturbations. Note, the background velocity profile U 0 has to be time independent and divergence free. Thus, inserting the above decomposition into the momentum equation (2) we obtain since we assume that our background profile does not vary with time and in the flow direction, the time derivative of U 0 and the term U 0 ·∇U 0 vanish. In order to ensure that the background velocity is not dissipated by viscosity the viscous term associated with it is dropped. The equation takes the form such that effectively the differential equation for the velocity perturbations only is solved. Here, the velocity perturbations,ũ, are initially zero, but the background velocity is incorporated through the two advective terms involving U 0 . By the same procedure the equations (1) and (3) become: where the effect of the background velocity, U 0 , on the density and temperature is taken into account. Mathematically, the equations for U 0 +ũ are the same as for u in the viscous method. Therefore, Table 1. Comparing the linear eigenvalue-solver results with those from non-linear calculations during the linear phase. For case I the shear amplitude is U 0 = 0.08 and 1/L u = 118 such that Ri = 8 × 10 −4 . Taking C k = 8 × 10 −5 results in a Pe = 34. For case II U 0 = 0.041, 1/L u = 20 and C k = 1 × 10 −4 such that Ri = 0.1 and Pe = 82.

Relaxation method
In the relaxation method, an external force ensures that the flow relaxes towards the initial profile on an arbitrary time-scale τ 0 . We define any quantityf as where the overbar denotes that the quantity f is horizontally averaged, and N x and N y are the resolutions in x-direction and in ydirection respectively. The force for the relaxation method in the momentum equation depends on the horizontally averaged velocitȳ u x (z) and is given by whereê x is the unit vector in x-direction. For this method it is crucial to ensure that the forcing is aligned with the flow direction and does not depend on the velocity variation along this horizontal direction, which would otherwise correspond to a local small-scale force. A local forcing should be avoided, because it will suppress any instability and can lead to non-physical behaviour. The relaxation method was used by Prat & Lignières (2013) to model shear driven turbulence and provides the advantage of an adjustable back reaction on the actual mean flow. It is a global forcing (due to the averaged operator) and a dynamical forcing, because it depends on the actual flow.

COMPARISON OF THE FORCING METHODS
Investigating saturated flows on long time-scales requires an external force to sustain a target velocity profile. Here we compare the effects of the different forcings presented in Section 2.2 on the development of the shear instability. Our investigation is divided into two parts. We compare the linear regime of a shear instability in a two-dimensional framework in Section 3.1 and the nonlinear regime is investigated in Section 3.2 using three-dimensional calculations. Here we first qualitatively examine two-dimensional slices and three-dimensional renderings of key quantities for different forcing methods in Section 3.2.1, before we calculate the horizontally averaged velocity, turbulent Reynolds numbers and turbulent length in Section 3.2.2. A detailed analysis of the energy budgets is provided in Section 3.2.4, where the theoretical framework for the energy budgets is introduced in Section 3.2.3. The relation between the external work done by the forcing and the amount of dissipated energy by viscosity is analysed in detail in Section 3.2.5 and finally we summarize our findings for the saturated regime in Section 3.2.6.

Linear regime
Since the initial linear phase of shear flow instabilities is purely two-dimensional, which becomes evident from Squire's theorem (Squire 1933), we focus here on calculations in a two-dimensional domain, which has a spatial resolution of N x = 512 and N z = 480. The stability of a shear flow in a stratified atmosphere is characterized by the non-dimensional Richardson number, Ri. Using the general definition of the Brunt-Väisälä frequency given by whereT = T (P t /P) 1−1/γ is the potential temperature, the minimum value of the Richardson number, Ri, across the layer is defined as where the derivative of the background velocity profile, defined in equation (9), with respect to z corresponds to a local turnover rate of the shear. In most cases the minimum Ri value is at z = 0.5, but for some parameter choices with large temperature gradient, θ, and broad shear width the minimum is shifted towards greater z.
Here we consider unstable shear flows with a Richardson number less than 1/4 at a point in the domain. We do not consider shear instabilities triggered by thermal diffusion for which larger values of Ri can be used (Dudis 1974;Zahn 1974;Lignières et al. 1999). Because the 1/4 criterion is a necessary, but not sufficient, requirement for instability we also solve the corresponding linear stability problem based on the approach used in Witzke et al. (2015) in addition to conducting the non-linear calculations. For simplicity, the Prandtl number is fixed to be unity whereas the dimensionless thermal diffusivity C k is varied from 10 −4 to 10 −5 . Taking the previous linear study by Witzke et al. (2015) into account, our parameters satisfy the following requirements: To ensure a stable stratification the polytropic index is set to be m = 1.6, the amplitude U 0 of the shear flow is chosen such that the Mach number in the middle of the domain remains less than 0.08, which avoids additional stabilisation by compressible effects. Furthermore, we take the initial Péclet number, which we define as at City University, London on August 16, 2016 http://mnras.oxfordjournals.org/ Downloaded from to be much greater than unity to avoid a destabilizing effect caused by thermal diffusion. Note, that due to the definition of U 0 and L u in equation (9) a factor of 4 is needed to be consistent with the generally used definition.
Here we want to compare two linearly unstable cases with distinct behaviour when no external forcing is included. An unforced case will decay if the instability grows on a larger time-scale than the viscous dissipation time-scale. Therefore, we consider two different cases corresponding to two different minimum Ri values. In case I we consider a very small value of Ri = 8 × 10 −4 , such that the system is unstable even without external forcing. Case II has a greater Ri = 0.1 and greater C k so that the system is unstable only if we introduce an external forcing. Without forcing, the initial shear flow diffuses quickly such that the Richardson number increases rapidly above 1/4 and no shear instability can be sustained. By considering these two different cases, we investigate how the forcing method used affects the development of a shear instability during the exponential growth phase. Unforced calculations provide a reference for the system's evolution without any external forces. All parameters for case I and case II, the resulting linear growth rates, and the wave number for the most unstable mode for each method are summarised in Table 1 using both an eigenvalue solver and direct numerical calculations of unforced and forced cases. The growth rates for the non-linear calculations are obtained by calculating where the angle brackets denotes that any quantity f is volume averaged as follows In addition, we fit horizontally-averaged velocity profiles in the shear direction with a hyperbolic tangent profile as in equation (9) to estimate the shear width during the exponential regime. In order to find the most unstable wave number, which corresponds to the wave number with the most energy, the kinetic energy spectrum is calculated as where the hat symbol denotes the Fourier transform of the corresponding quantity and the star symbol denotes the complex conjugate. The growth rates for the viscous forcing, the perturbation method and the relaxation method (where τ 0 = 0.01) are almost identical for both case I and II. The growth rates achieved by these methods in case II correspond to the growth rate calculated by using the EV-solver with a 12% relative error when taking the growth rate from the EV-solver as reference. For case I the error is 8%. The most unstable mode, k max , in non-linear calculations is the same for the viscous-and the perturbation method. Note that the most unstable wave number obtained by DNS is always slightly smaller than the one obtained from the eigenvalue solver, and more so for case I than for case II. This is due to a thinner shear width in case I, which is more affected by viscous dissipation, such that the instability is triggered when the shear profile is significantly broader. Therefore, k max deviates for case I more than for case II, where the initial shear width is broader. We now look at the effect of varying the relaxation time-scale τ 0 in the relaxation method. The instability develops always in the mid-dle of the domain and is visually similar to the instability observed by using either the viscous-or the perturbation method. The evolution of < w >, for different τ 0 parameters, for case I and case II is shown in Fig. 2 and the growth rates of these runs are summarised in Table 1. For all cases the onset of the instability is not sensible to the chosen relaxation time-scale τ 0 , but the growth rate decreases with increasing τ 0 . This is expected since a smaller τ 0 implies a larger restoring force as soon as the averaged velocity profile differs from the target equation (9). Therefore, for relaxation times that are larger the viscous dissipation might be unbalanced such that the initial state changes before an instability is triggered. The relaxation method leads either to the same instability or sustains an instability triggered by a smoother velocity profile, when the relaxation time τ 0 is increased. To put τ 0 in relation with typical dynamical times the initial turnover time of the shear flow t s = L u /U 0 is calculated, which is t s ≈ 0.1 for case I and t s ≈ 1.2 for case II. Furthermore, the initial viscous time-scale, t µ = L 2 u /µ, that accounts for the time on which the initially shear width is dissipated, gives another reference time. For case I the viscous timescale is t µ = 0.9 and for case II it is t µ = 25. Note, that both timescales represent initial time-scales given by the initial configuration and will change with time as the system evolves, especially in the saturated regime these times might be significant different from the initial values. So that a relaxation time greater than both the viscous and dynamical time-scale corresponds to a very weak back-reaction of the forcing to the change of the averaged velocity. Note that τ 0 < 1 means that the forcing relaxation is quicker than a sound crossing time which is unlikely to occur in physical system. By varying the relaxation times, τ 0 , we investigate how to chose an appropriate τ 0 with respect to the initial time-scales in order to recover the linear instability dictated by the initial state and at the same time a saturated regime that evolves towards a quasi-static state. As expected the numerical calculations using the viscous method and perturbation method lead to exactly the same result, because both methods are mathematically equivalent. Note however that the computational cost is slightly greater when using the perturbation method. In conclusion the viscous and perturbation method are the same, but the relaxation method shows different dynamics depending on the relaxation time. If we are to understand which of these forcing methods is most suitable to model shear flows in stellar interiors it is essential to conduct a comparison for the non-linear evolution.

Non-linear phase
In order to investigate the non-linear evolution of a stratified shear flow a three-dimensional domain is crucial (Thorpe 1987), because after a Kelvin-Helmholtz instability three-dimensional instabilities are triggered (Peltier & Caulfield 2003). These secondary instabilities lead to a turbulent collapse of horizontal vortices, which is suppressed in a two-dimensional setup as discussed by Scinocca (1995). Therefore, to capture the effect of different forcing methods on the entire dynamics we now focus on three-dimensional calculations, which have a spatial resolution of N x = 256, N y = 256 and N z = 360 . To avoid confinement effects associated with the upper and lower boundaries, we focus on a case with a temperature gradient θ = 5. In this case the instability is more likely to remain confined in a narrow central region as to minimize the importance of our particular choice of boundaries. The polytropic index is kept the same as in the previous section to ensure a stable stratification. A  (21), for the two-dimensional calculations for case I in (a) and for case II in (b). Different τ 0 parameters for the relaxation method are used and compared to the viscous method and unforced calculations. The vertical velocity is displayed in logarithmic scale and t is given soundcrossing time.
parameter search was conducted to find combinations of the other parameters that lead to a finite spread of the unstable region. The dynamical viscosity is taken of order 10 −4 and the Prandtl number, σ = 0.1, because low Prandtl number flows are more relevant for stellar interiors. The shear flow amplitude is set to U 0 = 0.2 and we take 1/L u = 80 such that Ri min = 0.003. As discussed above, the viscous method and the perturbation method are equivalent, such that in the following, only results obtained with the viscous method are shown but we have checked that the calculations are indeed the same when using the perturbation method. Our study compares the resulting non-linear dynamics obtained by using the viscous method and the relaxation method. Furthermore, we discuss the effect of varying τ 0 in the relaxation method and how does it compare with the viscous method, taken as reference.

Visualisation
In order to compare the flow evolution we start with a visualization of the vorticity at three different stages during the non-linear saturation. The first stage is the exponential growth phase, during which all forcing methods show a similar evolution, where the layer with non-zero vorticity spreads vertically. For the viscous forcing  and the relaxation method with τ 0 1, this layer remains significantly smaller than for calculations where the relaxation method with τ 0 ∼ O(1) is used. Note that the the viscous time-scale for this case is t µ ≈ 1.5. The second stage is chosen at a point when the instability starts to saturate and the fluid parcels overturn. In Fig. 3 the vorticity component perpendicular to the x-z plane for the viscous forcing, the relaxation method with τ 0 = 0.1 and τ 0 = 10 for the second and third stage are plotted in the x-z plane at y=0.8. For the second stage the dynamics differ between the different forcing methods, which can be seen in Fig. 3 (a), (c) and (e). In Fig. 3 (a) small patches of strong positive vorticity are merging together into each other along a thin horizontal layer and a few small negative vorticity patches are present. In comparison, when using the relaxation method with τ 0 = 10 large billows of smaller positive vorticity occupy a horizon-tal layer which is more extended in the vertical direction, see Fig.  3 (c). This can be explained by the smoother shear width, which is a consequence of the slow back-reaction of the forcing. Using a smaller τ 0 leads to a greater vorticity amplitude than achieved by the viscous method (see Fig. 3 (a) and (e), note the different color scales) while the spread of the instability remains similar. The third stage for the different methods is several sound crossing times after saturation, where the flow is evolving towards a quasistatic state. Comparing Fig. 3 (b), where the viscous method is used with Fig. 3 (d), and Fig. 3 (e), where the relaxation method with τ 0 = 10 and τ 0 = 0.1 is used respectively, the main differences are the vertical extent of the overturning region and the amplitude of the vorticity. Using a larger τ 0 = 10 leads to a similar situation as for the viscous method which becomes evident when comparing Fig. 3 (b) and (d). In Fig. 3 (b) the layer is thin and shows elongated at City University, London on August 16, 2016 http://mnras.oxfordjournals.org/ Downloaded from regions of strong positive vorticity, whereas in Fig. 3 (d) the region is significantly extended with small patches of strong positive and negative vorticity. A few larger regions are present further away from the middle of the domain. For the relaxation method with τ 0 = 0.1 a drastically different behaviour is observed, see Fig. 3 (f), where the vorticity amplitude is much greater and the region where overturning is present is taking up almost the entire domain. In addition, much more small scale turbulence occurs around z = 0.5. This can be explained by the form of the forcing used: The viscous method acts with a force, that is confined in a narrow region around the middle of the domain such that the instability can develop freely further away from z = 0.5. The relaxation method drives the fluid towards the target profile even far away from the middle of the domain. This supports additional spread of the shear instability and triggers more turbulent motion. While it is worth mentioning that for the relaxation method with τ 0 = 0.1 due to strong viscous heating convective motion is present in the upper half of the domain (which we checked by calculating the Brunt-Väisälä frequency), this dynamics will not be further discussed, as it is the result of the artificially low value of τ 0 used in that case. Furthermore, we expect this unrealistic instability to disappear at lower Prandtl numbers. Fig. 4 shows contour plots of the vertical velocity in three dimensions at approximately 40 sound-crossing times, which is well after the non-linear saturation. For the viscous method the patches of downwards and upwards motion form an alternating pattern along the x-direction where regions of the same velocity are arranged in small tubes that are extended along the y-direction, see Fig. 4 (a). Such a pattern is not present in Fig. 4 (b), where the relaxation method with τ 0 = 10 is used. Here, the regions of the same velocity form larger patches which are extended along the x-direction and changes sign along the y-direction. This indicates that a secondary instability which forms overturning billows along the y-direction is more dominating when the relaxation method with larger values of τ 0 is used. For the relaxation method with τ 0 = 0.1, displayed in Fig. 4 (c), the artificially strong forcing leads to an intense forward energy cascade and associated small-scale turbulent motions in the middle of the layer. The large-scale structures observed in the upper part of the domain are convective cells resulting from the large central viscous heating. When looking at the time evolution along Fig. 3 (a), (b) and the evolution along Fig. 3 (e), (f) a significant difference in the amplitude of the vorticity can be noticed. While for the viscous method the amplitude increases towards a peak during the saturation, Fig.  3 (a), and start to decrease afterwards, for the relaxation method the amplitude of the vorticity constantly increases and reaches a maximum after saturation, see Fig. 3 (b) and (f). This might be explained by the form of the forcing: Because the relaxation method adjusts the magnitude of the force according to the mean flow, the strength of the force increases constantly and lead to more overturning with time, whereas the force remains constant when using the viscous method, such that the overturning settles down after saturation. However, in order to check if this is indeed the case a more detailed analysis on the work done by the force and the total viscous dissipation is required, which is discussed in Section 3.2.4. Having compared the non-linear evolution for different forcing methods qualitatively we can conclude the following. The viscous method and the relaxation method with τ 0 ∼ O(t µ ) or larger result in similar, but still distinct, evolutions. Using a relaxation time τ 0 t µ but still greater than the dynamic time-scale t s ≈ 0.06 leads to a very different non-linear dynamics with great mixing and possible non-physical behaviour leading to convection. Therefore, we conclude that, for the saturated regime, the initial viscous time- scale, here t µ ≈ 1.5, gives a reference time for the choice of τ 0 and the case with τ 0 = 0.1 will be excluded from further analysis.

Horizontally averaged profiles
The system under consideration is stratified, such that most quantities will change with depth, z, throughout the domain. Therefore, investigating horizontally averaged profiles with depth provides further insight in the system dynamics during the saturated regime. Let us first define a local turbulent Reynolds number as follows whereρ(z) is the horizontally averaged density. Here u rms (z) is the root mean square of the fluctuating velocity averaged over the horizontal layers calculated as follows where U 0 (z) is the target velocity profile as defined in equation (9). The turbulent length-scale is taken as where k 2 = k 2 x + k 2 y is the horizontal wave number and the energy spectrum E(k, z) is averaged over horizontal layers such that it takes the form E(k, z) = 1 4 kû (k x , k y , z)ρu * (k x , k y , z)+ρu(k x , k y , z)û * (k x , k y , z). The resulting Reynolds number variations with depth for the four different runs after 60 sound-crossing times is shown in Fig. 5 (a). It becomes evident that the flow for all four cases is very different at the late stage of the calculation. The viscous method is characterised by two regions above and below the middle of the domain, where the flow has high Reynolds numbers about 850 and 1300 respectively. The asymmetry is due to significantly denser fluid at the bottom of the domain. For the relaxation method with τ 0 = 10, these two regions are narrower and the maximal Reynold numbers are approximately 800 and 900. For τ 0 = 1.0 the region of high Re t is spread with two peaks very close to z = 0.5. For these three methods a moderate turbulent flow confined in the middle is obtained. The corresponding turbulent length-scales, shown in Fig. 5 (b), reveal that for the three methods used the length-scales become less than 0.5 around z = 0.5 and increase towrads the boundaries. Using the viscous method leads for this particular case to smaller turbulent length scales than using a relaxation method. We now consider the horizontally averaged velocity profileū x (z) shown in Fig. 6 at a time,t, where the system evolves towards a quasi-static state. For all methods the shear profile is smoothed out. For the relaxation methodū x (z) remains a hyperbolic tangent profile, but for the viscous method a steep transition occurs around z = 0.5, which merges into a smoother region at the boundaries. This is caused by the different type of forcing: The force applied in the viscous method is solely defined by the shape of the target profile. Since the initial target profile has a thin width, the second derivative is large in this region, which causes a stronger forcing, compared to the outer parts of the domain. In contrast the relaxation method applies a force that depends on the deviation of the actual mean profile from the target profile, such a back reaction ensures the preservation of a hyperbolic tangent profile.

Theoretical framework for energy budgets
To get a more comprehensive insight on the processes involved when a linear shear flow instability saturates it is useful to track the evolution of the standard forms of energy during the transition phase and beyond. Following the same procedure as used in Landau & Lifshitz (1987) and Griffies (2004) we decompose the total energy into the kinetic energy, E kin , internal energy, I, and gravitational potential energy, E pot , which in our case are given as where the volume integral is taken over the whole domain and the internal energy for an ideal gas is considered. Then, in order to understand the energy evolution with time, and to get more detailed insights into energy budgets, it is useful to derive the energy changes in our system using equations (1) -(3). The time derivative of the kinetic energy becomes where S denotes the volume surface, ε is the positive viscous dissipation rate, H a is the change rate due to advection, H p is the rate of work done by expansion and contraction, H ρ denotes the exchange rate with the potential energy due to density flux and W is the work done by external forcing. For the rate of change in the internal energy we get where q = C k σ(γ − 1)|τ| 2 /2ρ such that ε is due to viscous heating and Φ temp corresponds to heat loss or gain at the surface, S . In the standard form, the changes in the gravitational potential energy are only due to the exchange of density flux as can be seen by taking the time derivative Summing equations (30) -(32) yields the total energy change of the system  Method with τ 0 = 1 with τ 0 = 10 before instability 0 < t < 2.5 0 < t < 3.5 0 < t < 4.5 exponential growth 2.5 < t < 4.5 3.5 < t < 6.5 4.5 < t < 10 saturation 4.5 < t < 23 6.5 < t < 30 10 < t < 30 quasi-steady state t > 23 t > 30 t > 30 which is only due to external forces, heat loss or gain at the surfaces, and advection. Note, that H a is negligible in our case, because our system is closed and mass is conserved. It can be seen that viscous dissipation, ε, density flux,H ρ , and work done by expansion and contraction, H p are exchanged between the kinetic energy and the potential energies. However, using the standard decomposition of energies it remains impossible to distinguish between reversible and irreversible energy exchange between these three energy budgets. In order to resolve this issue Winters et al. (1995) introduced a method to analyse the mixing behaviour of a turbulent flow, which can be used to track reversible and irreversible changes of different potential energies. This framework was further extended to compressible fluids by Tailleux (2013). For this method we decompose the gravitational potential energy of the system defined in equation (29) into two parts. One part is the so-called background potential energy defined as where the ρ is the adiabatically redistributed density. This definition is also appropriate for a compressible fluid. Another part is the available potential energy which corresponds to the difference between the actual potential energy E pot and E back . The available potential energy can be transformed into other types of energies, whereas the background energy can not be accessed and transformed in other types of energies. Therefore, the background potential energy can be interpreted as the part of the total gravitational potential energy that corresponds to the lowest energy level that can be reached if the system is adiabatically redistributed. While a system is evolving the background potential energy can be only changed by irreversible processes. In our numerical calculations the background potential energy is obtained by taking the actual density distribution and sorting it in an ascending order, such that the highest density is at the bottom of the domain. In a similar procedure the internal energy can be decomposed into a background internal energy budget and an available internal energy budget, for a more detailed discussion see Tailleux (2013). However, for our purpose it is sufficient to analyse only the budgets for the gravitational potential energy in order to understand the mixing behaviour of the system.

Energy budgets from numerical calculations
The Kelvin-Helmholtz instability converts the kinetic energy that is available to the base horizontal shear flow into vertical fluctuations that need to overcome the stably-stratified atmosphere. The gravitational potential energy of the system is changed during this process. In addition, after saturation the fluid starts to overturn, and is mixed, where irreversible processes change the potential energy. Here we want to investigate how the forcing contributes to the different forms of energy in the system. Using the separate components responsible for the change in different energy budgets presented above and tracking the changes of the kinetic energy, internal energy and different gravitational potential energy budgets, we will discuss how the system behaves for different forcings. Below we distinguish between four stages of the system's evolution that are: the time interval before the exponential growth of an instability, the exponential growth phase, the onset of saturation, and a fully saturated quasi-steady state. These stages are at different times for the three calculations that are discussed and can be found in Table 2.
Significant differences between the forcing methods become evident when following different energy budgets normalized by the initial value of the system's total energy with time as shown in Fig. 7. In general the sum of the three energies is increasing due to the external work done by the forcing. For the viscous forcing the kinetic energy remains almost constant until the instability starts to grow att ≈ 3, at that time E kin begins to decrease. The internal energy increases at an almost constant rate from the beginning of the calculation, which is due to viscous heating that extracts kinetic energy via ε in equation (30). As the kinetic energy remains almost constant in the beginning, it can be concluded that the amount of energy dissipated is fed back into the system due to external work, W. The decrease in the kinetic energy of the system during the exponential growth of the instability is a direct consequence of the Kelvin-Helmholtz instability extracting energy from the large-scale shear flow in order to overcome the potential energy associated with the stably-stratified atmosphere. This amount of energy is partly converted into vertical motion, which contributes to the kinetic energy, and partly exchanged into gravitational potential energy. The term in the energy change rate associated with this exchange is H ρ , which leads to a slight increase in E pot . Because the conversion of the mean-flow kinetic energy to vertical motion retains the energy in the kinetic energy budget, the decrease is small. During the saturation phase the negative rate of change in the kinetic energy grows. Whereas for the internal energy a plateau is present just after the exponential growth phase and this is followed by a steeper increase during the end of the saturation phase. The gravitational potential energy starts to increase faster from the beginning of the saturation phase. Comparing the gravitational potential energy and the internal energy shows that the changes of the internal energy are similar to the changes in the potential energy but with a time shift and a greater amplitude. The time delay is a direct consequence of the irreversible processes during the feedback of gravitational potential energy into kinetic energy, where previously kinetic energy is transformed reversibly and irreversibly into gravitational potential energy due to the term H ρ .
Eventually a quasi-static state is reached, where the kinetic energy eventually will decrease very little, but the potential energies will constantly increase due to viscous heating and irreversible mixing processes. Both processes extract kinetic energy due to ε, H ρ , and H ρ respectively, but only a part of ε is added to the system by the external force. Therefore, the system will always evolve very slowly, but remain statistically similar for a very long time. We now focus on the calculations where the relaxation method was used to sustain a shear flow. The time evolution of E kin for the relaxation method with τ 0 = 10 is similar to the viscous forcing. However, the early evolution is different because E kin decreases even before the instability starts to develop. This is expected since the kinetic energy initially contained in the initial shear flow is dissipated by viscosity over a short time-scale. Therefore, the exponential growth regime is shifted to later times, where a similar drop in E kin as was found in the viscous forcing case. In both cases, this reduction in kinetic energy corresponds to an increase in the potential energy in the system (see Fig. 7 (c)). In the non-linear regime the system also tends towards a quasi-steady state. Similar to the viscous forcing, for the relaxation method with τ 0 = 10 the background potential energy and internal energy increase slowly until the system start to saturate. During the saturation phase a steeper increase is present. Then, after several sound crossing times, both potential energies start to converge towards a constant small growth after the system saturated. This behaviour reveals that the energy induced by the forcing principally transfers into internal energy due to dissipation, but does not contribute to an increase in either kinetic energy or available potential energy for late time evolution. A calculation with a shorter relaxation time, τ 0 = 1.0, shows a different behaviour, where E kin increases with the onset of instability. This growth is due to the very intense external forcing present as soon as the averaged velocity profile deviates from the target profile, which is the case when the instability starts growing. For τ 0 = 1.0 this corresponds to a growth in the total kinetic energy over approximately 30 sound crossing times whereby E kin slowly oscillates around a fixed value. The case with τ 0 = 1.0 shows that a constantly large increase of the potential energy is present even for the saturated stage. Looking back at Fig. 7 (a) the kinetic energy converges towards a constant value at large times. This indicates again that the kinetic energy pumped into the system by external forcing, W, is used to balance viscous dissipation and partly converted into gravitational potential and internal energy. Here the amount of externally added energy is significantly greater than for the other two calculations. In contrast for the calculation with τ 0 = 0.1, which is not displayed here, the kinetic energy grows during the whole duration of the calculation. As discussed earlier, this growth in the kinetic energy of the system is in that particular case associated with a transition between a stably stratified atmosphere (the polytropic index is initially m = 1.6) and a convectively unstable atmosphere where large convective cells appear in the upper part of the domain (see Fig. 4 (c)). This transition is driven by the large viscous heating in the central shear region modifying the temperature profile and changing the sign of the Brunt-Väisälä frequency. While the interaction between a large-scale shear flow and thermal at City University, London on August 16, 2016 http://mnras.oxfordjournals.org/ Downloaded from convection is of interest (see for example Guerrero & Käpylä 2011;), this is beyond the scope of the current study. In Fig. 8 the evolution of the available gravitational potential energy is plotted for the three-dimensional calculations. This part of energy is due to the reversible part of H ρ in the energy equations. When looking at the available potential energy in Fig. 8, no available potential energy is present before the saturation of the shear instability for all cases. Such that the system is in a state of lowest possible potential energy. The E avail for the viscous method and relaxation method with τ 0 = 10 increases similarly, although the arch of E avail is shifted towards later times in the calculation with τ 0 = 10. During saturation the fluid is mixed mostly, which means that due to the onset of overturning the background density is modified. This is evident from the increase of E pot and E avail , which indicates reversible and irreversible mixing processes see Fig. 7 (b) and (c). After saturation less mixing occurs in the system. For the viscous method E avail converges towards zero for late times, whereas in the calculation using the relaxation method with τ 0 = 10 the available potential energy oscillates around a small value. Since available potential energy is directly related to mixing (Peltier & Caulfield 2003), the system evolves towards a state with little mixing. This means that, after a certain modification of the density profile, the overturning settles down and persists at a low level over a long period of time. In agreement with kinetic energy evolution for the relaxation method with τ 0 = 1.0, the available potential energy starts to growth more rapidly and the system seems to reach a type of quasi-static state at very late times. However, for both cases τ 0 = 1.0 and τ 0 = 10, with current limitation on numerical resources, it remains unclear if the available potential energy is saturated or will eventually decay. In order to clarify this, the calculations need to be evolved further. However, for the purposes of comparing the forcing methods in this paper it is immaterial. By using the relaxation method we can reach a long-lived state and different mixing behaviours exist depending on the relaxation time τ 0 , which persist sufficiently long to study long-time evolution of the generated turbulence.

Comparing total viscous dissipation and external work
To investigate how much of the energy induced into the system by forcing balances the viscous dissipation, which part remains as kinetic energy and what converts into potential energy, it is useful to study the work done by the forcing, given by W as well as the total viscous dissipation rate, ε, with time. These quantities can be found for all three calculations in Fig. 9. At the start of the calculation the viscous forcing will always almost exactly balance the viscous dissipation, because the velocity profile does not deviate from the target velocity such that the viscous force cancel the viscous dissipation exactly (see equation (10)). This is true until approximatelyt ≈ 5 when the instability starts to saturate. After saturation the work done by the viscous forcing is not sufficient to balance the additional dissipation associated with small-scale fluctuations in the system in that case. This is associated with a decrease in the total kinetic energy as already discussed previously. At late times the amount of ε converges towards the work done (see Fig. 9 (a)) and so the system evolves towards a quasi-static state, where a sustained turbulent flow is present. At the beginning of each calculation, using the relaxation method the work done by the forcing, W, is initially zero since the velocity profile exactly matches the target profile (see equation (16)). Therefore, depending on τ 0 , viscous dissipation is initially not balanced,  Figure 9. The evolution of total viscous dissipation rate of momentum, ε, and the work done by the forcing, W for the viscous method and the relaxation method are shown. In (a) the viscous forcing is used. The relaxation method with τ 0 = 10 is used in (b), and with τ 0 = 1.0 in (c).
as can be seen in Fig. 9 (b, c). As the initial shear flow diffuses, the associated dissipation decreases until it becomes equal to the external forcing leading to a quasi-steady state. For the case with τ 0 = 1.0 the force increases shortly after the start of the calculation such that a phase where the viscous dissipation is balanced is present before the shear flow instability start to saturate. After saturation the work done by the forcing is greater than ε, which explains the increase in kinetic energy noticed previously. Due to the long relaxation time, the case with τ 0 = 10 reveals a distinct behaviour, where W remains less than ε throughout the exponential growth phase after which it matches for a few soundcrossing times ε, see Fig. 9 (b). When the total viscous dissipation reaches a peak the work done remains insufficient to balance for the viscous dissipation. At large times when the system is evolving towards a quasi-static state the total viscous dissipation remains less than the energy input such that turbulence can be sustained.

Discussion
Our research has revealed several characteristics of the different forcing methods: The viscous method provides a well defined localised force, but without control on the resulting velocity profile of the saturated flow. The shear instability can freely develop further away from the middle of the domain, but no turbulent motion is sustained there. Therefore, modification of the background profiles are solely due to non-linear dynamics of the instability. This results in less control on the resulting averaged velocity profile further away from the middle layer. From energetic considerations it can be concluded that the additional energy, that is added to the system during the late time evolution by external forcing, approaches a constant value. This initially kinetic energy is mostly converted into potential energy via dissipation. Only a small part contributes to the turbulent dynamics by mixing processes.
On the contrary to the viscous method using the relaxation method provides control on the resulting velocity profile, because the force is proportional to the deviation of the horizontally averaged velocity. Therefore, the target profile is controlled even far away from z = 0.5. This corresponds to a global forcing, which can suppresses changes of the background velocity throughout the domain such that modifications of the background profile due to the shear instability are suppressed. This also can induce more mixing, which is not initially caused by the instability of the localised shear layer and therefore leads to non-physical behaviour further away from the middle plane. An additional parameter, the relaxation time τ 0 , provides control on the strength of the forcing. Investigating the energy evolution revealed that decreasing τ 0 results in significantly more kinetic energy induced into the system by the external force than viscous dissipation converts into internal energy. Therefore, turbulent motion is supported by the external forcing. For the choice of τ 0 two typical time-scales are of interest, the turnover time t s = L u /U 0 and the viscous time-scale t µ = L 2 u /µ. For τ 0 < t s an instability with almost the same properties as obtained by the EV-solver is triggered (see Section 3.1). However, from the energy evolution of the saturated flow it becomes obvious that only a τ 0 greater than or of the same order as the viscous time-scale, which is t µ ≈ 1.6 for this case, leads to a system which can reach a quasi-steady state. Therefore, having compared the non-linear evolution for the viscous forcing and the relaxation method, we conclude that both methods can be used depending on the properties of the dynamics that needs to be modelled. For investigations with focus on the resulting velocity profile during the saturated regime the relaxation method is more appropriate, where a careful choice of the relaxation time has to provide that no significant effects from the forcing can induce unphysical behaviour. This should be provided if τ 0 ∼ O(t µ ) or greater. On the other hand, if the non-linear evolution of a shear unstable flow is of interest, where the the mixing behaviour induced by the shear instability is the main focus, the viscous method provides a more appropriate forcing. Since the viscous method does not significantly affect the turbulent dynamics further away from the shear region, the turbulence induced by the instability can evolve freely.

CONCLUSIONS
Turbulent motions driven by a shear flow instability are subject to occur in a wide range of physical systems, where numerical calculations can provide a comprehensive insight to the physical processes. Here direct numerical calculations in two-and threedimensional Cartesian domains are used to analyse different forcing methods, which were exploited in the past to maintain a background shear flow. In order to determine how different forcing methods affect a saturated flow and to get a more thorough understanding on possible unphysical behaviour it is essential to investigate different forcings qualitatively, e.g. by tracking global quantities and horizontally averaged profiles of velocities. Testing three methods in the linear and non-linear regime reveals that two conceptually different methods, the viscous forcing and perturbation method, result in exactly the same solutions. For both of these methods the force term remains as initially set, such that there is no back reaction on the forcing by the velocity changes associated with the flow. Furthermore, in a few cases with a weak instability and moderate viscosity the unstable flow decays after saturation. The third method uses horizontally averaged velocity profiles of the actual flow to formulate a force that drives the shear flow back to the target profile after a relaxation time that can be chosen. Such a method provides a self regulating force and more control on the strength of the forcing due to different relaxation times. Comparing the exponential growth phase with solutions from a linear stability analysis shows that the growth rates achieved by all methods used are close to the predicted value, if for the relaxation method a sufficiently small relaxation time τ 0 is chosen. Cases with larger τ 0 lead to a slightly different shear instability, since the viscous dissipation is initially not balanced and the initial state evolves before an instability can occur. However, focusing on the non-linear evolution a significant difference in the system dynamics is revealed when using the relaxation method with different τ 0 . When choosing τ 0 greater than the viscous time-scale a non-linear evolution like for the viscous forcing method is achieved, where the long time evolution converges towards a quasi-static state. Energy induced into the system by the force balances the loss by viscous dissipation, but little additional kinetic and potential energy is obtained. In contrary a relaxation time τ 0 less than t µ leads to a system which is constantly forced and develops a turbulent region which spreads across a larger region in the vertical direction. Such cases do not tend to evolve towards a quasi-static state, which becomes evident due to their energy evolution. Moreover, the energy induced into the system is significant greater than the loss by dissipation such that the energy overrun is converted into kinetic energy of the disturbances and available potential energy of the system. Analysing the turbulent Reynolds number for late times shows that when decreasing τ 0 the horizontal layer of turbulent flow reaches a larger vertical extend and very high Re numbers. However, greater τ 0 and the viscous method develop a small region confined around the middle of the domain with moderate Re numbers. Thus the strength of the forcing has a strong impact on the spread of the resulting turbulent region. Interestingly, the mean flow resulting from viscous forcing develops a peculiar form around the middle plane, where a steep slope is present, while the relaxation method leads to a horizontally averaged velocity profile that generally preserves a hyperbolic tangent profile. Since the physical mechanism driving shear flows in different objects are not known in detail, the relaxation method provides a tool to adjust the force such that a more suitable flow can be achieved. Therefore, we conclude that the relaxation method provides a more suitable method to sustain a velocity profile when modelling stellar interior as for example the tachocline in our Sun or shear regions in more massive mainsequence stars. However, in order to study the non-linear evolution of a shear driven turbulent flow the viscous method or the equivalent perturbation method suit better, as no artificial dynamics due to the forcing affects the modification of the background profiles. The physical mechanism for the generation and maintenance of the differential rotation in the solar interior and especially the tachocline is not well understood (Gough & McIntyre 1998;). It is widely believed that external processes such as Reynold stresses, which originated in the convection zone, drives the shear flow in the tachocline (Miesch et al. 2008). However, we do not know what form the resulting force has that drives the shear flow in the tachocline. The viscous method correspond to a forcing confined within the shear region whereas the relaxation method corresponds to a bulk forcing. Prospective global-scale in-at City University, London on August 16, 2016 http://mnras.oxfordjournals.org/ Downloaded from vestigations might reveal which of these two forcings is more relevant when modelling the tachocline. Having established a detailed analysis of possible numerical methods to sustain a localised shear flow with minimised effect on the boundaries, possible applications of shear driven turbulence in stellar interiors have to be considered and are currently underway.