{'fraction_non_alphanumeric': 0.04015875151352079, 'fraction_numerical': 0.015897573882236872, 'mean_word_length': 4.182314664187776, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Response maxima in modulated turbulence Part II: Numerical simulations

Numerical simulations of fully developed turbulence driven by a modulated energy input rate or driving force are performed within two dynamical cascade models, the GOY shell model and a reduced wave vector set approximation of the Navier-Stokes equation (REWA). The frequency behavior of the system response is studied and compared with predictions from a variable range mean-field theory, which excludes turbulent fluctuations. In agreement with the mean-field approach we find a constant response amplitude for low driving frequencies and a 1/w-decay of the amplitude for high frequencies. In the mean-field theory, the finite cascade time scale had lead to an oscillating behavior of the response amplitude as a function of the driving frequency. In the simulations of both models we observe the main maximum. The higher maxima and minima are completely washed out by fluctuations.

Numerical simulations of fully developed turbulence driven by a modulated energy input rate or driving force are performed within two dynamical cascade models, the GOY shell model and a reduced wave vector set approximation of the Navier-Stokes equation (REWA). The frequency behavior of the system response is studied and compared with predictions from a variable range mean-field theory, which excludes turbulent fluctuations. In agreement with the mean-field approach we find a constant response amplitude for low driving frequencies and a 1/ω-decay of the amplitude for high frequencies. In the mean-field theory, the finite cascade time scale had lead to an oscillating behavior of the response amplitude as a function of the driving frequency. In the simulations of both models we observe the main maximum. The higher maxima and minima are completely washed out by fluctuations.

I. MOTIVATION
Many realistic turbulent flows are subject to modulated driving forces, as e.g. the atmosphere of the earth driven by the periodic heating of the sun or the pulsed flow through a pipeline. Three dimensional turbulence is characterized by an energy cascade from the outer length scale, where the forcing acts, to the dissipative scale, where most of energy is dissipated, see e.g. [1,2]. The down-cascading of energy from large to small scales takes a characteristic time τ . In a statistically stationary flow the energy dissipation rate equals the energy input rate. In a situation with time dependent energy input, on the other hand, this statement will only hold on average, whereas the energy dissipation at a certain time t is expected to depend on the energy input at an earlier time due to the finite time delay of the energy transfer.
In a previous work [3] the effect of an energy input rate modulated in time, e in (t) = ǫ 0 (1 + e sin ωt), with a modulation amplitude e ≪ 1 and a modulation frequency ω, has been studied within a variable range mean-field theory [4]. The response of the system can be observed in the second order velocity structure function of the flow field at the outer length scale L, D L (t) = (u(x + L, t) − u(x, t)) 2 = 6u 1,rms , which is equivalent to the Reynolds number Re(t) = u 1,rms (t)L/ν of the flow and the total energy E(t) = u 2 /2 of the system. Here, u 1,rms is the rms of one velocity component and ν is the viscosity. The response follows the oscillation of the energy input rate with almost constant modulation amplitude at low frequencies ω of the energy input rate, whereas the response amplitude strongly decreases (∝ 1/ω) at higher frequencies. The finite energy transfer time τ plays a crucial role in this theory. This time τ is the average time the energy stays within the system while it is transported by the interaction cascade from the large eddies towards the small eddies, where it is finally dissipated. This intrinsic time scale of the system is a multiple a of order 1 of the large eddy turnover time τ L , corresponding to the sum over the eddy turnover times on all scales. τ −1 determines the frequency at which the crossover takes place between the regime of constant response amplitude and decreasing amplitude. In addition, it leads to an oscillating behavior of the system response with driving frequency ω, where the maxima and minima are at frequencies connected to the inverse of the energy transfer time τ . In the limit of large frequencies ω, the extrema of the response can be estimated to be at frequencies ω r ≃ n π τ , n = 1, 2, 3, .... Recent experiments on modulated turbulence in a cylinder between two counter rotating disks [5] revealed evidence for the proposed response maxima. In accordance with the predictions from the mean-field theory [3], for small frequencies a constant response amplitude was measured. For large driving frequencies a 1/ω-decay of the velocity response amplitude was observed, again in agreement with the prediction from our mean-field approach. Note here, that both the velocity response as well as the energy response are the same up to a factor of 2, in linear order, cf. section VI of reference [3].
In the experiments the amplitude of the driving force rather than that of the energy input rate is modulated. Since the energy input rate is not a controlled quantity any more it can serve to measure the response of the system. Of course, also within the mean-field theory we can apply a modulated driving force, see [3]. The main features, the 1/ω-decay of the energy response amplitude for high frequencies and the constant response amplitude for low frequencies, pertain. The response maxima are only slightly shifted in comparison to the case of a controlled and modulated energy input rate. In the case of a modulated driving force, as in the experiments, the energy input rate as a response of the system, also shows maxima in addition to the mentioned mean features. These are at the same frequencies as the maxima of the total energy response amplitude.
In the mean-field approach, the (intermittent) fluctuations of the energy and, in particular, of the cascade time τ are not present. In experiments and numerical simulations these fluctuations are of course present, and  they may lead to broader and less pronounced response maxima and minima. Therefore, in this paper we shall study the frequency dependence of the response to a modulated energy input rate into a system where turbulent fluctuations are included. In particular, we shall address the question whether the response maxima and minima in REWA. The result is compared with the response amplitude as calculated within the mean-field model with the same e, τ , and b (dashed lines). (a) Log-log plot of the amplitude A versus frequency. The dotted line is ∝ 1/ω. The arrow indicates the mean-field crossover frequency ω M F cross τL = (6/b) 3/2 = 0.019. Inset: linear scale plot of the response amplitude. (b) Log-log plot of the compensated amplitude, i.e., A/(ωτL) −1 versus frequency. A clear maximum is observed in REWA at a frequency near to the first maximum of the mean-field amplitude. The arrow indicates the height of the maximum, i.e., a deviation from the 1/ω-decay by a factor of 1.4 in REWA.
can still be well identified in the presence of fluctuations. Furthermore, we not only consider a modulated energy input rate, but also discuss the slightly different case of a modulated driving force in order to compare with the above mentioned experiments.
An appropriate way to numerically study the prob- lem of modulated turbulence would be a direct numerical simulation of the Navier-Stokes equation for this specific time-dependent energy input rate. However, as we need high Reynolds numbers to achieve fully developed, isotropic, and homogeneous turbulence and, in addition, need the response of the system as a function of time for a wide range of driving frequencies, the computational demands would be too high. Therefore, we first study the problem within a dynamical cascade model of turbulence, the Gledzer-Ohkitani-Yamada (GOY) shell model [6,7,8,9,10,11,12,13]. With this model large Reynolds numbers and enough statistics within a reasonable computing time for each driving frequency can be achieved. The GOY model has been successfully used in a study about decaying and kicked turbulence [14]. In addition, to be even closer to a numerical Navier-Stokes simulation and to distinguish between real effects and artifacts of the turbulence model, we follow another approach. We calculate the response of the system to a modulated energy input rate within a reduced wave vector set approximation (REWA), [15,16,17], where the Navier-Stokes equation is solved on a reduced, geometrically scaling subset of wave vectors. This method is much closer to the Navier-Stokes dynamics than the GOY-model, as it contains (i) much more modes than GOY, (ii) it solves the Navier-Stokes equation for those modes and not only a model equation, and (iii) it is three dimensional.
Our main results are summarized in Figs. 1, 2, and 3. In Figs.1 and 2, the amplitude A of the energy response is shown as a function of the driving frequency for both the GOY model ( Fig.1) and the REWA simulation ( Fig.2) with a modulated energy input rate. This is compared with the results of the mean-field model with the corresponding parameters, i.e, the same modulation amplitude e and time scale τ . In Fig.3 the results from the GOY model solutions are shown for a modulated driving force and compared with the mean-field model. In all cases we observe a constant amplitude for low driving frequencies and a 1/ω-decay for high frequencies. This can in particular be observed in the compensated plots (parts b of all three figures), where A, compensated by its asymptotic amplitude, i.e., A/(ωτ L ) −1 , is plotted versus frequency. The 1/ω-decay of the response means, that for fast modulation no response is detectable any more. The remaining dissipation rate is that of the stationary system itself.
In the mean-field approach, a sequence of response maxima is present for both types of forcing, starting at a frequency ω ∝ 1/τ . In the simulations, this main maximum can also be observed, although it is weaker and broader, i.e., it is "washed out" by fluctuations. The higher order maxima and minima are not visible in the simulations, but are completely washed out by fluctuations. On the other hand, we emphasize that the turbulent fluctuations in the GOY model are strongly overestimated due to the extreme mode reduction in this model. In the REWA simulation, an artificially large Kolmogorov constant b is found. Using such large b in the mean-field approach also leads to a considerable weakening of the first maximum and a shrinking of the higher order maxima and minima towards very small amplitudes.
These results will be explained and discussed in detail in this paper, which is organized as follows. In the next section we study the modulated turbulence within the GOY shell model. Before calculating the response of the system to a modulated energy input rate as well as a modulated forcing in section II B and II C we briefly introduce the model and study its stationary properties in section II A. In section III we present our findings on modulated turbulence within the reduced wave vector set approximation. We summarize our results in section IV.

A. Stationary properties
The GOY shell model consists of a set of coupled ODEs for one-dimensional complex velocity modes u n [6,7,8,9,10,11,12,13]. These modes u n correspond to velocity differences |u(x + r n ) − u(x)| on scale r n . N modes are taken into account, n = 1, 2, ..., N , one complex velocity mode per cascade level n, defined by the wave numbers k n = λ n k 0 which are equally spaced on a logarithmic scale, here, λ = 2. The model equations read: where n = 1, ..., N , a = 1, b = −1/4, and c = −1/2. These are the traditional parameters. We impose boundary conditions on the u n , i.e., u n = 0 for n < 1 or n > N . We use N = 14 shells, a viscosity of ν = 1.01875 · 10 −4 , and k 0 = 2 −4 . The forcing acts on the largest scale, i.e., the first shell, n = 1. F is constant, Together with ν and k 0 this sets the time and length units as well as the Reynolds number. Equations (2) are integrated using a fourth order Runge-Kutta scheme with adaptive step size [18]. With the above chosen parameters the GOY dynamics is chaotic [13]. The system is forced on large scales while most of the energy is dissipated on small scales. It reaches a steady state, in which the velocities are stochastically fluctuating. In this sense the system has similar properties as three dimensional Navier-Stokes turbulence. The scaling behavior of structure functions and dissipation has been extensively studied in [10,11,12,13,19,20,21]. The deviations from K41 scaling due to intermittency observed in the GOY model are very similar to experimental values. Most of the previous studies have been done with 22 or more shells. As we use here only 14 shells in order to reduce the computational effort we explicitly check some scaling properties of, e.g. the structure functions and the energy spectrum in a simulation with constant, stationary forcing F = F 0 in Equations (2). N = 14 then turns out to be sufficient.
The Reynolds number of the system can be defined as follows. An outer length scale L is given by the smallest wave number k 1 , L = 1/k 1 . A typical velocity U is the velocity on that scale, |u 1 | 2 1/2 t . The average ... t is taken over time. With these length and velocity scales the Reynolds number of the present simulation is Re 0 = UL ν = 8.6 · 10 4 . The simulated time interval is several hundreds of large eddy turnover times |k 1 u 1 | −1 . For the second order structure function we use the following method. In [10] it has been suggested to study the scaling of Σ n,q = |ℑ(u n u n+1 u n+2 + 1−ǫ λ u n−1 u n u n+1 )| q/3 instead of the pure moments of the velocity s n,q = |u n | q in order to eliminate the period 2 and period 3 oscillations which are an artifact of the GOY model. The second order quantity Σ n,2 , corresponding to the second order structure function, is shown in Fig.4 as a function of wave number index n = log 2 (k n /k 0 ). One clearly observes an inertial subrange (ISR) between shell 2 and 9, where the second order structure function shows scaling with an exponent near to the K41 value 2/3, i.e., Σ n,2 ∝ k −0.71 n corresponding to ∝ r 0.71 for the structure function. The scaling exponent is not equal to the K41 value, because the model shows intermittency corrections. The higher wave numbers, i.e., smaller scales, n = 12 − 14 belong to the viscous subrange (VSR), where dissipation takes place. In this range the viscosity term is dominant and the velocity decays rapidly with k n , in fact The k-space the scaling regime of the spectrum is more extended, i.e., about three decades, as has been already shown in Fig.4.
more than exponentially [20,21]. The external forcing of the flow acts on shell n = 1, therefore the stirring subrange (SSR) contains the first shell only. The spectrum, which is obtained by fast Fourier transforming (FFT) a time series of the velocity u(t) = n Re(u n (t)) and raising it to the power two, can be compared to the energy spectrum under the assumption of the Taylor hypothesis. This spectrum is shown in Fig.5 for N = 14. We observe about one frequency decade of (nearly) Kolmogorov scaling, where the energy decays as ∝ f −5/3 with frequency f . For N = 14 the Kolmogorov scaling range is not yet well developed. For less shells this region becomes even narrower and the spectrum shows strong peaks at some intrinsic frequencies. When reducing the number of shells even more, the velocities relax to a stationary value without any fluctuations, i.e., the chaotic behavior of the system is lost.
In conclusion, the GOY model with N = 14 shells exhibits an inertial subrange scaling, although the ISR for the frequency spectrum is only narrow, but the spectrum in k-space has a scaling range of about three orders of magnitude. This seems acceptable for our goal to study the response of the system to a modulated driving, because we are only interested in global quantities like the total energy E(t) but do not need information on scale resolved quantities.
The time scales in the model have been determined as follows. For each shell n an eddy turnover time τ n is defined by [10]: This is also considered as the time scale for the turbulent energy transfer through the nth level. The time scale relevant for the energy loss on level n due to viscosity is defined as Both time scales are shown in Fig.6.
In the ISR between shell 2 and 9, where the energy transfer times are the relevant time scales for the dynamics, the decrease of the τ n with n is near to τ n ∝ 2 −2n/3 (dashed line in Fig.6) as expected for the turnover times of eddies of sizes r n /L ∝ (1/2) n . In this range the dissipation time scales τ d n are much larger than the τ n , meaning that the turbulent energy transfer is much faster than the viscous dissipation, and therefore the dominant process. In the VSR instead τ d n < τ n , i.e., on average the energy is dissipated by viscosity before it can be transferred to smaller scales.
The largest eddy's turnover time is in general defined by the velocity on the outer length scale L, i.e., on the length scale of the forcing, which in this case is 1/k 1 . However, in this model, the time τ 1 is disturbed due to finite size effects. Therefore, we extrapolate from the turnover times of the other shells. A linear fit in Fig.6 for log 2 τ n with n = 2, ..., 11 leads to τ f it n = 22.46 · 2 −0.63n . Without intermittency one would have τ n ∝ 2 −2n/3 ; the small deviation corresponds to the intermittent scaling of u n ∝ k −0.37 n or an intermittency correction of δξ 1 = ξ 1 − 1/3 ≃ 0.04. The extrapolation for n = 1 yields for the large eddy turnover time τ L = τ f it 1 = 14.5. The time scale corresponding to the energy transfer time τ used in the mean-field model [3] is the sum over the eddy turnover times of all energy-input and inertial-range shells, i.e., here τ ≃ 11 n=1 τ n ≃ 35.9. The factor τ /τ L = a between the transfer time and the large eddy turnover time is then a = 2.48.

Constant energy input rate
Until now, we have considered a constant forcing F = F 0 . The resulting energy input rate e in (t) = u * 1 (t)F 0 then fluctuates around its mean value because of the u * 1 (t)-fluctuations. In the mean-field theory [4] the energy input rate e in is constant instead. For closer comparison we also consider another type of forcing in the GOY-model: F (t) = ǫ 0 u 1 (t)/|u 1 (t)| 2 . This forcing F (t) fluctuates as u 1 (t). Then the energy input rate is e in = u * 1 F = ǫ 0 = const by definition. The ISR-scaling behavior as well as the energy spectrum then turn out to be similar to the previously discussed ones with the constant forcing F = F 0 . The energy transfer time τ is slightly larger in this case, namely τ = 39.5, and the large eddy turnover time is τ L = τ f it 1 = 15.57. Again, the large eddy turnover time is extrapolated from τ 2 , ..., τ 11 . This leads to the factor τ /τ L = a = 2.54 between the total time delay of the energy cascade and the large eddy turnover time. In the following sections we will study the timedependent cases where either the energy input rate e in , i.e., F = ǫ 0 u 1 /|u 1 | 2 (section II B) or the forcing F = F 0 is modulated (section II C).

B. Modulated energy input rate
In this section we apply a modulated energy input rate to the GOY model, i.e., we set the forcing F = F (t) in equations (2) as with a modulation amplitude e = 0.2. Then, the resulting energy input rate e in is = ǫ 0 (1 + e sin ωt), and has a prescribed modulation amplitude eǫ 0 by definition. The total energy of the system is calculated for a wide range of driving frequencies ω in order to study the frequency behavior of the response. The brackets ... denote the ensemble average. This ensemble average is performed as follows. From a long stationary simulation we collect an ensemble of 1500 starting configurations which we then let evolve according to Eqs.(2) but now including the modulation of the forcing F (t), Eq.(5), and average over these 1500 time series. To ensure that the different realizations can be considered as statistically independent, the time delay between the successive starting configurations for the different realizations is chosen to be about 100 large eddy turnover times. The adaptive step size routine controlling the Runge-Kutta integration algorithm for ODEs does not produce the same time steps for all time series. To overcome this, we have calculated also the same number of equidistant time steps by spline interpolation for all time series.
The oscillating response of the system ∆(t) is then studied in terms of the ratio between the energy E(t) with modulated energy input and the energy E 0 (t) without modulation, namely E and E 0 both are averaged over 1500 realizations. In spite of the averaging not only E but also E 0 still contains (weak) fluctuations. Therefore, we write E 0 (t), as E 0 is still slightly fluctuating around its mean value. E and E 0 contain about the same size of fluctuations. Accordingly, the energy input rate e in is normalized by its stationary value, e norm in (t) = e in (t)/ǫ 0 = 1 + e sin ωt. In Fig.7 the input rate e norm in (t) and the energy E norm (t) are plotted for four different driving frequencies. For the two low frequencies where ωτ L ≪ τ L /τ = 1/a ≃ 0.39, the energy follows the oscillation of the energy input rate with almost constant, but smaller amplitude. For higher frequencies the amplitude of the deviations of the normalized energy from its stationary value 1 strongly decrease, and a phase shift with respect to the energy input becomes visible. The same behavior of the energy has been observed in the mean-field theory [3].
To quantitatively access the frequency behavior of the response amplitude, we calculated time series of the total energy E(t) for 85 different driving frequencies varying over almost 3 decades between 0.012 ≤ ωτ L ≤ 9.3. The chosen frequencies are approximately equally spaced on a logarithmic scale. The normalized energy E norm (t) is fitted by a function of the form with three free parameters: E const , the amplitude A, and the phase shift Φ. E const is near to 1 for all frequencies, i.e., E const = 1.0022±0.0032. The fits (9) are included in Fig.7 as dashed lines but they are mostly indistinguishable from the solid lines for the energy itself. Fig.1 of section I shows the amplitude A, resulting from the fitting procedure, as a function of the dimensionless frequency ωτ L . A is almost constant for low frequencies and has a value of about 2/3. For higher frequencies the amplitude decreases as ∝ 1/ω. The same features have been observed in the mean-field calculations, see dashed lines. The long-dashed line in Fig.1 represents the low frequency limit of the mean-field theory, A ≃ 2/3, and the dotted line the high frequency limit, A ∝ 2/(3ωτ L ).
The crossover between the regime of constant amplitude and the one of 1/ω-decay of the energy response takes place at ω cross τ L ≃ 1/a ≃ 0.39, which is indicated by the arrow in Fig.1. In the mean-field approach this crossover is always at ωτ L = 1, independent of the factor a between the large eddy turnover time and the total time scale of the energy transfer. In experiments [5] the crossover frequency has been used to measure the energy cascade time scale. The present simulations confirm that this frequency gives the correct order of magnitude for the energy transfer time.
Response maxima, as observed in the mean-field model at frequencies connected with the inverse energy transfer time, are difficult to be identified in Fig.1a. There is some structure visible at ωτ L ≃ 0.31 and ωτ L ≃ 1.57. In Fig.1b, where the amplitude A, compensated by the asymptotic amplitude (ωτ L ) −1 , i.e., A/(ωτ L ) −1 is plotted versus frequency, this structure becomes more evident, and we see a clear maximum at a frequency of about ωτ L ≃ 1.57. This maximum probably corresponds to the mean-field maximum. Of course, the maximum in GOY is broadened and weakened due to the large fluc- tuations, and the higher order maxima and minima are apparently washed out completely. As in the mean-field theory no fluctuations are included, the energy cascade time τ is considered to be constant. However, in the GOY model this assumption is not true, as can clearly be seen in Fig.8. Here, a time series of the cascade time τ (t) = 11 n=1 τ n (t) is plotted, computed within the GOYshell model with non-modulated forcing F = ǫ 0 u 1 /|u 1 | 2 . The inset shows the probability distribution of τ /τ L . This distribution has its maximum at τ /τ L = 2.39, almost at the mean cascade time a = 2.54, and a width (FWHM/2) of about 0.46a. The width is almost half the size of the mean which indicates that the transfer time fluctuates strongly and therefore we have to expect that the response maxima are more or less washed out. However, these strong fluctuations are considered as an artifact of the GOY-model and not as a feature of real turbulence. The GOY-model contains only one velocity mode per cascade level instead of infinitely many modes in real turbulence. This one-mode approximation leads to an overestimation of the fluctuation strength. In order to confirm this, we performed another simulation with more modes per level within the reduced wave vector set approximation (REWA) of the Navier-Stokes equation. This will be presented in section III.

C. Modulated driving force
In this section we present further results within the GOY model based on a non-fluctuating driving force F 0 which is regularly modulated as was ǫ 0 in the previous section. This case may be more comparable to the exper-imental method in reference [5], because there the driving force is modulated. What cannot be modeled with GOY is the spatial inhomogeneity in the experiments [5].
In equations (2) we now apply a forcing with a modulation amplitude of e f = 0.2. As in section II B we calculate the ensemble averaged time series of the energy input rate e in (t), see first line of equation (6), and the total energy of the system E(t) cf. equation (7), for 89 different frequencies between 0.0144 ≤ ωτ L ≤ 3.04, again logarithmically equally distributed. The normalized energy E norm (t) and energy input rate e norm in (t) are then fitted by a function according to equation (9), with the parameters E const , A, Φ, and e in,const , A ein , Φ ein , respectively.
The amplitudes A and A ein are plotted in Fig.3 of section I as a function of the dimensionless frequency ωτ L . Also in this case of a modulated force, the response amplitude is almost constant for small frequencies, namely A ≃ 1, and decreases as 1/ω for high frequencies, see Fig.3a (full dots). Again, the long-dashed line represents the low frequency limit of the mean-field theory for A (which is A ≃ 1 in this case) as well as the dotted line the high frequency limit. As in the previous section, the crossover frequency between the two regimes is determined by the energy transfer time, i.e., ω cross τ L ≃ 1/a ≃ 0.40 (as a = 2.48 in this case), which is marked by the small arrow in Fig.3a. The amplitude of the energy input rate A ein starts with a value of about 1.5 for low frequencies and merges towards 1 for high frequencies, see Fig.3c. This indicates that at very large frequencies the velocity is not oscillating any more as it only feels a mean constant force. The oscillations of the energy input rate are then only a consequence of the oscillation of the driving force F . In the mean-field theory we have observed the same trend for both amplitudes. The corresponding mean-field results are included as dashed lines.
Both amplitudes A and A ein show a maximum at a frequency near to the crossover frequency ωτ L ≃ τ L /τ = a −1 = 0.40. In the compensated plot Fig.3b, where A/(ωτ L ) −1 is plotted as a function of frequency, a clear deviation from the dotted line representing A ∝ 1/ω can be observed. At this frequency, the mean-field theory predicts a first maximum for the energy response amplitude (Fig.3a) and a maximum directly followed by a minimum for the energy input rate, see Fig.3c. This frequency is connected to the energy transfer time τ . In the mean-field model further maxima of the energy response and wiggles of the energy input rate are observed at multiples of this frequency. However, in the GOY model all further maxima and minima are washed out presumably because of the strong fluctuations. The fluctuations of the energy cascade time scale are found to be similar as in the case of a fluctuating force shown in Fig.8 in the previous section.

III. MODULATED TURBULENCE IN A REDUCED WAVE VECTOR SET APPROXIMATION
As was pointed out in the introduction of this chapter a full numerical simulation of the Navier-Stokes equation for modulated turbulence at high Reynolds numbers is still not possible or requires low Reynolds numbers. Therefore, we first have considered the GOY shell model. This model correctly describes many features of turbulence, however, due to the one-mode approximation in each cascade level, it contains various artifacts. Namely, it strongly overestimates the strength of the fluctuations. The aim of this section is to study the problem of modulated turbulence within another model, the reduced wave vector set approximation (REWA) [15,16,17], which is much closer to the Navier-Stokes equation than the GOY model and contains much more modes per cascade level. Of course, as compared to full numerical simulations of the Navier-Stokes equation, it still contains a mode reduction in order to make the computational effort reasonable for the desired high Reynolds numbers. The present approximation has been introduced and extensively studied in [15,16,17] . Here, we use it together with a time-dependent driving. For completeness we briefly explain the approximation before we present the results with modulated driving.
To maintain the turbulent flow we apply a forcing as in [15,16,17]: The subset K in of K 0 by choice contains the wave vectors with the three smallest lengths. K in contains 14 vectors, namely K in = {±(2, 2, 2), ±(−1, 2, 2), ±(−2, 1, 1), +permutations}. In reference [16] it has been shown that the statistics of the solutions of the equations of motion do not depend on the particular choice of K 0 . This forcing corresponds to the same type of forcing, which has already been applied to the GOY model in section II B. It enforces the energy input rate to be modulated: n ) · f (p (l) n ) = ǫ 0 (1 + e sin ωt). (15) Eq.(11) is a set of 3N (l max + 1) coupled ODEs for the complex mode amplitudes u i (p (l) n ) which is numerically solved within the Burlisch-Stoer integration scheme with adaptive step size [18]. Length scales are measured in units of L 0 and time scales in units of L The main features of fully developed turbulence as irregular velocity signals, characteristic scaling of structure functions, etc. are well described within this approximation, as has been shown in [15,16,22]. The REWA solutions show small scale intermittency, which is produced by the competition between down-scale energy transport and viscous dissipation on the small scales [16,22]. Other mechanisms leading to intermittency in turbulence as e.g. nonlocal interactions between wave vectors are underestimated in this approximation [23]. The down-scale energy transport in the REWA fluid is less effective than in real turbulence, because in this approximation the larger wave vectors are more and more thinned out [24]. This is in contrast to the case of the complete set of wave vectors (e.g. in full grid simulations) where the density of states increases ∝ p 2 , whereas in the reduced wave vector set K the number of admitted wave vectors decreases as 1/p [17]. In reference [24] it has been shown that this reduced energy transport leads to an overestimation of the Taylor Reynolds number of the system as well as the Kolmogorov constant b, defined by D(r) = b(ǫr) 2/3 , by roughly one order of magnitude for our choice of N . In the present simulation we obtain b = 83.5 instead of b = 6 − 9 as in experiments [25,26,27]. Since D(L) is the energy density ∝ u 2 of the fluctuations in the fluid system, the large b value indicates that even in the REWA approximation the strength of the fluctuations is highly overestimated. The large Kolmogorov constant will change the relevant time scales in the system, as will be shown in section III B.
The characteristic time scale for the turbulent energy transfer on scale l can be estimated as: where p (l) denotes the mean wave number on scale l, i.e., it is the mean inverse eddy size in K l . As in the GOY model, Eq.(4), the time scale of viscous dissipation is τ d = [ν(p (l) ) 2 ] −1 . Again, in the ISR τ (l) > τ d , whereas in the VSR τ d > τ (l). From a simulation with stationary forcing, i.e., e = 0 in Eq. (14), the time delay of the energy down-transport τ sum is then estimated by the sum of all τ (l) in the ISR, τ sum = l∈ISR τ (l) ≃ 0.186. The largest of these τ (l), on the largest scale, τ (0) = 0.0632 can be regarded as a large eddy turnover time τ L . Thus, τ sum = 2.94τ L , and the factor between the cascade time scale and the large eddy turnover time is τ sum /τ L = a = 2.94. As we have seen in the GOY model, Fig.8, the energy transfer time is strongly fluctuating. We attributed these strong fluctuations to the one-mode per level approximation of the GOY model. Fig.9 shows a time series of the energy transfer time τ sum (t) (use Eq.(16) with u (l) (t)) in the present reduced wave vector set approximation (REWA), and, in the inset, the distribution of this time scale. Clearly, the fluctuations are much weaker than in the GOY model; they are zero in the mean-field approximation. The distribution is centered around τ sum /τ L = a with a width (FWHM/2) of about 0.02a.

