Nonlinear Processes in Geophysics A weakly nonlinear model for multi-modal evolution of wind-generated long internal waves in a closed basin

A weakly nonlinear evolution model that accounts for multi-modal interaction in a small, continuously stratified lake of variable depth is derived. In particular, an evolution model for the first two vertical modes in a lake that is subject to wind stress forcing is numerically simulated. Defining modal energies, energy transfer between the first and the second vertical modes is calculated for several different forms of the density stratification. Modal energy transfer mainly occurs during reflection of mode-one waves at the vertical end walls, and it is shown that the amount of energy transfer from the first to the second mode is greatly dependent on the shape of the stratification profile. Also, the initial modal energy partition at the wind setup is shown to depend significantly on the penetration depth of the internal shear stress induced by the wind stress, especially if the stress distribution extends into the upper levels of the metalimnion.


Introduction
Thermally stratified lakes are often subject to wind stress forcing, generating basin scale internal waves that are the primary energy source for driving material transport in lakes.For modeling of such long internal waves in lakes, a simple two-layer stratification model has been preferably used since its establishment in early 20th century.The model is a reasonable approximation as long as the lake is strongly stratified (e.g., during summer) and the density stratification is confined to a thin layer between a homogeneous epilimnion (upper warm, mixed layer) and a homogeneous hypolimnion (lower cold, stagnant layer).The stratification is generally continuous, and its structure varies seasonally, with consequent seasonal effects on the evolution of internal waves (An-Correspondence to: T. Sakai (tsakai@usc.edu)tenucci et al., 2000).In a continuously stratified fluid, as is well-known from linear analysis, the vertical distribution of fluid velocities and displacement of isopyicnal surfaces possess multi-modal structures.The two-layer model accounts only for the first baroclinic mode, and it implicitly neglects all the other baroclinic modes in the field.
In the first vertical mode, the isopycnal displacement is only pronounced near the thermocline, and the entire fluid column moves vertically in the same direction.In the second mode, the isopycnal displacement is pronounced at the upper and lower parts of the thermocline, and the vertical fluid motions there take place in opposite directions, stretching and compressing the metalimnion (the thermoclinic layer between the epilimnion and the hypolimnion).The second vertical velocity mode has a single vertical node at the thermocline.There are increasing numbers of nodes for higher vertical modes, and their locations are spread increasingly toward the upper and bottom surfaces, providing even shorter vertical length scales.Many field observations have been made which capture the multi-modal nature.Early discovery of multi-modal response was achieved by Mortimer (1952).He identified the second vertical mode from a vertical temperature record in Windermere, and applied a model having three homogeneous layers to calculate the frequency of the second mode seiche.Years later, field observations in several lakes revealed internal responses dominated by the first and second modes (Wiegand and Chamberlain, 1987;Münnich et al., 1992;Roget and Zamboni, 1997;Boehrer, 2000).In small lakes, when the frequency of the second vertical mode is near that of diurnal wind forcing, a resonant response of the second mode may occur.Observations have been reported of apparent resonance in actual lakes (Wiegand and Chamberlain, 1987;Münnich et al., 1992).Responses higher than the second mode have been also reported in small lakes (seventh to tenth mode dominated response in Frains lake by LaZerte, 1980; third mode dominated response in lake Banyoles by Roget and Zamboni, 1997).The multi-modal feature is not T. Sakai and L. G. Redekopp: Long internal waves in lakes an isolated phenomenon in lakes, but has also been reported in oceanography (e.g., Bogucki et al., 2005;Gerkema, 2003).
It is routine to solve the vertical normal mode equation for a given stratification profile to obtain eigenfrequencies and eigenfunctions for particular modes of interest.However, such normal mode equations are based on linear theory, and the solutions therefore are linearly independent, yielding no information about modal interactions.Gerkema (2001Gerkema ( , 2003) ) adopted a multi-modal approach, and formulated a weakly-non1inear and weakly-dispersive multi-modal evolution system that was successfully applied to the oceanography problem.Utilizing a stratification profile in the Bay of Biscay in which the third vertical mode is dominant, he demonstrated that energy may leak (not dissipation) from the third mode to lower modes.Also, the amount of energy leakage was shown to be highly dependent on the strength of the stratification.Hüttemann and Hutter (2001) observed emergence of solitary waves of the second vertical mode when a mode-one soliton ran over a sill in a long laboratory channel.At the same time, Vlasenko and Hutter (2001) simulated the twodimensional Navier-Stokes equations in the same configuration as Hüttemann's experiment, and detailed structure of the flow field was obtained, confirming that both the first and the second mode solitons are very close to those obtained by Korteweg-de Vries (KdV) theory.
To understand the full basin energetics in lakes, it is essential to determine the modal energy distribution among dominant vertical modes in a system allowing full bi-directional propagation of the linearly independent modes.Nonlinear models are the essential for capturing inter-mode energy transfer.In this paper, we derive a weakly-nonlinear, wind-forced evolution model by applying the multi-modal approach that yields an evolution equation for each vertical mode with inter-modal interactions through nonlinear terms.For fundamental study, we limit the vertical modes in the model to the first two modes which are energetically dominant in many cases.The model is numerically solved, and we study the modal energetics for various parameters of the modeled stratification and the wind forcing.