B. Modulated energy input rate
The response of the system to a modulated driving force, cf. Eq. (14), is calculated now in terms of the total energy of the system The modulation amplitude of the energy input rate e in (Eq. (15)) is chosen as e = 0.3. The average ... is performed as follows. We average over 25 to 30 realizations, which are obtained from Eqs. (11) with different starting values. The delay between the different starting values is one driving period. For the higher frequencies the period 2π/ω becomes too small to ensure statistical independence of the different realizations. Then, we chose the delay between the successive starting values of the different realizations to be several driving periods such that it is at least 60τ sum . The response is calculated for 150 (approximately equally spaced on a logarithmic scale) frequencies between 0.00016 ≤ ωτ L ≤ 3.0. The energy is normalized by E 0 , calculated from a stationarily forced solution with e = 0 and averaged as E. The oscillating response ∆(t) of the system is defined in the same way as for the GOY model, see Eq. (8). Then, as for the GOY calculations, the averaged and normalized signals E norm (t) are fitted with Eq. (9). The fit parameter E const is again near to 1 for all frequencies, E const = 1.0064 ± 0.0065. In Fig.10 the time averaged responses and the normalized energy input rates are plotted for four different driving frequencies. Also the fits according to Eq.(9) have been included = 1 + e sin ωt (dotted lines) and energy E norm (solid lines) for four different modulation frequencies ω as calculated in the REWA simulation. The energy input rate is modulated with a modulation amplitude of 30% of the constant energy input rate, i.e., e = 0.3 in Eq. (14). Also included is the fit to the energy data cf. Eq.(9) as dashed lines but these are in all cases indistinguishable from the solid lines. The averaged time series of E norm are repeated once for better visibility. (a) ωτL = 3.16 · 10 −4 , (b) ωτL = 3.16 · 10 −3 , (c) ωτL = 3.16 · 10 −2 , (d) ωτL = 0.316. For larger ωτL the energy is indistinguishable from 1 on this scale. The crossover to the 1/ω-decay regime is in this simulation at ωτL ≃ 1.1 · 10 −2 , between the frequencies of b) and c).
as dashed lines but are indistinguishable from the solid lines for the energy signal itself.
We observe in Fig.10 for REWA the same features as in Fig.7 for the GOY model and in Fig.1 of reference [3] for the mean-field model. For the two lower frequencies the response amplitude remains almost constant and is about 2/3 of the amplitude of the energy input rate, whereas for the two higher frequencies the response amplitude strongly decreases. This trend becomes more clear in Fig.2a in section I, where the amplitude of the response -determined from the fit (9) -is shown as a function of the driving frequency (full dots). For low driving frequencies A ≃ 2/3, whereas for high frequencies the amplitude decreases as 1/ω. The crossover between the regime of constant amplitude and the one of 1/ω-decay takes place at ω cross τ L ≃ 0.011, i.e., at a much smaller frequency than expected from the original case of the mean-field theory in which the crossover was at ω MF cross τ L = 1 with b = 6. We understand this as follows. In section III A it was mentioned that the Kolmogorov constant in the REWA simulation is b = 83.5 instead of b = 6 − 9 as in experiments. In the figures of the mean-field approach [3] we have set b = 6. The mean-field solution for a general b revealed that the crossover frequency decreases with increasing b while the positions of the response maxima are left unchanged. For b = 83.5 the mean-field crossover fre-quency is at ω MF cross τ L = (6/b) 3/2 ≃ 0.019 in close agreement to what we observe in the REWA simulations. The response amplitude calculated from the mean field model with b = 83.5 is included in Fig.2a as dashed line. Apart from the changed crossover frequency we observe that, in the mean-field calculations, the first response maximum at ωτ L ≃ 0.1 is considerably smaller and broadened as compared to the case with b = 6, given as the dashed line in Fig.1a. In agreement with this, our REWA simulations (with a value b = 83.5) show a broad maximum in the response amplitude at ωτ L ≃ 0.028. This means it occurs at a similar frequency as the mean field model. The maximum becomes more clear in the compensated plot, Fig.2b, where A/(ωτ L ) −1 is shown as a function of frequency. There, we observe a deviation from the 1/ωdecay of the amplitude by a factor 1.4 at the maximum in the REWA simulations. The mean field maximum has a height of 2.8. The subsequent maxima and minima in the mean field model occur at frequencies where the amplitude is already very small (A ≤ 10 −2 ) because the crossover to the 1/ω-regime takes place at a much smaller frequency whereas the response maxima stay at the same frequencies as for a smaller b. Therefore, the higher order maxima are not visible in the REWA simulations. The cascade time shows about 2% fluctuations in REWA as shown in Fig. 9. However, at small response amplitudes these fluctuations are already large enough to wash out the higher order response maxima.
In conclusion, the REWA system reproduces qualitatively the features of modulated turbulence as predicted by the mean-field model including the first response maximum. The latter is considerably weakened due to the large Kolmogorov constant in REWA. Another consequence of the large b is, that the crossover between constant amplitude and 1/ω-decay is shifted towards smaller frequencies, and therefore the higher order maxima and minima are already at very small amplitudes where the fluctuations in the cascade time scale are finally large enough to wash them out. We cannot clarify, at present, how close the response in direct numerical simulations -which lead to an order of magnitude smaller b and thus have much smaller fluctuations -will come to the mean-field features, but we expect a clearly visible first maximum at least.

IV. CONCLUSIONS
We have simulated the response of modulated turbulence within two numerical models. Namely, we have used the GOY shell model and the reduced wave vector set approximation of the Navier-Stokes equation (REWA). The results are compared with predictions from a mean-field theory. For a modulated energy input rate this mean-field theory had predicted a constant response amplitude for low frequencies and a 1/ω-decay for high driving frequencies. In addition, at certain frequencies connected with the energy cascade time scale, a sequence of maxima and minima of the response amplitude is observed.
Both numerical models well reproduce the basic trend, i.e., the constant amplitude for small ω and the 1/ωdecrease for large ω. The main response maximum can be observed in both numerical models, although it is weakened due to fluctuations. The higher order maxima and minima as predicted by the mean-field theory cannot be identified in the simulations. Obviously, they are strongly washed out by fluctuations. In the GOY model the large fluctuations are explicitly visible in a broad cascade time distribution. We attribute these fluctuations to the strong mode reduction in the model, i.e., they are an artifact of the model properties and not a feature of real turbulence. In the reduced wave vector set approximation of the Navier-Stokes equation these fluctuations are much weaker, and, we believe, more realistic for real turbulence, as more modes are taken into account. However, due to an overestimated Kolmogorov constant b in REWA the higher order maxima are considerably reduced, and therefore washed out by the fluctuations although being smaller than the GOY-fluctuations. Therefore, we believe, that in real turbulence with a realistic Kolmogorov constant the first maximum should be clearly observable and possibly also the higher order extrema in the response. Thus, the predictions of the meanfield model, which excludes all fluctuations, might be quite reasonable for real turbulence. To further study the response maxima numerically, it is necessary to perform full numerical simulations of the Navier-Stokes equation, as then all relevant time scales including their fluctuations are reproduced realistically, which turned out to be essential for the observation of the response maxima.
Recent experiments on modulated turbulence [5] revealed evidence for the response maxima. These experiments may be more comparable to the case of a modulated force instead of a modulated energy input rate. We have studied this case also within the mean-field model and have found basically the same behavior of the energy response as for a modulated energy input rate. In addition, the amplitude of the energy input rate showed "wiggles" at the same frequencies where the energy response had maxima. In the experiments the response maxima were measured in the energy input rate, which can be regarded as a response of the system as well in this case. Also the constant amplitude for low driving frequencies and the 1/ω-decay of the velocity responsewhich in leading order is corresponding to a 1/ω-decay of the energy response, as well -have been observed in the experiments. Here, we have studied the case of a modulated driving force within the GOY shell model. Also in this simulation, the response amplitude behaves basically as in the mean-field model, i.e., it decreases as 1/ω. The energy response amplitude as well as the amplitude of the energy input rate show the main maximum. Due to the above mentioned large fluctuations all higher order maxima are washed out in the GOY model.
There are two regimes in the frequency behavior of the response amplitude, namely a constant amplitude at low ω and a decreasing amplitude at high ω. The present simulations give further confidence that the crossover frequency between these two regimes gives the correct order of magnitude of the cascade time scale, i.e., in experiments it can be used to measure this time scale as suggested in reference [5].
Both models in the present study were able to reproduce the main features of the frequency behavior of the response amplitude in modulated turbulence as predicted by the mean-field model, however, both also have their shortcomings, which prevent us from correctly predicting the behavior of real turbulence in all quantitative details. Therefore, we believe, that it is worth to further study modulated turbulence numerically as well as experimentally.