Derivation of a nonlinear, multi-modal system
We consider a closed basin containing continuously stratified fluid.To isolate the physics from the effect of the earth's rota-tion, we neglect Coriolis acceleration of fluid elements.Assuming width of the lake is sufficiently narrow and uniform, we start with the two-dimensional, incompressible, Boussinesq approximated, inviscid equations of motion that are perturbed from the basic state of hydrostatic equilibrium: (1) In these equations, subscripts denote partial derivatives, p is the density normalized pressure, τ is the density normalized horizontal stress to account for wind forcing and bottom friction along the horizontal (x) axis, σ is the perturbed buoyancy (σ =ρg/ρ 0 ), and N is buoyancy frequency defined by using a reference density ρ 0 , static density ρ s (z) and gravity g as We assume the field domain that is enclosed by a nondeformable upper surface, vertical end walls and nonuniform (variable depth) bottom surface (Fig. 1).
Solutions of Eq. ( 1) are dictated by slip-free, impermeable and non-deformable boundary conditions: (3) The length of the lake is denoted by L, and h b is the variable depth of the fluid.We introduce a long wave scaling to horizontal space and time coordinates, and also define a slow space coordinate ξ , where µ∼h/ l 1 is the long wave scaling parameter with h a depth scale (e.g., thickness of surface mixed layer) and l is typical (long) wave length scale.We assume that the topography varies slowly in space (i.e., h b is a function of ξ only) so that topographic interaction terms will appear in the second-order approximation.We then expand the dependent variables in an asymptotic series of the form where ∼a/ h 1 is the amplitude parameter, a representing a typical amplitude of long internal waves.Scaling of the stress τ is intentionally taken to be second-order in the Nonlin.Processes Geophys., 16, 487-502, 2009 www.nonlin-processes-geophys.net/16/487/2009/ amplitude parameter so that wind forcing or benthic frictional dissipation appears first in the second order balance (i.e., weak forcing and weak damping), as do the leading effects of variable depth.We introduce the familiar KdV scaling (i.e., =µ 2 ), in order to balance weak nonlinearity and leading-order non-hydrostatic correction at the same level of approximation.
Transforming independent variables using Eq. ( 4), and substituting Eq. ( 5) into Eq.( 1), the leading-order balance gives the linear equation set T − N 2 w (1) = 0. (6) From Eq. ( 6) one can obtain a single equation in favor of w (1) : One can seek, in the sense of a consistent asymptotic approximation, a slowly-varying normal mode solution in the form where φ=0 at z=0 (upper surface) and z=−h b (bottom surface) to satisfy leading order boundary conditions.The dependence of the eigenfunction on the slow longitudinal coordinate ξ arises from the inhomogeneity of the wave guide in the propagation direction (i.e., variable depth h(ξ )).We note that the analysis can be readily extended to include a slowlyvarying wave guide width, as in the single-mode model presented in our unpublished report, but we choose to not include that further complication in the present study.Substituting Eq. (8) into Eq.( 7) one can construct a standard boundary value problem along the vertical line for every ξ : where c n is an eigenvalue, φ n is the corresponding eigenfunction, and primes denote partial derivatives with respect to z.
The corresponding orthogonality relation is 0 and where δ mn is Kroneker's delta.
All dependent variables are now expanded using the consistency implied by Eq. ( 6): Substituting Eq. ( 11) into Eq.( 6), employing Eq. ( 10), and eliminating W (1) n and P (1) n , we obtain the coupled pair of evolution equations As evident in Eq. ( 11) Z n is the modal isopycnal amplitude and U n is the modal amplitude for the horizontal velocity.Eq. ( 12) defines a set of independent, linear, bi-directional waves propagating with their respective eigenspeeds c n .
Proceeding to the next order balance using Eq. ( 5) leads to the inhomogeneous set = −{(u (1) σ (1) ) X + (w (1) σ ( 1) ) z }. (13) Note that the leading-order stress term τ (1) appears and that the boundary condition in the vertical direction gives The term u (1) ξ in the first equation in Eq. ( 13) is the leading-order effect of slowly varying depth; the bracketed terms in the second equation (x-momentum) contain the leading-order nonlinear acceleration; the term w (1) T in the third equation (z-momentum) is the leading-order nonhydrostatic correction; and the bracketed terms in the last equation (continuity) define the leading-order buoyancy flux correction.
We expand the second-order variables u (2) , p (2) and σ (2) in the same manner as in Eq. ( 11), albeit w (2) can not be expanded by φ n because φ n does not satisfy Eq. ( 14).It is not necessary to expand w (2) to define the evolution at this level of approximation.Substituting the expansion of the dependent variables into the conservation of mass equation (the T. Sakai and L. G. Redekopp: Long internal waves in lakes first equation in Eq. 13), multiplying by φ n , then integrating over the physical z domain and using Eq. ( 10), one obtains Using Eq. ( 14), the integral term on the left side of Eq. ( 15) can be evaluated via integration by parts 0 and w (2) can be eliminated by using the last equation in Eq. ( 13).Through this procedure Eq. ( 15) yields an evolution equation for the isopycnal amplitude function Z n .The evolution equation for U n can be obtained by substituting the expansions into momentum equations in Eq. ( 13) and using Eq. ( 10).After some algebraic manipulation, the evolution of the leading-order field variables is obtained in the form: (1) iX U (1) j and We point out that the first-order stress τ (1) has been divided into two parts: τ (1) s is the wind stress at the lake surface, and τ (1) b is the bottom shear stress.The coefficients appearing in Eqs.(17,18) are defined by the following set of relations: where The term U (1) iX Z (1) j appearing in Eq. ( 18) (i.e., continuity equation) derives from the vertical buoyancy flux term, and one observes that the equation cannot be integrated with respect to X because of the presence of this term.Hence the velocity and isopycnal amplitudes can not be decoupled into a single, second-order-in-time wave equation.Combining Eq. ( 12) and Eqs.(17,18), and transforming the independent variables back to their non-scaled form, we obtain the weakly-nonlinear evolution equation set: and where The r ni coefficient, containing effects of variable depth (alt., spatially varying eigenvalue), can be further evaluated by using Eq. ( 9) and Eq. ( 10).The final expression is given here without derivation: In evolution equations Eqs.(20, 21) there are essentially three kinds of modal-coupling terms: nonlinear ones, nonhydrostatic ones, and topography induced ones.The nonhydrostatic coupling in Eqs.(20, 21) derives from the fact that the hydrostatic representation of the modal structure is no longer valid at second order.The linear terms scaled by the coefficients r ij and s ij appearing on the right hand side of the equations capture the effect of self-modal distortion and cross-modal transfer resulting from variable depth effects.That these terms appear at this order is a direct consequence of the slowly-varying depth assumption, which in turn is required in order to derive a rationally-based evolution system where the underlying dynamics is represented by the linear modes of the system.A more general derivation of the topographic coupling terms for uni-directional wave propagation over topographies that vary both along and transverse to the propagation direction has been developed by Griffiths and Grimshaw (2007) in their study of internal tide at continental shelf regions.
The wind stress τ s can be expressed in terms of a friction velocity u * 0 and a prescribed, dimensionless stress distribution function F (x, t).
The vertical distribution of the horizontal stress induced by the wind, varying from its surface value τ s , is also needed to determine the coefficient N 2 φ n τ in Eq. ( 20).Assuming only the wind stress contributes to the integral, we model τ using a static, vertical stress distribution function τ (z) where τ is dimensionless, and τ (0)=1.The bottom stress τ b can be modeled assuming that the boundary layer is turbulent and using a friction coefficient C f : The value of C f for shallow water flows is typically quite small, of order 10 −3 (cf., Baines, 1995), with recommended value C f =0.0025 for weakly nonlinear long wave theory (Grimshaw, 2002).We fixed C f at this value throughout this study, and the value of u b in Eq. ( 26) is the inviscid, waveinduced velocity at the bottom surface.We take the contribution of bottom friction to a particular mode in the form 3 The two-mode evolution model In this study we limit the number of active vertical modes to the lowest two (V1: mode-one; V2: mode-two).This restriction is made in order to reduce the complexity of the evolution model while retaining the energetically-dominant modes in the system.To reduce the number of free parameters in the evolution equation, we introduce non-dimensional variables by use of the following scales: x by L; Z n and z by the surface mixing layer thickness h 1 ; t by 2L/c 0 , where c 0 is a reference phase speed taken as a spatial average of V1 phase speed; c n (or c 0 ) by N 0 h 1 , where N 0 is the maximum buoyancy frequency; and U n by N 0 h 2 1 .After recasting the Eqs.(20, 21) in a dimensionless form, we have an evolution equation set for V1 in the form: and and corresponding set for V2 in the form: and The variables U and Y in these equations are the velocity amplitude and the isopycnal amplitude of V1, respectively, and V and Z are the V2 counterparts.The terms on the right hand side, as for example in Eq. ( 28), represent nonlinear self-modal and cross-modal interactions, non-hydrostatic (dispersive) effects, variable depth effects, followed by wind forcing and bottom frictional damping, respectively.
T. Sakai and L. G. Redekopp: Long internal waves in lakes The shallow water parameter S, defined as S=h 1 /L, enters as a quadratic scaling factor in the dispersive terms appearing in Eq. ( 28) and Eq. ( 30), and as an inverse scaling factor multiplying the bottom friction terms in the same equations.Also appearing in Eq. ( 28) and Eq. ( 30) is the Wedderburn number W .This parameter is inversely proportional to the magnitude of the wind stress (cf., Imberger andPatterson, 1990 andHorn et al., 2001), is defined by the relation The Wedderburn number measures the baroclinic pressure gradient relative to the vertical gradient of the imposed wind stress represented in terms of the water friction velocity u * 0 .Using measured values of the drag coefficient induced by wind blowing over a wavy free surface, as summarized by Phillips (1966), typical values of the Wedderburn number can be roughly estimated by the relation W ≈10 5 SU −2 10 , where U 10 is the wind speed in ms −1 measured at an elevation of 10 m above the mean water surface.Data from both laboratory experiments and field measurements reveal that significant internal wave dynamics are observed when the Wedderburn number lies in the range 1<W <5 (cf., Horn et al., 2001).The quantity ksn in Eq. ( 28) and Eq. ( 30) is the modified wind forcing factor defined as With this evolution model in hand, we set up the vertical structure which qualitatively represents typical stratification profiles in lakes.We adopt a three-layer, continuously varying density structure comprising a well-mixed layer (epilimnion) of thickness h 1 , a thermoclinic layer (metalimnion) of thickness h 2 having uniform density gradient, and a weakly-stratified deep water (hypolimnion) of thickness h 3 .The square of the corresponding buoyancy frequency is expressed by a simple formula in a dimensionless form: where δ is a smoothing parameter across the interface between the metalimnion and either the epilimnion or hypolimnion, and N 2 h is the buoyancy frequency in the hypolimnion relative to the value in the metalimnion.
We present in Fig. 2 selected profiles of N 2 (z) for different h 2 , N 2 h and h 3 , and their corresponding eigenfunctions of V1 (φ 1 ) and V2 (φ 2 ).
The eigenfunctions are normalized by their maximum values.If there is no stratification in the hypolimnion (Fig. 2a, c), φ 1 attains the maximum value within the metalimnion, and φ 2 has extremal values near the top and the bottom portions of the metalimnion, implying that isopycnal displacements of V1 and V2 are both pronounced in the metalimnion, which is stretched and squeezed by V2.A slight increase in the stratification of the hypolimnion leaves the shape of φ 1 nearly unchanged, but φ 2 is significantly altered having its maximum value shifted downward into the hypolimnion (Fig. 2b).Stratification in the hypolimnion enhances vertical displacements in the hypolimnion via V2, and also enhances the horizontal motions at the lake bottom where the gradient of φ 2 is maximum.It is evident, therefore, that weak stratification in the hypolimnion can significantly enhance benthic stimulation from wind-forced V2 internal waves.
The most notable result is the sensitivity of the coefficients with respect to the thickness of the metalimnion (see Table 1).Some of the coefficients (µ 222 , ν 211 , µ 122 ) become very large as the metalimnion thickness decreases.Selfnonlinear coefficients of V2 are larger than their V1 counterparts by roughly an order of magnitude.This is due to the fact that the gradient of φ 2 in the middle of the metalimnion be- comes larger for thinner metalimnion (see Fig. 2a).Depending on the wave amplitudes, the appearance of these large coefficients can cause the corresponding nonlinear terms to be larger than linear terms.The asymptotic assumption that was used to derive the evolution model then becomes disordered, necessitating that the model be restricted to wind-forcings that yield smaller amplitudes.In fact, when the evolution model was simulated (numerical method and run configuration are briefly described in Sect.4), numerical instability was encountered as the strength of the wind stress forcing was increased.In Fig. 3, for example, the threshold Wedderburn number for achieving stable numerical integration up to two V1 seiche periods as a function of the metalimnion thickness for a fixed total depth is plotted.For thinner metalimnion, our numerical code is not capable of performing long-time integration for strong wind forcing.We also found during numerical testing that the instability is pronounced for higher numerical resolution.When all the nonlinear coupling terms between the two modes are turned off, numerical integration becomes stable even for strong wind forcing.The precise mechanism of how the instability is triggered is not straight forward.We conjecture that the large nonlinear coefficients enhance excessive energy transfer between V1 and V2 (i.e., widely disparate length scales), and also generate excessive energy levels at high wave numbers.
Based on the data presented by Horn et al. ( 2001), showing that the regime of active, wind-driven, V1 nonlinear internal wave motion occurs when 0.2<W −1 <1, and the representation of the Wedderburn number in terms of the shallow water parameter S and the 10 m wind speed U 10 given in the discussion following Eq.( 32  wave dynamics provided the metalimnion is not too thin (i.e., h 2 ≥0.8, say).Adding stratification to the hypolimnion decreases the magnitude of the nonlinear coupling coefficients (see Table 2), making the model less nonlinear and, hence, allowing larger wind energy input to the system.Reducing the thickness of the hypolimnion decreases the magnitude of V1 nonlinear coefficients, while V2 counterparts change only slightly.It is well known from KdV theory that the weak self-nonlinearity is identically zero for a completely symmetric vertical structure.This is realized in Table 3 for h 3 =1 (h 1 =h 2 =1), albeit values of the cross nonlinear terms (µ 112 , µ 121 , • • •), which are not included in the table, are non-zero for this configuration.

Wind forced response of two-mode model
When the wind stress is applied over the lake surface, a horizontal shear stress progressively penetrates across the epilimnion.Heaps and Ramsbottom (1966) found an analytical solution for the shear stress distribution by solving the twodimensional, linear hydrostatic equations.In this case, the wind induced shear stress decreases to zero linearly across the epilimnion.For a continuously stratified field, however, the shear stress distribution is not well known.The most commonly assumed form takes the shear stress diminishing linearly to a zero value at the base of the epilimnion, and the stress is zero beyond (e.g., see Monismith, 1987).With this assumption, the stress term N 2 φ n τ in the wind stress factor ks is identically zero.In this study, we adopt the linear stress function, but allow the stress to penetrate the metalimnion (see Fig. 4).
The stress function is expressed by introducing a stress penetration depth h s as a parameter: Furthermore, in this study the wind stress function F (x, t) is uniform over the surface of the lake and, is switched on and off in time.Lake response to variations in the spatiotemporal character of the wind forcing was studied in some detail in the uni-modal study by authors (unpublished report).
Equations (28-31) are solved numerically using the 4thorder compact finite difference scheme (Lele, 1992) for spatial discretization, and a forward-in-time 3rd-order Adams-Bashforth scheme.The 4th-order compact filter (Lele, 1992;Slinn and Riley, 1998) is applied every 10 time steps for dealiasing and stabilization.Numerical resolution that follow was chosen using 1025 points for spatial domain and a time step t=5×10 −5 .This high resolution configuration, in conjunction with the spatial descretization scheme with "spectral-like" resolution, sufficiently resolves steep nonlinear fronts and oscillatory waves.
Figure 5 shows the evolution of isopycnal amplitudes (Y and Z) at several selected times following the initiation of a wind event in a lake of uniform depth with h 2 =1, h 3 =3.
In what follows, a smoothing parameter δ=0.1 used in Eq. ( 34), and a shallow water parameter of S=1/500, defining the basin length in terms of the upper mixed layer depth, are used except as otherwise noted.Uniform, rightward wind stress of W =1.5 is applied for the first quarter V1 seiche period (t=1/4), and the wind is turned off thereafter.The wind stress penetration depth is chosen as h s =1 (no penetration to the metalimnion).The model was integrated from an initial condition at rest.Looking at V1, at the end of the wind setup time (t=1/4), the isopycnal surface is tilted almost linearly across the domain.After the wind setup, the surface tends to return dynamically to its initial equilibrium state, forming a basin scale seiche.After one-seiche period, a nonlinear wave front develops (as indicated by arrow in the figure), and the front steepens as it propagates.When the front becomes steeper (about t=1.75), an oscillatory wave packet forms behind the front, owing to the realization of an approximate balance between the weak nonlinearity and the leading non-hydrostatic effect.The oscillatory waves spread as the degenerating front moves back and forth in the domain.
The initial response of V2 appears near the end walls.It should be noted that the V2 response immediately following the wind setup has an oppositely-signed displacement of the reference isopycnal relative to the V1 response, but of comparable magnitude.The shorter horizontal scale of the the V2 displacement at t=1/4, as compared to the V1 displacement, is a direct consequence of their different long wave phase speeds.
The negative volume on the left side is steepened, and it evolves into a high wave number oscillatory wave packet as it propagates toward the basin interior.The wave phase speed of V1 is about one-third of V2 (c 1 =0.939; c 2 =0.284).A distinct wave packet appears at t≈2 for V1, and t≈0.75 for V2.Since the self nonlinearity of V2 is much stronger than that of V1 (see Table 1, µ 111 =0.554; µ 222 =3.312), wave lengths in the V2 wave packet are much shorter than their V1 counterparts.More interesting, when V1 waves reflect from the end walls, footprints of V1 waves are evident in the V2 domain just during the V1 wave reflection process (e.g., see V2 panels at t=1.75, 2.25, 2.75 in Fig. 5), which implies that energy is transferred from V1 to V2 during V1 wave reflections.Footprints of V2 waves can also be seen in the V1 domain, but their amplitudes are very small and not so significant energetically.
The space-time dynamics associated with the reflection of V1 serves as an effective generator of V2.This arises because the "stagnation" of V1 induces a bulging of the metalimnion, which at leading-order resembles a V2 modal distortion of the isopycnals.There may well be a generation of higher modes in this reflection precess, but the reflection of V1 (alt., the collision of V1 waves) will principally induce an energy transfer to V2 so long as the peak of the V1 eigenfunction and the nodal point of the V2 eigenfunction are positioned near the mid-point of the metalimnion.The symmetry/anti-symmetry of these modes will be shifted by the presence of stronger stratification in the hypolimnion, whereupon one expects greater energy flow to higher modes (V3, V4, etc.) during reflection.
The numerical result discussed above, however, should be looked at with caution.The field response for both modes depends quite sensitively on the wind stress penetration depth h s .Figure 6 shows several corresponding fields for evolution when the wind stress penetrates down to the mid-level of the metalimnion (h s =1.5).
T. Sakai and L. G. Redekopp: Long internal waves in lakes The response of V1 is qualitatively similar to that in Fig. 5, but the response of V2 is quite different, exhibiting a substantially diminished wind energy input to V2.Hence, there is no generation of a V2 nonlinear front.Energy transfer from V1 to V2 during wave reflection is still clearly observed.
It is instructive to define the modal energy and to quantify aspects of the modal energetics.From the Euler and continuity equations in the Boussinesq limit, one can show (e.g., see Gerkema, 2003) that the conserved energy density is where The energy density given in Eq. ( 36) has the familiar structure as a sum of kinetic energy and potential energy.However, the potential energy is expressed by an infinite series.For a reasonable calculation of the energy, we define the total energy in the system with a variable buoyancy frequency by the integral relation where only a leading order correction of the potential energy for non-uniform N 2 (z) is included.Field variables are expressed by using two vertical modes: Substituting Eq. ( 39) into Eq.( 38), and evaluating the integral assuming the lake depth is uniform, one obtains: where • • • denotes integration over the horizontal domain, and all coefficients are related to coefficients in Eqs.(28)(29)(30)(31).Terms in Eq. ( 40) are selectively grouped into V1 or V2.
In each group, the first bracketed term represents the kinetic energy and the last bracketed term represents the potential energy.
Using Eq. ( 40), modal energies were calculated at the end of wind forcing (t=1/4) as a function of h s , and results are exhibited in Fig. 7.
Energies are normalized by the total energy for evolution with a stress penetration depth corresponding to h s =0.Energy input to V2 dramatically decreases as the wind stress penetrates the metalimnion.The V1 energy also decreases as the stress penetrates down, but the change in energy is much less than that in V2.In the same figure, values of the wind forcing factor ksn are shown for the same range of h s .The forcing factor and modal energies exhibit the same trend.The forcing factor of V2 decreases to zero as the wind stress penetrates down to the half depth of the metalimnion.In this case,  the wind energy is not injected to V2 directly.Energy of V2 does not vanish completely, however, because energy is still transferred from V1 to V2 during the wind forcing through nonlinear coupling.In Fig. 7, corresponding results for thick metalimnion (h s =2) are also presented.Their trends are the same as those of former case, but the decrease in modal energies as h s increases is slower.Although it is not shown in the figure, the forcing factor | ks | vanishes as the wind stress penetration increases to near the half-depth of the metalimnion.Modal energy partition becomes more sensitive as the metalimnion becomes thinner.
In Eq. ( 40), the most interesting terms are ZY 2 in V1 and Y Z 2 in V2.These terms are the correlation between the potential energy of one mode and isopycnal amplitude of another mode.We call these terms the energy transfer terms.Of course, energy transfer is processed through 'all' dependent variables which are governed by the evolution equation set, but the energy transfer terms solely provide explicit modal energy exchange among all the other energy terms in Eq. ( 40).Another type of the modal interaction U x V x with non-hydrostatic coupling coefficients (i.e., d 12 and d 21 ) can be equally distributed to both modes (note the identity (c 2 1 /I 1 )d 12 =(c 2 2 /I 2 )d 21 in Eq. 40), hence there is no explicit energy transfer through this term.Here we define, for convenient quantification purposes, the amount of explicit modal energy transfer E tr as a difference of the energy transfer terms where α 1 and α 2 are abbreviated representations of the energy transfer coefficients We interpret that if E tr >0, the amount of energy |E tr | is transferred from V2 to V1, and vice versa for E tr <0.
Figure 8 shows modal energies and E tr as a function of time for h 2 =1 and h 2 =2 with h s =1.
Energies are normalized by the total energy at the end of forcing (t=1/4).The total energy is not necessarily conserved after the wind forcing, because Eq. ( 40) is still an approximation, and the evolution model includes bottom friction damping.Fluctuation amplitudes of the total energy for h 2 =1 is larger than that for h 2 =2.This quite probably occurs because the nonlinearity of the evolution model for a thinner metalimnion is larger, requiring a higher-order correction in the energy expression.Looking at the h 2 =2 case, energy damping due to the bottom friction is negligible.Total energy for the h 2 =1 case seems slightly damped due to the numerical filtering to suppress high wave number noises that arise from larger nonlinearity of V2.In all cases, the modal energies oscillate in time, and they are out of phase.The amount of energy transfer also oscillates in every half V1 seiche period.Furthermore, the energy is transferred from V1 to V2 for most of the time (E tr <0).Comparing with Fig. 5, the energy transfer occurs when V1 waves reflect against the end walls, leaving their footprints in V2 domain.When V1 waves leave the wall after the reflection, the energy transferred into V2 during reflection is returned to V1, with no permanent energy transfer between the modes.Figure 9 shows vertically integrated potential and kinetic energy densities of each mode during V1 wave reflection at the right end wall (x=1).We chose h s =1.5 with h 2 =1 to focus more particularly on the energy transfer from V1 to V2 during reflection by suppressing initial energy input to V2.In the V1 packet, the potential energy is larger than the kinetic energy by more  than a factor of two.In the V2 packet, the potential energy is much larger than the kinetic energy which is almost negligible.During V1 wave reflection, the potential energy of V1 dominates near the end wall, because the isopycnal amplitude increases due to superposition of incident and reflected waves and, also, the horizontal velocity (kinetic energy) approaches zero at the end wall.The dominance of potential energy in V1 is preferably transferred to V2.
In Table 4, the maximum amount of the energy transfer from V1 to V2 and values of the modal energy transfer The energy transfer becomes larger for smaller h 2 or N 2 h , as indicated by the fact that α 1 becomes much larger than α 2 .Values of α 1 are much larger than α 2 for all choices of h 3 , indicating a dominance of energy transfer from V1 to V2.The coefficients ν 211 and σ 211 in α 1 arise from the vertical integral of the buoyancy flux terms (σ u and σ w) of V1, which are coupled in the V2 evolution equation (Eq.30).ν 211 arises from the vertical buoyancy flux (σ w), and σ 211 arises form the horizontal buoyancy flux (σ u).As seen in Table 1-3, ν 211 is much larger than σ 211 .This implies that the vertical buoyancy flux of V1 plays an important role in transferring energy from V1 to V2. Especially, a thinner metalimnion corresponds to larger ν 211 , resulting in more pronounced energy transfer to V2.It can also be observed from Table 4 that wind stress penetration into the metalimnion gives little change in the amount of energy transfer.In Fig. 10 we show the amount of energy transfer from V1 to V2 as a function of the Wedderburn number for selected profiles of buoyancy frequency.The fractional amount of energy transfer increases quadratically with respect to the inverse of the Wedderburn number for all cases.

Discussions
In order to provide some insight as to nominal time scales one might encounter in applications, we present in Table 5 data for the V1 and V2 seiche periods in lakes with a Brunt-Väisälä frequency of 0.02 [s −1 ] in the metalimnion, and with respective depth and length scales as given.The seiche period for V1 is the normalizing time scale for dynamics shown in Fig. 5, and in all following figures.It is evident that the V2 seiche period is typically a factor of three or four times longer  than that of V1.However, the nonlinear steepening of the V2 front arising from a wind event clearly occurs much earlier in time than does the steepening of the V1 front.As a consequence, the self-modal nonlinear interaction of V2 is continually propagating through a variable background formed by the larger-scale V1 field.Sustained wind for periods on the order of one-eight to one-half V1 seiche period leads to significant energy deposited into the basin-scale internal wave field, which subsequently cascades to smaller scales through process studied here.
Examination of our evolution model shows that the modal energy partition to V1 and V2 depends greatly on the penetration depth of the wind stress into the metalimnion.Monismith (1987) observed in his laboratory experiment the dominant response of V1 without appearance of V2 for strong wind stress forcing.He also observed that V1 and V2 coexist in response to weak wind stress forcing.This suggests that the wind stress penetration depth may depend on the strength of the wind stress (Wedderburn number), and such vertical penetration of the wind stress in conjunction with the vertical stratification profile determines the initial energy input to each vertical mode.
Recently, Stashchuk et al. (2005) conducted numerical experiments to simulate the fully-nonlinear, baroclinic response of a near two-layer fluid in a rectangular, long, narrow tank, an experimental configuration similar to that employed by Horn et al. (2001).Stashchuk et al. (2005) discovered that the metalimnion (interfacial) layer thickens behind the nonlinear wave front (see Figs. 6 and 7 in their paper).They concluded that the layer thickening is attributed to the interaction of the wave front with the vertical side walls (see Fig. 8 in their paper).They did not conclude that the widening is attributed to V2, due to the fact that the field velocity profile did not comply with typical V2 eigenprofiles.However, V1 and V2, or even higher vertical modes, can coexist in the field, and they are superimposed, at least in a linear sense.This requires a decomposition of the disturbance field into vertical modes, which indeed the present multi-modal approach does, to determine the modal energy spectrum.Their results exhibited that the widening part of the density layer persistently follows the tail of the nonlinear wave train.If this layer widening is attributed to V2, this phenomenon implies the permanent energy transfer from V1 to V2, which is not observed in the model simulations presented here.Their density layer is very thin (about 0.2h 1 ), and the initial V1 wave amplitude is large (0.9h 1 ).This configuration, however, is outside the operational range of the present weakly-nonlinear model.
Although the variable depth terms in the evolution model admit modal energy transfer, our model simulations exhibited that inter-mode transfer is negligible in variable depth cases.For example, Fig. 12 shows evolution of isopycnal amplitudes in a lake having a sloping topography as depicted in Fig. 11.
Wind stress penetration depth was chosen to be h s =1.5 to suppress energy input to V2 in order to see modal interaction during V1 waves propagating over the sloping bottom.At t=3, V1 waves are about to pass over the slope, but at this time no significant interaction is observed in the V2 packet.Modal interaction during V1 wave reflection is rather outstanding as observed at t=2.75 and t=3.25.In Fig. 13, corresponding pictures for the uniform depth configuration are shown.Wave amplitudes are larger than the former case because the wind forcing factor k sn and nonlinearity are larger in the uniform depth case.It can be observed that there are a series of small footprints of the V1 oscillatory waves in the V2 domain, but the V1 wave reflection contributes a much larger footprint to the V2 domain.Qualitative structures of the modal components of the wave field in the variable depth case and the uniform counterpart are still very similar.
The weak cross-modal transfer due to varying depth was explored further by examining the coefficients κ ij and λ ij appearing in Eqs.(28-31) for a range of topographic slopes.Using a linear slope prescribed by h 3 =h 30 −(h 3s / l s )x, where h 30 , h 3s and l s are the maximum depth of the hypolimnion, the change in hypolimnion depth, and the horizontal length of the slope, respectively, coefficient values were computed and presented in Table 6.The coefficient values shown pertain to a configuration with h 30 =5h 1 , l s =L/3, h 2 =h 1 , and N 2 h =0.Furthermore, coefficient values are reported at the mid-slope depth position (i.e., at x s =l s /2).One notes that the coefficient values in Table 6 are significantly less than unity, except for the self modal coefficient κ 11 appearing in the V1 velocity evolution equation Eq. ( 28), and for the cross modal coefficient λ 12 appearing in the V1 isopycnal evolution equation Eq. ( 29).Corresponding ratios of the cross-modal to the self-modal coefficients are shown, and the ratios are typically order one except for those corresponding to the V1 amplitude evolutions.This suggests that in the present slow-topography formulation the energy transfer due to variable depth, even though relatively small, can occur in the V1 field principally through the V1 self-modal, topographic interaction, and through the V2-induced, crossmodal interaction if the amplitudes of the V2 field are sufficiently large.

Fig. 4 .
Fig. 4. Penetration of the wind stress into the metalimnion.

Fig. 5 .
Fig. 5. Evolution of the isopycnal amplitude for (a) mode-1 and (b) mode-2 with the wind stress penetration depth h s =1.The abscissa covers the full length (0, 1) of the basin in the scaled x-coordinate.

Fig. 6 .
Fig. 6.Evolution of the isopycnal amplitude for (a) mode-1 and (b) mode-2 with the wind stress penetration depth h s =1.5.The abscissa covers the full length (0, 1) of the basin in the scaled x-coordinate.

Table 1 .
Coefficients of selected nonlinear terms for different metalimnion thickness h 2 .The thickness and the stratification of the hypolimnion are fixed as h 3 =3 and N 2 h =0, respectively.The smoothing parameter is chosen as δ=0.1 for h 2 =1 and 2, δ=0.05 for h 2 =0.5, and δ=0.025 for h 2 =0.25.

Table 2 .
Coefficients of selected nonlinear terms for different stratification N 2h in the hypolimnion.The thicknesses of the metalimnion and the hypolimnion are fixed as h 2 =1 and h 3 =1, respectively.The smoothing parameter is fixed as δ=0.1.

Table 3 .
Coefficients of selected nonlinear terms for different hypolimnion thickness h 3 .The thickness of the metalimnion, the stratification of the hypolimnion and the smoothing paramter are fixed as h 2 =1, N 2 h =0 and δ=0.1.
), it is clear that the present model is capable of simulating stable, two-mode internal Model stability limit in inverse Wedderburn number W −1 as a function of the metalimnion thickness h 2 for a fixed total depth (h 1 +h 2 +h 3 =5); (h s =1.0).

Table 4 .
The modal energy E tr transferred from vertical mode-1 to mode-2 and the energy transfer coefficients α 1 and α 2 for various h 2 , N 2 h and h 3 (W =1.5.)

Table 6 .
Coefficients of variable topography terms for different values of the slope height h 3s .