@prefix vivo: .
@prefix edm: .
@prefix ns0: .
@prefix dcterms: .
@prefix skos: .
vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Civil Engineering, Department of"@en ;
edm:dataProvider "DSpace"@en ;
ns0:degreeCampus "UBCV"@en ;
dcterms:creator "Imam, Yehya Emad"@en ;
dcterms:issued "2012-01-09T18:16:49Z"@en, "2011"@en ;
vivo:relatedDegree "Doctor of Philosophy - PhD"@en ;
ns0:degreeGrantor "University of British Columbia"@en ;
dcterms:description """The baroclinic response to wind is examined in narrow elongated lakes. The main objective is to link the excitation and modulation of baroclinic modes to lake bathymetry and stratification, temporal and spatial characteristics of wind forcing, and damping. Three lake bathymetries are examined which represent 1) variable-depth single basins with straight thalweg, 2) two-basin lakes with straight thalweg, and 3) two-basin lakes with arms at which the thalweg bends sharply. The bathymetries are examined using idealized lake forms as well as data from two Canadian lakes. Modal analysis and a three-dimensional hydrodynamic numerical model are used to analyze the baroclinic response. Also, a modal-based forced model is developed to simulate the decoupled modal responses to wind and provide a direct link between lake and forcing characteristics and modal composition of the response.
This study shows that coupling between wind-forcing spatial structure and bathymetry determines which modes are excited, while wind-forcing temporal patterns modulate the magnitude of excited modes. It is found that, when wind is near-uniform, the first horizontal mode, H1, dominates the response regardless of bathymetry, because the near-uniform wind couples with the spatial distribution of layer flow for H1. In sub-basins separated by geometric constrictions (sills and contractions) or sharp bends in the thalweg relative to wind direction, the wind induces local metalimnetic tilts that are superimposed on the domain-wide H1 tilt. The sub-basin tilts are attributable to higher horizontal modes which are equivalent to the H1 modes of the decoupled sub-basins.
Furthermore, this study demonstrates the following: 1) interbasin exchange due to H1 shifts from two to more layers due to interaction of vertical modes, 2) geometric constrictions result in strong damping of H1 which causes high forcing-response coherence and broadens the resonance bandwidth, and 3) along-thalweg depth variability in single basins increases the number of excited modes and localizes the interface shear for asymmetric basins and causes the opposite effects for symmetric basins. The findings of this study contribute to understanding the baroclinic response to wind in lakes of complex bathymetry."""@en ;
edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/39953?expand=metadata"@en ;
skos:note "MODAL DECOMPOSITION OF THE BAROCLINIC RESPONSE TO WIND IN ELONGATED LAKES by Yehya Emad Imam B.Sc., Cairo University, 2002 M.Sc., Cairo University, 2005 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2011 © Yehya Emad Imam, 2011 ii Abstract The baroclinic response to wind is examined in narrow elongated lakes. The main objec- tive is to link the excitation and modulation of baroclinic modes to lake bathymetry and stratification, temporal and spatial characteristics of wind forcing, and damping. Three lake bathymetries are examined which represent 1) variable-depth single basins with straight thalweg, 2) two-basin lakes with straight thalweg, and 3) two-basin lakes with arms at which the thalweg bends sharply. The bathymetries are examined using idealized lake forms as well as data from two Canadian lakes. Modal analysis and a three-dimensional hydrodynamic numerical model are used to analyze the baroclinic response. Also, a modal-based forced model is developed to simulate the decoupled modal responses to wind and provide a direct link between lake and forcing characteristics and modal composition of the response. This study shows that coupling between wind-forcing spatial structure and bathymetry determines which modes are excited, while wind-forcing temporal patterns modulate the magnitude of excited modes. It is found that, when wind is near-uniform, the first horizon- tal mode, H1, dominates the response regardless of bathymetry, because the near-uniform wind couples with the spatial distribution of layer flow for H1. In sub-basins separated by geometric constrictions (sills and contractions) or sharp bends in the thalweg relative to wind direction, the wind induces local metalimnetic tilts that are superimposed on the domain-wide H1 tilt. The sub-basin tilts are attributable to higher horizontal modes which are equivalent to the H1 modes of the decoupled sub-basins. Furthermore, this study demonstrates the following: 1) interbasin exchange due to H1 shifts from two to more layers due to interaction of vertical modes, 2) geometric constrictions result in strong damping of H1 which causes high forcing-response coherence and broadens the resonance bandwidth, and 3) along-thalweg depth variability in single basins increases the number of excited modes and localizes the interface shear for asymmetric basins and causes the opposite effects for symmetric basins. The findings of this study contribute to understanding the baroclinic response to wind in lakes of complex bathymetry. iii Preface For all the chapters of this thesis, I have written and modified the original manuscripts, reaching the final versions presented in this dissertation. The research work was done by myself; I have analyzed the field data, developed (or prepared) and applied the models used in this study, and obtained the results. For each chapter, members of my PhD committee, Dr. Bernard Laval, Dr. Roger Pieters, and Dr. Gregory Lawrence, provided guidance in focusing the scope, improving the structure, and edited the original and revised manuscripts. Dr. Eddy Carmack, from my supervisory committee, provided feedback on preliminary results and contributed to the final review of the dissertation. Chapter 3 was submitted for journal publication. The submission is titled “The baroclinic response to wind in a small two-basin lake,” by Y. Imam, B. Laval, and G. Lawrence. In preparing Chapter 4, I have benefited from comments by two anonymous reviewers; they suggested adding model validation which was done using data from another year, reducing the assumptions in the model which was done by improving it (e.g., incorporating cross-sectional variability), and discussing implications for the selective withdrawal facility. Preliminary results in Chapter 5 were presented in Y. Imam, B. Laval, R. Pieters, and G. Lawrence. “Characterizing the internal wave field in a large multi-basin reservoir,” Proceedings of the Fifth International Symposium on Environmental Hydraulics, Tempe, Arizona, 2007. iv Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix List of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Objectives and Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 11 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.1 Family of Curved-bottom Basins . . . . . . . . . . . . . . . . . . . 14 2.2.2 Model Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Table of Contents v 2.3.1 Baroclinic Modes . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3.2 Response to Impulsive Wind . . . . . . . . . . . . . . . . . . . . . 24 2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.4.1 Continuum between r = 2 and r = ∞ . . . . . . . . . . . . . . . . 33 2.4.2 Parameterization of Velocity Shear . . . . . . . . . . . . . . . . . . 35 2.4.3 Relevance to Laboratory Experiments . . . . . . . . . . . . . . . . 35 2.4.4 Depth of Flat-bottom Basin . . . . . . . . . . . . . . . . . . . . . 36 2.4.5 General Application . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3 The Baroclinic Response to Wind in a Small Two-basin Lake . . . . . . . . . 40 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2.1 Field Site and Data . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2.2 Two-layer Variable Cross-section Model (TVC) . . . . . . . . . . . 44 3.2.3 Multi-layer Uniform-depth Model (MLM) . . . . . . . . . . . . . . 49 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3.1 Field Observations . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3.2 Spectral Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3.3 Two-layer Variable Cross-section Model (TVC) Results . . . . . . 54 3.3.4 Multi-layer Model (MLM) Results . . . . . . . . . . . . . . . . . . 60 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.4.1 Resonance with Wind . . . . . . . . . . . . . . . . . . . . . . . . 64 3.4.2 Excitation of Higher Horizontal Modes . . . . . . . . . . . . . . . 65 3.4.3 Higher Vertical Modes and Exchange Flow . . . . . . . . . . . . . 68 3.5 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4 The Damped Baroclinic Response to Wind in a Reservoir . . . . . . . . . . . 72 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2.1 Site Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2.2 Field Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.2.3 Modal Analysis Model . . . . . . . . . . . . . . . . . . . . . . . . 76 Table of Contents vi 4.2.4 Modal-based Forced Model . . . . . . . . . . . . . . . . . . . . . 77 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.3.1 Field Observations . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.3.2 Modeled Baroclinic Response . . . . . . . . . . . . . . . . . . . . 87 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.4.1 Effect of Bathymetric Idealizations on Damping . . . . . . . . . . 92 4.4.2 Coherence and Damping . . . . . . . . . . . . . . . . . . . . . . . 93 4.4.3 Comparison to Other Lakes . . . . . . . . . . . . . . . . . . . . . 96 4.4.4 Source of Damping . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.4.5 Implications to CWRF . . . . . . . . . . . . . . . . . . . . . . . . 101 4.5 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5 The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir . . . 104 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.2.1 Field Site and Data . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.2.2 Hydrodynamic Modeling . . . . . . . . . . . . . . . . . . . . . . . 107 5.2.3 Modal Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.2.4 Modal Decomposition . . . . . . . . . . . . . . . . . . . . . . . . 112 5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.3.1 Observed Baroclinic Response . . . . . . . . . . . . . . . . . . . . 115 5.3.2 Simulated Baroclinic Response . . . . . . . . . . . . . . . . . . . 117 5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 6.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 A Morphometric Parameters of Nechako Reservoir . . . . . . . . . . . . . . . . 155 Table of Contents vii B Coherence and Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 C Baroclinic Modal Analysis in Lakes . . . . . . . . . . . . . . . . . . . . . . . 159 viii List of Tables 1.1 Comparison of modal-based forced models (MFMs). . . . . . . . . . . . . 8 1.2 Main features of models used in this study. . . . . . . . . . . . . . . . . . . 10 3.1 Parameters of the two-layer idealization. . . . . . . . . . . . . . . . . . . . 47 3.2 Parameters for simulations using the MLM. . . . . . . . . . . . . . . . . . 61 4.1 Parameters of the two-layer (N = 2) and three-layer (N = 3) stratification for modal analysis using the VDC model. For each segment of the record, the top to bottom density difference is taken to be the same for N = 2 and 3 and is evenly distributed over N− 1 interfaces. The last row in the table gives the parameters used for the 1994 TVC simulation. . . . . . . . . . . . 78 4.2 Damping parameters in several basins, as reported in or estimated from the literature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.1 Characteristics of baroclinic modes which dominated the simulated baro- clinic response. The listed modes have average and maximum modal ener- gies > 0.25 GJ and > 0.5 GJ, respectively. . . . . . . . . . . . . . . . . . . 123 A.1 Morphometric parameters of Nechako Reservoir. Shoreline development index (SD) is defined as SD = P(4piA)−0.5 where A and P are the surface area and shoreline length of the reservoir, respectively. . . . . . . . . . . . 155 C.1 Model parameters in studies of baroclinic modes in stratified lakes. Shallow zones refer to partitions of the domain with a number of density layers that is less than the maximum. . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 ix List of Figures 1.1 (a) Typical thermal stratification in summer. (b) Baroclinic response of two- layer stratification showing tilting of the density interface. (c) Baroclinic response of three-layer stratification. The tilt of the upper and lower in- terfaces changes the thickness of the intermediate layer which represents the metalimnion. The vertical scale in (b,c) is exaggerated relative to the horizontal scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Features of models used in modal analysis (MAMs). . . . . . . . . . . . . 7 2.1 Bottom profiles for (a) Baldegg Lake, (b) Wood Lake, (c) the north basin of Windermere, and (d) Mundaring Weir reservoir. Also shown are (a-c) CB and (d) WB profiles fitted to the maximum lake depth. Lake bottom-profiles were digitized from (a,c) Lemmin and Mortimer (1986), (b) Wiegand and Chamberlain (1987), and (d) Shintani et al. (2010). Interface depths were taken from the same sources. . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Schematic of the curved-bottom basin (CB) illustrating geometrical dimen- sions and density stratification. Note that the vertical scale is exaggerated. . 14 2.3 The effect of the parameters m and r on the bed profile of the curved-bottom basin. Inflection points of the profiles are indicated by dots. The profiles in panel (c) have no inflection points. . . . . . . . . . . . . . . . . . . . . . . 15 2.4 Spatial structure of (a-d) the lower-layer horizontal flow and (e-h) the in- terface deflection for the four lowest modes. Solid lines are for CB2 and dashed lines are for the flat-bottom basin CB∞. For display purposes, the distributions Q2,n (x) and Z2,n (x) are scaled to a maximum of unity. . . . . . 23 2.5 Shapes of along-thalweg wind components for (a) CB2, (b) WB2, and (c) the flat bottom CB∞ and WB∞. The shapes are scaled to a maximum of unity. . 23 List of Figures x 2.6 Comparison of the baroclinic responses to impulsive wind in CB2 and CB∞. Panels (a-d) show the interface deflection at the upwind endwall ( xL−1 =−0.5). Panels (e-h) show the flow in the lower layer at the center of the basin( xL−1 = 0 ) . Each column of panels corresponds to a damping ratio ζ1 as noted on the top. Interface deflection is normalized by the dynamic equilibrium deflection at xL−1 = −0.5. Layer flow is normalized by the magnitude of the undamped maximum layer flow at xL−1 = 0 in CB∞. Time is normalized by the period of the fundamental mode in CB∞. Note the different vertical scales for the panels corresponding to ζ1 = 1.0 and 2.0. . . 26 2.7 Ratio between maximum layer flows in CB2 and CB∞ as a function of damping. 27 2.8 Snapshots along the thalweg of CB2 and CB∞ showing (a) the interface deflection and (b) the flow in the lower layer. Snapshots are for ζ1 = 0.1. For each row of panels, snapshots in the two basins correspond to the same phase of the baroclinic response tT−11 as noted on the far left. Normalization of the interface deflection and the layer flow follow Figure 2.6. . . . . . . . 28 2.9 Along-thalweg distribution of velocity in the (a,d) upper and (b,e) lower layers and (c,f) velocity shear in (a-c) the CB∞ and (d-f) the CB2 for m= 0.5. The distributions are shown for tT−11 = 1/16, 1/8, and 1/4 (as noted on panel f) and for no damping (ζn = 0). The velocities U1 and U2 and the velocity shear ∆U are normalized by the maximum value of U1 in CB∞. . . . . . . . 29 2.10 The maximum and the spatial average velocity shear in CB2, WB2, and CB∞ (WB∞) for ζ1 = 0.1. The velocity shear is normalized by its peak undamped magnitude in CB∞, ∆Um = 4ρ1h1 (τsT1)−1 (Spigel, 1978). The normalized velocity shear in WB∞ and CB∞ are identical. Time is normalized by the period of the fundamental mode, T1, of the respective basin. . . . . . . . . . 30 2.11 (a-e) Interface deflection and (f-j) lower-layer flow at different locations in WB2 ( 0≤ xL−1 ≤ 0.5). Wind is uniform, steady, and blows in the negative x direction (i.e., the downwind end is xL−1 = 0). The deflections and flow are shown for ζ1 = 0.1. Interface deflection is normalized by the steady state deflection. Layer flow is normalized by the magnitude of the undamped maximum layer flow at xL−1 = 0.25 in WB∞. Time is normalized by the period of the fundamental mode in WB2. . . . . . . . . . . . . . . . . . . . 32 List of Figures xi 2.12 Along-thalweg distribution of velocity shear in WB2 for the lowest three baroclinic modes. Distributions are scaled to a maximum of unity. Higher modes n> 3 also have their peak ∆U at xL−1 = 0.5 (not shown). . . . . . . 33 2.13 The ratio of the periods of baroclinic modes in CB2 and a flat-bottom basin of depth H1+h2. The ratio is shown for the four lowest modes. . . . . . . . 37 3.1 Map of Amisk Lake. Isobaths are shown every 10 m from a water surface of ∼ 612 m above sea level. The 20 m and 40 m isobaths are shown in solid black. The dashed line shows the thalweg of Amisk Lake. Along-thalweg distance, x, is positive to the north with the origin, x = 0, at the south end where the thermocline intersects the bottom. Lower right inset shows the location of Amisk Lake in the province of Alberta, Canada. . . . . . . . . . 43 3.2 Temperature profiles collected on 30-June and on 28-July, 1991. Also shown is the average of the profiles and its two-layer idealization. Note that the thermocline was shallower than the channel crest. . . . . . . . . . . . . . . 44 3.3 (a) Variation of the width of Amisk Lake in the vertical and along the thal- weg. The vertically-hatched region on top represents the distribution of wind forcing along the thalweg. (b) Along-thalweg distribution of the cross- sectional areas (S1 and S2) and the interfacial width (b2). . . . . . . . . . . 48 3.4 Observations at Amisk Lake from 28 June to 30 July, 1991 (days 179-211). (a) Wind speed and direction. Wind speed is smoothed with a 1 h moving average. Wind direction is displayed for speeds greater than 2 m s-1. Dotted lines indicate the average orientation of the thalweg. (b) Water temperature at the mooring. Depths of thermistors are indicated on the records. (c) Along- thalweg water currents recorded by the ADCP. Positive velocities are from south to north. The black contour lines indicate the depth of flow reversal. Instances of two-layer flow are indicated by (+). . . . . . . . . . . . . . . 52 3.5 Histogram of total volumetric transport through the cross-section at the ADCP site. Grey bars indicate two-layer exchange (52%) and white bars stacked on top indicate exchange with three or more layers (48%). . . . . . . . . . . . 53 List of Figures xii 3.6 Power spectra of (a) the along-thalweg component of wind stress, (b) the fluctuations of water temperature at 5 m and 8 m at the mooring, and (c) the along-thalweg component of velocity in the surface and bottom bins of the ADCP. Spectra were generated after removing linear trends and were smoothed in the frequency domain (Saggio and Imberger, 1998). Dashed lines show the 90% confidence interval based on a chi-squared distribution. Arrows in panel (a) indicate the estimated periods of the H1, H2, and H3 modes. Shown in (d,e) are the magnitude-squared coherence (solid) and phase (boxes) between (d) the along-thalweg wind stress and the water tem- perature at 5 m and (e) the along-thalweg wind stress and the velocity in the bottom bin. Dashed lines show coherence at the 90% confidence level (Thompson, 1979). In all panels, Grey bars indicate the ∼ 1 d and 7.3 h periods described in the text. . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.7 (a) Wind stress over Amisk Lake. Observed and simulated (H1-H9) flow in the (b) upper and (c) lower layers at the ADCP. Observed flow corresponds to two-layer exchange. Also shown is the contribution of the H1 mode to the simulated flow. (d) Observed and simulated interface deflection at the thermistor mooring. In (a-c), stress and flow are positive from south to north. 57 3.8 Amplitude spectra of simulated steady-state oscillations for (a) flow at the ADCP, |Dn,m|Qn, and (b) interface deflection at the mooring, |Dn,m|Znω−1f ,m. Also shown are the phase differences between the wind components F(n,m)s and the associated (a) flow at the ADCP,∠F(n,m)s −∠q(n,m)2 =∠ ( Bm (Dn,mQn) −1), and (b) the interface deflection at the mooring, ∠F(n,m)s −∠q(n,m)2 − pi2 . The amplitude spectra are symmetric about ω f ,m = 0 while the phase spectra are anti-symmetric (not shown). . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.9 (a) Components n= 1−3 of wind forcing corresponding to a steady spatially- uniform wind stress of 0.02 N m-2. Shown for the H1, H2, and H3 modes are (b) the steady-state interface deflection, (c) the peak flow in the lower layer, and (d) the peak depth-averaged velocity in the lower layer. The markers q and s indicate the locations of the ADCP and the thermistor mooring, respectively. Grey shading indicates the approximate extent of the channel connecting the north and south basins. In (a,c,d), direction of forcing, flow, and velocity is positive from south to north. . . . . . . . . . . . . . . . . . 60 List of Figures xiii 3.10 (a) Three-layer stratification for simulation s2 in the box basin. (b,c) Two- layer exchange. (d) Three-layer exchange. . . . . . . . . . . . . . . . . . . 62 3.11 (a) Surface, (b) intermediate, and (c) bottom flows from the MLM simulation s2 (three-layered) at the ADCP. Surface and bottom flows from the MLM simulation s1 (two-layered) are shown in panels (a) and (c), respectively. Panel (b) also shows the observed flow in the intermediate layer when three- layer exchange occurred. Flow is positive from south to north. . . . . . . . 63 3.12 Amplitude of steady flow oscillations for the H1 mode as a function of the ratio of forcing to modal frequency, ω? = ω f ,mω−1u,n . Two sets of curves are shown, one set for variable ω f ,m and the second for variable ωu,1 (see legend). Each set shows curves for ζ1 = 1.0, 1.3, and 1.7. The flow amplitude is normalized by wind forcing and modal structure, ∣∣∣D?1,m∣∣∣. Note the flatness of the solid curves over a wide range around the modal frequency, ω? = 1. . 66 3.13 Baroclinic flows simulated using the MLM based on the stratifications of s2 (a-c), s3 (d-f), and s4 (g-i) and forced by step wind. The flows are shown for modes V1H1 (top panels), V2H1 (middle panels), and the summation of V1H1 and V2H1 (bottom panels). Flows are normalized by the maximum flow from simulation s3. Arrows in the bottom panels indicate the shift from two- to three- layer exchange. . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.1 Map of Nechako Reservoir showing the locations of Kenney Dam, the spill- way, and the powerhouse offtake. The focus of this study is the two eastern basins, Knewstubb and Natalkuz Lakes (shaded). The TVC model is applied along the thalweg of the two basins (dashed-dotted line). Marked are the locations of thermistor moorings and wind stations (KN, NA, SK). Top left insert shows the location of Nechako Reservoir within the province of British Columbia, Canada. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 List of Figures xiv 4.2 (a) Width of cross-section versus depth and distance along the thalweg. (b,c) Bathymetric parameters for the TVC model: (b) cross-sectional area above and below the density interface and (c) cross-sectional width at the level of the interface. (d) Wind forcing on day 234.5, 1994, when the average wind speed at the moorings was 8.4 m s-1 blowing from the WSW. Positive wind forcing is in the direction of increasing x. In all panels, the thalweg origin, x= 0, is at Kenney Dam. . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.3 Wedderburn number for Knewstubb and Natalkuz Lakes for observations in 1994. For calculating the Wedderburn number, the average along-thalweg wind stress τs was low-pass filtered (Stevens and Lawrence, 1997) using a fifth-order Butterworth filter with T1/4 = 1.4 d cutoff. Also shown is the Wedderburn number below which non-linear steepening becomes pro- nounced, Wn.s. (Horn et al., 2001). Wn.s. is estimated for record segments, S1-S4, using the two-layer parameters (ρ2− ρ1 and h1) in Table 4.1, the average depth below the density interface h2 = 44 m−h1, and ζ1 = 0.1. . . 82 4.4 Histograms of wind direction at (a) KN, (b) NA, and (c) SK during 24-Jun to 12-Oct (days 175-285), 1994. Data shown are for wind speeds greater than 2 m s-1. Sum of occurrences relative to record length are noted on each panel. 83 4.5 (a) Power spectra of wind speed (solid lines). (b,c) Power spectra of the deflections of the 10°C isotherm (solid lines) at the (b) KN and (c) NA moorings for record segments S1-S4 (Table 4.1). For the record segments, periods of horizontal modes (V1Hn) are indicated by circles with each mode connected by a grey line. Periods of H1 and H2 are from the VDC with three- layers and higher modes (H3-H8, increasing from left to right) are from the VDC with two-layers. (d) Magnitude-squared coherence between the W to E wind component and the level of the 10°C isotherm at mooring KN for segment S3 of the record (19-Aug to 15-Sept). For clarity, the power spectra in panels (a-c) are vertically shifted. Spectra were generated after removal of the linear trend and were subsequently smoothed to improve confidence (Boegman et al., 2003). Dotted lines indicate the 90% confidence levels. . . 84 List of Figures xv 4.6 (a) Wind components along the W to E (black) and N to S (grey) directions at platform KN, 1994. (b,c) Isotherm levels at (b) KN and (c) NA. Isotherms are spaced at 2°C intervals. The 10°C isotherm is marked by a black line. Wind and water temperature data were low-pass filtered using a 5th order Butterworth filter with 12 h cutoff. . . . . . . . . . . . . . . . . . . . . . . 86 4.7 The average relative weights, Cn, of wind components (black line). Vertical bars indicate one standard deviation on each side of the average Cn. . . . . . 88 4.8 Comparison of the observed (thin, black) and simulated (thick, grey) deflec- tions of the interface at the (a) KN and (b) NA moorings, 1994. . . . . . . . 89 4.9 Histograms of the simulated deflections for the H1 and H2 modes at the (a) KN and (b) NA moorings. The histograms are for the hourly deflections between day 175 and 285, 1994. . . . . . . . . . . . . . . . . . . . . . . . 90 4.10 (a) Wind speed components along the W to E (black) and N to S (grey) directions at platform KN, 2005. (b) Isotherm levels at mooring KN. (c) Comparison of the simulated (thick, grey) and observed (thin, black) deflec- tions of the interface at mooring KN. In panel (b), isotherms start from 6°C at the bottom and increase upwards at 2°C intervals. The 10°C isotherm is marked by a black line. Wind and water temperature data in panels (a) and (b) were low-pass filtered using a 5th order Butterworth filter with 12 h cutoff. 91 4.11 Comparison of interface deflection, Z2, along the thalweg between the TVC and the VDC with N = 2 for the (a) H1 and (b) H2 modes. The interface deflection is normalized by its maximum absolute magnitude. Shading in- dicates approximate extent of the Narrows. The thalweg origin, x = 0, is at Kenney Dam. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.12 Damping ratio ζn from VDC versus the percent change in AnZ (n) 2 from its TVC value. The damping ratios give the least RMS error compared to the calibrated TVC response for mode H1 and H2. . . . . . . . . . . . . . . . 94 List of Figures xvi 4.13 Wavelet power spectrum of (a) W to E wind speed and (b) 10°C isotherm level at mooring KN in 1994. Wavelet power is normalized by the variance. (c) Magnitude-squared coherence between the W to E wind speed and the 10°C isotherm level. In all panels, hatching masks the cone of influence (Torrence and Compo, 1998). Black contours bound zones exceeding the 95% confidence level calculated following Grinsted et al. (2004). In panel (c), arrows indicate the phase relationship between the wind and the 10°C isotherm level; in-phase deepening, _ deepening with quarter period lag, and in-phase rise in isotherm level with increasing W to E wind speed. Dashed white lines in panels (b,c) indicate the period of the V1H1 mode estimated from the VDC with N = 3. . . . . . . . . . . . . . . . . . . . . 97 4.14 Across-thalweg echo-sounding transects (a) near Kenney Dam (x= 4km) and (b) near the west end of the Narrows (x= 33km), October 2005. Note the trees protruding into the metalimnion (s), the small width of the Nechako Canyon below 40 m (a), and the irregular bathymetry within the Narrows (b). Image contrast has been improved to better reveal the trees and the metalimnion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.1 (a) Map of Nechako Reservoir. (b) Map of the eastern part of Nechako Reservoir showing Knewstubb and Natalkuz Lakes, shaded in (a). Reservoir shoreline corresponds to the high water mark, at 853.4 m above mean sea level. The solid triangles indicate the location of the thermistor moorings and floating weather stations deployed in 1994. Locations of CTD profiles, also conducted in 1994, are indicated by numbered circles. . . . . . . . . . 106 5.2 The numerical bathymetry of Knewstubb and Natalkuz Lakes which is used to ELCOM. The insets show the longitudinal profiles along the thalweg of Knewstubb and Big Bend arms (C1) and along the thalweg of Knewstubb and Natalkuz Lakes (C2). The origin of the profiles is at the north east end of Knewstubb Arm from which distance is marked at 10 km intervals (hollow circles). The locations of mooring KN and NA are indicated by hollow triangles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 List of Figures xvii 5.3 (a) Average temperature profiles on day 225 from observations at mooring KN and NA. (b) Background stratification from the mooring data during days 225.5-240 and a three-layer idealization of the stratification. . . . . . . . . 109 5.4 Spatial structures Zi,r of the upper (i= 2) and lower (i= 3) density inter- faces for the dominant baroclinic modes. . . . . . . . . . . . . . . . . . . . 112 5.5 Wind forcing and thermocline response between day 225.5 and 240 in 1994 at mooring KN (left panels) and NA (right panels): (a,b) observed wind direction, (c,d) observed wind speed, (e,f) observed isotherm depths, and (g,h) simulated isotherm depths. Isotherms in (e-h) are displayed at 1°C intervals with solid black lines denoting the 7°C and 16°C isotherms. In panels (c-h), the dotted vertical lines, P1, P2, and P3, indicate features which are described in the text. Note that peak P3 in (c) is slightly shifted to the left. Arrows in (g,h) denote the times of snapshots S1-S3 at which the simulated metalimnetic response is investigated. Dashed lines indicate the domain- wide average depths of the simulated 7°C and 16°C isotherms. . . . . . . . 116 5.6 (a) Wind direction and (b) speed during CTD profiling on 8 June 1994 (day 159). (c) Temperature profiles from the CTD downcasts. Locations of the casts in Knewstubb and Natalkuz Lakes are marked on Figure 5.1. Note the colder surface temperature and shallower metalimnion in the west side of Natalkuz Lake (casts 4 and 5). In (a,b), grey hatching indicates the duration of profiling. Wind conditions are from station SK. . . . . . . . . . . . . . . 118 5.7 Snapshots of simulated isotherms along curtains (a,c,e) C1 and (b,d,f) C2. Times of snapshots (a,b) S1, (c,d) S2, and (e,f) S3 are marked on Figure 5.5g,h. Solid black lines denote the 7°C and 16°C isotherms at the time of the snapshot. Dotted lines indicate the domain-wide average depths of the 7°C and 16°C isotherms. . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.8 Deflections of the upper (η2) and lower (η3) interfaces at mooring KN (a,c) and NA (b,d). The superposition of deflections from the dominant modes listed in Table 5.1 (black lines) are compared to deflections simulated with ELCOM (grey lines). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.9 Snapshot S2 (day 235.4) of metalimnetic deflections 0.5(η2,r+η3,r) of dom- inant modes. Black contours indicate zero deflections (i.e. nodal lines). . . 124 List of Figures xviii 5.10 Normalized modal amplitudes of interface deflection between day 225.5 and 240 for modes which dominate the baroclinic response of Knewstubb and Natalkuz Lakes. Vertical dotted lines indicate the time of snapshot S2. Modal amplitudes are normalized by the amplitude at S2. . . . . . . . . . . 125 5.11 Snapshot S2 (day 235.4) of variations in metalimnetic thickness η2,r−η3,r for dominant modes. Black lines indicate no changes in thickness. The modal structure in (g) was classified as H5 according to the number of nodal lines in the metalimnetic deflection (Figure 5.9g). . . . . . . . . . . . . . . 127 5.12 Cumulative energy input by wind (solid) and dissipated energy (dashed) for the dominant modes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.13 Damping ratios of the dominant baroclinic modes corresponding to the 10th, 50th (median), and 90th percentiles. For each mode, damping ratios esti- mated when energy dissipation rate Dr was less than 10% of its maximum value were excluded from calculating the percentiles. The solid line repre- sents a linear best fit to the median ζr as explained in the text. . . . . . . . . 128 5.14 Schematic of typical features of the baroclinic response to wind and modal contributions. Black lines indicate the deflected upper and lower bounds of the metalimnion (solid) and their initial levels (dashed). Modes which to contribute to the observed metalimnetic tilts along the sub-basins are denoted in ellipses while modes which contribute to the observed metalimnetic com- pression or expansion at sub-basin ends are denoted in rectangles. Arrows indicate along-thalweg wind direction. . . . . . . . . . . . . . . . . . . . . 130 B.1 Power spectra of (a) the n = 1 wind component integrated along the thal- weg of Knewstubb and Natalkuz Lakes, ´ L 0 F (1) s dx, and (b) the interface deflection simulated by the TVC at mooring KN due to the H1 mode for two values of the damping ratio ζ1. For generating the spectra, linear trends were removed, a cosine-taper window was used, and the spectra were smoothed in the frequency domain (Saggio and Imberger, 1998). Dotted lines show the 90% confidence interval based on a chi-squared distribution. Grey hatching denotes the period of the H1 mode. . . . . . . . . . . . . . . . . . . . . . . 157 List of Figures xix B.2 (a) Magnitude-squared spectral coherence between wind forcing and the interface defection at KN due to H1 for two values of the damping ratio ζ1. Dashed lines show coherence at the 90% confidence level (Thompson, 1979). Shown in (b,c) is the corresponding magnitude-squared coherence using wavelet transform for the two values of ζ1. White hatching masks the cone of influence (Torrence and Compo, 1998). Black contours bound zones exceeding the 95% confidence level (Grinsted et al., 2004). In all panels, grey hatching indicates the period of the H1 mode. . . . . . . . . . . . . . 158 xx List of Acronyms Bathymetry CB Curved-bottom basin WB Wedge-shaped basin FB Flat-bottom basin Modal analysis models MAM Modal analysis model TDP Two-layer Defant procedure TVD Two-layer variable-depth model VDC Variable-depth complete model with N layers TVDC Two-layer variable-depth complete model ThVDC Three-layer variable-depth complete model Modal-based forced models MFM Modal-based forced model TVC Two-layer variable cross-section model xxi Acknowledgments First and foremost, I thank Allah1 for enabling me to go this far and accomplish this work. I’m thankful to Him for blessing me with family, supervisors, colleagues, and friends who have supported me during the long path of PhD. Since he who does not thank people, hasn’t thanked God, thanks are due to all those whose path have crossed mine during the course of my PhD program and have left a positive mark on me personally or academically. In particular, I’m grateful to Dr. Laval and Dr. Lawrence for academic and financial support. Dr. Laval gave me the opportunity and space to explore and experiment with different research projects. He was patient with my writing style and provided valuable advice. I’m also grateful to Dr. Laval for introducing me to field work on lakes. I’m thankful to Dr. Lawrence for never conveying his critical comments on my manuscripts and presentations except coupled with encouragement. To Dr. Pieters, I owe what I know about laboratory work related to instrument calibration and maintenance as well as a large part of my knowledge about field work and instruments. His decency in criticizing my work was unmatched except by his careful review of my manuscripts. I’ve tremendously benefited from both qualities. Thanks are also due to all my colleagues in the EFM group for sitting through the informal EFM lunches, listening to my long presentations, and providing valuable feedback. My sincerest thanks go to my wife and colleague, Yasmin, for her encouragement and patience. I also wish to acknowledge the care and favors of my parents. If it was not for the grace of God and them, I would have accomplished nothing. Finally, I acknowledge using the Matlab wavelet transform software developed by Aslak Grinsted and the ELCOM model developed by CWR, University of Western Australia. 1Allah is the Arabic word for God. xxii To my family 1Chapter 1 Introduction 1.1 Motivation Motions induced by wind play an important role in the transport of nutrients and contami- nants within lakes and, consequently, affect water quality and lake ecology (MacIntyre and Melack, 1995; Evans et al., 2008; Pannard et al., 2011). It is therefore crucial to characterize the temporal patterns and spatial structures associated with wind-induced motions. In density-stratified lakes, wind-induced motions below the surface mixed-layer are typi- cally dominated by baroclinic basin-scale wave modes (Wüest and Lorke, 2003). Excited and energized by wind, the baroclinic modes constantly adapt to changing wind conditions and usually persist as free oscillatory motions when wind subsides (Mortimer, 1952; Monismith, 1986; Wiegand and Chamberlain, 1987). While the study of the baroclinic response to wind has a long history (Watson, 1903, 1904; Wedderburn, 1912; Mortimer, 1953), the emphasis has been on free oscillation of baroclinic modes (internal seiches) during calm periods be- tween wind storms. This emphasis was motivated by frequent observations of these oscilla- tions in many lakes (Mortimer, 1953; Lemmin, 1987; Wiegand and Chamberlain, 1987) and by the relative ease of explaining oscillations in the absence of external forcing (Mortimer, 1952). In contrast, the theory for the evolution of basin modes subjected to wind forcing was only recently extended from narrow, flat-bottom, uniform-width basins (Heaps and Ramsbottom, 1966; Spigel, 1978; Monismith, 1985, 1987) to arbitrary bathymetry (Shimizu et al., 2007; Shimizu and Imberger, 2008). Prior to this extension, several studies examined the effect of bathymetry and stratifica- tion on the periods and structures of baroclinic modes in real and idealized lakes (Hutter et al., 1983; Lemmin and Mortimer, 1986; Münnich, 1996; Roget et al., 1997; Fricker and Nepf, 2000). However, with the exception of inferences about resonance with diurnal wind in some lakes (Hutter et al., 1983; Wiegand and Chamberlain, 1987; Münnich et al., 1992; Antenucci and Imberger, 2003; Vidal et al., 2007), it was seldom explained why Chapter 1. Introduction 2 specific baroclinic modes dominated the response to a given wind in a given bathymetry and stratification. Recently, this was addressed by Shimizu et al. (2007) and Shimizu and Imberger (2008) for lakes characterized by relatively simple geometry, namely, Lake Tiberias (Kinneret) which is an elliptic single-basin lake, and Lake Biwa which is composed of a simple-shaped main basin with small bays and a baroclinically-decoupled shallow shelf. For other geometries, however, there is still a need to examine the link between the modal composition of the response and the effect of lake bathymetry and stratification, temporal and spatial characteristics of wind forcing, and damping. Here, this link is examined for elongated lakes with non-uniform bottom profiles, multiple basins, and multiple arms. Before listing the objectives of this work, the next section provides background material. Detailed review of relevant literature is presented in the introduction to each chapter (Chap- ter 2 to Chapter 5). Section 1.3 lists the study objectives and presents an overview of the structure of the thesis. 1.2 Background 1.2.1 Processes Freshwater lakes are typically temperature stratified during summer. Vertical temperature profiles usually reveal three distinct zones; the epilimnion which is often well mixed, the metalimnion which is characterized by a sharp thermal gradient, and the hypolimnion where thermal gradients are weak (Mortimer, 1952; Boehrer and Schultze, 2008) (Figure 1.1a). The main source of energy which induces basin-wide motions in stratified lakes is wind (Imbo- den, 2005). When the metalimnion is thin, wind-induced motions are well approximated by a two-layer idealization of the stratification (Figure 1.1b). However, when the metalimnion is thick, three or more density layers are needed to adequately represent the stratification and capture changes in the thickness of the metalimnion (Figure 1.1c). In a narrow, elongated lake with two-layer stratification, the lake response to wind is composed of vertical and planimetric circulatory motions that are mainly confined to the surface layer (Mortimer, 1952; Heaps and Ramsbottom, 1966; Toffolon and Rizzi, 2009). Superimposed on the circulatory motions, the free surface and density interface tilt in op- posite directions to balance the surface wind stress (Heaps and Ramsbottom, 1966; Shintani et al., 2010). This adjustment of the surface and interface to dynamic equilibrium is often Chapter 1. Introduction 3 Figure 1.1: (a) Typical thermal stratification in summer. (b) Baroclinic response of two- layer stratification showing tilting of the density interface. (c) Baroclinic response of three- layer stratification. The tilt of the upper and lower interfaces changes the thickness of the intermediate layer which represents the metalimnion. The vertical scale in (b,c) is exaggerated relative to the horizontal scale. accompanied by underdamped oscillations as observed in many lakes (Mortimer, 1953; Hutter et al., 1983; Lemmin, 1987; Roget et al., 1997; Antenucci et al., 2000). For two-layer stratification, oscillatory deflections arise from two types of standing grav- ity waves, the barotropic (external) and baroclinic (internal) waves (Gill, 1982). The period and spatial structure of baroclinic standing waves are determined by lake bathymetry and stratification (Münnich, 1996; Fricker and Nepf, 2000) and, hence, are referred to as the natural (eigen) baroclinic modes of the lake. Amplitudes of metalimnetic oscillations asso- ciated with baroclinic modes are typically three orders of magnitude larger than those due to barotropic modes (Mortimer, 1952; Gill, 1982). Moreover, the baroclinic modal periods are much longer (Hutter et al., 1991). For these reasons, baroclinic modes are more important to the relatively slow metalimnetic dynamics induced by wind. In lakes with a large lateral extent, the rotation of Earth (Coriolis force) plays a role in the baroclinic response to wind, particularly in lakes further from the equator. Rotation affects the periods and spatial structures of baroclinic gravity modes and admits another class of waves, the vorticity modes (Hutter et al., 2011). The importance of rotation can be assessed using the internal Rossby radius of deformation, R= cm f (1.1) where cm is the linear phase speed for the vertical mode of interest and f = 2Ωsinφ is the Chapter 1. Introduction 4 Coriolis parameter in which Ω = 7.292×10−5rad s-1 is the rotational speed of Earth and φ is the central latitude of the lake. When R is smaller than the lateral horizontal dimension of the lake B, rotation is important to the baroclinic modes and response to wind. Frequently, R and B are combined in a dimensionless ratio, the Burger number S= 2RB−1 (Antenucci and Imberger, 2001; Appt et al., 2004) such that, if S 1, rotation has a negligible effect on the periods and structures of baroclinic modes. As will be shown, rotation does not play a major role in the lakes considered here (Chapters 2 to 5). In many cases, the metalimnion possesses considerable thickness that gives rise to a second vertical baroclinic mode in which the upper and lower bounds of the metalimnion oscillate in anti-phase (Mortimer, 1952; Wiegand and Chamberlain, 1987; Münnich et al., 1992; Roget et al., 1997). Higher vertical modes are also possible. In reservoirs characterized by near-linear stratification, higher vertical modes, up to the fifth, have been observed (Pérez- Losada et al., 2003; Vidal et al., 2005, 2007). As wind forcing intensifies, the amplitudes of baroclinic modes increase and the baro- clinic response starts to deviate from linear behavior. Manifestations of this deviation include transfer of energy between vertical modes (Sakai and Redekopp, 2009) and transfer of modal energy into non-linear internal surges (Thorpe et al., 1972; Hunkins and Fliegel, 1973; Horn et al., 2001). In flat-bottom basins stratified into two-layers, Horn et al. (2001) classified the evolution of baroclinic modes based on the initial endwall displacement ηo of the tilted density interface. When the ratio of ηo to the upper layer thickness h1 is small, the evolution of the interface is well described by linear superposition of oscillations of baroclinic modes. For large ηo/h1 combined with a small ratio of h1 to total depth, baroclinic modes steepen and degenerate into non-linear surges which ultimately disperse into internal solitary waves (Boegman et al., 2005b). The energy contained in the internal solitary waves is dissipated as they break over the basin end slopes (Boegman et al., 2005a). However, the lake systems investigated here display primarily linear behavior and are well described by linear models. Another important path for dissipation of energy from baroclinic modes is through turbu- lence in the bottom (benthic) boundary layer (Gloor et al., 2000; Marti and Imberger, 2006) where turbulence and mixing are stronger than within the lake interior by one to two orders of magnitude (Goudsmit et al., 1997; Wüest and Lorke, 2003). The damping of baroclinic modes that results from energy dissipation is frequently parameterized by linear drag (Heaps and Ramsbottom, 1966; Bengtsson, 1973; Platzman, 1978; Bäuerle, 1998; Shimizu and Imberger, 2008). In this study, modal damping is parameterized by linear drag and the Chapter 1. Introduction 5 estimated damping parameters are compared for different modes and lakes. Having described processes relevant to the baroclinic response to wind, next the methods used to analyze the baroclinic response are described. 1.2.2 Methods In analyzing the baroclinic response in lakes, modal analysis models (MAMs) are com- monly used to estimate the periods and spatial structures of baroclinic modes based on lake bathymetry and stratification (Figure 1.2). The most basic MAM is the layered flat-bottom 1D model which gives the natural periods Tn of longitudinal baroclinic modes as (Münnich et al., 1992), Tn = 2L ncm (1.2) where L is the basin length and n ∈ Z+ is the rank of the longitudinal mode corresponding to the number of nodes in free surface or interface deflection. For two-layer stratification, the wave celerity of the baroclinic mode is given by, cm = √ g ( 1− ρ1 ρ2 ) h1h2 h1+h2 (1.3) where g is the gravitational acceleration; ρ1 and ρ2 are the densities of the upper and lower layers, respectively; and h1 and h2 are the thicknesses of the upper and lower layers, respectively. Equation 1.2 with cm given by equation 1.3 is known as the two-layer Merian formula (Lemmin and Mortimer, 1986), and is used frequently because of its simplicity (e.g., Mortimer, 1953; Stevens and Lawrence, 1997; Lawrence et al., 1997; Laval et al., 2008). More sophisticated MAMs allow for better representation of lake bathymetry. For ex- ample, the two-layer Defant procedure (TDP, Lemmin and Mortimer, 1986) accounts for cross-sectional variability along the thalweg of elongated lakes. Even better, variable depth complete models (VDC) such as the two- and three- layer VDC models (TVDC and ThVDC, Salvadé et al., 1992) provide full representation of lake bathymetry (i.e. depth in the two horizontal dimensions) and account for shallow zones which contain less than the maximum number of specified density layers. However, the TDP and VDC-like models do not include rotation and, hence, are only suitable for narrow lakes. For modal analysis in lakes of large lateral dimensions, MAMs with rotation were de- veloped by, for example, Bäuerle (1994), Lemmin et al. (2005), and Shimizu et al. (2007). Chapter 1. Introduction 6 More recently, a model for dissipative baroclinic modes was added to the range of available MAMs (Shimizu and Imberger, 2008). In contrast to layered MAMs, modal analysis with continuous (non-layered) stratification is currently confined to longitudinal direction (1D) without width variability, rotation, or dissipation (e.g., Münnich, 1996; Fricker and Nepf, 2000; Vidal et al., 2007; Pannard et al., 2011). The periods and structures of baroclinic modes obtained from MAMs can be used in analyzing the baroclinic response in two ways, namely, spectral and spatial decomposition. Of the two ways, it is more common to decompose the observed baroclinic response spec- trally by using Fourier transform (e.g., (Mortimer, 1953; Hutter et al., 1983; Lemmin, 1987; Roget et al., 1997)) or continuous wavelet transform (e.g., Stevens et al., 1996; Stevens, 1999; Antenucci and Imberger, 2003; Naithani et al., 2003; Vidal et al., 2007) to identify periods (or equivalently frequencies) that correspond to high spectral or wavelet power. These periods are matched against periods estimated using MAMs to identify the excited modes. Traditionally, partial confirmation of the excited modes is obtained by comparing the phase of the baroclinic response at few observation stations to the spatial modal structures also obtained from MAMs (e.g., Wiegand and Chamberlain, 1987; Lemmin et al., 2005). More recently, spatial decomposition of the response over the whole domain has been proposed by Shimizu et al. (2007). This spatial decomposition requires that the baroclinic response be known over the entire domain and at different times so that the temporal evolu- tion of baroclinic modes can be inferred. However, because field observations are typically collected at only a few locations, time-stepping hydrodynamic models are needed to simulate the three-dimensional (3D) thermal and velocity fields under wind forcing. These simulated fields are then decomposed using the spatial structures of baroclinic modes obtained from MAMs. The output of MAMs can also be used in dynamic forced models to simulate the baro- clinic response to wind. These dynamic models are referred to here as modal-based forced models (MFMs) (Table 1.1). These MFMs take the baroclinic modes from MAMs, decouple and simulate the individual modal responses, and then superimpose them to obtain the total response. With this method of simulating the response, MFMs provide a natural way to explore the link between excitation and evolution of baroclinic modes and lake and forcing characteristics. The most basic MFM is that of Heaps and Ramsbottom (1966) which simulates the evolution of longitudinal baroclinic modes in narrow flat-bottom basins with two-layer strat- Chapter 1. Introduction 7 Stratification Horizontal Dimensions and Bathymetry Rotation Energy Dissipation Modal Analysis Models Layered Continuous 1D 2D Variable depth Constant depth e.g. Munnich et al. (1992) Variable cross-section e.g. TDP, Lemmin and Mortimer (1986) Rectangular cross-section Neglects rotation e.g. TVDC Salvade et al. (1992) Includes rotation e.g. Lemmin et al. (2005) 1D 2D Variable width Constant width e.g. Fricker and Nepf (2000) Variable depth Constant depth e.g. Wiegand and Chamberlain (1987) Variable depth Constant depth Variable depth Constant depth Includes shallow zones e.g. ThVDC Salvade et al. (1992) Neglects shallow zones e.g. TVD, Hutter et al. (1983) Non-dissipative e.g. Lemmin et al. (2005) Dissipative e.g. Shimizu and Imberger (2008) Notes: - 1D: one-dimensional analysis along the thalweg of elongated lakes. - 2D: two-dimensional analysis in the horizontal. - 2D analysis with constant depth accounts for lake surface geometry. - Dissipative models have only been recently applied (Shimizu and Imberger, 2008). Dissipation is typically small and its effect on modal periods and structures is secondary. - Models using 2D analysis with continuous stratification or 1D analysis with continuous stratification and variable depth and width have not been applied yet. - Shallow zones refer to lake parts with depth shallower than the inerface(s) between the layers. - See List of Acronyms for meanings of TDP, TVD, TVDC, and ThVDC. Figure 1.2: Features of models used in modal analysis (MAMs). Chapter 1. Introduction 8 Table 1.1: Comparison of modal-based forced models (MFMs). Feature Model Heaps and Ramsbottom (1966) Monismith (1985) Monismith (1987) Shimizu and Imberger (2008)a Stratification two-layer multi-layer continuous multi-layer Dimension 1D 1D 1D 2D Bathymetry uniform depth, uniform width uniform depth, uniform width uniform depth, uniform width variable depth Rotation no no no yes Decoupled baroclinic modes first vertical longitudinal modes vertical modes vertical modes horizontal and vertical modes Damping yes no no yes Nature of model Analytical numerical numerical numerical aSame as Shimizu et al. (2007) with the addition of damping. ification. By decoupling vertical modes, this flat-bottom MFM was extended by Monismith (1985, 1987) to multi-layered and continuous stratification. Only recently were MFMs extended to basins of variable bathymetry and layered stratification using decoupled rotating dissipative modes (Shimizu and Imberger, 2008). However, due to the very recent develop- ment of this model, it has not been applied to many lakes. In fact, it seems to have only been applied to Lake Tiberias (Kinneret) (Shimizu and Imberger, 2008, 2010). In this study, periods and structures of baroclinic modes are obtained using the TDP and VDC models. Rotation is neglected because the lakes considered here have small lateral extents. Modes obtained from the TDP are used in an MFM that is developed here and is referred to as the two-layer variable cross-section model (TVC). Since it uses the TDP modes, the TVC is used here to simulate the baroclinic response in uniform-width basins of variable along thalweg depth (Chapter 2) and in variable-width variable-depth basins (Chapter 3 and 4). The temporal dependence in the TVC is derived for three different forcing functions in time: step, periodic, and impulsive wind (Chapters 2 to 4, respectively). For examining the effect of the metalimnion thickness in Chapter 3, the two-layer stratification in the TVC is inadequate and a three-layer MFM with a flat-bottom is used. Finally, a three- dimensional (3D) hydrodynamic numerical model is used to simulate the response in one of Chapter 1. Introduction 9 the examined lakes (Chapter 5). The output of this 3D model is then spatially-decomposed using modes obtained from ThVDC. The main features of the models used in this study are summarized in Table 1.2. 1.3 Objectives and Overview The study examines the modal composition of the baroclinic response to wind in narrow, elongated lakes with one or more basins and multiple arms. The main objective is to link the modal composition of the response to lake characteristics; (a) bathymetry, (b) stratification, and (c) damping, and to forcing characteristics including (d) the spatial wind distribution and (e) temporal wind pattern. To accomplish the stated objective, this study combines field observations, modal-based forced modeling, hydrodynamic numerical modeling, and modal analysis and decomposi- tion. Chapter 2 addresses the baroclinic response to wind in elongated single-basin lakes of non-uniform depth. Analytical solutions are used to examine the response of particular mathematical, curved- and wedge-shaped bottom profiles and compare to the response in flat-bottom basins. Chapters 3 to 5 use field observations from two water bodies located in western Canada, namely, Amisk Lake which is a small (5 km2) two-basin lake, and Knewstubb and Natalkuz Lakes which are large connected lakes (210 km2) comprising the eastern basins of Nechako Reservoir (1135 km2). The baroclinic response of Amisk Lake is analyzed in Chapter 3. Modal-based forced models (MFMs) are used to investigate the excitation of horizontal and vertical baroclinic modes as well as resonance with wind. Chapters 4 and 5 examine the baroclinic response of Knewstubb and Natalkuz Lakes. In Chapter 4, bathymetric idealizations are invoked to iden- tify the dominant horizontal (longitudinal) modes in the two lakes and to illustrate the strong damping that characterizes the response. To obtain a more complete picture of the baroclinic response, a 3D hydrodynamic numerical model is applied in Chapter 5 to Knewstubb and Natalkuz Lakes without idealizations in bathymetry or stratification. Three-layered modal decomposition of the simulated response is used to infer the dominant horizontal and vertical modes. Chapter 1. Introduction 10 Table 1.2: Main features of models used in this study. Model Features and application Modal analysis models (MAMs) TDP • Two-layer Defant procedure (Lemmin and Mortimer, 1986) • Applicable to narrow, elongated lakes with two-layer stratification • Accounts for cross-sectional variability • Assumes linearity, hydrostatic pressure, no mixing, and no rotation • Used in Chapters 3 and 4 TVDC • Two-layer Variable Depth Complete model • Applicable to narrow, elongated lakes with two-layer stratification • Accounts for full bathymetry including shallow zones and side arms (unlike the TDP) • Assumes linearity, hydrostatic pressure, no mixing, and no rotation • Used in Chapters 4 and 5 ThVDC • Same as the TVDC but uses three density layers instead of two. The third layer is added to account for the thickness of the metalimnion. • Used in Chapters 4 and 5 Modal-based forced models (MFMs) TVC • Two-layer Variable Cross-section model (developed in this study) • Uses baroclinic modes obtained from the TDP • Simulates individual responses of baroclinic modes • Analytical temporal dependence to step, periodic, and impulsive wind • Response to general wind patterns obtained using convolution or superposition • Used in Chapters 2 to 4 • In Chapter 2, uniform width is assumed and the TVC is referred to there as the two- layer uniform-width variable-depth model (TUV) Box MFM • Adapted from Heaps and Ramsbottom (1966) and Monismith (1985) • Simulates individual responses of horizontal and vertical baroclinic modes in a flat- bottom basin • Analytical temporal dependence to step wind • Response to general wind patterns obtained using convolution • Used in Chapter 3 3D hydrodynamic numerical model ELCOM • Estuary Lake and Coastal Ocean Model (Hodges et al., 2000) • Solves for the free surface and 3D fields of density and velocity at every time step • Utilizes a structured finite difference grid • Includes non-linearity and mixing • Used in Chapter 5 11 Chapter 2 The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 2.1 Introduction The baroclinic response to wind is an important driver of water motion in stratified lakes and reservoirs (Imboden, 2005) and can have a pronounced impact on the advection and mixing of nutrients and contaminants. Examples include the transport associated with the baroclinic response in Lake Lucerne (Van Senden and Imboden, 1989), in Amisk Lake (Lawrence et al., 1997), and in Lake Geneva (Umlauf and Lemmin, 2005). The baroclinic response can also lead to upwelling of anoxic hypolimnetic water which may cause fish kills (Marti-Cardona et al., 2008) or lead to downwelling of warm surface water which could affect cold water releases from hypolimnetic outlets in reservoirs (Stevens and Lawrence, 1997). The summertime thermal stratification in temperate lakes can often be idealized into two layers (Thorpe, 1974; Csanady, 1975; Hutter et al., 2011). For this idealized stratification, Heaps and Ramsbottom (1966), referred to hereafter as HR, derived an analytical dynamic model for the baroclinic response to wind in a two-layered, narrow, flat-bottom basin. For the same basin shape, Monismith developed a numerical model for multiple layers (Monismith, 1985) and continuous stratification (Monismith, 1987). More recently, Shimizu et al. (2007) and Shimizu and Imberger (2008) formulated a numerical model for the baroclinic response of layered stratification in basins of irregular bathymetry. This numerical model was used to investigate the baroclinic response of Lake Biwa (Shimizu et al., 2007) and Lake Tiberias (Shimizu and Imberger, 2008). The above-mentioned models can be classified as modal-based forced models (MFMs) since they construct the baroclinic response to wind from natural basin modes which are dictated by stratification and bathymetry. In contrast to using MFMs, analysis of the baro- clinic response is often performed by detecting the spectral signatures of baroclinic modes Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 12 in the observed response (Lemmin, 1987; Roget et al., 1997; Lemmin et al., 2005) or in the response simulated using non-modal numerical models (Hodges et al., 2000; Appt et al., 2004; Lemmin et al., 2005). In this paper, we develop and apply an MFM because, unlike spectral analysis, MFMs reveal the coupling between the spatial and temporal characteristics of the wind field and the excited baroclinic modes. Moreover, MFMs are more useful than spectral analysis in identifying the excited modes for strongly damped systems where the response follows the forcing with little free baroclinic oscillation (Chapter 4). The linear, hydrostatic, and no-mixing assumptions are common to MFMs. Linear drag is included in some MFMs to allow for energy dissipation (Heaps and Ramsbottom, 1966; Shimizu and Imberger, 2008) and thus provide more realistic modeling of the baroclinic response (Shimizu et al., 2007). Despite these assumptions, MFMs have been successfully used in analyzing the baroclinic response of several lakes and reservoirs (Heaps and Rams- bottom, 1966; Monismith, 1985; Stevens, 1999; Shimizu and Imberger, 2008). The applica- tion of MFMs is most suitable during episodes of weak wind forcing or strong stratification which are associated with small deflections of the thermocline relative to the thickness of the epilimnion (Horn et al., 2001). During such episodes, internal surges, solitary waves, and shear-induced instabilities are insignificant and the baroclinic response is viscously damped (Horn et al., 2001). Lakes and reservoirs typically have non-uniform bottom profiles. Single-basin lakes may have curved bottoms which are deepest near the center and shoal towards the endwalls. Examples include Wood Lake, Baldegg Lake, and the north basin of Windermere (see bottom profiles in Figure 2.1a-c). For reservoirs impounded in canyons and valleys, a wedge shape is common with a vertical dam at one end and a bottom shoaling towards the other end (Boegman et al., 2005a; Vidal et al., 2007; Vidal and Casamitjana, 2008). An example of such bathymetry is the Mundaring Weir reservoir (Figure 2.1) (Okely and Imberger, 2007; Shintani et al., 2010). The effect of bottom profile on natural baroclinic modes has been investigated in many lakes (Lemmin and Mortimer, 1986; Salvadé et al., 1992; Lemmin et al., 2005) as well as for a few idealized basin shapes (Münnich, 1996). However, the effect of bottom profile on the baroclinic response to wind has been examined for only a few specific sites (e.g., Lake Biwa, Shimizu et al., 2007; Lake Tiberias, Shimizu and Imberger, 2008). The need for a more general understanding has motivated this paper, the objective of which is to characterize the differences in the baroclinic response to wind between flat-bottom, curved-bottom, and Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 13 Figure 2.1: Bottom profiles for (a) Baldegg Lake, (b) Wood Lake, (c) the north basin of Windermere, and (d) Mundaring Weir reservoir. Also shown are (a-c) CB and (d) WB profiles fitted to the maximum lake depth. Lake bottom-profiles were digitized from (a,c) Lemmin and Mortimer (1986), (b) Wiegand and Chamberlain (1987), and (d) Shintani et al. (2010). Interface depths were taken from the same sources. wedge-shaped basins. We start with a description of an analytical family of curved-bottom (CB) and wedge- shaped (WB) basins, then present the derivation of the MFM used here. In the results section, we examine the differences in the natural baroclinic modes and response to hypothetical wind between flat, curved, and wedge basins. In the discussion, we look at changes in the baroclinic modes and response as the bottom deviates from flat. We also discuss the application of the MFM to elongated basins of arbitrary bathymetry. Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 14 ρ2 ρ1 L H1 H2 h2 (x) x η1 η2 q1 q2 τs Figure 2.2: Schematic of the curved-bottom basin (CB) illustrating geometrical dimensions and density stratification. Note that the vertical scale is exaggerated. 2.2 Methods 2.2.1 Family of Curved-bottom Basins The CB is elongated with a length L and a uniform width. The width is smaller than the internal Rossby radius such that, for a first approximation, the effect of rotation of the earth can be neglected. Density stratification in the CB is two-layered with an upper layer of uniform thickness H1 and a lower layer of variable thickness h2 described by: h2 (x) = H1 { m ( 1− ∣∣2xL−1∣∣r) 1−m(1−|2xL−1|r) } (2.1) where x is the along-thalweg distance measured from the center of the basin (−0.5≤ xL−1 ≤ 0.5); r is an exponent which determines the bottom profile; and m = H2 (H1+H2) −1 is the ratio between the maximum thickness of the lower layer, H2, and the total depth at the same location (0 < m< 1) (Figure 2.2). As defined by equation 2.1, CB is symmetric, deepest at the center, and shoals towards the endwalls. For a fixed m, CB approaches a flat-bottom basin as r→ ∞ (Figure 2.3). Setting r = ∞ in equation 2.1 gives a lower layer of uniform thickness, h2 (x) = H2. For a fixed upper layer thickness, CB becomes infinitely deep (H2→ ∞) as m approaches unity (Figure 2.3). The bottom is concave upward everywhere along the thalweg for m≤ 0.5(1− r−1). Above Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 15 −0.5 0 0.50 0.5 1 1.5 2 h 2 H −1 1 xL−1 (a) m= 2/3 −0.5 0 0.5 xL−1 (b) m= 1/2 −0.5 0 0.5 xL−1 (c) m= 1/5 r = 2 r = 4 r = 6 r =∞ Figure 2.3: The effect of the parameters m and r on the bed profile of the curved-bottom basin. Inflection points of the profiles are indicated by dots. The profiles in panel (c) have no inflection points. this limit ( m> 0.5 ( 1− r−1)), the bottom adjacent to the endwalls is concave downward. For m→ 1, the inflection points in the bottom curvature shift from the endwalls towards the center (Figure 2.3). In the range 0 ≤ xL−1 ≤ 0.5 (or alternatively −0.5 ≤ xL−1 ≤ 0), equation 2.1 describes the bottom of wedge-shaped basins (WB) with a vertical dam at xL−1 = 0 and a slope towards the other endwall. As such, the bottom profile of WB has the same characteristics as half a CB with similar r, and tends to a flat-bottom basin, WB∞, as r→ ∞. To examine the effect of bottom profile on the baroclinic response to wind, we derive an analytical solution for the baroclinic responses in curved-bottom and wedge-shaped basins of r = 2, CB2 and WB2, respectively, and compare them to the responses in the flat-bottom basins CB∞ and WB∞ (r = ∞). Comparison of the baroclinic responses includes the ther- mocline deflections, the associated flows in the epilimnion and the hypolimnion, and the velocity difference across the interface (velocity shear). Noteworthy, comparison to same- depth flat-bottom basins (CB∞ and WB∞) was motivated by the desire to bound the responses of CB and WB profiles with r > 2. The effect of the depth of the flat-bottom basin (FB) on the comparison will be discussed later. In the next section, we first derive the temporal evolution of baroclinic modes under wind forcing for arbitrary along-thalweg bottom profiles. Afterwards, we solve analytically for the periods and spatial structures of the baroclinic modes for the particular profiles CB2, WB2, and CB∞ . The temporal evolution of the modes with the modal periods and spatial structures provide a complete model of the response to wind forcing used later in the results section. Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 16 2.2.2 Model Derivation Governing equations Wind stress is applied to the thalweg of a variable-depth basin (Figure 2.2). Neglecting cross- thalweg motions, the linearized hydrostatic equations of motion for two-layer stratification are given by (Heaps and Ramsbottom, 1966): ∂q1 ∂ t =−gH1∂η1∂x + 1 ρ1 (τs− τi− τd1) (2.2) ∂q2 ∂ t =−gh2ρ1ρ2 ∂η1 ∂x −gh2 ( 1− ρ1 ρ2 ) ∂η2 ∂x + 1 ρ2 (τi− τd2) (2.3) ∂q1 ∂x + ∂η1 ∂ t − ∂η2 ∂ t = 0 (2.4) ∂q2 ∂x + ∂η2 ∂ t = 0 (2.5) where η1 (x, t) and η2 (x, t) are respectively the surface and interface deflections (positive upwards); q j (x, t) is the horizontal volumetric flow per unit width in layer j ( j = 1 denotes the upper layer and j = 2 denotes the lower layer); t ≥ 0 is the time coordinate; g is the gravitational acceleration; ρ j is the density of layer j (Figure 2.2). In the momentum equations (equations 2.2 and 2.3), τs (x, t) is the wind stress, τi (x, t) is the interfacial shear stress, and τd j (x, t) is the boundary drag on each layer. Analogous to HR, Bengtsson (1973), Platzman (1978), Bäuerle (1998), and Shimizu and Imberger (2008), we adopt a Guldberg-Mohn parameterization for the drag. With this parameterization, the rate of energy dissipation is a constant fraction of the kinetic energy in the basin (Spigel, 1978; Shimizu and Imberger, 2008). This leads to exponential decay of the baroclinic response resulting from a wind pulse (Spigel, 1978). The boundary drag is parameterized as τd j (x, t) = 2κbρ jq j (2.6) and, similarly, the interfacial shear is parametrized as τi (x, t) = κi (ρ1q1−ρ2q2) (2.7) where κi and κb ≥ 0 are the boundary and interfacial drag coefficients, respectively. Both coefficients have dimensions of inverse time. Noteworthy, setting κi = 0 corresponds to a Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 17 frictionless interface which is typically assumed in layered models (Heaps and Ramsbottom, 1966; Shimizu and Imberger, 2008). However, for generality, we maintain κi 6= 0. In the final mathematical model, κb and κi are collected into a single damping parameter. The governing equations are simplified by applying the rigid lid and Boussinesq approxi- mations (Gill, 1982). The rigid lid approximation, q1 =−q2, filters out the barotropic modes and retains the baroclinic modes (Gill, 1982; Monismith, 1987; Wang et al., 2001). With these approximations and the linear parameterizations for the drag, equation 2.2 and 2.3 reduce to, ∂q2 ∂ t +2(κb+κi)q2 =−g ( 1− ρ1 ρ2 ) H1h2 H1+h2 ∂η2 ∂x − 1 ρ1 h2 H1+h2 τs (2.8) The governing equation for the baroclinic flow in the lower layer is obtained by using equation 2.5 with 2.8 to eliminate η2, ∂ 2q2 ∂ t2 +2(κb+κi) ∂q2 ∂ t = g ( 1− ρ1 ρ2 ) H1h2 H1+h2 ∂ 2q2 ∂x2 − 1 ρ1 h2 H1+h2 ∂τs ∂ t (2.9) Response to wind component To evaluate the solution of equation 2.9 under general wind conditions, we first consider the case where, from a state of rest, steady wind forcing is applied at t = 0 to a basin of arbitrary bottom profile. Let the solution of equation 2.9 be separable, q2 (x, t) =M (t)Q2 (x) (2.10) where M (t) and Q2 (x) are the temporal and spatial dependencies, respectively. Accordingly, equation 2.9 reduces to the ordinary differential equation of a damped harmonic oscillator (De Silva, 2007), d2M dt2 +2ζω dM dt +ω2M = 0 (2.11) where ζ = (κb+κi)ω−1 is the damping ratio and ω is the natural (undamped) angular frequency given by the eigenvalue problem, ω2 =−g ( 1− ρ1 ρ2 ) H1h2 H1+h2 1 Q2 d2Q2 dx2 (2.12) For a basin of arbitrary bottom profile, the numerical or analytical solution of equation 2.12 yields the natural baroclinic modes which are controlled by bottom profile (h2 (x)) and stratification (H1, ρ1, and ρ2). To designate the quantities associated with a single mode, Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 18 we use the subscript n ∈ Z+. For a given mode, the natural frequency, ωn, constitutes the eigenvalue, and the corresponding spatial structure, Q2,n (x), constitutes the eigenvector. Before finding ωn and Q2,n (x) for CB2 and WB2, we examine the temporal evolution Mn (t) of the baroclinic modes. To solve equation 2.11 for Mn (t), two initial conditions are imposed. At t = 0, the lake is forced from a state of rest which necessitates that Mn (0) = 0 and that the interface is level, [∂η2,n/∂x]t=0 = 0. By substituting the second initial condition in equation 2.8, we find that Mn (t) is independent of x when the spatial distribution of the wind forcing is related to that of the baroclinic flow by τs,n (x, t) =−H (t) ( H1+h2 (x) h2 (x) ) AnQ2,n (2.13) whereH (t) is the Heaviside function and An is the constant magnitude of the spatial wind distribution of mode n. The wind component, τs,n, triggers only the baroclinic mode n with the modal structure Q2,n. The layer flow associated with this mode is given by q2,n (x, t) = Mn (t)Q2,n (x) (equation 2.10) where Mn has the well-known form (De Silva, 2007), Mn(t) = An ρ1ωn √ 1−ζ 2n sin (√ 1−ζ 2n ωnt ) exp(−ζnωnt) 0≤ ζn < 1 An ρ1 t exp(−ωnt) ζn = 1 An ρ1ωn √ ζ 2n−1 sinh (√ ζ 2n −1ωnt ) exp(−ζnωnt) ζn > 1 (2.14) For a given mode n, equation 2.14 indicates that the baroclinic flow, q2,n, depends on the value of ζn. In the undamped response (ζn = 0), the baroclinic flow in each layer undergoes indefinite oscillations after application of the steady wind. As the damping ratio increases (0 < ζn < 1), the response to wind is increasingly damped; the flow oscillations diminish in magnitude after fewer cycles asymptoting to a quiescent state until, at ζn = 1, the flow oscillations are totally inhibited. At this point, the response to wind is critically-damped and the layer flows increase to a peak and then decrease to quiescence without subsequent oscillations. With further increases in the damping ratio (ζn> 1), the response becomes over- damped; the rise to peak flow and then the relaxation towards the quiescent state become progressively slower. The interface deflection η2,n is obtained by substituting q2,n =MnQ2,n into equation 2.5 Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 19 and integrating with respect to time to get, η2,n (x, t) = An ρ1ω2n {1−Nn (t)exp(−ζnωnt)}Z2,n (x) (2.15) where Z2,n (x) =−dQ2,n/dx is the spatial structure of the interface deflection for mode n and Nn (t) = cos (√ 1−ζ 2n ωnt ) + ζn√ 1−ζ 2n sin (√ 1−ζ 2n ωnt ) 0≤ ζn < 1 1+ωnt ζn = 1 cosh (√ ζ 2n −1ωnt ) + ζn√ ζ 2n−1 sinh (√ ζ 2n −1ωnt ) ζn > 1 So far, we have derived the analytical solution for the response to a sudden imposition of steady wind of a spatial distribution that is related to Q2,n (equations 2.10, 2.14, and 2.15). This solution is used to construct the response to wind forcing of general spatial and temporal characteristics as illustrated in the next section. Response to general wind forcing Following HR, the response to impulsive wind of arbitrary spatial distribution, τs=H (t)V (x), can be reproduced by superposition of the individual modal responses to wind components of the form given by equation 2.13. The along-thalweg wind distribution V (x) is decomposed in terms of the eigenvectors Q2,n (which form a linearly independent set), V (x) = ( H1+h2 h2 ) ∞ ∑ n=1 AnQ2,n . (2.16) In this decomposition, the modal magnitudes An are given by, An =− ´ L/2 −L/2Q2,nV (x)dx´ L/2 −L/2 h −1 2 (H1+h2)Q 2 2,ndx , (2.17) which was derived by multiplying both sides of equation 2.16 by Q2,p and using the orthog- onality property, ´ L/2 −L/2 h −1 2 (H1+h2)Q2,pQ2,ndx = 0 if p 6= n (Shimizu et al., 2007). The modal magnitudes An are used to evaluate the temporal dependencies of individual modes, Mn (t), and, subsequently, the total lower-layer flow, q2 (x, t) = ΣnMnQ2,n. The interface deflection is obtained in a similar way, by summing equation 2.15 over the excited modes (i.e., modes with An 6= 0). Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 20 The preceding analysis is for a sustained constant wind forcing of arbitrary spatial dis- tribution. Treatment of temporally-varying wind is possible by means of convolution. Con- ditional upon a spatial distribution not changing in time, τs = Fs(t)V (x), the response to time-varying wind is given by the convolution, ∂Fs ∂ t ?q2 = ˆ t 0 ∂Fs(t ′) ∂ t ′ ?q2 ( x, t− t ′)dt ′, (2.18) where q2 (x, t) is the response to the step wind of general spatial distribution V (x) (as sug- gested by HR for a flat-bottom basin). To be able to model the dynamics under wind forcing, the baroclinic modal structures and their associated frequencies (or equivalently periods) have to be evaluated first. As such, we turn our attention back to the eigenvalue problem described by equation 2.12 and solve it analytically for the particular basins CB2, WB2, and CB∞. Baroclinic modes of particular basins The key to arriving at the baroclinic modal structures and periods for CB2 was to recast equation 2.12 into a Legendre differential equation in terms of Z2,n. This was done by multiplying both sides of equation 2.12 by Q2,n, differentiating with respect to x, replacing −dQ2,n/dx with Z2,n, substituting for h2 using equation 2.1 with r = 2, using the transfor- mation x′ = 2xL−1, and re-arranging to get, d dx′ [( 1− (x′)2) dZ2,n dx′ ] +ω2nL 2 ( 4g ( 1− ρ1 ρ2 ) mH1 )−1 Z2,n = 0 . (2.19) This equation is a Legendre differential equation for which the solutions are Legendre poly- nomials given by (Tang, 2007), Z2,n ( x′ ) = x′ n= 1 1 2 ( 3(x′)2−1 ) n= 2 2n−1 n x ′Z2,n−1 (x′)− n−1n Z2,n−2 (x′) n> 2 . (2.20) These solutions for Z2,n are polynomials of degree n, which are valid when the coefficient of Z2,n in equation 2.19 equates to n(n+1) (Tang, 2007). From this equality, the undamped period, Tn = 2piω−1n , of an arbitrary mode n is given by, Tn|CB2 = piL ( n(n+1)g(1− ρ1 ρ2 )mH1 )−1/2 . (2.21) Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 21 Using the solutions for Z2,n, the modal structures for layer flow are obtained from, Q2,n = −´ x−L/2Z2,ndx. These Q2,n structures are polynomials of degree n+ 1 which satisfy the no-through flow boundary condition at xL−1 =±0.5. Equations 2.20 and 2.21 describe the structures and periods of baroclinic modes in the CB2 basin. In contrast, the baroclinic modes in the flat-bottom basin, CB∞, have the sinu- soidal forms Q2,n ∣∣ CB∞ = sin ( npi ( xL−1−0.5)) and Z2,n∣∣CB∞ =−npiL−1 cos(npi (xL−1−0.5)) (Heaps and Ramsbottom, 1966). The corresponding modal periods are given by the two- layer Merian formula (Lemmin and Mortimer, 1986), Tn|CB∞ = 2L ( n2g(1− ρ1 ρ2 )mH1 )−1/2 . (2.22) Noteworthy, by substituting these modal structures, Q2,n ∣∣ CB∞ and Z2,n ∣∣ CB∞ , and periods, Tu,n|CB∞ , in equation 2.10 and 2.15, we obtain the HR model for the baroclinic response to wind in a flat-bottom basin (Heaps and Ramsbottom, 1966, equation 85 and related) - noting that the modal magnitudes of HR equate to m−1An. The baroclinic modes of the wedge-shaped basin WB2 are described by the even modes (n= 2,4,6 . . .) of CB2. Because of the symmetry of CB, the Q2,n structures of the even CB modes are anti-symmetric and thus are characterized by a node in the layer flow at xL−1 = 0. These even modes satisfy the eigenvalue problem (equation 2.12) for WB including the no- through flow at xL−1 = 0 and xL−1 = 0.5. Consequently, the baroclinic modes of WB2 are described by one-half of the even CB2 modes. For example, the fundamental mode of WB2 corresponds to the one-half 0≤ xL−1 ≤ 0.5 of the second mode (n= 2) of CB2. 2.3 Results 2.3.1 Baroclinic Modes Deviation from a flat-bottom affects the modal periods. The ratio of the modal periods for CB2 and CB∞ is given by βTn = 0.5pi √ n(n+1)−1 (equation 2.21 and 2.22). The ratio is smallest for the fundamental mode; βT1 = 1.11, and increases to 1.28 for the second mode. For higher n, βTn asymptotes to 0.5pi = 1.57 as n→ ∞. Deviation from a flat bottom also affects the modal structures. Figure 2.4 compares Q2,n and Z2,n in CB2 and CB∞ for the four lowest modes (n= 1−4). The nodes and anti-nodes of the interface deflection (or equivalently the anti-nodes and nodes of the layer flow) shift Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 22 towards the endwalls in CB2. For the second mode, the nodes of the interface shift from xL−1 = ±0.25 in CB∞ to ±0.29 in the CB2 (Figure 2.4f). For the third mode, the nodes (other than the one at the basin’s center) shift from xL−1 = ±0.33 in CB∞ to ±0.39 in CB2 (Figure 2.4g). Another deviation in the modal structures is that, for a given mode, the interface deflec- tion at the endwalls becomes larger than at the intermediate anti-nodes and, concurrently, the anti-nodes in the layer flow decrease in amplitude from the center towards the endwalls. For the second mode of CB2, the interface deflection at the endwalls is double that at the center of the basin (Figure 2.4f) while for the third mode, the end wall deflections are 2.2 times the deflections at intermediate anti-nodes (Figure 2.4g). This ratio increases for higher modes. In comparison, in a flat-bottom basin, the interface deflection at the endwalls is of the same magnitude as at the intermediate anti-nodes and the anti-nodes in the layer flow all have the same amplitude. As previously noted, the modes of the wedge-shaped basin WB2 correspond to one- half the even modes of CB2 (Figure 2.4b,d,f,h). As such, the features of modal structures are similar in WB2 and CB2. The sloping bottom in WB2 amplifies the modal interface deflection at the shoaling end relative to the deflection at the vertical wall ( xL−1 = 0 ) (Figure 2.4b,d). The peaks in modal layer flow are greater near the dam and decreases towards the opposite endwall. Differences in bottom profiles affect the shapes of along-thalweg wind components through the Q2,n structures and the ratio H1h−12 (x) (equation 2.16). For example, the n = 1 wind component has a uniform distribution in CB2, a linear distribution with no magnitude at the vertical endwall in WB2, and a half-sine distribution in CB∞ and WB∞ (Figure 2.5). The shape of the wind components control which baroclinic modes are excited by a given along- thalweg wind distribution V (x) (i.e. n for which An 6= 0, equation 2.17). As shown in the next section, the shapes of wind components are important to understanding the differences in the excitation of modes by a hypothetical wind between CB2, WB2, and the flat bottom CB∞ and WB∞. Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 23 −1 0 1 (a) Q 2 ,n n = 1 (b) n = 2 (c) n = 3 (d) n = 4 −0.5 0 0.5 −1 0 1 (e) Z 2 ,n xL−1 −0.5 0 0.5 (f ) xL−1 −0.5 0 0.5 (g) xL−1 −0.5 0 0.5 (h) xL−1 Figure 2.4: Spatial structure of (a-d) the lower-layer horizontal flow and (e-h) the interface deflection for the four lowest modes. Solid lines are for CB2 and dashed lines are for the flat- bottom basin CB∞. For display purposes, the distributions Q2,n (x) and Z2,n (x) are scaled to a maximum of unity. −0.5 0 0.5 −1 −0.5 0 0.5 1 xL−1 ( H 1h −1 2 + 1 ) Q 2, n (a) CB2 −0.5 0 0.5 xL−1 (c) CB∞, WB∞ 0 0.25 0.5 xL−1 (b) WB2 n = 1 n = 2 n = 3 Figure 2.5: Shapes of along-thalweg wind components for (a) CB2, (b) WB2, and (c) the flat bottom CB∞ and WB∞. The shapes are scaled to a maximum of unity. Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 24 2.3.2 Response to Impulsive Wind Curved-bottom basin We examined the effect of bottom profile on the forced baroclinic response of CB2 using a hypothetical case. A uniform wind stress, τs, was applied along the thalweg of CB2 and CB∞. The stress was initiated at t = 0 and sustained thereafter. Our choice of a uniform wind distribution is adequate for elongated basins that are confined in valleys (Vidal et al., 2007). The choice was also motivated by the desire to assess the parameterization of the velocity shear in the mixed-layer model of Spigel and Imberger (1980) which is based on uniform wind. For CB2, the shapes of wind components in Figure 2.5a clearly show that the n = 1 component fully represents the applied uniform wind. Accordingly, only the fundamental baroclinic mode (corresponding to n = 1) is triggered and higher modes are not excited as An = 0 for n > 1. For the fundamental mode, the modal magnitude obtained from equation 2.17 is A1 =−mτs. For CB∞, decomposition of the applied uniform wind using equation 2.17 reveals an infi- nite number of sinusoidal wind components which excite the corresponding sinusoidal baro- clinic modes. Because of the symmetry in forcing and basin shape, only odd modes, which have symmetric forms, are triggered by the uniform wind distribution. The magnitudes of odd modes are given by 4(npi)−1mτs (adapted from HR) where the factor n−1 reflects the diminishing modal magnitudes for higher modes. In superimposing the individual modal responses to get the total baroclinic response, contributions from modes higher than n = 99 were truncated. In both CB2 and CB∞, the deflection of the interface and flow in the lower layer were sim- ulated for different modal damping ratios ζn. The damping ratio for the fundamental mode, ζ1, was set to the same value in both basins. In CB∞, the damping ratios for higher modes were set to ζn = (ζ1ω1)ω−1n = ζ1n−1 which maintains mode-independent drag coefficients (i.e., constant κi and κb).1 Under uniform wind, maximum interface deflections from superimposing all excited modes occurred at the endwalls. Figure 2.6a-d shows the evolution of the total interface deflection at the upwind endwall in CB2 and CB∞ for different damping regimes (ζ1 = 0, 0.1, 1.0, and 2.0). The effect of ζn on the interface deflection is similar to that on the layer 1This assumption is justified in Chapter 5. See Figure 5.13 and associated discussion. Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 25 flow which has been described in section 2.2.2. As ζn was increased from ζn = 0 to 1, the interface oscillations decayed more rapidly towards the dynamic equilibrium position. For ζn≥ 1, the oscillations of the interface were completely damped and the interface asymptoted to the equilibrium position. Under the identical forcing and independent of damping, the interface in CB2 and CB∞ oscillated about or asymptoted to the same dynamic equilibrium tilt; ∂η2/∂x=−τs [(ρ2−ρ1)gH1]−1. At the steady state, this interfacial tilt balances the pressure gradient associated with the surface tilt which was caused by the sustained wind stress (Spigel and Imberger, 1980). However, the evolution of the response towards steady state was different between the basins. Primarily, the dynamic response evolved over different time-scales. Under the spatially-uniform wind, these time-scales were dictated by the period of the fundamental baroclinic mode which dominated the response. For overdamped responses (ζ1 > 1), the longer modal period in CB2 (βT1 = 1.11) led to a slower evolution towards steady state than in CB∞. For the underdamped responses (ζ1 < 1), the difference in the modal periods led to the oscillatory responses moving out of step as time progressed. Aside from the differences in the periods, the oscillatory patterns were different. Most notable, for ζn = 0, the interface in CB2 oscillated in a pure sinusoidal pattern while the interface in CB∞ oscillated in a piecewise linear pattern which, for xL−1 =±0.5, reduces to a triangular wave form (Spigel, 1978). Because the dynamic equilibrium tilt of the interface, ∂η2/∂x, is the same in CB2 and CB∞, the water volumes ultimately transported, above and below the interface, from one- half of each basin to the other half has to be equal in both basins. However, because the baroclinic response in the two basins has different T1 periods and temporal patterns, the layer flows differ in their peak magnitudes between the two basins. This is shown in Figure 2.6e-h which displays the evolution of the layer flow at the center of the basins where the maximum flow occurs. The difference in the peak magnitudes is easily explained for the undamped case (ζ1 = 0) in which the interface attains the same tilt in the two basins every quarter of the fundamental period of the respective basin. Thus, over the different time intervals (1 4T1 ) , the volume transported across any section has to be identical in both basins. At xL−1 = 0, the longer T1 period in CB2 combined with the sinusoidal pattern leads to a smaller peak flow than in CB∞. The ratio of the peak flow at x = 0 in CB2 to that in CB∞ is given by β q (ζ1 = 0) = 1/ √ 2 = 0.71. Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 26 Figure 2.6: Comparison of the baroclinic responses to impulsive wind in CB2 and CB∞. Panels (a-d) show the interface deflection at the upwind endwall ( xL−1 =−0.5). Panels (e-h) show the flow in the lower layer at the center of the basin ( xL−1 = 0 ) . Each column of panels corresponds to a damping ratio ζ1 as noted on the top. Interface deflection is normalized by the dynamic equilibrium deflection at xL−1 =−0.5. Layer flow is normalized by the magnitude of the undamped maximum layer flow at xL−1 = 0 in CB∞. Time is normalized by the period of the fundamental mode in CB∞. Note the different vertical scales for the panels corresponding to ζ1 = 1.0 and 2.0. Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 27 0 0.5 1 1.5 2 0.7 0.8 0.9 1 ζ1 β q 1√ 2 Figure 2.7: Ratio between maximum layer flows in CB2 and CB∞ as a function of damping. As damping becomes much stronger (ζ1 1), the interface attains the equilibrium tilt at t→ ∞ (Figure 2.6). The same water volumes previously transported over 14T1 for the case of no damping are now transported over much longer times. Consequently, the peak layer flows in the two basins are greatly reduced. The difference between the peak flows diminishes and the ratio β q approaches one. Hence, β q has a minimum value of 0.71 at ζ1 = 0 and tends to one as ζ1 increases (Figure 2.7). To further highlight the differences in the responses of the two basins, we examined the spatial structure of the response. Figure 2.8 shows snapshots of the interface deflection and the layer flow along the thalweg of the two basins. The snapshots are shown for the case of weak damping, ζ1 = 0.1. The phase difference between the responses of the two basins due to T1 was removed by comparing corresponding phases of the baroclinic response; i.e. similar tT−11 where T1 is for the respective basin. In CB2, because of the bottom curvature, horizontal pressure gradients developed around xL−1 = 0 as soon as wind forcing was applied (Figure 2.8a). The interface tilted with a uniform slope associated with the fundamental mode (Figure 2.4d). In contrast, the interface in CB∞ did not tilt around xL−1 = 0 until pressure gradients built up at endwalls and propagated towards the center of the basin. For example, the interface remained flat between xL−1 = −0.2 and 0.2 till tT−11 = 0.2. In the central part of CB∞, the spatially-uniform wind forced a uniform flow which increased linearly with time up to tT−11 = 0.25 (Figure 2.8b). Because there is no flow through the endwalls, flow near the endwalls was non-uniform. With time, the non- uniformity propagated towards the center, creating a triangular distribution at tT−11 = 0.25. In contrast, the layer flow in CB2 maintained the constant shape of Q2,1 (Figure 2.4a) as the peak flow increased with time. The biggest difference between the flows in the two basins Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 28 Figure 2.8: Snapshots along the thalweg of CB2 and CB∞ showing (a) the interface deflection and (b) the flow in the lower layer. Snapshots are for ζ1 = 0.1. For each row of panels, snapshots in the two basins correspond to the same phase of the baroclinic response tT−11 as noted on the far left. Normalization of the interface deflection and the layer flow follow Figure 2.6. occurred near the endwalls (Figure 2.8b). At xL−1 = ±0.4, the flow in CB2 was 60% less than in CB∞ at tT−11 = 0.05. By tT −1 1 = 0.25, the flow in CB2 was 27% greater. In contrast, for xL−1 =−0.28 and 0.28, the flow in CB2 remained between -29% and +10% of the flow in CB∞. The velocity shear across the interface is ∆U (x, t) = |U1−U2| where U1 (x, t) = q1H−11 and U2 (x, t) = q2h−12 are the depth-averaged velocities in the lower and upper layers, re- spectively. The velocity shear is an important parameter in modeling the deepening of the surface mixed-layer as a result of shear-induced entrainment (Spigel and Imberger, 1980; Spigel et al., 1986; Yeates and Imberger, 2003). In CB∞, U1, U2, and ∆U follow the spatial Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 29 CB2 (d) (e) −0.5 0 0.5 xL−1 1/16 1/8 1/4 tT−11(f) 0 0.5 1 1.5 CB∞ U 1 (a) −1.5 −1 −0.5 0 U 2 (b) −0.5 0 0.5 0 1 2 xL−1 ∆ U (c) Figure 2.9: Along-thalweg distribution of velocity in the (a,d) upper and (b,e) lower layers and (c,f) velocity shear in (a-c) the CB∞ and (d-f) the CB2 for m= 0.5. The distributions are shown for tT−11 = 1/16, 1/8, and 1/4 (as noted on panel f) and for no damping (ζn = 0). The velocities U1 and U2 and the velocity shear ∆U are normalized by the maximum value of U1 in CB∞. distribution of the layer flows because the thicknesses of the layers are constant (Figure 2.9a- c). In CB2, the variable thickness of the lower layer leads to an increase inU2 away from the center (x= 0). In contrast, the velocity in the upper layer U1 has its peak at the center and vanishes towards endwalls (Figure 2.9d,e). Because of the particular form of CB2 (equation 2.1), the difference between U1 and U2 results in a uniform distribution of ∆U (Figure 2.9f). In CB2, the velocity shear evolves as a sinusoidal function in time (Figure 2.10). In comparison, the spatial-averaged ∆U in CB∞ is a quadratic function in time and is smaller in magnitude than ∆U |CB2 for corresponding tT−11 (Figure 2.10). The maximum ∆U in CB∞ occurs at the center of the basin ( xL−1 = 0 ) and evolves as a triangular wave form in time (Figure 2.10). The uniform ∆U in the CB2 is approximately bounded between the maximum and average ∆U of CB∞ (Figure 2.10). Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 30 Figure 2.10: The maximum and the spatial average velocity shear in CB2, WB2, and CB∞ (WB∞) for ζ1 = 0.1. The velocity shear is normalized by its peak undamped magnitude in CB∞, ∆Um = 4ρ1h1 (τsT1)−1 (Spigel, 1978). The normalized velocity shear in WB∞ and CB∞ are identical. Time is normalized by the period of the fundamental mode, T1, of the respective basin. Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 31 Wedge-shaped basin A similar simulation to the above was conducted for WB2 ( 0≤ xL−1 ≤ 0.5). Starting from a state of rest, uniform steady wind was applied to WB2. Unlike in CB2, wind decomposition using equation 2.17 gave An 6= 0 for all n. Because wind components for WB2 taper towards the vertical endwall (e.g., n= 1−3 in Figure 2.5b), an infinite number of wind components are needed to fully represent the uniform wind near the vertical endwall as in the case of the flat-bottom basin CB∞ (and similarly WB∞). Further, because the applied wind distribution is symmetric, both odd and even wind components which are asymmetric resulted from the decomposition of wind over WB2. While each wind component excited the corresponding baroclinic mode, the contribution of higher modes to the baroclinic response diminished as n increased and, overall, the fundamental mode (n= 1) dominated the response. Figure 2.11 shows the evolution of the interface deflection and layer flow at different locations along the thalweg of WB2. With little contributions from modes n > 10, the superposition of deflections and flows of the modes n= 1 to 10 satisfactorily represents the total response (Figure 2.11, grey lines). This can be observed in the decay of flow oscillations along thalweg towards quiescence (Figure 2.11f-j) and the convergence of the interface deflections at xL−1 = 0 and xL−1 = 0.5 towards their steady state values as transient oscillations decay (normalized deflections of -1 and 1, respectively, Figure 2.11a,e). Comparison of the n = 1 response (dashed line) to the total response (n= 1−10) indicates the dominance of the fundamental mode at most locations along the thalweg (Figure 2.11). Relative to n = 1, the second and third modes caused pronounced interface deflections and layer flows towards the shoaling end of WB2 (Figure 2.11d,e,f,j). The simulation reveals an asymmetry in the magnitude of transient interface oscillations between the ends of WB2. At tT−11 = 0.58, the interface at the upwind end xL −1 = 0.5 rose to ∼2.5 (normalized by the steady state deflection) while, at the other end (xL−1 = 0), the interface depressed to only 1.3 at tT−11 = 0.51 (Figure 2.11a,e). This asymmetry in the transient deflections resulted from the asymmetry of the modal Z2,n structures as previously described (Figure 2.4f,h). Further, because the modal frequencies ωu,n for n > 1 are not harmonics (i.e. whole multiples) of ωu,1, the transient oscillations attained irregular shapes. For tT−1 < 1, peaked crests and broad troughs are observed in the interface deflection at xL−1 = 0.5 (Figure 2.11e) Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 32 Figure 2.11: (a-e) Interface deflection and (f-j) lower-layer flow at different locations in WB2 ( 0≤ xL−1 ≤ 0.5). Wind is uniform, steady, and blows in the negative x direction (i.e., the downwind end is xL−1 = 0). The deflections and flow are shown for ζ1 = 0.1. Interface deflection is normalized by the steady state deflection. Layer flow is normalized by the magnitude of the undamped maximum layer flow at xL−1 = 0.25 in WB∞. Time is normalized by the period of the fundamental mode in WB2. while broad crests and narrow troughs are observed at xl−1 = 0.375 (Figure 2.11d). In Baldegg Lake which has a similar bottom profile to CB2 (Figure 2.1a) and in which higher horizontal modes were detected (Lemmin, 1987), asymmetric interface oscillations were observed. As simulated for WB2, the asymmetry of interface oscillations in Baldegg Lake may result from the non-harmonic modal frequencies combined with the relative magnitudes of the modal oscillations. Unlike in CB2, the distribution of velocity shear in WB2 was non-uniform and, unlike in CB∞ or WB∞, the magnitude of ∆U was not greatest at the basin center. For most of the time, ∆U in WB2 was maximum at the shoaling end ( xL−1 = 0.5 ) and was always zero at the vertical endwall ( xL−1 = 0 ) . This followed from the distribution of ∆U for individual modes which have their individual peak values at the shoaling end while, towards the other end, U1, U2, and consequently ∆U tend to zero for all modes (Figure 2.12). The spatial-average ∆U in WB2 evolved in a relatively regular temporal pattern, was Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 33 Figure 2.12: Along-thalweg distribution of velocity shear in WB2 for the lowest three baroclinic modes. Distributions are scaled to a maximum of unity. Higher modes n > 3 also have their peak ∆U at xL−1 = 0.5 (not shown). always greater than the average ∆U in WB∞, and occasionally exceeded the maximum ∆U in WB∞ (e.g., tT−11 < 0.14 and 0.33 < tT −1 1 < 0.7) (Figure 2.10). On the contrary, the maximum ∆U in WB2 had an irregular temporal pattern and persistently and considerably exceeded the maximum ∆U in WB∞. When the layer velocities due to different modes were aligned, ∆U at the shoaling end of WB2 reached high values which were up to several times the average ∆U (e.g., peak at tT−11 = 0.64. Figure 2.10). This relatively localized high ∆U was associated with the modal ∆U distributions being progressively peaked adjacent to the shoaling end for higher modes (Figure 2.12). 2.4 Discussion 2.4.1 Continuum between r = 2 and r = ∞ For CB (WB) profiles with r > 2, the structures and periods of baroclinic modes and the forced response are bounded by those in CB2 (WB2) and CB∞ (WB∞). This bounding is useful as the lower-layer thickness along the thalweg of many lakes can satisfactorily be described by CB profiles with r ≥ 2 (Figure 2.1a-c). Similarly, reservoirs impounded in canyons and valleys can be described by WB profiles with r ≥ 2 (Figure 2.1d). As r→ 2+, the modal structures, Q2,n and Z2,n, smoothly shift from sinusoidal forms in CB∞ (WB∞) to polynomial forms in CB2 (WB2) (Figure 2.4). For the associated modal Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 34 periods, a greater portion of the bottom deviates from the maximum depth, H1 +H2, as r→ 2+ and, as internal waves traverse the length of the basin, the waves slow down over the shallower depths. Consequently, as r→ 2+, the modal periods increase over the values from the two-layer Merian formula (equation 2.22) and asymptote to values from equation 2.21. Under a spatially-uniform wind, the response in the flat-bottom basin is composed of an infinite number of symmetric baroclinic modes. The fundamental mode dominates. Con- tributions of higher modes to interface deflection and layer flow diminish with n−1. In CB2, the response to a spatially-uniform wind is comprised solely of the fundamental mode. Because of partial topography sheltering or a curved thalweg, irregular along-thalweg wind distributions may arise. Since uniform wind triggers the fundamental mode of CB2, irregular wind distributions will trigger higher modes in CB2. The excited modes can be evaluated by decomposing the wind distribution (equation 2.17). Similarly, the modes excited in CB∞ can be evaluated. For intermediate CB profiles (r > 2) subjected to a given wind distribution, we expect a smooth transition in the number and relative magnitude of baroclinic modes that contribute to the response. In transitioning from r = ∞ to r = 2, the velocity shear ∆U shifts towards a uniform distribution. This is accomplished by increased magnitudes of U2, and consequently ∆U , as the bottom shoals towards the ends of the basin. Velocity observations from Lake Alpnach Gloor et al. (1994) seem to qualitatively support our analytical solution, U2 may have ele- vated magnitudes as the bottom shoals towards basin ends even though the baroclinic flow vanishes towards the basin ends. For all WB profiles (except WB∞), because of the asymmetry of the bottom profile, spatially-uniform wind triggers odd and even baroclinic modes which are all asymmetric. For r = ∞ (WB∞), the even modes transform into anti-symmetric modes and are suppressed because of the symmetry of the uniform wind. The fundamental mode which dominates the response of WB2 and WB∞ is expected to dominate the response for WB with r> 2. As r→ 2+, asymmetry gradually increases between the magnitudes of transient interface oscillations at the ends of WB. The magnitudes of oscillations decrease at the vertical endwall and increase at the opposite end. The location of maximum velocity shear ∆U shifts from the center of the basin towards the shoaling end. The elevated velocity shear and interface deflection at the shoaling end may have implications on the entrainment of stream flows which typically enter wedge-shaped reservoirs from the reservoir end opposite to the dam. Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 35 2.4.2 Parameterization of Velocity Shear For CB profiles, the tendency towards a uniform along-thalweg distribution of ∆U as r→ 2+ leads to more horizontal uniformity in entrainment and deepening of the interface. Conse- quently, as r→ 2+, one would expect one-dimensional mixed-layer models to better simulate the horizontally-uniform mixing. In the mixed-layer models of Spigel and Imberger (1980) and Spigel et al. (1986), the parameterization of ∆U is based on its value at the center of CB∞ for t ≤ 0.25 T1|CB∞ , before horizontal pressure gradients are established at that location. The rationale for this parameterization with a time cutoff is that it maintained the equivalence with the velocity shear induced by wind in a horizontally unbounded domain where no horizontal pressure gradients build up. In the unbounded domain, the assumption of horizontal uniformity is valid. However, in CB2, WB2, and generally variable-depth basins with bottom curvature, the horizontal pressure gradients will build at locations along the thalweg as soon as uniform wind forcing is applied. Consequently, the rationale for the parameterization of ∆U by Spigel and Imberger (1980) is not generally inline with the baroclinic response of variable-depth basins. It might be more useful to think of ∆U as a spatial average rather than the velocity shear at a location with no horizontal pressure gradients. Numerically, changing the parameterization of ∆U to reflect the spatial average in CB∞ might be significant. The entrainment induced by shear and the interfacial thickening asso- ciated with billowing depend on the square of ∆U (Spigel and Imberger, 1980; Spigel et al., 1986). However, given the high level of idealization in mixed-layer models, modification to the parameterization of ∆U might not be warranted. In particular, the application of these models is often extended to periods of strong forcing and weak stratification which leads to large interfacial deflections. For such conditions, the linear assumption is violated. Internal surges may form and, in curved-bottom basins, the surges would shoal over the slope, breaking and causing mixing which subsequently leads to entrainment from the lower layer and thickening of the interface (Boegman et al., 2005a). 2.4.3 Relevance to Laboratory Experiments As noted by Shintani et al. (2010), almost all laboratory experiments on entrainment and deepening of the surface layer were conducted in flat-bottom tanks (Wu, 1973; Kranenburg, 1985; Niño et al., 2003). For comparison, it would be useful to conduct similar experiments Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 36 in tanks fitted with bottom profiles of CB2 in which ∆U is theoretically uniform under uniform surface shear. Further, laboratory experiments investigating the degeneration of the baroclinic response in basins with shoaling bottoms can make use of our analytical solution for the baroclinic modes of WB2 and CB2. Boegman et al. (2005a) conducted such experiments in a tilting flume with two-layer stratification using relatively steep slopes at one end of the flume. Boegman et al. (2005a) noted that their flume did not allow experimenting with very mild slopes. However, if the flume allowed for very mild slopes, the horizontal extent of the slopes will shift them to basin-scale features rather than local topographic features. Consequently, the slopes will affect the structures and periods of baroclinic modes and the energy released to each mode as the initially tilted interface relaxes. For the experiments of Boegman et al. (2005a), these parameters were calculated assuming a flat-bottom basin. As an alternative, the experimental flume can be fitted with bottom profiles emulating WB2 or CB2. Our solution for the baroclinic modes can then be utilized. The experimentation with the wedge shaped-basin and the symmetric curved-bottom basin can help isolate the effect of the direction of the interface deflection on interacting with the shoaling bottom (Boegman et al., 2005a). 2.4.4 Depth of Flat-bottom Basin In our comparison of the baroclinic responses in CB2 (or WB2) and the flat-bottom basin, we have chosen to match the maximum depth of the basins, H1+H2. This choice was motivated by the desire to bound the responses of CB and WB profiles with r> 2. This choice of depth also provides a simple comparison between CB2 (WB2) and CB∞ (WB∞) that is independent of m (see Figure 2.6). The effect of the depth of the flat-bottom basin (FB) on its baroclinic response deserves some discussion. In regard to the modal periods, the choice of a depth smaller than H1+H2 for FB generally reduces the difference in Tn between CB (WB) and FB. The depth of FB also influences the modal magnitudes An (equation 2.17). This, however, does not change the contribution of the different modes to the magnitude of the deflections of the interface. For FB, the factor ω−2n An giving the magnitude of the deflections for a mode n (equation 2.15) is independent of the thickness of the lower layer. On the contrary, for the flow in the upper and lower layers, the depth of FB affects the contribution of each mode to the amplitude of the flow as given by the factor ω−1n An (equation 2.14). Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 37 0 0.2 0.4 0.6 0.8 1 0.9 1 1.1 1.2 1.3 1.4 m β T n ,m n = 1 n = 2 n = 3 n = 4 Figure 2.13: The ratio of the periods of baroclinic modes in CB2 and a flat-bottom basin of depth H1+h2. The ratio is shown for the four lowest modes. To illustrate the effect of the depth of FB on the modal periods and the baroclinic flows, consider a flat-bottom basin, FBh2 , having the average depth of CB2, H1+h2. The choice of FBh2 is motivated by its common use in estimating modal periods in lakes (Mortimer, 1953; Stevens and Lawrence, 1997; Lawrence et al., 1997; Laval et al., 2008). The ratio between the periods of baroclinic modes in CB2 and FBh2 depends on m and is given by, βTn,m = β T n √√√√√ 1 m − √ 1 m −1 arctan (√ m 1−m ) . (2.23) As given by this equation, values of βTn,m are bounded between √ 2 3β T n and βTn (Figure 2.13). For example, βT1,m is bound between 0.91 and 1.11 while β T 2,m is bound between 1.05 and 1.28. For m = 0.83, the period of the fundamental mode in CB2 and FBh2 is identical, βT1,m = 1 (Figure 2.13). For any other mode or value of m, there is a difference between the modal periods in CB2 and FBh2 . This difference is however less than between CB2 and CB∞ ( βTn,m ≤ βTn ) . For the layer flow, the ratio between the flow in FBh2 and CB∞ is given by the factor βTn −1βTn,m ≤ 1 which decreases with the decrease of m (equation 2.23). If the damping ratio is low (ζ1 < 0.8), the difference in peak flows between CB2 and FBh2 is less than between CB2 and CB∞. Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 38 2.4.5 General Application Equations 2.12, 2.14, and 2.15 describe the baroclinic response to wind in elongated basins of arbitrary bottom profiles and constitute a two-layer uniform-width variable-depth model (TUV). The TUV is suitable for simulating the baroclinic response to wind when the re- sponse tends to be linear. This occurs when the forcing is weak relative to the stratification; i.e., when the Wedderburn numberW = (ρ2−ρ1)gH21 (τsL)−1 1 (Monismith, 1986; Horn et al., 2001). Being a modal-based forced model (MFM), the TUV separates the effects of the tem- poral pattern and the spatial distribution of the wind forcing from each other and from the effects of stratification and bathymetry. The baroclinic response is synthesized from 1) the structures and periods of baroclinic modes which are determined by basin bathymetry and stratification (equation 2.12), 2) the modal magnitudes An determined by the along-thalweg wind distribution (equation 2.17), and 3) a temporal variation term dictated by a damping parameter and the wind time-series. The TUV is useful in revealing the building blocks of the response and their relative importance which, for the purpose of analysis, is more useful than the direct numerical solution of the governing equation (equation 2.9) (Shimizu et al., 2007). Compared to the recent two-dimensional numerical model of Shimizu and Imberger (2008), the TUV provides for simpler analysis of the baroclinic response in elongated lakes. The simplification is particularly useful in preliminary investigations. For application of the TUV to a given basin, the eigenvalue problem (equation 2.12) is solved numerically to find the modal periods and structures. The along-thalweg wind distribution is decomposed into components which trigger corresponding baroclinic modes (equation 2.17). The interface deflection and layer flow are computed for each mode (equa- tions 2.14 and 2.15), superimposed, and convolved with the time derivative of the wind stress to obtain the total baroclinic response of the basin. By fitting the response to observed deflections or flow, the damping ratios can be calibrated (Shimizu and Imberger, 2008). For basins with non-rectangular cross-sections, which is the more common case, the TUV can still be applied with some approximation. The thickness of the lower layer, h2 (x), is taken to represent the average cross-sectional depth at a distance x along the thalweg. The thickness h2 (x) is computed from the cross-sectional area below the interface S (x), and the average width at the level of the interface B, h2 (x) = B−1S (x). Chapter 2. The Baroclinic Response to Wind in Elongated Basins of Non-uniform Depth 39 2.5 Conclusions We have derived a modal-based forced model (MFM) for the baroclinic response to wind in narrow uniform-width, variable-depth basins with two-layer stratification. The model provides a simple tool for simulating wind-forced thermocline deflections and baroclinic flow without the complexity introduced by the lateral variability. For a particular curved-bottom basin, CB2, the analytical solution for the baroclinic modes was derived. The even baroclinic modes (n= 2,4,6 . . .) of CB2 fully describe the modes of a wedge-shaped basin (WB2) representing a half of CB2. Combined with the modal-based forced model, this solution for the modes provides a complete analytical de- scription of the baroclinic response to wind in CB2 and WB2. Comparison of the baroclinic responses of CB2, WB2, and a flat-bottom basin FB re- vealed the effect of depth variability. The important similarity in the response of these basins to spatially-uniform wind is that the response is controlled by the period of the fundamental mode, T1. For the same basin length and stratification, the difference in T1 between the basins depends on the relative depth of the basins. For similar depths, T1 in WB2 is greater by 15% and 28% than T1 in CB2 and FB, respectively. The main difference in the baroclinic responses of CB2, WB2, and FB is that pressure gradients form in CB2 and WB2 as soon as uniform wind is applied while in FB pressure gradients propagate from the endwalls to the center. Further, the baroclinic response of WB2 revealed that asymmetry in the basin shape alters the symmetry of the transient interface deflections and layer flows between the ends of the basin. The velocity shear in CB2 was found to be uniformly distributed along the thalweg. This suggests that shear-induced mixing occurs more uniformly along the thalweg of variable- depth basins where the bed shoals towards the endwalls. The horizontal uniformity of mixing supports application of one-dimensional vertical mixing models. In contrast, velocity shear in WB2 had a non-uniform distribution with high localized shear at the shoaling end of the basin. This has implications on the suitability of vertical mixing models for wedge-shaped basins and on the entrainment of stream inflows. 40 Chapter 3 The Baroclinic Response to Wind in a Small Two-basin Lake 3.1 Introduction Baroclinic motions forced by surface wind are an important dynamical feature driving trans- port in stratified lakes (Ostrovsky et al., 1996; Roget et al., 1997). When stream flow is small relative to lake volume, in-lake water quality may primarily be affected by the baroclinic response to wind (Fischer et al., 1979; Umlauf and Lemmin, 2005). The characterization of the baroclinic response in terms of the excitation and amplification of horizontal and vertical baroclinic modes is important for understanding the associated transport and its implication on water quality. The baroclinic response to wind is typically dominated by the first vertical and hori- zontal baroclinic mode, V1H1 (Mortimer, 1953). Higher horizontal modes, however, are also commonly observed (e.g., Mortimer, 1953; Hutter et al., 1983; Lemmin, 1987; Roget et al., 1997). The excitation of higher horizontal modes has been attributed to complicated bathymetry and spatial wind distribution (Mortimer, 1952, 1953), wind direction relative to surface geometry (Roget et al., 1997; Shimizu et al., 2007), or transfer of modal energy from the degeneration of the first mode (Lemmin, 1987). Some excited modes may be amplified due to resonance with wind forcing (Antenucci et al., 2000; Antenucci and Imberger, 2003; Vidal and Casamitjana, 2008). The difference between the forcing frequency and the frequency of an excited mode determines the strength of resonance (Thorpe, 1974; Wiegand and Carmack, 1986). Further, tuning into resonance with wind forcing is affected by modal damping (Shimizu and Imberger, 2008), strength of wind forcing (Boegman and Ivey, 2007; Wake et al., 2007), and structure of stratification (Wake et al., 2007). Higher vertical modes are typically excited when the metalimnion is broad or the modes Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 41 come into resonance with wind (Wiegand and Chamberlain, 1987; Münnich et al., 1992; Roget et al., 1997; Vidal et al., 2007; Pannard et al., 2011). In particular, the second-vertical first-horizontal mode, V2H1, has frequently been observed to compress and expand the met- alimnion thickness at basin endwalls (e.g., Saggio and Imberger, 1998; Okely and Imberger, 2007). While this feature of the V2H1 mode has been extensively studied, the effect of V2H1 on the vertical flow structure of inter-basin exchange has rarely been investigated in multi-basin lakes (Boehrer et al., 2000). Our objective in this paper is to characterize the baroclinic response to wind in a small two-basin lake. In particular, we examine the excitation of horizontal and vertical modes and the effect of damping on the resonance of the V1H1 mode with wind forcing. To simulate the excitation of horizontal modes, we use a dynamic two-layer model which incorporates the cross-sectional variability of the lake. The model also incorporates forcing by periodic wind components which are used to decompose the lake response over the observed frequencies in the wind. We also report higher vertical modes that are inferred from the vertical structure of the exchange between the two basins of the lake. For a given vertical structure, the exchange may be associated with multiple vertical modes. To examine the evolution and interaction of these modes, we use a multi-layered box model which is adapted from Heaps and Ramsbottom (1966) and Monismith (1985). We start with a description of the field site and data. The two models mentioned above are then described. After that, we present field observations and results from the models. Finally, we discuss the link between the excitation and modulation of baroclinic modes and lake bathymetry, stratification, and the spatial and temporal wind characteristics. 3.2 Methods 3.2.1 Field Site and Data Amisk Lake, Alberta, Canada, is 8.5 km long, with a mean surface width of 0.6 km (Prepas, 1990) (Figure 3.1). The lake is divided into two basins by a narrow and elevated channel. While the south basin has a slightly greater surface area than the north basin, it is almost twice as deep (Figure 3.1). Stream inflow is low with bulk residence time of ∼ 8 years (Prepas, 1990) and will be neglected. The baroclinic response of Amisk Lake to wind forcing plays an important role in in- Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 42 terbasin oxygen transport (Lawrence et al., 1997; Rueda et al., 2010). Between 1988 and 1993, oxygen was artificially injected in the hypolimnion of the north basin. However, hypolimnetic oxygen in the south basin also increased as a result of oxygen transport from the north basin (Lawrence et al., 1997). Lawrence et al. (1997) and Rueda et al. (2010) estimated that 30-80% of the increase in oxygen was due to the baroclinic response to wind which establishes the importance of the baroclinic transport in Amisk Lake and motivates this investigation of the detailed baroclinic response. To simulate the baroclinic response of Amisk Lake to wind forcing, we use field data from 1991 (see Lawrence et al., 1997, for detailed description). The data are from a wind station, a thermistor mooring, and an ADCP, marked on Figure 3.1. For our purposes, we focus on the period when the ADCP data was available, between day 179 and 211 (28 June to 30 July, 1991). During this period, the internal Rossby radius is ∼ 1.9 km, about twice the maximum width of the lake (∼ 1 km) and, hence, the effect of rotation on the baroclinic response will be neglected. In 1991, temperature profiles were collected along the thalweg of Amisk Lake on 30-June (day 181) and on 28-July (day 209). Between these surveys, the thermocline deepened but remained above the 12 m crest of the elevated channel, maintaining the baroclinic connection between the two basins (Figure 3.2). Wind forcing per unit length of the thalweg of the lake, Fs, was estimated from the recorded wind speed and direction as follows. The wind speed was adjusted from the sensor height of 4 m to 10 m by a factor of 1.09 (Charnock, 1955) and transformed to wind stress by a standard bulk formulation (Fischer et al., 1979). This wind stress was projected along the average thalweg orientation and multiplied by the cross-sectional width at the surface, b1, to get, Fs (x, t) = τs (t)b1 (x), where τs is the along-thalweg wind stress (τs > 0 if directed from south to north), t is time, and 0 ≤ x ≤ 8.13km is the along-thalweg distance at the depth of the thermocline (5 m) with the origin at the south end of the lake (Figure 3.1). Though the wind stress τs is assumed to be spatially uniform, the resulting along-thalweg distribution of the wind forcing is irregular, mirroring the variability in the surface width of the lake (see Fs in Figure 3.3a). Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 43 0 1 20.5 km Amisk Lake ADCP Wind station and thermistor mooring Albert a 112°39’W 5 4 °3 7 ’N 5 4 °3 4 ’N x = 0 Figure 3.1: Map of Amisk Lake. Isobaths are shown every 10 m from a water surface of ∼ 612 m above sea level. The 20 m and 40 m isobaths are shown in solid black. The dashed line shows the thalweg of Amisk Lake. Along-thalweg distance, x, is positive to the north with the origin, x = 0, at the south end where the thermocline intersects the bottom. Lower right inset shows the location of Amisk Lake in the province of Alberta, Canada. Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 44 4 8 12 16 200 5 10 15 20 Temperature (◦C) D ep th (m ) Channel crest 30-June 28-July Average Two-layer Figure 3.2: Temperature profiles collected on 30-June and on 28-July, 1991. Also shown is the average of the profiles and its two-layer idealization. Note that the thermocline was shallower than the channel crest. 3.2.2 Two-layer Variable Cross-section Model (TVC) In this section, a model is developed to simulate the baroclinic response to wind in Amisk Lake by superimposing the decoupled responses of horizontal baroclinic modes. The spatial structures and associated periods of the modes are obtained for a two-layer stratification from the two-layer Defant procedure (TDP) (Mortimer, 1979) which is applicable along the thalweg of elongated lakes and accounts for cross-sectional variability (e.g., Imboden et al., 1983; Lemmin and Mortimer, 1986; Roget et al., 1997). The TDP neglects mixing, non- linear effects, cross-thalweg motions, and deviation from hydrostatic pressure distribution and uses the rigid lid and the Boussinesq assumptions. We apply the TDP by solving the eigenvalue problem1, ω2u,nQn =−g ( 1− ρ1 ρ2 ) S1S2 S1+S2 d dx ( 1 b2 dQn dx ) (3.1) where ωu,n is the natural (undamped) frequency and Qn (x) is the spatial structure of the horizontal volumetric flow in the lower layer for the baroclinic mode n ∈ Z+; g is the gravitational acceleration; ρ j is the density of layer j (where j = 1 is the upper layer and j = 2 is the lower layer); S j (x) is the cross-sectional area of layer j; and b2 (x) is the basin width at the equilibrium level of the density interface. The baroclinic modes obtained from 1Equation 3.1 is not given by Lemmin and Mortimer (1986) but derived independently. Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 45 equation 3.1 constitute the basis that will be used to simulate the baroclinic response to wind forcing. To derive the equation governing the baroclinic response, we laterally integrate the mo- mentum and continuity equations for the upper and lower layers. Then, by adopting the assumptions of the TDP and by applying linear parameterizations for drag forces (see Chap- ter 2, for a similar treatment), we arrive at the governing equation for the baroclinic flow, ∂ 2q2 ∂ t2 +2κn ∂q2 ∂ t −g ( 1− ρ1 ρ2 ) S1S2 S1+S2 ∂ ∂x ( 1 b2 ∂q2 ∂x ) =− 1 ρ1 S2 S1+S2 ∂Fs ∂ t (3.2) where q2 (x, t) is the horizontal volumetric flow in the lower layer which is equal and opposite to that in the upper layer; and κn is the mode-dependent drag coefficient. We seek solutions of equation 3.2 that have separable spatial and temporal dependencies for q2. This requires restricting the along-thalweg distribution of wind forcing to be time- independent, Fs (x, t) =V (x)τs (t) , (3.3) where V (x) is the along-thalweg wind distribution and τs (t) is the temporal pattern. With this restriction, the spatial wind distribution,V (x), is decomposed using the modal structures, Qn (x), as basis functions. Similarly, the temporal pattern of the wind forcing, τs (t), is decomposed using the complex exponential as basis (i.e. using Fourier transform). This spatial and temporal decomposition of the wind forcing give components of the form, F(n,m)s =− ( S1+S2 S2 ) (AnQn) ( Bm exp ( iω f ,mt )) (3.4) where An ∈R is the magnitude of the spatial mode n, Bm ∈C is the magnitude of the temporal mode m, ω f ,m is the angular frequency of mode m, and i is the imaginary unit. In equation 3.4, the factor (S1+S2)S−12 was necessary for the separation of the spatial and temporal dependencies of the baroclinic response as follows. For a single wind component F(n,m)s , the solution for the layer flow is q (n,m) 2 =Mn,m (t)Qn (x) where Mn,m is the temporal dependency. Subject to the initial conditions of no motion and a flat density interface, the form of Mn,m (t), corresponding to a spatial mode, n, and a forcing mode, m, is given by the real part of Mn,m = Dn,m exp ( iω f ,mt )−Dn,m exp(−ζnωu,nt)Nq (t) (3.5) Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 46 where ζn = κnω−1u,n is the damping ratio for mode n, Dn,m = AnBm ρ1 ω f ,m( ω2f ,m−ω2u,n ) i+2ζnωu,nω f ,m , and (3.6) Nq (t) = cos (√ 1−ζ 2n ωu,nt ) + 1√ 1−ζ 2n ( ωu,n ω f ,m i−ζn ) sin (√ 1−ζ 2n ωu,nt ) 0≤ ζn < 1 1+ ( ωu,n ω f ,m i−1 ) ωu,nt ζn = 1 cosh (√ ζ 2n −1ωu,nt ) + 1√ ζ 2n−1 ( ωu,n ω f ,m i−ζn ) sinh (√ ζ 2n −1ωu,nt ) ζn > 1. (3.7) The first term on the right hand side of equation 3.5 represents the final stage of the response of mode n in which the direction of flow in each layer oscillates at the frequency of the forcing, ω f ,m, with an amplitude2, |Dn,m|, that decreases with increasing the damping ratio ζn. The second term on the right hand side of equation 3.5 describes transients that eventually decay leaving only the steady-state oscillations described by the first term. The time-scale of the transients is dictated by the undamped frequency ωu,n while the damping ratio ζn determines the nature of the transients which can be underdamped (0≤ ζn < 1), critically-damped (ζn = 1), or over-damped (ζn > 1) (Chapter 2). Noteworthy, for ω f ,m = 0; i.e. a steady wind applied at t = 0, the steady-state term in equation 3.5 vanishes leaving only the second term which describes the flow transients triggered by the onset of the wind (Chapter 2). Using the solution for the layer flow q(n,m)2 = Mn,mQn, the interface deflection, η , is derived analytically from conservation of mass in the lower layer, η(n,m) (x, t) = Dn,m ω f ,m { exp ( i ( ω f ,mt− pi2 )) + iexp(−ζnωu,nt)Nη (t) } Zn (x) (3.8) where Zn (x) =−b−12 dQndx is the spatial structure of the interface deflection for mode n, and Nη (t) = cos (√ 1−ζ 2n ωu,nt ) + 1√ 1−ζ 2n ( ω f ,m ωu,n i+ζn ) sin (√ 1−ζ 2n ωu,nt ) 0≤ ζn < 1 1+ ( ωu,n+ iω f ,m ) t ζn = 1 cosh (√ ζ 2n −1ωu,nt ) + 1√ ζ 2n−1 ( ω f ,m ωu,n i+ζn ) sinh (√ ζ 2n −1ωu,nt ) ζn > 1 (3.9) 2|Dn,m| is the complex magnitude of Dn,m, given by |Dn,m| = √ [ℜ(Dn,m)]2+[ℑ(Dn,m)]2 where ℜ(Dn,m) and ℑ(Dn,m) are the real and imaginary parts of Dn,m, respectively. Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 47 Table 3.1: Parameters of the two-layer idealization. Date h1 Θ1 Θ2 ρ2−ρ1 (m) (°C) (°C) (kg m-3) 30-June 3.6 19.8 8.3 1.579 28-July 5.6 19.9 7.7 1.638 30-June to 28-July 4.5 19.9 8.1 1.618 Equations 3.5 and 3.8 represent the response to a single combination of the temporal components (m-modes) and the spatial components (n-modes) of the wind forcing, F(n,m)s (equation 3.4). The total layer flow driven by Fs is given by the superposition of responses to individual combinations of n and m; q2 = Σn,mq (n,m) 2 which can be further simplified by collecting the temporal contributions for a given baroclinic mode; i.e. q2 = ΣnQn (ΣmMn,m). The total interface deflection is treated similarly. Hereafter, the above equations are referred to as the two-layer variable cross-section model (TVC). As a first step to finding the response of Amisk Lake to wind using the TVC, we evaluated the horizontal baroclinic modes of the lake. Equation 3.1 was discretized using a conservative finite difference scheme. The discretized equation was solved using the Matlab eigenvalue solver, EIGS, to determine the spatial structure of the layer flow, Qn (x); the corresponding period, Tu,n = 2piω−1u,n ; and the spatial structure of the interface deflection, Zn (x). To obtain the stratification parameters for the TVC, the average of the temperature pro- files of 30-June and 28-July was idealized into two layers (Figure 3.2). The density interface was placed at the depth of the maximum temperature gradient, h1 = 4.5 m. The layer temperatures, Θ j, were calculated from the volumetric mean temperatures above and below the interface, respectively. Densities, ρ j, were derived from temperatures using Chen and Millero (1986) and assuming a constant salinity of 0.15 (PSU) based on observed specific conductivity of 300 µS/cm (Prepas, 1990). Table 3.2 lists the parameters of the two-layer idealization for the average stratification of the entire period and for the individual average stratifications of 30-June and 28-July. To obtain the bathymetric parameters for TVC, the variation of the width with depth was extracted from 95 cross-sections (Figure 3.3a). The widths were vertically integrated in each layer to obtain S1 and S2. The bathymetric parameters, S1, S2, and b2, were interpolated at Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 48 South Basin North Basin Wind station, Thermistor Mooring ADCP Width (m) 0 500 1000 Fs (x, t) (a) 0 10 20 A re a (× 1 0 3 m 2 ) S1 S2 0 2 4 6 8 0 200 400 600 800 Distance, x (km) W id th (m )b2 (b) Figure 3.3: (a) Variation of the width of Amisk Lake in the vertical and along the thalweg. The vertically-hatched region on top represents the distribution of wind forcing along the thalweg. (b) Along-thalweg distribution of the cross-sectional areas (S1 and S2) and the interfacial width (b2). regular 25 m intervals for input to the finite difference scheme (Figure 3.3b). Damping ratios for the baroclinic modes can be estimated from the relaxation of each mode under calm wind conditions (Heaps and Ramsbottom, 1966; Arneborg and Liljebladh, 2001). For the n = 1 mode, the damping could not be characterized using this method because calm episodes during the record (day 179-211) were not long enough compared to the modal period of n = 1, Tu,1 =21.4 h. Instead, we calibrated the simulated flow at the ADCP site against the observed flow to infer ζ1, similar to the way modal damping in Lake Tiberias was inferred from observed thermocline deflections (Gomez-Giraldo et al., 2006; Shimizu and Imberger, 2008). For modes n> 1, the modal periods were short (. 9 h) relative to the calm episodes. Using the logarithmic decrement method (Mortimer, 1952; De Silva, 2007), ζn for n> 1 was estimated to be 0.05 from the decay of the ∼ 7 h thermocline oscillations during day 191 when calm wind conditions prevailed (see Figure 3.4a,b). Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 49 3.2.3 Multi-layer Uniform-depth Model (MLM) Vertical baroclinic modes in Amisk Lake were investigated using a multi-layered model (referred to hereafter as MLM) which is adapted from the models of Heaps and Ramsbottom (1966) and Monismith (1985). In a basin of uniform depth and width containing N layers of uniform thicknesses h j ( j = 1,2 . . .N), the vertical wave modes are determined from the eigenvalue problem (Monismith, 1985; Münnich et al., 1992), β 2r R (r) = αR(r) (3.10) where βr is the celerity of the vertical mode r; α is an N×N square matrix whose elements are given by, α jk = gh j ρk ρ j k < j gh j k ≥ j ; (3.11) and R(r) is a vector of dimension N representing the amplitudes of the layer flows associated with the mode r. Vertical baroclinic modes3 represent r = 1 to N− 1 and, combined with horizontal modes, will be referred to using the notation of Münnich et al. (1992). The VrHn mode represents the rth vertical mode combined with the nth horizontal mode. The flow in layer j associated with the VrHn mode is given by, q(n,r)j =Mn,r (t)R (r) j sin (npi L x ) , j = 1,2 . . .N (3.12) where Mn,r (t) is the temporal dependence of the flow which, in addition to βr and a damping parameter, depends on the temporal pattern of the wind forcing and the initial conditions. For a steady wind imposed impulsively, F(n)s =An sin ( npixL−1 ) , and for at-rest initial conditions, the temporal dependence Mn,r is given by, Mn,r (t) = Cn,r sin (√ 1−ζ 2n,rωn,rt ) exp(−ζnωn,rt) 0≤ ζn < 1 Cn,rωnpt exp(−ωn,rt) ζn = 1 Cn,r sinh (√ ζ 2n,r−1ωn,rt ) exp(−ζnωn,rt) ζn > 1 (3.13) where ωn,r = npiβrL−1 and ζn,r = κnω−1n,r are the undamped frequency and the damping ratio associated with the VrHn mode, respectively, and Cn,r is a coefficient obtained by solving ΣN−1r=0 R (r) j βrφn,rCn,r = δ1 j AnL npiρ1 j = 1,2 . . .N 3r = 0 represents the barotropic mode. Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 50 in which δ is the Kronecker delta and φn,r = √ 1−ζ 2n,r 0≤ ζn < 1 1 ζn = 1√ ζ 2n,r−1 ζn > 1 . (3.14) For a general wind distribution Fs = ΣnAn sin ( npixL−1 ) , the total flow in each layer is obtained by superposition of modal contributions, q j = Σn,rq (n,r) j , as done for the TVC (Heaps and Ramsbottom, 1966). For general temporal wind patterns, the resultant response is given by the convolution of the response to steady wind with the temporal derivative of the wind forcing (Chapter 2). To apply the MLM, Amisk Lake was idealized as a box basin (BB) which has a flat bottom and a uniform width. Compared to Amisk Lake, BB has the same thalweg length along the interface, L = 8.13 km, and the same stratification parameters (ρ j for j ≤ N and h j for j < N). The thickness of the bottom layer ( j = N) was taken as the average depth of Amisk Lake below the Nth interface. The width of BB was taken as the average width (460 m) at the static equilibrium level of the thermocline. The BB width is irrelevant to the periods or structures of the baroclinic modes, but is used to compute the total flows. 3.3 Results 3.3.1 Field Observations Between 28-June and 30-July (day 179-211), the wind direction over Amisk Lake was variable, blowing at times from the south along the thalweg of the lake and at other times across the thalweg, particularly from the west (Figure 3.4a). Wind speeds were mild to moderate with magnitudes smaller than 5.5 m s-1 for 90% of the time (Figure 3.4a). The winds over the rest of the summer of 1991 and over the summer of 1990 were similar (not shown) which suggests that the observed wind over day 179-211 in 1991 is representative of the typical summertime forcing over Amisk Lake. Water temperature recorded at the mooring (x = 3.5 km) fluctuated with magnitudes dependent on the depth (Figure 3.4b). The largest fluctuations were exhibited at 5 m depth with a range of 8.5°C. These large fluctuations relative to the shallower and deeper thermis- tors confirm the placement of the density interface at 4.5 m as representative of the average Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 51 thermocline. The thermocline was shallower at the start (day 179) and steadily deepened thereafter (Table 3.1 and Figure 3.4b). The ADCP, at the north end of the channel (x= 5.1 km), shows a variable velocity struc- ture for the along-thalweg currents (Figure 3.4c). Two-layer exchange dominated the record, occurring during 52% of the time, while three-layer exchange occurred during 25% and four-layer exchange occurred during 17% of the record. Besides occurring more frequently, two-layer exchange was associated with higher velocities and flows. Figure 3.5 shows the histogram of total volumetric transport through the cross-section at the ADCP site. For two- layer exchange, the average transport was 50 m3s-1 with a standard deviation of 27 m3s-1. For higher vertical structures, the average transport was 25 m3s-1 with a standard deviation of 11 m3/s. It should be noted that the total volumetric transport at the ADCP due to the operation of the oxygen diffuser in the North basin was estimated to be 2 m3s-1 (Lawrence et al., 1997) which is negligible compared to the observed transport (Figure 3.5). Only two-layer exchange can be resolved by the TVC. The occurrence of three-layered exchange will be discussed at a later stage using the MLM. For the two-layer exchange, the depth of flow reversal occurred about an average depth of 5 m with a standard deviation of 1.3 m. This is in good agreement with the at-rest depth of the density interface and justifies the comparison of the observed two-layer flows with the flows simulated by the TVC. 3.3.2 Spectral Analysis The time-scales characterizing water currents and temperatures are related to the temporal scales of wind forcing and to the periods of the baroclinic modes (Wiegand and Carmack, 1986; Antenucci and Imberger, 2003). To examine this relationship in Amisk Lake, the ob- served wind, water temperature, and currents were analyzed using the conventional Fourier transform (Bendat and Piersol, 2010). For the along-thalweg component of the wind stress, spectral analysis shows a significant diurnal periodicity (∼ 1 d) (Figure 3.6a). This ∼ 1 d period was also common to the power spectra of the 5 m and 8 m temperature records at the mooring (Figure 3.6b) and to the currents in the surface and bottom bins at the ADCP (Figure 3.6c). Another prominent spectral peak at a period of 7.3 h appears in the power spectra of the 5 m and 8 m temperature records and the currents in the bottom ADCP bin (Figure 3.6b,c), while the power spectra of the wind and the currents in the top ADCP bin do not show this 7.3 h period (Figure 3.6b,c). Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 52 0 5 10 W in d S p ee d (m s− 1 ) (a) E S W N E W in d D ir ec ti o n 10 15 20 T em p er a tu re (◦ C ) 3m 5m 8m 10m 12m (b) Julian Day in 1991 D ep th (m ) (c) 180 185 190 195 200 205 210 0 5 10 Velocity (m s−1) −0.13 0 0.13 Figure 3.4: Observations at Amisk Lake from 28 June to 30 July, 1991 (days 179-211). (a) Wind speed and direction. Wind speed is smoothed with a 1 h moving average. Wind direction is displayed for speeds greater than 2 m s-1. Dotted lines indicate the average orientation of the thalweg. (b) Water temperature at the mooring. Depths of thermistors are indicated on the records. (c) Along-thalweg water currents recorded by the ADCP. Positive velocities are from south to north. The black contour lines indicate the depth of flow reversal. Instances of two-layer flow are indicated by (+). Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 53 0 20 40 60 80 100 120 140 0 5 10 15 Flow (m3s−1) O cc u rr en ce (% ) Figure 3.5: Histogram of total volumetric transport through the cross-section at the ADCP site. Grey bars indicate two-layer exchange (52%) and white bars stacked on top indicate exchange with three or more layers (48%). Spectral coherence between the along-thalweg wind and the measurements at the moor- ing shows that, at the ∼ 1 d period, the wind and the fluctuations of water temperature at 5 m and 8 m were highly coherent (Figure 3.6d). The wind and the currents at the ADCP were also highly coherent at ∼ 1 d (Figure 3.6e). At the 7.3 h period, coherence between the wind and the water temperatures and velocities were low but statistically significant (Figure 3.6d,e). To identify the excited baroclinic modes from the power spectra of water temperature and currents, we estimated the modal periods of Amisk Lake using equation 3.1 based on the average stratification of the June and July profiles (Table 3.1 and Figure 3.2). The period of the H1 mode, Tu,1 = 21.4 h, was close to the prominent∼ 1 d spectral peak while the periods of the H2 (8.3 h) and H3 (6.4 h) modes were within 14% of the 7.3 h spectral peak (Figure 3.6). Noteworthy, for the June or July stratification in Table 3.1, the estimated periods for the H1, H2, and H3 modes varied by less than 8% from the periods corresponding to the average stratification. In summary, there is a strong periodicity in the forcing and the response at a period of ∼ 1 d with high spectral coherence. It is likely that the H1 mode which has a comparable period is following the ∼ 1 d periodicity in the forcing (Lawrence et al., 1997) as has been observed in Lake Biwa (Saggio and Imberger, 1998) and Lake Tiberias (Lemckert et al., 2004). Besides the ∼ 1 d signal, a less energetic 7.3 h signal is manifested in the response, Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 54 in particular, the temperature fluctuations. The near absence of the 7.3 h periodicity from the forcing suggests the excitation of an underdamped baroclinic mode. The excited mode is either H2 or H3, which will be investigated using the TVC. For later comparison with the TVC results, we examined the phase relationships between the observed wind and the observed temperature fluctuations and water currents. At the∼ 1 d period, the current in the bottom bin of the ADCP was ' 180◦ out of phase with the wind (Figure 3.6e) while the current in the surface bin was in-phase with the wind (not shown). At the thermistor mooring, the ∼ 1 d temperature fluctuations at the 5 m thermistor were out of phase with the wind by about -90° (Figure 3.6d) which means that the underlying thermocline deflections lagged the wind by ∼6 h. The negative phase difference indicates that wind blowing from the south towards the north (in the positive x direction) gradually raises the thermocline at the mooring (positive deflection) which registers as a drop in temperature at the 5 m thermistor. At the 7.3 h period, the phase between the wind and the temperature fluctuations at the 5 m thermistor shifted abruptly between 0° and 180° (Figure 3.6d), while the phase between the wind and the bottom currents at the ADCP shifted from -90° to +90° (Figure 3.6e). These phase relationships will also be seen in the TVC results. 3.3.3 Two-layer Variable Cross-section Model (TVC) Results For day 179 to 211, the baroclinic flow at the ADCP was simulated using the TVC with different values of the damping ratio for the first horizontal mode, ζ1. For the times when two-layer exchange was observed (and excluding day 186, 187, and 188 as explained below), the root mean square (RMS) error was calculated between the simulated and observed flows in the upper and lower layers. The value of ζ1 = 1.3 gave the least RMS error, 15 m3s-1. For a 5% increase in this RMS error, the value of ζ1 ranged between 1.0 and 1.7 which corresponds to a change of -23% and +31% from the optimum ζ1 = 1.3, respectively. For comparison, in Lake Tiberias, a 5% increase in the RMS error resulted in a range of ζ1 between -25% and +62% of the optimum ζ1 (Shimizu and Imberger, 2008) which shows a higher sensitivity of the RMS error to ζ1 than for Amisk Lake. This is attributable to the stronger damping in Amisk Lake (ζ1 = 1.3) than in Lake Tiberias (ζ1 = 0.08) (Shimizu and Imberger, 2008). In the preceding calibration of ζ1, we excluded days 186-188 which contained only ∼ 14% of the occurrences of two-layer exchange between day 179 and 211. During days 186- 188, non-linearities seemed to affect the observed flows and, because of its linearity, the TVC Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 55 14d 4d 2d 1d 12h 6h Period S p ec tr a l D en si ty (N 2 m −4 H z− 1 ) (a)10−1 101 S p ec tr a l D en si ty (◦ C 2 H z− 1 ) 5m 8m (b)102 104 S p ec tr a l D en si ty (m 2 s− 2 H z− 1 ) Top Bott om (c)10−2 101 0 0.5 1 C o h er en ce (d) −180 −90 0 90 180 P h a se (◦ ) 10−6 10−5 10−4 0 0.5 1 Frequency (Hz) C o h er en ce (e) −180 −90 0 90 180 P h a se (◦ ) Figure 3.6: Power spectra of (a) the along-thalweg component of wind stress, (b) the fluctuations of water temperature at 5 m and 8 m at the mooring, and (c) the along-thalweg component of velocity in the surface and bottom bins of the ADCP. Spectra were generated after removing linear trends and were smoothed in the frequency domain (Saggio and Imberger, 1998). Dashed lines show the 90% confidence interval based on a chi-squared distribution. Arrows in panel (a) indicate the estimated periods of the H1, H2, and H3 modes. Shown in (d,e) are the magnitude-squared coherence (solid) and phase (boxes) between (d) the along-thalweg wind stress and the water temperature at 5 m and (e) the along-thalweg wind stress and the velocity in the bottom bin. Dashed lines show coherence at the 90% confidence level (Thompson, 1979). In all panels, Grey bars indicate the ∼ 1 d and 7.3 h periods described in the text. Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 56 significantly overestimated the peak flows. If the flows from days 186-188 were accounted for, the minimum RMS error increases to 23 m3s-1 corresponding to ζ1 = 1.75. In summary, the value of ζ1 & 1 which indicates that the response of the H1 mode is overdamped, directly following the wind pattern without free oscillations. Using the value of ζ1 = 1.3 shows good agreement between the observed and simulated flows in the upper and lower layers (Figure 3.7b,c). The similarity between the baroclinic flows and the wind pattern can be seen by comparing panels (a) and (b,c) in Figure 3.7. In addition to the flow, the TVC was used to simulate the deflections of the thermocline at the location of the mooring. In Figure 3.7d, the simulated deflection is compared to the average deflection of the 12° and 16°C isotherms obtained from the temperature records at the mooring by linear interpolation and detrending. The simulated deflection reproduces the same order of magnitude of the observed deflections, ∼ 1 m, with an overall RMS error of 0.4 m. The discrepancy between the simulated and observed deflections may be attributed to three factors. First, the deflections inferred from the thermistors (. 1 m) are small compared to the spacing between the thermistors (2−3 m) which introduces errors in the interpolation of isotherms. Second, the wind stress may not be uniform which affects the along-thalweg distribution of Fs and, subsequently, the amplitude An of the baroclinic modes (equation 3.4). Third, the modal periods of the H2 and H3 modes differed from the 7.3 h period detected in the power spectra of the temperature records (Figure 3.6b). Because the H2 and H3 modes are weakly-damped, their modal responses are shifted in time from the observed 7.3 h deflections. Nevertheless, the TVC reproduced the modal components of the deflections as will be shown below. We now examine the modal composition of the simulated response at the ADCP and the mooring. Because of the strong damping of H1 and the short modal periods of higher modes, the transient response to each forcing frequency ω f ,m decays rapidly (< 2 d) relative to the duration of the simulation (32 d). Accordingly, to characterize the simulated response, it is sufficient to analyze the steady-state forced oscillations whose amplitudes are given by |Dn,m|Qn (x) for layer flow and by |Dn,m|Zn (x)ω−1f ,m for interface deflection (equation 3.5 and 3.8). The phase difference between the wind and the steady oscillations of the flow are given by4 ∠F(n,m)s −∠q(n,m)2 = ∠ ( Bm (Dn,mQn) −1). The interface deflections lag the flow by, pi/2, a quarter of a period (compare equation 3.5 and 3.8). For the flow at the ADCP, the H1 mode dominates the distribution of |Dn,m|Qn over 4The symbol ∠ denotes the complex phase. Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 57 −0.09 0.09 W in d S tr es s, τ s (N m −2 ) (a) −80 0 80 F lo w (m 3 s− 1 ) (b) −80 0 80 F lo w (m 3 s− 1 ) (c) Simulated, H1-H9 Simulated, H1 Observed 180 185 190 195 200 205 210 −1 0 1 Julian Day in 1991 D efl ec ti o n (m ) (d) Simulated Observed Figure 3.7: (a) Wind stress over Amisk Lake. Observed and simulated (H1-H9) flow in the (b) upper and (c) lower layers at the ADCP. Observed flow corresponds to two-layer exchange. Also shown is the contribution of the H1 mode to the simulated flow. (d) Observed and simulated interface deflection at the thermistor mooring. In (a-c), stress and flow are positive from south to north. Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 58 most of the forcing frequencies ω f ,m, in particular, at the ∼1 d period (Figure 3.8a). In comparison, the amplitudes of |Dn,m|Qn are insignificant for H2 and H3 except for high frequencies corresponding to the natural frequencies of these modes. The relative amplitudes agree with the spectral analysis of the ADCP data for the surface and bottom currents. There is also agreement in the phase of the baroclinic response relative to the wind (Figure 3.8a). Around the period of H1 (∼1 d), the TVC shows that the wind (F(n,m)s ) and the flow in the bottom-layer (q(n,m)2 ) are out of phase. In the vicinity of the periods of H2 and H3, phase differences vary abruptly between−pi/2 and pi/2 as was also seen in the ADCP currents (Figure 3.6e). The abrupt phase shift for H2 and H3 relative to H1 is attributable to the weak damping of the higher modes, since for linear oscillators the abruptness of the phase shift at the modal period increases for weaker damping (Humar, 2002; Bendat and Piersol, 2010). For the interface deflections at the mooring, the H1 mode dominates the distribution of |Dn,m|Zn (x)ω−1f ,m for low frequencies (Figure 3.8b). Each of the H1-H3 modes exhibit a peak in |Dn,m|Zn (x)ω−1f ,m corresponding to their respective modal periods. The peak for H2 is larger than for H1 and H3. The H2 mode is likely responsible for the 7.3 h spectral peak in the temperature fluctuations as H2 is more manifested than H3 in the simulated flow at the ADCP and interface deflections at the mooring. The phase differences between the wind and the interface deflections in the TVC agree well with the observed phase differences (compare Figure 3.8b to Figure 3.6d) noting that thermocline deflections are in anti-phase with temperature fluctuations. To examine the spatial variability of the baroclinic response of Amisk Lake, we used the TVC to simulate the response of the H1-H3 modes to steady wind forcing by setting ω f ,0 = 0 in equation 3.4, 3.5, and 3.8. Figure 3.9a shows the along-thalweg distribution of wind components F(n,0)s for n= 1−3 based on a spatially-uniform stress of 0.02 N m-2 which is the average stress between day 205.4 and 205.9 (Figure 3.7). In response to the sustained wind forcing, the interface asymptoted to the dynamic equilibrium deflections shown in Figure 3.9b. The prominent feature is that the H1 mode dominates the interface deflection at the endwalls of Amisk lake. The deflections associated with H2 and H3 become more significant relative to H1 in the vicinity of the connecting channel. Within the channel, the peak layer flow is dominated by the H1 mode (Figure 3.9c) as found from analyzing the flow at the ADCP. The H2 and H3 modes become significant for the flow in the south and north basins, respectively (Figure 3.9c). Finally, the along-thalweg distribution of average velocity in the lower layer shows that the highest velocities occur in the channel and are associated with Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 59 14d 4d 1d 8h 2h Period, 2piω−1f,m F lo w A m p li tu d e (m 3 s− 1 ) (a) 0 2 4 -pi 0 pi P h a se (r a d ) Frequency, (2pi)−1ωf,m (Hz) D efl ec ti o n A m p li tu d e (× 1 0 −2 m ) (b) 10−6 10−5 10−4 0 3 6 -pi 0 pi P h a se (r a d ) H1 H2 H3 Figure 3.8: Amplitude spectra of simulated steady-state oscillations for (a) flow at the ADCP, |Dn,m|Qn, and (b) interface deflection at the mooring, |Dn,m|Znω−1f ,m. Also shown are the phase differences between the wind components F(n,m)s and the associated (a) flow at the ADCP, ∠F(n,m)s −∠q(n,m)2 = ∠ ( Bm (Dn,mQn) −1), and (b) the interface deflection at the mooring, ∠F(n,m)s −∠q(n,m)2 − pi2 . The amplitude spectra are symmetric about ω f ,m = 0 while the phase spectra are anti-symmetric (not shown). Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 60 −10 0 10 20 F (n ,0 ) s (N m −1 ) n = 1 n = 2n = 3(a) In te rf a ce D efl ec ti o n (m ) (b) −1 0 1 H1 H2 H3 F lo w (m 3 s− 1 ) (c) −40 −20 0 20 x (km) V el o ci ty (× 1 0 −2 m s− 1 ) (d) 0 2 4 6 8 −4 −2 0 2 Figure 3.9: (a) Components n = 1− 3 of wind forcing corresponding to a steady spatially- uniform wind stress of 0.02 N m-2. Shown for the H1, H2, and H3 modes are (b) the steady- state interface deflection, (c) the peak flow in the lower layer, and (d) the peak depth-averaged velocity in the lower layer. The markersq ands indicate the locations of the ADCP and the thermistor mooring, respectively. Grey shading indicates the approximate extent of the channel connecting the north and south basins. In (a,c,d), direction of forcing, flow, and velocity is positive from south to north. H1 (Figure 3.9d). This high velocity above the channel bed is perhaps why the H1 mode is strongly damped relative to the H2 and H3 modes. 3.3.4 Multi-layer Model (MLM) Results The velocity profiles at the ADCP site showed that the exchange through the connecting channel occurred over three or more layers for about 50% of the time between day 179 and 211. These observations indicate vertical mode two and higher. To model these modes, we applied the MLM to the box basin (BB) idealization of Amisk Lake. Four simulations were conducted of which two are discussed in this section. Because of the secondary Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 61 Table 3.2: Parameters for simulations using the MLM. Simulation Thickness (m) Density (kg m-3) Period (h) Damping h1 h2 h3 ρ2−ρ1 ρ3−ρ2 V1H1 V2H1 ζ1,1 ζ1,2 s1 4.5 15.0 - 1.618 - 19.3 - 1.3 - s2 3.0 4.5 14.3 0.965 0.784 20.4 39.4 2.5 s3 3.0 1.0 17.8 0.965 0.784 21.1 72.3 4.4 s4 1.5 4.5 15.8 0.965 0.784 23.5 46.0 2.6 contribution of higher horizontal modes to the exchange flow at the ADCP site, only the first horizontal mode, H1, was considered in the MLM simulations. The damping ratio of the V1H1 mode was fixed at the TVC value, ζ1,1 = 1.3, and the ratios for higher modes were set to ζ1,r = ζ1,1 ( ω1,1ω−11,r ) (Table 3.2), assuming that κn = ζ1,rω1,r depends on the horizontal velocity structure (n) but not the vertical velocity structure (r). The empirical formula by Spigel and Imberger (1980) indicates that the e-folding time-scale of baroclinic modes, Te,n = κ−1n , is proportional to √ωn,r (Chapter 4) which suggests higher damping ratios, ζ1,r = ζ1,1 ( ω1,1ω−11,r )1.5 , than assumed. Stronger damping does not affect our results on the exchange structure and, hence, we chose the smaller damping given by ζ1,r= ζ1,1 ( ω1,1ω−11,r ) (also, see Chapter 5). The box basin idealization did not distort the simulated flow at the ADCP. Based on the two-layer stratification used in the TVC (Table 3.2), the first simulation, s1, adequately reproduced the flows simulated by the TVC at the ADCP (not shown). The linear correlation coefficient between the s1 and TVC flows was 0.97 and the RMS error was 8 m3/s which is smaller than between the TVC and the observed two-layer flows (15 m3s-1). The second simulation, s2, used a three-layer stratification to account for the thickness of the metalimnion (Figure 3.10a). The three-layered stratification slightly increased the period of the first vertical mode, V1H1, from 19.3 h in s1 to 20.4 h in s2 (Table 3.2). More importantly, the additional density interface introduced a second vertical baroclinic mode, V2H1, with a period of 39.4 h. Noteworthy, close to the period of V2H1, a spectral peak at a period of 43 h is observed in the power spectrum of the temperature fluctuations at the 8 m thermistor (Figure 3.6b). For three-layer stratification, the exchange flow at the ADCP was classified into two distinct cases. Two-layer exchange occurred when two adjacent layers flowed in the same Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 62 998 10000 3 7.5 21.8 Density (kg m−3) D ep th (m ) (a) (b) (c) (d) Surface flow Intermediate flow Bottom flow Figure 3.10: (a) Three-layer stratification for simulation s2 in the box basin. (b,c) Two-layer exchange. (d) Three-layer exchange. direction opposite to the third layer (Figure 3.10b,c), while three-layer exchange occurred when the flow in the intermediate layer was opposite to the other two layers (Figure 3.10d). The flow in single or adjacent layers was labeled as surface, intermediate, and bottom as shown in Figure 3.10b-d. In simulation s2, three-layer exchange occurred 16% of the time (Figure 3.11). During three-layer exchange in s2, the magnitude of the intermediate flow was comparable to that of the observed flow (Figure 3.11b). The 10, 50, and 90th percentiles of the magnitude of intermediate flow in s2 were 3, 7, 14 m3s-1 which closely corresponds to 5, 9, and 14 m3s-1 for the observed flow. For the remaining 84% of the time, two-layer exchange occurred in s2 because the three-layer idealization of the stratification did not allow more than two degrees of freedom for the exchange structure (Figure 3.10b-d). By comparing the surface and bottom flows in s2 to those from simulation s1 (Figure 3.11a,c), we find that the flows corresponding to the two-layer exchange were only slightly influenced by the addition of the third intermediate layer in s2. Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 63 −80 0 80 F lo w (m 3 s− 1 ) (a) 2−layer 3−layer F lo w (m 3 s− 1 ) (b) −80 0 80 Observations 3−layer 180 185 190 195 200 205 210 −80 0 80 Julian Day in 1991 F lo w (m 3 s− 1 ) (c) 2−layer 3−layer Figure 3.11: (a) Surface, (b) intermediate, and (c) bottom flows from the MLM simulation s2 (three-layered) at the ADCP. Surface and bottom flows from the MLM simulation s1 (two- layered) are shown in panels (a) and (c), respectively. Panel (b) also shows the observed flow in the intermediate layer when three-layer exchange occurred. Flow is positive from south to north. Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 64 3.4 Discussion 3.4.1 Resonance with Wind Comparison of the TVC and MLM results to observations in Amisk Lake suggests that most of the baroclinic response to wind in lakes can be understood through linear oscillators with multiple internal degrees of freedom (multiple baroclinic modes; r= 1→N−1 and n= 1→ ∞). When lakes are forced with periodic wind, the TVC shows that the amplitude of steady modal oscillations in layer flow at any location along the thalweg, |Dn,m|Qn, is the product of two terms (equations 3.5-3.7). The first term AnQn (x) is a spatial component dictated by the coupling between the along-thalweg wind distribution and bathymetry. This coupling is manifested in using the Qn structures obtained for the lake bathymetry (equation 3.1) to decompose the along-thalweg wind distribution and obtain An (equation 3.4). The second term |Dn,m|A−1n is a temporal component, proportional to the wind stress |Bm| and inversely proportional to ωu,n and ζn. Additionally, this temporal component, |Dn,m|A−1n , depends on the ratio between the forcing and modal frequency ω? = ω f ,mω−1u,n such that |Dn,m|A−1n is maximum at ω? = 1 and decreases as ω? deviates from unity (Humar, 2002). Wind forcing most effectively energizes each baroclinic mode near its respective modal frequency ωu,n. If the spectral energy of the wind is concentrated (i.e. |Bm| is highest) near a particular ωu,n, the transfer of wind energy to the corresponding baroclinic mode is maximized and the flow oscillations for that mode are amplified as observed in many lakes (e.g., Thorpe, 1974; Wiegand and Carmack, 1986; Antenucci and Imberger, 2003; Vidal et al., 2007). The temporal component |Dn,m|A−1n normalized by wind stress, ∣∣D?n,m∣∣= ∣∣Dn,mB−1m ∣∣A−1n , is inversely proportional to ωu,n. For Amisk Lake, the peak ∣∣D?n,m∣∣ for the H1 mode is ideally ωu,2ω−1u,1 = 2.6 times the peak ∣∣D?n,m∣∣ for the H2 mode if the damping of the two modes is the same, ζ1 = ζ2. However, because the H1 mode is more strongly damped, the peak of ∣∣D?n,m∣∣ for H1 is only 0.1 times that for H2. This shows the importance of modal damping in the relative amplification of baroclinic modes. Finally, because the spectral content of the wind is concentrated near ωu,1 (Figure 3.6a), the maximum of the component |Dn,m|A−1n for H1 rises to 0.4 relative to that of H2. This in turn shows the effect of the spectral distribution of the wind relative to the various modal frequencies. The dependence of modal resonance on ζn and ω? deserves further discussion. Consider first the normalized amplitude of flow ∣∣∣D?1,m∣∣∣ for a fixed modal frequency ωu,1 and a variable Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 65 forcing frequency ω f ,m (Figure 3.12, solid lines). As noted previously, the maximum ∣∣D?n,m∣∣ corresponds to ω? = 1. For strong damping, the ∣∣D?n,m∣∣ curves flatten up considerably in the vicinity of ω? = 1. For Amisk Lake, because of the strong damping of the H1 mode, the flow associated with H1 responds with similar magnitude to a broad bandwidth of forcing frequencies aroundωu,1. For 0.5≤ω?≤ 2, ∣∣∣D?1,m∣∣∣ varies by less than 14% (Figure 3.12, ζ1 = 1.3). Unfortunately, direct comparison to resonance bandwidth from experimental results was not possible as available results under periodic forcing are for weak damping (ζ1 < 0.1) (e.g., Boegman and Ivey, 2007; Wake et al., 2007). It is also important to examine the change in the normalized flow amplitude ∣∣∣D?1,m∣∣∣ due to a variable ωu,1 relative to a fixed dominant forcing frequency ω f ,m (Figure 3.12, dashed lines). For example, Lawrence et al. (1997) assumed the baroclinic exchange between the basins of Amisk Lake was stationary over a period of three months during the summer although ωu,1 changes with the background stratification (within about ±25% of the average ωu,1 value). For ζ1 = 1.3, a 20% decrease in ωu,1 from ω? = 1.1 to 0.9 leads to an 18% decrease in ∣∣∣D?1,m∣∣∣. Assuming that wind forcing is stationary and modal structures Qn do not vary considerably, the flow amplitude ∣∣D1,m∣∣Qn decreases as stratification weakens (ωu,1 decreases). The exchange between the basins of Amisk Lake is probably higher during the middle of the summer than in early- and late- summer when stratification is weaker. Velocity measurements, however, are not available to validate this hypothesis. 3.4.2 Excitation of Higher Horizontal Modes As has been noted, the record of water temperature from the 5 m deep thermistor at the mooring suggests higher horizontal baroclinic modes. The 7.3 h peak in the power spectrum of the record is attributable to the second horizontal mode, H2, as simulated with the TVC. Higher horizontal modes have been observed in other lakes. In the south basin of Winder- mere, Mortimer (1953) identified late-summer deflections of the thermocline as V1H2 mode and suggested that the cause of excitation was the irregular bathymetry or irregular wind distribution. In Lake Biwa, Shimizu et al. (2007) showed that the wind direction relative to the surface geometry influenced which horizontal modes were excited. The TVC model is linear and up to this point the explanation of higher modes is based on linear mechanisms. Note that a different explanation involving a non-linear mechanism has been proposed for the excitation of higher horizontal modes. Lemmin (1987) reported Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 66 0 0.5 1 1.5 2 0 2 4 6 8 10 12 14 16 Normalized frequency, ω⋆ N o rm a li ze d F lo w A m p li tu d e (m 3 s− 1 [N m −1 ]− 1 ) ζ1 = 1.0 1.0 1.3 1.3 1.7 ζ1 = 1.7 2piω−1u,1 =21.4 h, variable ωf,m 2piω−1f,m =24.0 h, variable ωu,1 Figure 3.12: Amplitude of steady flow oscillations for the H1 mode as a function of the ratio of forcing to modal frequency, ω? = ω f ,mω−1u,n . Two sets of curves are shown, one set for variable ω f ,m and the second for variable ωu,1 (see legend). Each set shows curves for ζ1 = 1.0, 1.3, and 1.7. The flow amplitude is normalized by wind forcing and modal structure, ∣∣∣D?1,m∣∣∣. Note the flatness of the solid curves over a wide range around the modal frequency, ω? = 1. Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 67 observations of the H2 to H5 modes in Baldegg Lake along with a dominant H1 mode. Lem- min (1987) hypothesized that these higher modes were excited through an energy cascade whereby the wind imparts energy into the H1 mode and, through non-linear interaction, the energy in the H1 mode is transferred into higher modes. The cascade of energy contin- ues into higher modes of smaller length scales down to the scales of viscous dissipation. Lemmin (1987) partly supported his hypothesis by showing the consistency between the power spectra of isothermal deflections and the Garrett-Munk spectrum. The explanation by Lemmin (1987) led Lemmin et al. (2005) to distinguish between large and small lakes in regards to the excitation of horizontal modes. In lakes as large as Lake Geneva, Lemmin et al. (2005) linked the excitation of horizontal modes to the spatial wind distribution while, in small lakes, Lemmin et al. (2005) maintained that other factors were responsible for the excitation. It should be noted that Lemmin (1987) did not pursue non-linear analysis to support his hypothesis and also did not examine if the spatial wind distribution over Baldegg Lake contributed to the excitation of the higher modes. Note that in laboratory experiments an energy flux from lower to higher horizontal modes did not occur. For example, Boegman et al. (2005b) observed the evolution of an initially sloped density interface in a two-layer stratified laboratory flume. Almost all of the potential energy in the initial interfacial slope (∼ 99%)was contained in the H1 mode. As the interface relaxed, this energy was partly transferred to a progressive internal surge and partly lost by viscous dissipation. There was no excitation of higher linear horizontal modes other than the ones that were triggered by the initial conditions which is equivalent to wind forcing in our case. The TVC model does not include non-linear coupling between modes, instead it is the coupling between the spatial distribution of wind forcing and basin bathymetry that determines which baroclinic modes are excited (equation 3.1 and 3.4). For Amisk Lake, this coupling is responsible for the excitation of the observed H2 mode and possibly H3. The along-thalweg distribution of wind forcing corresponding to n = 1, F(1,m)s , is particularly peaked over the connecting channel (Figure 3.9a). This peaked shape results from the quantity ( S1S−12 +1 ) Q1 in equation 3.4 being largest for the connecting channel. To have any considerable forcing over the north and south basins (see Fs in Figure 3.3a), higher wind components are required (Figure 3.9a, F(2,m)s and F (3,m) s ). These components trigger corresponding higher baroclinic modes (Figure 3.9b,c). Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 68 3.4.3 Higher Vertical Modes and Exchange Flow As previously noted, three-layer exchange at the ADCP site is associated with smaller total volumetric transport (Figure 3.5). The MLM simulation s2 reveals a very similar pattern. For two-layer exchange, the average transport in s2 was 57 m3s-1 with a standard deviation of 54 m3s-1. For the three-layer exchange, the average transport was 23 m3s-1 with a standard deviation of 16 m3s-1. To explain this feature, we examine the baroclinic flows driven by steady wind. As previously explained, these flows were used to construct the s2 response between day 179 and 211 using convolution. For the V1H1 mode, the velocity structure is two-layered with flow in the upper and intermediate layers being opposite to that in the lower layer (Figure 3.13a). For the V2H1 mode, the exchange is three-layered. The flow in the intermediate layer is opposite to the flow in the overlying and underlying layers (Figure 3.13b). The V1 mode dominates the response in the lower layer while the V2 dominates the flow in the intermediate layer. The summation of flows from V1 and V2 results in the intermediate and lower layers flowing in the same direction opposite to the upper layer. This two-layer exchange initially dominates. However, because the V1 mode has a shorter period and decays faster than V2, the V2 mode subsequently dominates leading to a three-layer exchange (Figure 3.13c). During this phase (t & 0.5 d), the flows are considerably diminished compared to the flows in the initial mixed V1-V2 phase (t . 0.5 d) (Figure 3.13c). Changes in the thicknesses of the layers do not change the qualitative picture of how the exchange at the ADCP evolves with time. To show this, we ran two additional MLM simulations, s3 and s4 (Table 3.2). Compared to s2, simulation s3 used a thicker intermediate layer while simulation s4 used a thinner upper layer. The total depth was kept constant on the expense of the lower layer. The density differences across the two interfaces were also kept constant. In s3, the small thickness of the intermediate layer reduced the importance of the V2 mode relative to V1 (Figure 3.13d,e). In s4, the shallower upper layer increased the flows associated with the V2 mode, but the V1 flows in the intermediate and lower layers were also increased (Figure 3.13g,h). Despite the changes in the simulations s3 and s4, two-layer exchange initially dominated with the intermediate and lower layers flowing in the same direction. Subsequently, the V2 mode dominated with a three-layer exchange as the flow in the lower layer reversed to the direction of the upper layer (Figure 3.13f,i). In summary, the thickness of the metalimnion mainly affects the time at which the exchange shifts from two- Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 69 −1 0 1 s2 V 1H 1 (a) −1 0 1 N o rm a li ze d F lo w V 2H 1 (b) 0 0.5 1 1.5 −1 0 1 Time (day) V 1H 1+ V 2H 1 (c) s3 (d) Upper Intermediate Lower (e) 0 0.5 1 1.5 Time (day) (f ) s4 (g) (h) 0 0.5 1 1.5 Time (day) (i) Figure 3.13: Baroclinic flows simulated using the MLM based on the stratifications of s2 (a-c), s3 (d-f), and s4 (g-i) and forced by step wind. The flows are shown for modes V1H1 (top panels), V2H1 (middle panels), and the summation of V1H1 and V2H1 (bottom panels). Flows are normalized by the maximum flow from simulation s3. Arrows in the bottom panels indicate the shift from two- to three- layer exchange. to three- layer; the smaller the thickness, the later the shift (Figure 3.13c,f). For the same metalimnetic thickness, a shallower epilimnion slightly delays the shift and causes larger three-layer exchange flows due to V2H1 (Figure 3.13c,i). The response to the idealized steady wind explains the observed exchange structure in Amisk Lake. At the beginning of a wind storm, the V1 and V2 modes are both excited with a two-layer exchange due to the dominance of V1. Wind events over Amisk Lake do not last for long and when the wind subsides, the flows due to the V1 and V2 modes are readjusted in the opposite direction to the flows during the wind event. As the V1 mode decays, the V2 takes over albeit with diminishing flow magnitudes. Our results from the MLM simulations are in agreement with observations by Wiegand and Chamberlain (1987) in Wood Lake, a small single-basin lake. Wiegand and Chamber- lain (1987) observed that the V1 mode initially dominated the deflections of metalimnetic Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 70 isotherms. As the V1 deflections decayed, the V2 mode appeared prominently in the de- flections of the upper and lower metalimnetic isotherms. Similar fast decay of the V1 mode relative to higher modes (V2 and V3) was also observed in Lake Bromont (Pannard et al., 2011). 3.5 Summary and Conclusions To characterize the baroclinic response of Amisk Lake to wind forcing, we have formulated and applied the TVC model which generalizes the linear dynamic model of Heaps and Ramsbottom (1966) and the two-layer Defant procedure of Lemmin and Mortimer (1986). The TVC describes the baroclinic response to wind forcing which is spatially and temporally variable. The response is synthesized from 1) the structures and periods of baroclinic modes which are determined by lake bathymetry and stratification, 2) the spectral components of the wind time-series, 3) the modal magnitudes An determined by the along-thalweg wind distribution, and 4) the damping parameters of the modes (ζn). For the purpose of analy- sis, the TVC is advantageous over hydrodynamic numerical models which do not directly reveal the modal composition of the baroclinic response or provide basic explanation of the coupling between the response and wind characteristics, bathymetry, and stratification. As opposed to post-decomposition of the response simulated from numerical models, the baroclinic response in the TVC is constructed from the elementary spatial and temporal components. By applying the TVC to Amisk Lake, we have found that the excitation of higher hori- zontal modes results from the coupling between the along-thalweg wind distribution and the bathymetry of the lake. Higher modes (H2 and H3) contribute to the flow within the basins of Amisk Lake but the exchange between the basins is dominated by the lowest horizontal mode, H1. Despite the deviation of the period of H1 (21.4 h) from diurnal periodicity, the flow associated with H1 resonates with the strong diurnal periodicity of the wind over Amisk Lake. This results from the strong damping of the H1 mode which broadens the range of resonant frequencies. In contrast, damping of the second horizontal mode, H2, was found to be much weaker. This is explained by noting that H1 is associated with the highest velocities within the connecting channel (Figure 3.9c). The corresponding bottom drag would be stronger than for higher horizontal modes. Because of the thickness of the metalimnion, higher vertical modes were excited in Chapter 3. The Baroclinic Response to Wind in a Small Two-basin Lake 71 Amisk Lake. Unlike most previous studies which detect vertical modes from records of water temperature at basin endwalls (Wiegand and Chamberlain, 1987; Roget et al., 1997; Vidal et al., 2007), we have detected vertical modes from observed water currents in the channel that connects the basins of Amisk Lake. The observations and a three-layered box model show that a two-layer flow structure dominates the baroclinic exchange between the basins. For the initial phases of the response to wind, the superposition of the first (V1H1) and second (V2H1) vertical modes favors two-layer exchange. With similar damping time- scales (κ−1n ), the V1H1 mode which has a smaller modal period is damped faster than V2H1. The exchange shifts to three layers, but with smaller magnitudes. This may explain similar observations of the relative timing of the V1 and V2 modes in other lakes; e.g., Wood Lake (Wiegand and Chamberlain, 1987). For Amisk Lake, implications of our study include the importance of the first horizontal mode and two-layer exchange to the transport of oxygen, and other elements, between the north and south basins. The H1 mode likely dominates the interbasin exchange over the stratified season. Because of its strong damping, the H1 mode would remain in near- resonance with the diurnal wind despite seasonal shifts in the modal period. 72 Chapter 4 The Damped Baroclinic Response to Wind in a Reservoir 4.1 Introduction Since the work of Mortimer (1952; 1953), weakly-damped metalimnetic oscillations have been emphasized as the “normal” baroclinic response of stratified lakes to wind forcing (e.g., Wiegand and Carmack 1986; Lemmin 1987; Roget et al. 1997; Shimizu and Imberger 2008). Metalimnetic tilts are set up during wind storms and, after the wind subsides, the tilts relax triggering free oscillations of the metalimnion. The oscillations persist for several cycles until they decay or are disturbed by new wind storms (Mortimer, 1952). In contrast, there are few observations in lakes of over-damped baroclinic responses for which wind-forced metalimnetic deflections slowly rebound to static equilibrium after the wind ceases (e.g., Imberger 1985). The forced deflections and free oscillations of the metalimnion are composed of basin- scale baroclinic modes. These modes can have complicated horizontal and vertical nodal structures dictated by bathymetry and stratification (Wiegand and Chamberlain 1987; Lem- min et al. 2005; Vidal et al. 2007). Identification of these excited modes is important as they contain most of the wind energy transferred to the interior of stratified lakes (Roget et al. (1997); Horn et al. (2001); Wüest and Lorke (2003)). For weakly-damped baroclinic responses, the excited baroclinic modes are identified from their spectral signatures in the records of water temperature or currents (Hutter et al. 1983; Münnich 1996; Fricker and Nepf 2000; Lemmin et al. 2005). For this identification, the periods of the baroclinic modes are estimated using modal-analysis models (MAMs) which typically involve some idealization of the stratification, bathymetry, or both. For over- damped baroclinic responses, identification of excited modes from their spectral signatures is not possible. Modal-based forced models (MFMs) such as those of Heaps and Ramsbottom Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 73 (1966), Monismith (1985; 1987), and Shimizu and Imberger (2008) are needed to identify the excited modes. MFMs use periods and spatial structures of baroclinic modes obtained from modal analysis models to simulate the evolution of individual modes under wind forcing. Calibration of the MFM response against field data indicates the strength of damping for each mode (Heaps and Ramsbottom, 1966; Shimizu and Imberger, 2008). In this paper, we analyze the baroclinic response to wind in a complex reservoir where the main excited baroclinic mode is shown to be over-damped. Identification of the excited modes is first attempted using spectral analysis of temperature data. For this identifica- tion, the periods of baroclinic modes are obtained using a MAM which accounts for depth variability in the two horizontal dimensions and uses multi-layer stratification. Because the baroclinic response is found to be overdamped, an MFM is used to decompose and simulate the response. This MFM uses two density layers and accounts for cross-sectional variability along the thalweg of the domain. The paper is structured as follows. In the next section, we describe the study site, the field dataset, and details of the used MAM and MFM. In results, we present the field observations and the outcome of the models. The discussion section includes the effect of bathymetric idealizations on the results, the detection and source of the strong damping of the baroclinic response, and the implications for the field site. Also included in the discussion section is a comparison to damping in other lakes. 4.2 Methods 4.2.1 Site Description Nechako Reservoir, located in the Interior Plateau of British Columbia, Canada, is composed of a series of basins which are connected through narrow and shallow flooded river reaches (Figure 4.1). The basins of Nechako Reservoir form a ring shape which extends 200 km east- west and 75 km north-south. Contrasting with this large geographical extent (15×103 km2), the surface area of the reservoir1 is only 1135 km2 reflecting the narrowness of the basins. Nechako Reservoir was impounded by the Kenney Dam in the early 1950s. There is currently no outlet at Kenney Dam. Water is withdrawn at the west end of the reservoir for hydroelectric power generation, and at the Skins Lake spillway, 80 km west of Kenney Dam 1See Appendix A Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 74 125°W126°W127°W 5 3 °4 5 ’N 5 3 °1 5 ’N Kenney Dam Nechako River Spillway Tunnel to Powerhouse Narrows Natalkuz Lake Nechako Reservoir Southern arc Northern arc SK KN NA 0 20 4010 km Tetachuck River Knewstubb Lake Figure 4.1: Map of Nechako Reservoir showing the locations of Kenney Dam, the spillway, and the powerhouse offtake. The focus of this study is the two eastern basins, Knewstubb and Natalkuz Lakes (shaded). The TVC model is applied along the thalweg of the two basins (dashed-dotted line). Marked are the locations of thermistor moorings and wind stations (KN, NA, SK). Top left insert shows the location of Nechako Reservoir within the province of British Columbia, Canada. Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 75 (Figure 4.1). A cold water release facility has been proposed at Kenney Dam to regulate summertime river temperature downstream in the Nechako River (Mitchell et al. 1995). In examining the baroclinic response of Nechako Reservoir to wind forcing, bathymetric idealizations were needed to deal with the complexity of the reservoir. The first level of idealization was decoupling the two basins upstream of Kenney Dam - Knewstubb and Natalkuz Lakes - from the rest of Nechako Reservoir. The two basins are baroclinically disconnected from the southern arc of the reservoir by the Tetachuck River (Figure 4.1). In contrast, the connection of Knewstubb and Natalkuz Lakes with the northern arc of the reservoir is deeper than the typical summertime metalimnion. However, as will be shown, the baroclinic dynamics observed at Kenney Dam can mostly be explained by considering only Knewstubb and Natalkuz Lakes which have a combined volume, surface area, and mean depth of 6.6 km3, 200 km2, and 34 m, respectively. The maximum recorded depth is 113 m in the middle of Natalkuz Lake while the depth at Kenney Dam, 84 m, is the maximum in Knewstubb Lake. Knewstubb and Natalkuz Lakes are connected through a contraction and sill, the Narrows (Figure 4.1). At the most constricted point, the Narrows is about 50 m deep and 500 m wide at the surface. 4.2.2 Field Data A comprehensive dataset was collected in 1994 (Triton 1995) and forms the basis of the work presented here. Between 24 Jun and 12 Oct (Julian days 175 to 285), weather stations were deployed on floating platforms in Knewstubb (KN) and Natalkuz (NA) Lakes 1.8 km and 38 km, respectively, along the thalweg (the deepest path) from Kenney dam (Figure 4.1). Each station logged hourly vector averages from a wind speed and direction sensor mounted 3 m above the water surface. Wind speed and direction were also recorded at a land-based station at Skins Lake spillway (Figure 4.1). Observations in 1994 included thermistor chains at platforms KN and NA. Thermistors accurate to 0.2°C were vertically spaced at 2 m intervals. The thermistors sampled and logged at 10 min intervals. Conductivity from CTD casts was low (35-65 mS cm-1 at 25°C) and the summertime density structure is controlled by temperature. To validate the model results, we used additional data collected in 2005 at platform KN. Between 20 Jul and 7 Oct (day 201 to 280), wind speed and direction were recorded at 1 min intervals. Moored temperature loggers at KN recorded water temperature with a similar Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 76 configuration to 1994. 4.2.3 Modal Analysis Model To identify the excited baroclinic modes, we calculated the natural (undamped) periods of the modes Tu using a layered, variable-depth MAM, the VDC (Salvadé et al. 1992). This model is based on the linear and hydrostatic approximations to the equations of motion and neglects mixing and rotational effects. For Knewstubb and Natalkuz Lakes, the internal Rossby radius from July to Sept 1994 ranged between 2.4 km and 3.3 km and is larger than the average surface width (1.9 km) which, to a first approximation, justifies neglecting the rotation of Earth in applying the VDC model. For a water body stratified into N-layers, the VDC gives the baroclinic modes by the eigen-value problem, −ω2uZi = gΣNj=1 { ∇ · ((ΣNk=1c jkdikhk)∇Z j)} i= 1,2 . . .N (4.1) where ωu = 2piT−1u is the undamped frequency of a given baroclinic mode; Zi (x,y) is the corresponding horizontal structure of the deflection of interface i from its equilibrium level (the surface corresponds to i = 1); hk is the equilibrium thickness of layer k (top layer corresponds to k = 1); g is the gravitational acceleration; ∇ = ( ∂ ∂x , ∂ ∂y ) ; and x and y are the horizontal coordinates. In equation 4.1, the coefficients c jk and dik are given by c jk = ρ j−ρ j−1 ρk k ≥ j 0 k < j , and dik = 1 i< k+10 else (4.2) in which ρk is the density of layer k (note that ρ0 = 0 for c1k). The baroclinic modes admitted by the VDC can be classified as VmHn where m ∈ Z+ and n ∈ Z+ denote the vertical and horizontal nodal structures, respectively (Münnich et al. 1992; Roget et al. 1997; Vidal et al. 2007). For N−layer stratification, there are N − 1 possible vertical nodal structures for baroclinic modes, m = 1,2 . . .N− 1 (Turner 1980). In this study, we focus on in-phase deflections of metalimnetic isotherms, i.e., V1Hn modes, since these have a larger effect on the cold water release facility at Kenney Dam. As will be discussed, higher vertical modes (m> 1) cause relatively small changes in the metalimnetic thickness at Kenney Dam. To apply the VDC to Knewstubb and Natalkuz Lakes, equation 4.1 was solved numeri- cally using a conservative finite difference scheme. A uniform horizontal grid of 100 m × 100 m Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 77 cells was used. In solving equation 4.1, a no-flow boundary condition was imposed at the perimeter of each density layer which means that the entire domain was accounted for including shallow zones (Salvadé et al., 1992; Shimizu et al., 2007). Because the evolution of the background stratification affects the periods of baroclinic modes (Wiegand and Carmack, 1986; Antenucci et al., 2000), we evaluated the modal peri- ods for four segments of the 1994 record. The segments, denoted as S1, S2, S3, and S4, are about four weeks each, centered around day 189 (8 Jul), 217 (5 Aug), 245 (2 Sept), and 272 (29 Sept), respectively. The parameters for the stratification of each segment were obtained by averaging the temperature profiles from the KN and NA moorings (Table 4.1). The modal periods of the V1H1 and V1H2 modes were first estimated using N = 3. For higher horizontal modes, the distinction between the V1 and V2 modes was ambiguous because the nodal lines of the surface and the density interfaces were not vertically aligned (Shimizu and Imberger, 2008; Vidal and Casamitjana, 2008). For this reason, we estimated the periods of higher horizontal modes H3 to H8 of the V1 structure using N = 2 which admitted only V1 modes. Using N = 2 eliminated the thickness of the metalimnion and likely led to underestima- tion of modal periods. For the four segments of the record, the V1H1 periods predicted by VDC ranged between 7.7-9.3 d for N = 3 and 5.4-8.0 d for N = 2. For the V1H2 mode, the VDC periods were 3.8-4.8 d for N = 3 and 2.7-3.8 d for N = 2. The greatest difference between using N = 2 and 3 occurred when the metalimnion was thickest, i.e. during record segment S2, and hence the lower bound of the modal periods was the most affected. 4.2.4 Modal-based Forced Model To simulate the baroclinic dynamics in Knewstubb and Natalkuz Lakes, we used the two- layer variable cross-section model (TVC, Chapter 3). The TVC is a modal-based forced model (MFM) which simulates the baroclinic response to wind by superposition of the separate responses of baroclinic modes. For our application here, we simplified the TVC by using impulsive wind forcing instead of periodic Fourier wind components. The TVC is applicable along the thalweg of elongated lakes. For Knewstubb and Natalkuz Lakes, we discarded Big Bend Arm, Intata Reach, and Euchu Arm and applied the TVC along the thalweg running from Kenney Dam to the west end of Chelaslie Arm (Figure 4.1). For this thalweg, we evaluated the periods Tn = 2piω−1n and spatial structures of interface Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 78 Table 4.1: Parameters of the two-layer (N = 2) and three-layer (N = 3) stratification for modal analysis using the VDC model. For each segment of the record, the top to bottom density difference is taken to be the same for N = 2 and 3 and is evenly distributed over N − 1 interfaces. The last row in the table gives the parameters used for the 1994 TVC simulation. Record segment Top to bottom N = 2 N = 3 Julian days density difference h1 h1 h2 (kg m-3) (m) (m) (m) S1 175-202 0.77 23.5 17 13 S2 203-230 1.40 20.5 11 19 S3 231-258 1.15 22.5 15 15 S4 259-285 0.72 24.5 19 11 S1-S4 175-285 1.01 22.7 - - deflections Z(n)2 for V1Hn modes by solving the eigenvalue problem 2, ω2Z2 =−g ( 1− ρ1 ρ2 ) 1 b2 d dx ( S1S2 S1+S2 dZ2 dx ) (4.3) where S1 (x) and S2 (x) are the cross-sectional areas of the upper and lower layers, respec- tively, and b2 (x) is the basin width at the equilibrium level of the interface between the two layers. Here, 0≤ x≤ 75.6 km is the distance along the thalweg (whereas, in equation 4.1, x and y were the horizontal coordinates). Under the action of wind, the deflection of the interface between the two layers η2 (x, t) is given by the superposition of the contributions from horizontal baroclinic modes, η2 (x, t) = Σ∞n=1η (n) 2 (x, t) where η (n) 2 is the deflection due to mode n, and t ≥ 0 is the time ordinate. The modal interfacial deflection η(n)2 is given by the convolution 3, η(n)2 (x, t) = An (t)?η (n) 2,imp (4.4) where An is the modal amplitude obtained by decomposing the wind forcing per unit length 2This form of the eigenvalue problem is equivalent to equation 3.1. It was more advantageous here to write the eigenvalue problem in terms of Z2 because we focus here on interface deflections and to provide a similar presentation to equation 4.1. 3The convolution integral is given by An (t)?η (n) 2,imp (x, t) = ´ t 0 An (t− t ′)η (n) 2,imp (x, t ′)dt ′. Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 79 of the thalweg, Fs (x, t), in terms of the baroclinic modal structures, Fs (x, t) = ( 1+ S1 S2 ) Σ∞n=1An (t)Q (n) 1 (x) (4.5) in which Q(n)1 (x) = ´ x 0 b2Z (n) 2 dx is the spatial structure of the upper layer flow. In equation 4.4, the interface deflection η(n)2,imp is the response to a unit impulse from a state of rest and is given by η(n)2,imp (x, t) = Nn (t) ρ1ωn exp(−ζnωnt)Z(n)2 (x) (4.6) where Nn (t) = 1√ 1−ζ 2n sin (√ 1−ζ 2n ωnt ) 0≤ ζn < 1 ωnt ζn = 1 1√ ζ 2n−1 sinh (√ ζ 2n −1ωnt ) ζn > 1 in which ζn is the ratio of the decay constant κn for a mode n to its undamped frequency ωn. Depending on the value of the damping ratio ζn, the baroclinic response can be underdamped (ζn < 1), critically-damped (ζn = 1), or over-damped (ζn > 1) (Arneborg and Liljebladh, 2001). The evolution of initial conditions should be superimposed on the forced response. For a given mode n, the interface deflection η(n)2,int relaxes according to (De Silva, 2007) η(n)2,int (x, t) = DnZ (n) 2 exp(−ζnωnt)× cos (√ 1−ζ 2n ωnt ) +ζnNn (t) 0≤ ζn < 1 1+Nn (t) ζn = 1 cosh (√ ζ 2n −1ωnt ) +ζnNn (t) ζn > 1 (4.7) where Dn is the scalar amplitude of the initial interface deflection associated with mode n, η(n)2,int (x, t = 0). To simulate the baroclinic response of Knewstubb and Natalkuz Lakes over the whole 1994 record (day 175-285), the upper layer depth and density difference were set to the average of the segments of the record S1-S4 (Table 4.1). In applying the TVC, the variation of width with depth was extracted from 70 cross-sections along the thalweg of Knewstubb and Natalkuz Lakes (Figure 4.2a). The cross-sectional widths were interpolated to regular intervals of 100 m along the thalweg and then used to compute the bathymetric parameters S1 (x), S2 (x), and b2 (x) (Figure 4.2b,c). The thalweg was forced by wind from platform KN over x< 28 km (east of the Narrows) and the rest was forced by NA wind. Derivation of wind stress from wind speed followed Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 80 the steps of Chapter 3 for Amisk Lake. Wind stress projected along the thalweg, τs, was assumed to be uniform across the channel and was laterally integrated to obtain the wind forcing Fs (x, t) (see example in Figure 4.2d). The along-thalweg distribution of Fs (x, t) was decomposed at 1 h intervals to obtain the modal amplitudes An (t) of the spatial wind com- ponents. This was done by multiplying equation 4.5 by Q(p)1 /S1 and using the orthogonality property, ´ L 0 ( S−11 +S −1 2 ) Q(n)1 Q (p) 1 dx= 0 if p 6= n (Shimizu et al., 2007), to get, An (t) = ´ L 0 S −1 1 Q (n) 1 Fsdx´ L 0 ( S−11 +S −1 2 ){ Q(n)1 }2 dx . (4.8) The relative weights of the wind components were assessed using the measure, Cn (t) = ´ L 0 Fs {( 1+S1S−12 ) AnQ1,n } dx´ L 0 F 2 s dx (4.9) which represents a single component n of the right hand side of equation 4.5 after multiplying both sides of that equation by Fs, integrating over x, and dividing by the left hand side. It follows from this derivation that Σ∞n=1Cn (t) = 1 and, if the distribution of wind forcing is composed of a single component, then Cn = 1 for that component and zero for all other components. Each wind component n with a modal amplitude An (t) 6= 0 triggers the corresponding horizontal baroclinic mode n (equation 4.4 and 4.5). To determine the relative importance of baroclinic modes at moorings KN and NA, we examined the equilibrium deflection, η2e, which is ultimately attained under steady wind. The deflection η2e is given by the convolution η(n)2,imp with a temporally constant wind forcing (equation 4.4). For a given spatial wind distribution, each mode attains an equilibrium deflection η(n)2e = T 2 n AnZ (n) 2 ( 4pi2ρ1 )−1 (Chapter 2). The deflection η(n)2e is proportional to the square of the undamped modal period as well as the product of the constant modal amplitude An and the magnitude of the deflection Z(n)2 evaluated at the moorings. In applying the TVC (which is a linear model) to Knewstubb and Natalkuz Lakes, we have assumed that the baroclinic response of the two lakes to wind forcing reflects a linear response in which baroclinic modes are damped faster than they can steepen and degenerate into internal surges and solitary waves. This assumption is justified by the linearity criteria of Horn et al. (2001) which indicates that a linear response is expected if the Wedderburn number W = (ρ2− ρ1)gh21/τsL is greater than a critical value Wn.s. = 3(h1−h2)/2piζ1h2. Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 81 F s (× 1 0 2 N m −1 ) (d) −5 0 5 x (km) D ep th (m ) (a) 0102030405060700 50 100 Width (km) 0 2 4 A re a (× 1 0 4 m 2 ) (b) 0 4 8 S1 S2 b 2 (k m ) (c) 0 2 4 Figure 4.2: (a) Width of cross-section versus depth and distance along the thalweg. (b,c) Bathymetric parameters for the TVC model: (b) cross-sectional area above and below the density interface and (c) cross-sectional width at the level of the interface. (d) Wind forcing on day 234.5, 1994, when the average wind speed at the moorings was 8.4 m s-1 blowing from the WSW. Positive wind forcing is in the direction of increasing x. In all panels, the thalweg origin, x= 0, is at Kenney Dam. Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 82 Julian Day (1994) W 180 200 220 240 260 280 0.1 1 10 100 Observed Wn.s. Figure 4.3: Wedderburn number for Knewstubb and Natalkuz Lakes for observations in 1994. For calculating the Wedderburn number, the average along-thalweg wind stress τs was low-pass filtered (Stevens and Lawrence, 1997) using a fifth-order Butterworth filter with T1/4 = 1.4 d cutoff. Also shown is the Wedderburn number below which non-linear steepening becomes pronounced, Wn.s. (Horn et al., 2001). Wn.s. is estimated for record segments, S1-S4, using the two-layer parameters (ρ2−ρ1 and h1) in Table 4.1, the average depth below the density interface h2 = 44 m−h1, and ζ1 = 0.1. For Knewstubb and Natalkuz Lakes, using a typical value of ζ1 = 0.1 gives W > Wn.s. for most of the time (∼ 2% of the record length, Figure 4.3). As will be shown, the calibrated value of ζ1 is an order of magnitude higher, which further reduces Wn.s. and justifies neglecting non-linearities and applying the TVC. Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 83 0 10 20 (a) KN: 61% O cc u rr en ce (% ) 0 10 20 (b) NA: 70% O cc u rr en ce (% ) N E S W N 0 10 20 (c) SK: 49% O cc u rr en ce (% ) Wind Direction O cc u rr en ce (% ) O cc u rr en ce (% ) O cc u rr en ce (% ) Figure 4.4: Histograms of wind direction at (a) KN, (b) NA, and (c) SK during 24-Jun to 12-Oct (days 175-285), 1994. Data shown are for wind speeds greater than 2 m s-1. Sum of occurrences relative to record length are noted on each panel. 4.3 Results 4.3.1 Field Observations Wind field From mid-summer to early fall of 1994, wind at the Knewstubb platform (KN) near Kenney Dam was similar to that at the Natalkuz platform (NA), 20 km to the south, and Skins Lake (SK), 80 km east of the dam. Wind direction histograms for days 175-285 show that wind was predominantly from the WSW at platforms KN and NA, veering to WNW at SK (Figure 4.4). For these dominant westerlies, the recorded temporal patterns were similar at the three stations. For the 110 d record, the W to E component of wind speed at platforms KN and NA had a linear correlation coefficient of 0.75 while correlation between NA and SK (65 km apart) was 0.66. Further, at all stations, spectral analysis of wind speed revealed a statistically significant narrow peak at a period of 1 d and a broad peak around 8 d (Figure 4.5a). Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 84 0.1110100 101 103 105 KN NA SKS p ec tr a l D en si ty ([ m s− 1 ]2 H z− 1 ) (a) 10−3 100 103 106 V 1 H 1 V 1 H 2 S1 S2 S3 S4 S p ec tr a l D en si ty (m 2 H z− 1 ) (b) 10−3 100 103 106 V 1 H 1 V 1 H 2 S1 S2 S3 S4 S p ec tr a l D en si ty (m 2 H z− 1 ) (c) 0.1110100 0 0.25 0.5 0.75 1 Period (d) C o h er en ce (d) Figure 4.5: (a) Power spectra of wind speed (solid lines). (b,c) Power spectra of the deflections of the 10°C isotherm (solid lines) at the (b) KN and (c) NA moorings for record segments S1-S4 (Table 4.1). For the record segments, periods of horizontal modes (V1Hn) are indicated by circles with each mode connected by a grey line. Periods of H1 and H2 are from the VDC with three-layers and higher modes (H3-H8, increasing from left to right) are from the VDC with two-layers. (d) Magnitude-squared coherence between the W to E wind component and the level of the 10°C isotherm at mooring KN for segment S3 of the record (19-Aug to 15-Sept). For clarity, the power spectra in panels (a-c) are vertically shifted. Spectra were generated after removal of the linear trend and were subsequently smoothed to improve confidence (Boegman et al., 2003). Dotted lines indicate the 90% confidence levels. Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 85 Baroclinic dynamics The baroclinic response to wind is affected by stratification and, in Knewstubb and Natalkuz Lakes, the stratification (the thickness of the epilimnion and metalimnion and the density differences) evolved considerably over the 110 d record (Figure 4.6). Prior to day ∼230 (18 Aug), warm isotherms were created at the surface and subsequently deepened particularly between day 190 (9 Jul) and 210 (29 Jul). This general intensification of the stratification led to thinning of the epilimnion to less than 10 m deep and led to a corresponding increase in the thickness of the metalimnion. After day ∼230, the epilimnion cooled and deepened thereby sharpening the metalimnion (Figure 4.6b,c). Baroclinic deflections within the metalimnion are best represented by the 10°C isotherm which persisted within the metalimnion throughout the record (Figure 4.6b,c). The metal- imnetic isotherms at KN appear to deflect in phase (Figure 4.6b). Visually, the deflections at KN were dominated by the first vertical mode, V1, which is the focus of this study. At NA, the metalimnetic deflections were in phase when the metalimnion was relatively thin (days 178-194 and days 256-280), but the metalimnion frequently varied in thickness when it was broad (days 195-250). In particular, the metalimnion contracted considerably during the strong storm of days 233-236.5 (Figure 4.6a,c). During this limited time, the second vertical mode affected the observed dynamics which will be addressed in Chapter 5. Power spectra of the 10°C isotherm over the entire record at both KN and NA (not shown) do not reveal statistically significant spectral peaks. Because changes in the periods of baroclinic modes over the record can blur the power spectra, we examined the spectra for each segment of the record, S1 to S4 (Figure 4.5b,c). No statistically significant peaks at the 90% confidence level appear in any of these spectra except during S3 at KN for which a relatively narrow peak is observed around 1.1 d (Figure 4.5b,c). At both moorings, broad peaks of low significance can be discerned around a period of 8 d during S3 and a period of 2-5 d during S4 (Figure 4.5b,c). To determine the excited baroclinic modes, modal periods estimated from the VDC model were overlaid on the power spectra of the 10°C isotherm level (Figure 4.5b,c). The periods of the few spectral peaks can be matched with some horizontal modes. The ∼1 d peak at KN during S3 can be interpreted as one or more of H6 to H8 V1 modes. The ∼8 d signal at both moorings during S3 matches the period of the V1H1 mode. The low- significance broad power distribution between 2-5 d at both moorings during S4 corresponds Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 86 −5 0 5 W in d S p ee d (m s− 1 ) (a) 6◦C 14◦C16◦C18◦C D ep th (m ) (b) 0 20 40 6◦C 14◦C16◦C18◦C D ep th (m ) Julain Day (1994) (c) 180 190 200 210 220 230 240 250 260 270 280 0 20 40 Figure 4.6: (a) Wind components along the W to E (black) and N to S (grey) directions at platform KN, 1994. (b,c) Isotherm levels at (b) KN and (c) NA. Isotherms are spaced at 2°C intervals. The 10°C isotherm is marked by a black line. Wind and water temperature data were low-pass filtered using a 5th order Butterworth filter with 12 h cutoff. Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 87 to the periods of the H2-H5 modes. However, this association between the spectral peaks and certain modal structures is ambiguous because the spectral peaks are also common to the power spectra of the wind, in particular, the 1 d and 8 d peaks (Figure 4.5a). Power is also observed in the spectra of the wind between 2 d and 5 d though at low significance. When the wind carried relatively high spectral power, the coherence between the wind and the metalimnion deflections was high. For example, at KN, during segment S3, the magnitude- squared coherence between the 10°C isotherm level and the W to E component of the wind was about 0.94 for periods greater than 5 d (Figure 4.5d). The coherence was also high (up to 0.88) for periods around 1.1 d (Figure 4.5d). This high coherence with the wind suggests that the deflections at these periods simply follow the temporal pattern of the wind forcing and makes the occurrence of free oscillations of the metalimnion unlikely as these oscillations would not be correlated with the wind pattern. Accordingly, the high coherence indicates strong damping of the response as will be discussed later. Identification of the excited modes from the spectral signature of the metalimnetic deflections becomes unfeasible. For identification of the excited modes and the associated damping, we resort to the TVC model. 4.3.2 Modeled Baroclinic Response In this section, we analyze the baroclinic response of Knewstubb and Natalkuz Lakes to wind forcing using the TVC model which simulates the evolution of individual baroclinic modes. The damping parameters of the excited modes are calibrated using the 1994 data at mooring KN and NA. The calibrated model is validated using the 2005 data at mooring KN. Decomposition of the along-thalweg wind forcing Fs (x, t) during the entire record, days 175-285, 1994, indicates the dominance of the n= 1 wind component followed by the n= 2 component (Figure 4.7). The relative weightCn of wind components for n= 1 has an average of 0.41 and, for n= 2, the averageCn=0.23 while the averageCn for higher modes (n> 2) is less than 0.07. Because of the high values ofCn for n= 1 and 2, the corresponding baroclinic modes H1 and H2 are expected to dominate the overall baroclinic response. To assess the relative importance of the baroclinic modes to the observed interface deflec- tions at moorings KN and NA, we examined the equilibrium deflection, η(n)2e , for n = 1−8 attained at the moorings under the average wind of days 233-236.5. During this period, strong winds of up to 11.4 m s-1 occurred from the WSW (255° at KN and 243° at NA) Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 88 Component, n C n 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 Figure 4.7: The average relative weights, Cn, of wind components (black line). Vertical bars indicate one standard deviation on each side of the average Cn. (Figure 4.6a). At mooring KN, the calculated η(n)2e indicates that the H1 mode dominated the observed metalimnetic deflections. The equilibrium deflection of the H1 mode was η(1)2e = −10.8 m. This downward deflection is consistent with strong WSW winds deepening the metalimnion down at Kenney Dam which is the downwind side and raising the metalimnion at the west end of Natalkuz Lake which is the upwind side. In comparison, the next mode in deflection amplitude was H2 with η(2)2e = +1.8 m. Higher modes had amplitudes less than 0.3 m. The observations at mooring KN show a depression of the 10°C isotherm of about 6.3 m between day 233 and 235.5 (Figure 4.6b) which is reasonably close to η(1)2e −η(2)2e = −9 m. At mooring NA, the picture was different with the H1 and H2 modes having comparable equilibrium deflections, η(1)2e = 2.3 m and η (2) 2e =−2.6 m. Higher modes had magnitudes less than 0.3 m. Because the analysis of η(n)2e indicated that the H1 and H2 modes were the most signifi- cant at the KN and NA moorings, only these two modes were included in the TVC simula- tions. The damping ratios, ζn, and initial conditions, Dn, for these two modes were calibrated by fitting the corresponding deflections simulated with the TVC to the fluctuations of the 10°C isotherm at the two moorings. Before fitting the deflections, the seasonal deepening of the metalimnion was eliminated by removing the average linear trend at moorings KN and NA (0.03 m d-1) from the 10°C isotherm depth. The detrended fluctuations of the 10°C isotherm at the KN and NA moorings were both used in calculating the root mean square Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 89 D efl ec ti o n (m ) (a) KN −8 −4 0 4 D efl ec ti o n (m ) Julian Day (1994) (b) NA 180 200 220 240 260 280 −8 −4 0 4 Figure 4.8: Comparison of the observed (thin, black) and simulated (thick, grey) deflections of the interface at the (a) KN and (b) NA moorings, 1994. error for deflections from trial simulations. From these simulations, the damping ratios which produced the least RMS error were ζ1 = 1.3 and ζ2 = 0.7. The corresponding RMS errors were 1.34 m and 1.56 m at moorings KN and NA, respectively, while the linear correlation coefficients were 0.81 and 0.59. Quali- tatively, there is good agreement between the best fit simulation and the observed fluctuations of the 10°C isotherm particularly at mooring KN (Figure 4.8). To test the sensitivity of the TVC towards the damping ratios, we determined the range of ζ1 and ζ2 which increased the RMS error by up to 5% of its minimum value (Shimizu and Imberger, 2008) and did not reduce the correlation by more than 0.05 from its value corresponding to the least RMS error. This range was 1.2 to 2.0 for ζ1 and 0.6 to 1.4 for ζ2. The narrower range for ζ1 reflected the overall dominance of the H1 mode at the two moorings. To characterize the composition of the interface deflections at the KN and NA moorings, we generated histograms for the individual simulated deflections of the H1 and H2 modes (Figure 4.9). The histograms show that each of the two modes was out of phase between mooring KN and NA. The two modes were also associated with out of phase motions at each mooring. At mooring KN, the histograms show that the H1 mode was responsible for the large amplitude depressions in the interface. In comparison, the H2 mode caused much Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 90 O cc u rr en ce (% ) (a) KN 0 25 50 H1 H2 O cc u rr en ce (% ) Deflection (m) (b) NA −10 −8 −6 −4 −2 0 2 4 0 25 50 Figure 4.9: Histograms of the simulated deflections for the H1 and H2 modes at the (a) KN and (b) NA moorings. The histograms are for the hourly deflections between day 175 and 285, 1994. smaller upward deflections. At mooring NA, the H1 mode led to small-amplitude upward deflections while the H2 mode led to downward deflections of slightly larger magnitudes (Figure 4.9). To validate the results of the TVC, we simulated the interface deflection at mooring KN in 2005 using modes H1 and H2 with the same calibrated ζn. Wind and water temperature observations at mooring KN are plotted in Figure 4.10. As in 1994, wind from the WSW was dominant as evident by the higher W to E wind speeds in Figure 4.10a. However, in 2005, the metalimnion was thinner than in 1994, particularly between day 201 and 250 (compare Figure 4.6b and 4.10b). This suggests that higher vertical modes were less significant in 2005 than in 1994. The TVC simulation extended over the length of the 2005 record, from day 201 to 280. The stratification parameters used in the simulation were h1 = 24 m and ρ2− ρ1 = 0.93 kg m-3. Because we lacked wind measurements at NA in 2005, wind from the KN mooring was used to force the entire thalweg. For comparison with the observed deflections at mooring KN, the linear trend, 0.09 m d-1, was removed from the observed level of 10°C isotherm. Overall, the TVC reproduced the observed deflections (Figure 4.10c). Discrepancy at the beginning of the simulation (till day ∼208) arises because we have assumed that the Nechako Reservoir was initially at rest. Compared to 1994, the RMS error in 2005 rose Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 91 −5 0 5 W in d S p ee d (m s− 1 ) (a) 6◦C 16◦C 14◦C D ep th (m ) (b) 0 20 40 Julian day (2005) D efl ec ti o n (m ) (c) 201 220 240 260 280 −8 −4 0 4 Figure 4.10: (a) Wind speed components along the W to E (black) and N to S (grey) directions at platform KN, 2005. (b) Isotherm levels at mooring KN. (c) Comparison of the simulated (thick, grey) and observed (thin, black) deflections of the interface at mooring KN. In panel (b), isotherms start from 6°C at the bottom and increase upwards at 2°C intervals. The 10°C isotherm is marked by a black line. Wind and water temperature data in panels (a) and (b) were low-pass filtered using a 5th order Butterworth filter with 12 h cutoff. to 1.50 m and the linear correlation coefficient dropped to 0.64. Nevertheless, the 2005 validation simulation still suggests the dominance of the H1 mode at mooring KN and indicates reasonable confidence in the calibrated damping ratio for this mode. Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 92 4.4 Discussion 4.4.1 Effect of Bathymetric Idealizations on Damping In applying the TVC to Knewstubb and Natalkuz Lakes, we have made idealizations pertain- ing to bathymetry, namely, the side branches were discarded and the lateral flow variability along the main branch was ignored. The effect of these bathymetric idealizations on the estimated damping of the baroclinic response was assessed by comparing the TVC results with estimates from the VDC model. As a first step, we examined the effect of bathymetric idealizations on the modal periods and structures. The bathymetric idealizations led to underestimation of the modal periods. While the TVC predicted a period of 5.3 d for the H1 mode and 2.5 d for the H2 mode, the VDC with the same stratification (h1 = 22.7 m and ρ2− ρ1 =1.01 kg m-3) predicted 6.8 d and 3.2 d, respectively. The VDC also estimated longer periods for the higher horizontal modes n= 3−8. Further, the bathymetric idealizations altered the modal structures, Z(n)2 , with greater effect on higher modes. For the H1 mode, the agreement of the along-thalweg modal structure between the TVC and the VDC is reasonable (Figure 4.11a). The VDC predicted that the nodal line is located within the Narrows. The TVC gave the same prediction with a difference of only -1.2 km. The TVC also reproduced the relative steepness of the interface. As predicted by the VDC, the steepness in the TVC was comparable within both Knewstubb and Natalkuz Lakes and was considerably less than the steepness within the Narrows. For the H2 mode, there was less agreement between the TVC and the VDC. While both models predicted that the two nodal lines are located closer to the Narrows than Kenney Dam or the west end of Natalkuz Lake, the TVC differed from the VDC by +5.6 km for the nodal line in Natalkuz Lake and -3.5 km for the nodal line in Knewstubb Lake. Hence, the intermediate region between the two nodal lines in the VDC is smaller by ∼9 km than predicted by the TVC. The next step is to examine the effect of the bathymetric idealizations on the baroclinic response to wind forcing. Without these idealizations, equation 4.4 and 4.6 are valid for simulating the response with ωn and Z (n) 2 given by the VDC (equation 4.1). The correspond- ing modal amplitudes An differ from those of the TVC reflecting the coupling between the modal transport in the upper layer and the two-dimensional wind field (Shimizu et al., 2007). Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 93 D efl ec ti o n (m ) (a) H1 −1 0 1 TVC VDC Distance (km) D efl ec ti o n (m ) (b) H2 01530456075 −1 0 1 Figure 4.11: Comparison of interface deflection, Z2, along the thalweg between the TVC and the VDC with N = 2 for the (a) H1 and (b) H2 modes. The interface deflection is normalized by its maximum absolute magnitude. Shading indicates approximate extent of the Narrows. The thalweg origin, x= 0, is at Kenney Dam. The changes in ωn, Z (n) 2 , and An for the realistic bathymetry can not exactly reproduce the calibrated response of the TVC. This is because the terms AnZ (n) 2 ,ωn √ ζ 2n −1 and ζnωn in the combined equations 4.4 and 4.6 can not be simultaneously matched with their corresponding TVC values. A rough estimate of the corresponding change in ζn can be obtained by using an empirical approach similar to that of Gomez-Giraldo et al. (2006). Different percent changes in AnZ (n) 2 from its TVC values were considered assuming each percent change was constant with time. The responses for the H1 and H2 modes were calculated using the VDC estimate of ωn and were separately fitted to the corresponding responses from the calibrated TVC results. The damping ratios which produced the least RMS error show that, for a relatively large change in AnZ (n) 2 , both modes H1 and H2 are strongly damped. Thus, the TVC results are not far removed from the baroclinic response in the actual bathymetry (Figure 4.12). 4.4.2 Coherence and Damping An important point that requires examination is the relation between the strength of damping and the coherence between wind forcing and metalimnetic deflections. At the Knewstubb and Natalkuz moorings, we have found that, for frequencies corresponding to high spectral Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 94 −40 −20 0 20 40 0 0.5 1 1.5 2 2.5 3 Change in AnZ (n) 2 (%) ζ n H1 H2 Figure 4.12: Damping ratio ζn from VDC versus the percent change in AnZ (n) 2 from its TVC value. The damping ratios give the least RMS error compared to the calibrated TVC response for mode H1 and H2. power in the 10°C isotherm level, the isotherm level was highly coherent with the wind. This was clearly shown for segment S3 of the record (Figure 4.5b,d). Previous investigators have inferred the persistence of unforced basin-scale metalimnetic oscillations from high spectral power in isotherm fluctuations at the periods of basin-scale baroclinic modes, combined with low coherence between wind forcing and the energetic isotherm fluctuations (e.g., Lake Lugano, Mysak et al., 1983; Kootenay Lake, Wiegand and Carmack, 1986). Maintaining that wind initiates these basin-scale metalimnetic oscillations, the low coherence indicates weak modal damping.4 However, while low coherence in these cases reflects weak damping, high coherence between the wind and metalimnetic deflections does not necessarily imply strong damping, in particular, if the wind is periodic. In this latter case, high coherence is typically attributed to steady oscillations forced by the wind periodicity as for example in Lake Banyoles (Roget et al., 1997), Lake Biwa (Saggio and Imberger, 1998), and Lake Tiberias (Lemckert et al., 2004). In these lakes, wind had a significant diurnal periodicity and the baroclinic response is inferred to resonate with this periodic forcing. Such interpretation of the high coherence is in line with the ideal picture of a linear oscillator with constant parameters (stratification in this case) being driven with periodic forcing. Regardless of the strength of damping, the coherence approaches unity provided that the record is long relative to the transients initiated by the onset of forcing 4See Appendix B for further discussion. Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 95 (De Silva, 2007). Thus, under periodic wind forcing, high coherence between the wind and metalimnetic deflections does not necessarily imply strong damping. The interpretation of coherence is further complicated by subjective distinction between high and low coherence. Saggio and Imberger (1998) referred to coherence at ∼0.6 as high. In addition, they tested significance against the 90% confidence level. In comparison, in Lake Lugano, Mysak et al. (1983) considered coherence at ∼0.6 to be low and to be statistically insignificant if less than the coherence at the 90% confidence level. In our case, however, characterization of coherence is not an issue as the coherence levels of interest reached up to 0.94. For Knewstubb and Natalkuz Lakes, we interpreted high spectral coherence (Figure 4.5d) as an indication of strong damping because the wind over the two lakes is aperiodic (Figure 4.6a and 4.10a) and hence triggers an aperiodic response. Aperiodicity is not well reflected in global power spectra (Figure 4.5a-c) which only show the spectral power distri- bution with frequency but hide the fluctuation of spectral power with time. To examine the spectral power in the time-frequency domain, we applied continuous wavelet transform to the wind forcing and the corresponding metalimnetic deflections at mooring KN (Figure 4.13). In applying the wavelet transform, a complex Morlet wavelet was used following similar studies (Stevens et al., 1996; Antenucci et al., 2000; Naithani et al., 2003; Vidal et al., 2007) and estimation of the cone of influence and confidence levels followed Torrence and Compo (1998) and Grinsted et al. (2004). The wavelet power spectrum (WPS) of the W to E wind speed at KN in 1994 shows regions of high spectral power at periods consistent with the global power spectrum in Figure 4.5a (e.g., 1 d period on day 249, 2-5 d period around day 272, and ∼ 8 d period around day 235). However, the WPS reveals that the wavelet power corresponding to these periods was not constant with time but fluctuated considerably which reflects the aperiodic pattern of wind forcing. Qualitative comparison of the WPS of W to E wind speed and 10°C isotherm level indicates similar regions of high power (Figure 4.13a,b). However, there are more extensive regions of significant wavelet power (above 95% confidence level) in the WPS of the 10°C isotherm level. To facilitate the comparison, we examined the wavelet transform coherence (WTC) which gives the distribution over time of the spectral correlation between the wind and the deflec- tions (Grinsted et al., 2004). The WTC also shows the change in the phase relationships in the time-period domain for regions with high wavelet coherence. The WTC spectrum of W to E wind component and metalimnetic deflections shows extensive regions of high coherence Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 96 (Figure 4.13c). These high-coherence regions fully overlap the high wavelet power in the WPS of the wind and the metalimnetic deflections (Figure 4.13a,b) indicating the absence of large free metalimnetic oscillations and implying strong damping. In contrast, if damping was weak, free oscillations of the V1H1 mode would be triggered and these oscillations would not be correlated with the wind pattern. This would lead to zones of low coherence in the WTC and highly variable phase relationships depending on the onset of wind storms relative to the phase of the free modal oscillations. Because of strong damping, the WTC spectrum shows only two phase relationships (indicated by arrows in Figure 4.13c). The first phase relationship (slanted right-pointing arrows) indicates a lagged deepening in the isotherm level following increase in wind speed from the west as expected for the H1 mode. This phase relationship dominates and persists with time for periods greater than ∼ 4 d. The second phase relationship (slightly slanted left-pointing arrows) is observed for shorter periods (. 4 d), in particular, towards the end of the records (day ∼270), and indicates a nearly in-phase rise in the isotherm level with increase in westerly winds. The shift in phase between short and long periods is due to the increasing contribution of the H2 mode at the short periods compared to the slower H1 mode, noting that the wavelet power carried by the significant regions at short periods is considerably lower than at longer periods (see color scale associated with Figure 4.13a,b). In summary, the regions of high coherence in the WTC spectrum overlapped the period of the V1H1 mode and the phase relationship associated with these regions was in agreement with the downward deflections obtained from the TVC (Figure 4.9a). The extensive regions of high coherence (for periods > 4 d) together with the aperiodic pattern of the wind suggest strong damping of the V1H1 mode. At shorter periods (. 2 d), regions of high coherence are not contiguous. This suggests that higher horizontal modes, which have short periods, are less damped than the V1H1 mode. 4.4.3 Comparison to Other Lakes The application of the TVC to Knewstubb and Natalkuz Lakes indicated that the H1 and H2 modes were strongly damped. The calibrated damping ratios, ζn, were 1.3 for H1 and 0.7 for H2. These damping ratios are considerably higher than those reported or deduced for other lakes for which ζ1 = κ1T1 (2pi)−1 is listed in Table 4.2 alongside with the e-folding time-scale Te = κ−11 (where κ1 is the decay constant of the H1 mode). The small values of Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 97 Figure 4.13: Wavelet power spectrum of (a) W to E wind speed and (b) 10°C isotherm level at mooring KN in 1994. Wavelet power is normalized by the variance. (c) Magnitude- squared coherence between the W to E wind speed and the 10°C isotherm level. In all panels, hatching masks the cone of influence (Torrence and Compo, 1998). Black contours bound zones exceeding the 95% confidence level calculated following Grinsted et al. (2004). In panel (c), arrows indicate the phase relationship between the wind and the 10°C isotherm level; in-phase deepening, _ deepening with quarter period lag, and in-phase rise in isotherm level with increasing W to E wind speed. Dashed white lines in panels (b,c) indicate the period of the V1H1 mode estimated from the VDC with N = 3. Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 98 ζ1 (< 0.2) in Table 4.2 explains why, in the corresponding lakes, identification of the excited modes from spectral analysis was successful and relatively straightforward. As a tangential point, the damping ratios in Table 4.2 suggests that the damping is higher in larger lakes (i.e. with greater length L). In the TVC and similar models with linear damping, the power dissipated Pd is equal to a constant fraction of the available kinetic energy Ek (Spigel, 1978; Shimizu and Im- berger, 2008). This fraction5 is given by PdE−1k = 4κ (Arneborg and Liljebladh, 2001). For Knewstubb and Natalkuz Lakes, the decay constants corresponding to the calibrated ζn are κ1 = 1.55 d-1 and κ2 = 1.76 d-1 which give PdE−1k ∼75 KW GJ-1. This value is considerably larger than PdE−1k = 4−26 KW GJ-1 for the lakes listed in Table 4.2 except for Sooke Lake Reservoir ( PdE−1k = 66 KW GJ −1). As PdE−1k in Nechako Reservoir is higher than in the other lakes, this suggests excessive damping in Nechako Reservoir. The only other example that we found in the literature that has a higher PdE−1k is Gullmar Fjord for which Arneborg and Liljebladh (2001) estimated PdE−1k = 240 KW GJ -1. However, this higher value of PdE−1k may be related to the density difference across the pycnocline in the Gullmar Fjord which was up to eight fold that across the metalimnion in Knewstubb and Natalkuz Lakes and the other lakes of Table 4.2. 4.4.4 Source of Damping Wüest et al. (2000) proposed an empirical relation which gives the e-folding time-scale Te as one day per 38 m of maximum water depth. For Knewstubb and Natalkuz Lakes, this relation gives ζ1 = 0.3 as opposed to the calibrated value of 1.3. Another formula by Spigel and Imberger (1980) gives the e-folding time-scale as Te = 1.9× 103ρ1h1h √ νT−1n (τse)−1 where τs is the wind stress, h is the average water depth (h =34 m for Knewstubb and Natalkuz Lakes), ν is the kinematic viscosity of water (ν = 1.308× 10−6m2s−1 for water at 10°C, Mays, 2009), and e is the bed roughness height. According to this formula, the e-folding time is inversely proportional to wind stress. For a given lake, the lower bound for Te corresponds to the maximum wind stress τs. Using the maximum value of τs at mooring NA ( 0.26N m−2 ) and e = 0.06 m as set by Spigel and Imberger (1980) for the north basin 5The ratio of the dissipated power Pd to the available kinetic energy Ek has a unit of inverse time similar to κ . However, I’ve retained units of power to energy in KW GJ-1 for the values of PdE−1k . The unit KW GJ -1 provides a quick idea of the dissipation rate when given the energy content in dominant modes which is typically reported in tens of joules (e.g., Antenucci and Imberger, 2001; Shimizu et al., 2007). Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 99 Table 4.2: Damping parameters in several basins, as reported in or estimated from the literature. Basin Length L (km) T1 (d) Te (d) ζ1 Time of observation Source Kootenay Lake 100 12 10.6 0.18 Oct - Nov Stevens et al. (1996) 5 5.7 0.14 Aug - Sept Untersee (Lake Zurich) 29 1.8 2.1 0.14 Sept Horn et al. (1986) Lake Tiberias 22 0.9 1.8 0.08 June Shimizu and Imberger (2008) Sooke Lake Reservoir 4 0.3 0.7 0.07 Sept, Oct Stevens (1999) Windermere north basin 7 0.8 3.1 0.04 June Mortimer (1952) 0.5-0.6 2.0-2.3 0.04 August, Sept Heaps and Ramsbottom (1966); Spigel and Imberger (1980) Wood Lake 7 0.5 2.0 0.04 August Wiegand and Chamberlain (1987) Baldegg Lake 4.5 0.3-0.4 > 7 < 0.009 July, August Lemmin (1987) of Windermere, the lower bound of the e-folding time for Knewstubb and Natalkuz Lakes is Te = 44 h. The corresponding damping ratio for the H1 mode is ζ1 = 0.46 which would be the upper bound as the average wind stress at mooring NA is only 0.04 N m-2. The difference between the above estimates (ζ1 < 0.5) and the calibrated (ζ1 = 1.3) damping ratios for the H1 mode as well as the much smaller damping ratios for other stratified water bodies (Table 4.2) suggest that there is excessive damping in Knewstubb and Natalkuz Lakes. Causes for this excessive damping may include the drag against submerged trees, the geometrical constriction within the Narrows, and the radiation of energy to the adjacent basins. The Nechako Reservoir was not logged prior to impoundment. Submerged trees pro- truding upwards in the water column act as roughness elements draining energy from the baroclinically-driven flow in the hypolimnion and, more importantly, on the side slopes intersecting the metalimnion. Examples of submerged trees are shown in Figure 4.14 from echo-sounding transects near Kenney Dam and across the west end of the Narrows. Overall, up to 65% of the bed of Knewstubb and Natalkuz Lakes may be covered by coniferous trees. For stands of coastal Douglas-fir, Weiskittel and Maguire (2006) estimated that the surface area of stems and branches range between 1 m2 and 5 m2 per square meter of ground. Carrying this to Nechako Reservoir, the submerged trees increase the surface area Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 100 Figure 4.14: Across-thalweg echo-sounding transects (a) near Kenney Dam (x= 4km) and (b) near the west end of the Narrows (x= 33km), October 2005. Note the trees protruding into the metalimnion (s), the small width of the Nechako Canyon below 40 m (a), and the irregular bathymetry within the Narrows (b). Image contrast has been improved to better reveal the trees and the metalimnion. of the bed by up to five times and can explain the four fold increase in the damping ratio over the value estimated from Wüest et al. (2000). We can back calculate the bed roughness height from the formula of Spigel and Imberger (1980) to obtain the calibrated damping ratio ζ1 = 1.3. Using τs = 0.04 N m-2, we obtain e = 1.11 m. This value of e lies just above the upper bound of the range given by Chow (1959) for roughness heights in natural river beds (0.03−0.91m). The case of Nechako Reservoir is not unique. Many reservoirs around the world were not cleared of trees prior to impoundment (Tenenbaum 2004). In such reservoirs, the impact of salvaging the underwater timber should be investigated carefully. Removal of the submerged trees could lead to more energetic metalimnetic oscillations with implications on selective withdrawal facilities (Stevens and Lawrence 1997; Anohin et al. 2006) and the general ecology of the reservoirs. In investigating the damping of the baroclinic response, the Narrows deserves attention due to its geometric constriction. The flow regime within the Narrows was characterized using the TVC model. Flow in the upper and lower layers (q1 and q2, respectively) were obtained from conservation of mass; q1 (x, t) = −q2 (x, t) = ´ x 0 b2 ∂η2 ∂ t dx. For the H1 mode, the maximum flow occurred at the node in the interface deflection (x= 32 km, Figure 4.11a) while, for the H2 mode, the maximum flow occurred at 53 km. Over the simulated period in 1994, the maximum flows for the H1 and H2 modes had average magnitudes of 470 m3s-1 Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 101 and 390 m3s-1, respectively. The superposition of the flow from the H1 and H2 modes gives the total flow. During the simulated period, the maximum total flow occurred between 24 km and 56 km, in particular, between 40 km and 56 km (73% of the time). During 80% of the simulated period, the maximum average velocity in the lower layer occurred at the most constricted section of the Narrows (x= 28km). The simulated velocity at this section had an average magnitude of 0.07 m/s. The maximum value reached 0.3 m s-1 on day 234.8 in response to strong winds that started on day 233.5. Combined with the sub- merged trees and the irregular bathymetry (Figure 4.14b), the elevated velocities simulated within the Narrows may significantly contribute to dissipation of the energy contained in the baroclinic flow. Field measurements are however needed to validate this hypothesis. The maximum composite Froude number within the Narrows was less than 0.35 for 99% of the time and only briefly reached 0.73 on day 234.8. Accordingly, the simulated flow was internally subcritical and internal dissipation within hydraulic jumps at the density interface is unlikely. Another cause for the strong damping is leakage of energy from the excited baroclinic modes in Knewstubb and Natalkuz Lakes to the northern arc of Nechako Reservoir. To study this energy leakage, field measurements would be needed within the connection between the two partitions. A good example of the needed measurements is the winter-time dataset col- lected by Arneborg and Liljebladh, 2001 at the entrance of the Gullmar Fjord. Alternatively, the coupling between Knewstubb and Natalkuz Lakes and the rest of the Nechako Reservoir can be accounted for by using the more sophisticated two-dimensional model of Shimizu and Imberger (2008). To achieve this coupling, surveying the bathymetry of the northern arc is needed as well as deployment of wind monitors and thermistor moorings there to provide a reservoir-wide picture of the baroclinic response. 4.4.5 Implications to CWRF The baroclinic response to wind can affect the temperature of outflow from water release facilities (Anohin et al. 2006; Stevens and Lawrence 1997) which, with the proposal to build a cold water release facility (CWRF) at Kenney Dam, provided the impetus for studying the baroclinic response of Nechako Reservoir. From 20 July to 20 August of each year, the CWRF would be required to release up to qw =170 m3s-1 at less than 10°C (Mitchell et al., 1995). To avoid releasing water warmer Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 102 than 10°C, the CWRF should be deeper than the 10°C isotherm level by more than the drawdown due to water withdrawal and the maximum depression due to the baroclinic response to wind. For point withdrawal and two-layer stratification, the drawdown of the density interface (corresponding to the 10°C isotherm) is given by 0.7gq 2/5 w ( 1− ρ1ρ2 )−1/5 (equation 10 in Wood, 2001) which equates to 14 m for the two-layer stratification of the TVC in 1994. For the baroclinic response to wind, the TVC predicts that the deflection at Kenney Dam is almost the same as that at mooring KN (Figure 4.11). Disregarding the time of occurrence relative to the operational period of the CWRF, the maximum downward deflection at the KN mooring in 1994 was ∼ 9 m. Downward interface deflections in 2005 were smaller (Figure 4.10c) but were compensated by the slightly deeper level of the interface. Based on these two years (1994 and 2005), the minimum depth to the center of the CWRF, hw, should be 45 m. Based on studies of summertime river temperatures downstream of Kenney Dam, the CWRF would be required to meet the target withdrawal temperature and discharge for a 200 year return period (Mitchell et al., 1995). For this return period, the records of wind and water temperature available in Knewstubb and Natalkuz Lakes are too short to conduct a reliable return period analysis of the wind-forced metalimnetic deflections. To deal with this shortage, long-term onsite monitoring of the wind is needed to infer the severity, frequency, and duration of wind storms. Additionally, long-term wind records from nearby stations that are correlated to the available onsite records can be transferred using regression analysis to represent the wind conditions on Knewstubb and Natalkuz Lakes. This data would then be used by the TVC model to determine maximum possible metalimnetic deflections and, hence, the safe withdrawal level to achieve the target withdrawal discharge and temperature for the 200 year return period. 4.5 Concluding Remarks We have investigated the baroclinic response to wind in Nechako Reservoir using field data and a forced one-dimensional two-layer variable cross-section (TVC) model. The wind field over the eastern basins of Nechako Reservoir, Knewstubb and Natalkuz Lakes, was found to be nearly-uniform under prevalent westerly winds. This approximate knowledge of the wind field was useful in determining which baroclinic modes were excited. Decomposition of the wind field using the TVC indicated that the first (H1) and second (H2) horizontal modes of Chapter 4. The Damped Baroclinic Response to Wind in a Reservoir 103 Knewstubb and Natalkuz Lakes were more strongly excited than other modes. Calibration against thermistor moorings confirmed that these two modes were sufficient to explain most of the variability in the metalimnetic deflections, in particular, at Kenney Dam. The calibration showed that the H1 and H2 modes were strongly damped. Using a two- dimensional model (VDC), modal periods were estimated without bathymetric idealizations. A simplified analysis using these modal periods showed that the strong damping predicted by the TVC was relevant to the real bathymetry. The strong damping was also inferred from field data. The temporal pattern of the wind speed and the isotherm deflections at the moorings were highly coherent. This has been shown to indicate strong damping since the wind did not have a dominant periodicity corresponding to the excited baroclinic modes. Additionally, the spectral peaks for the isotherm deflections were broad and insignificant. This was the case even though spectra were examined for segments of the record to accommodate changes in the modal periods with the evolution of the background stratification. Wavelet coherence provided a clearer picture of the time-dependent spectral correlation and phase relationship between the wind and the isotherm deflections. Damping in Nechako Reservoir was shown to be stronger than in many lakes. The excessive damping may be due to submerged trees that cover a large portion of the bed. The geometrical constriction in the Narrows between Knewstubb and Natalkuz Lakes leads to high velocities which may form a localized high dissipation zone. Additionally, radiation of the energy of the H1 and H2 modes out of Knewstubb and Natalkuz Lakes may contribute to the fast decay of the metalimnetic deflections. Although it idealizes the bathymetry, the TVC with the calibrated damping ratios re- produces the metalimnetic deflections at Kenney Dam with reasonable accuracy. This is valuable for assessing the performance of the proposed withdrawal facility provided longer wind records are collected. As a broader implication, the ability to simulate complex systems by means of relatively simple tools such as the TVC is valuable in studying the impacts of physical processes on the ecology and management of complex water bodies. 104 Chapter 5 The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 5.1 Introduction Wind forcing on stratified lakes energizes basin-scale baroclinic modes (Hutter et al., 2011). Excitation, amplification, and decay of baroclinic modes are influenced by lake bathymetry, the spatial and temporal characteristics of wind forcing, density stratification, and dissipa- tive processes (Lemmin et al., 2005; Shimizu et al., 2007). Identification of modes which dominate the baroclinic response to wind is important to understanding in-lake transport and its effect on water quality. The conventional method of identifying dominant modes is through spectral or time-frequency analysis (Mortimer, 1953; Hutter et al., 1983; Lemmin, 1987; Münnich et al., 1992; Lemmin et al., 2005); however, identification using this method becomes difficult when modal damping is strong. In this case, modal decomposition of the spatial structure of in-lake temperature and velocity fields at different times can provide useful information on the dominant modes (Shimizu et al., 2007). The baroclinic response to wind in single-basin lakes has been extensively studied by laboratory experiments in flat-bottom rectangular and circular tanks (e.g., Monismith, 1986; Stevens and Imberger, 1996; Boegman et al., 2005b; Wake et al., 2007) and field studies (e.g., Kootenay Lake, Wiegand and Carmack, 1986; Baldegg Lake, Lemmin, 1987; Wood Lake, Wiegand and Chamberlain, 1987; Lake Tiberias (Kinneret), Gomez-Giraldo et al., 2006). Studies on the baroclinic response of multi-basin lakes have been largely confined to lakes of relatively simple surface geometries (e.g., Lake Banyoles, Roget et al., 1997; Upper Lake Constance, Appt et al., 2004; Lake Geneva, Lemmin et al., 2005; Lake Biwa, Shimizu et al., 2007). Studies of the baroclinic response to wind in multi-arm multi-basin lakes and reservoirs have received less attention (e.g., Quesnel Lake, Laval et al., 2008). This study contributes to understanding the baroclinic response to wind in lakes of com- Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 105 plex bathymetry and surface geometry by characterizing the response in two connected lakes having different orientations and multiple arms. A previous investigation of the response in these two lakes focused on explaining the summertime thermistor records from two moor- ings (Chapter 4) which indicated that the first and second horizontal modes were strongly damped. Here, we focus on the response to a strong summer storm. The objective is to describe the basin-wide features of the baroclinic response and to identify the dominant horizontal and vertical baroclinic modes. Because strong damping was indicated in the previous investigation, conventional spectral decomposition of the baroclinic response is not useful and we resort to decomposition of the spatial structure. A three-dimensional hydrodynamic numerical model is used to simulate the baroclinic response to wind. After validation against the mooring data, the model output is spatially decomposed in terms of baroclinic mode structures obtained using modal analysis with a three-layer stratification and full representation of lake bathymetry. The field site and the tools used are described in the next section. In results, we present field observations with focus on data collected during a summer storm, the simulated re- sponse to the storm, and modal decomposition of the simulated response. In discussion, we describe the typical features of the baroclinic response under prevailing wind conditions and the contribution of dominant modes. 5.2 Methods 5.2.1 Field Site and Data Knewstubb and Natalkuz Lakes are the eastern basins of Nechako Reservoir which is located in the northern interior plateau of the Canadian Cordillera (Figure 5.1a). Knewstubb and Natalkuz Lakes have a combined surface area of about 210 km2, a maximum depth of 113 m, and an average depth of 31 m, respectively. Knewstubb Lake is composed of Knewstubb Arm, Big Bend Arm, and Knewstubb mid-reach which connects to the main basin of Natalkuz Lake through a constricted and elevated reach, the Narrows (Figure 5.1b). Besides the main basin, Natalkuz Lake includes Intata Reach, Chelaslie Arm, and Euchu Arm. Knewstubb and Natalkuz Lakes connect to the rest of Nechako Reservoir through Tetachuck River which restricts baroclinic exchange and through Intata Reach (Figure 5.1b). The baroclinic exchange through Intata Reach is not addressed in this study because previous Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 106 0 10 205 km Natalkuz Lake Knewstubb Mid-reach Knewstubb Arm Big Bend ArmIntata Reach Chelaslie Arm Euchu Arm Narrows Tetachuck River (a) (b) 1 2 3 4 5 Knewstubb Lake NA KN 126° W 53°30’ N SK Figure 5.1: (a) Map of Nechako Reservoir. (b) Map of the eastern part of Nechako Reservoir showing Knewstubb and Natalkuz Lakes, shaded in (a). Reservoir shoreline corresponds to the high water mark, at 853.4 m above mean sea level. The solid triangles indicate the location of the thermistor moorings and floating weather stations deployed in 1994. Locations of CTD profiles, also conducted in 1994, are indicated by numbered circles. investigation (Chapter 4) showed that the observed baroclinic motions in Knewstubb and Natalkuz Lakes can mostly be explained by considering the two lakes without including the exchange through the Intata Reach. Wind over Knewstubb and Natalkuz Lakes is predominantly from the WSW-WNW (Chapter 4). The wind field over the lakes is affected by synoptic scale winds which cause frequent summertime occurrences of W and NW winds over the east Pacific (Kendrew and Kerr, 1955). Westerlies were also found to prevail over Babine Lake which lies about 70 km to the north of Nechako Reservoir (Farmer, 1978). This study is primarily concerned with data from the summer of 1994 (Triton 1995). The dataset included two rafts, KN and NA (Figure 5.1b). Each raft was equipped with a thermistor mooring and a weather station. Thermistors were vertically spaced every 2 m and logged water temperature at 10 min intervals with an accuracy of 0.2°C. Each weather station included an RM Young 05103 wind monitor which was mounted at a height of 3 m above Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 107 the water surface. Vector-average wind speed and direction were logged on hourly intervals. Air temperature, relative humidity, and solar radiation were also measured at station NA. CTD casts were collected on 8 June 1994 (day 159) using a Seabird SBE 19 SEACAT profiler. Locations of the casts in Knewstubb and Natalkuz Lakes are marked in Figure 5.1b. During profiling, station KN and NA were not operational; however, wind speed and direc- tion were available from a land-based station at SK (Figure 5.1a). 5.2.2 Hydrodynamic Modeling Field observations in Knewstubb and Natalkuz Lakes were not sufficient for reliable identifi- cation of dominant baroclinic modes and calibration of their damping rates. For example, in Lake Tiberias (166 km2), data from six thermistor moorings were used to estimate the energy and damping rates of five selected baroclinic modes (Shimizu and Imberger, 2008). In con- trast, in Knewstubb and Natalkuz Lakes (210 km2), only two moorings were available and dominant modes were not known apriori as in Lake Tiberias. By using a three-dimensional hydrodynamic numerical model to expand the field data, we were able to infer the dominant modes and their damping rates. The Estuary Lake and Coastal Ocean Model (ELCOM) was selected for simulating the three-dimensional hydrodynamics in Nechako Reservoir. Without prior calibration, ELCOM has adequately reproduced the observed baroclinic response to wind in several stratified water bodies; e.g., Lake Tiberias (Hodges et al., 2000; Laval et al., 2003; Gomez-Giraldo et al., 2006); Lake Constance (Appt et al., 2004); Lake Biwa (Shimizu et al., 2007); and Valle de Bravo Reservoir (Okely et al., 2010). ELCOM solves the Reynolds-averaged Navier-Stokes and transport equations over a structured numerical grid (see Hodges et al., 2000, for a detailed description of ELCOM). For Knewstubb and Natalkuz Lakes, we constructed a non-uniform horizontal grid with relatively high resolution in Knewstubb Arm and within the Narrows. In these regions, 50×50 m cells were used to adequately resolve the narrow submerged Nechako Canyon that runs through the Narrows and along Knewstubb Arm (Figure 5.2)1. The cell width increased gradually (∼5% per grid cell) to a maximum of 200 m in Knewstubb mid-reach and 500 m in Chelaslie Arm. In the vertical, layer thickness was constant at 0.5 m for the top 30 m, increasing gradually by ∼10% to a maximum of 4.5 m below a depth of 70 m. 1Also, see Figure 4.14 in Chapter 4. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 108 320 330 340 350 360 370 380 390 5920 5930 5940 Easting (km) N or th in g (km ) NA KN Curtain C2 Distance (km) D ep th (m ) 0102030405060 100 50 0 Curtain C1 Distance (km) D ep th (m ) 0 5 10 15 100 50 0 Depth (m) 0 20 40 60 80 100 Figure 5.2: The numerical bathymetry of Knewstubb and Natalkuz Lakes which is used to ELCOM. The insets show the longitudinal profiles along the thalweg of Knewstubb and Big Bend arms (C1) and along the thalweg of Knewstubb and Natalkuz Lakes (C2). The origin of the profiles is at the north east end of Knewstubb Arm from which distance is marked at 10 km intervals (hollow circles). The locations of mooring KN and NA are indicated by hollow triangles. Bathymetry data from pre-impoundment areal surveys and from echosounding in 2005 were merged and interpolated to the centers of grid cells based on a Delaunay triangulation. Within the submerged Nechako Canyon in Knewstubb Arm, the numerical grid was manu- ally smoothed to allow a better connectivity between the cells covering the narrow canyon. In Chelaslie Arm, bottom depth was linearly interpolated along the thalweg of the western side where bathymetry data was lacking. The bathymetry along the thalweg of Knewstubb and Big Bend arms (curtain C1) and along the thalweg of Knewstubb and Natalkuz Lakes (curtain C2) are shown in Figure 5.2. The two curtains C1 and C2 are used to illustrate the isotherms simulated using ELCOM. The period covering 13-28 August 1994 (day 225.5 to 240) was simulated with a time step of 45 s. Only temperature was assumed to control the density field as specific con- ductance was low (averaging ∼60 mS cm-1) with little vertical structure. Simulations were initialized with a non-uniform thermal field synthesized by averaging temperature profiles Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 109 5 10 15 200 10 20 30 40 50 60 Temperature (°C) D ep th (m ) (a) KN NA 998.6 999.3 1000 Density (kg m−3) (b) Background Idealization Figure 5.3: (a) Average temperature profiles on day 225 from observations at mooring KN and NA. (b) Background stratification from the mooring data during days 225.5-240 and a three-layer idealization of the stratification. from day 225 at each of the two moorings (Figure 5.3a), applying the average KN profile to Knewstubb and Big Bend Arms, applying the average NA profile to Natalkuz Lake, and interpolating the average profiles at the two moorings over Knewstubb mid-reach using inverse weighted distance. At the surface, wind measured at KN was applied to Knewstubb and Big Bend arms while wind measured at NA was applied to the rest of the domain. Heat fluxes through the water-atmosphere interface were excluded from the simulation. Because the simulated duration was relatively short, surface heat fluxes had secondary effect on the metalimnetic response to wind which was confirmed by sensitivity simulations not reported here. On the bottom and sides of the domain, a free-slip condition was selected. Sensitivity simulations showed little difference in the results for free- and no-slip conditions. For simplicity, stream inflows and outflows were not included in the simulations. Be- tween day 225.5 and 240, the recorded water level dropped by only 0.09 m. Also, over the same period, the inflows from the major tributaries feeding Knewstubb and Natalkuz Lakes (including Tetachuck River, Figure 5.1b) amounted to an average of about 75 m3s-1. The inflows ranged in temperature between 12°C and 16°C. Based on the background thermal stratification at mooring KN, the volume contained between these isotherms within Knew- Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 110 stubb and Natalkuz Lakes was 0.8 km3. This gives a residence time in our domain of 124 d for the given stream flow. Over the short duration of the simulations (14.5 d), the stream inflows would have secondary effects on the baroclinic response to wind. In Intata Reach, fluxes between Natalkuz Lake and the northern basins of Nechako Reservoir were not simulated; i.e. Intata Reach was modeled with a no-flux boundary. As will be shown by the results of ELCOM, the domain including Knewstubb and Natalkuz Lakes reproduces most of the observed metalimnetic response to wind without the need to include either the stream inflows or the exchange fluxes with the rest of the reservoir. 5.2.3 Modal Analysis Layered modal analysis was used to obtain the undamped periods Tr of baroclinic modes and the corresponding spatial structures of the deflection Zi of density interface i= 1,2 . . .N where i = 1 denotes the free surface and N is the maximum number of layers used in the modal analysis. Because changes in metalimnetic thickness are observed at mooring KN and NA, three density layers with two internal interfaces (i.e. N = 3) were needed to capture the changes in the thickness of the metalimnion. The three-layer, two-dimensional, variable- depth, complete model (ThVDC) was used for the modal analysis (Salvadé et al., 1992). A summary of the ThVDC is given here for completeness and to introduce the notation used later in the modal decomposition. The ThVDC solves the eigenvalue problem, −ω2r Zi,r = N ∑ j=1 ∇ · {( N ∑ k=1 c jkdikhk ) ∇Z j,r } 1≤ i≤ N (5.1) where r ∈ Z+ is a mode index increasing with decrease in period Tr and with r = 1 denoting the mode with the longest period (V2H1); ωr = 2piT−1r is the undamped frequency of mode r with modal structures Zi,r (x,y); hk (x,y) is the at-rest thickness of layer k (top layer corresponds to k = 1); N ≤ 3 is the number of layers in a given partition of the domain (determined by layer thicknesses relative to local depth); ∇ = ( ∂ ∂x , ∂ ∂y ) is the del operator in the horizontal Cartesian coordinates x and y; and c jk and dik are respectively the reduced Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 111 gravity and a coefficient given by, c jk = gρ jρ−1k j = 1 g ( ρ j−ρ j−1 ) ρ−1k k ≥ j > 1 0 k < j and dik = 1 i< k+10 else (5.2) in which ρk is the density of layer k and g is the gravitational acceleration. Equation 5.1 is subject to no flow at the shoreline and conservation of flow between domain partitions which contain different numbers of layers. The flow associated with mode r is described by, ( Qk,r,Pk,r ) =− 1 ωr k ∑ j=1 c jkhk∇Z j,r 1≤ k ≤ N (5.3) where Qk (x,y) and Pk (x,y) are the spatial structures of depth-integrated flow in layer k in the x and y directions, respectively. To apply the ThVDC to Knewstubb and Natalkuz Lakes, the background stratification be- tween day 255.5 and 240 was idealized into three layers (Figure 5.3b). The layer thicknesses were taken as h1 = h2 = 12 m (in domain partitions with k < N, otherwise hk is variable) and the density differences between the three layers were taken as ρ2− ρ1 = ρ3− ρ2 = 0.67 kg m-3. The upper (i= 2) and lower (i= 3) interfaces between the layers correspond to the approximate levels of the 16°C and 7°C isotherms on day 225, respectively. For this reason, these two isotherms were used to compare and analyze the metalimnetic deflections in the field data and ELCOM results. In applying the ThVDC to Knewstubb and Natalkuz Lakes, equation 5.1 was discretized on a uniform Arakawa C-grid (Lemmin et al., 2005; Shimizu et al., 2007) with 55×55 m cells. The relatively fine grid resolution was needed to resolve the bathymetry within the Narrows. Partitions of the domain with less than three layers (N = 1 and 2) were included as they occupied 38% of the total surface area and 12% of the total volume. The baroclinic modes obtained from the ThVDC were classified using the convention for flat-bottom basins, VmHn where m = 1,2 and n are the number of vertical and horizontal nodal lines, respectively. From the many modes of Knewstubb and Natalkuz Lakes that can be excited by wind, seven modes will be shown to dominate the baroclinic response. For these dominant modes, the spatial structures Zi,r of the density interfaces and the corre- sponding VmHn classification are shown in Figure 5.4. For the low modes of the Knewstubb and Natalkuz Lakes, the classification was straightforward (e.g., V1H1 and V2H1, Figure Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 112 Figure 5.4: Spatial structures Zi,r of the upper (i= 2) and lower (i= 3) density interfaces for the dominant baroclinic modes. 5.4a,b). However, for higher modes, the classification was more subjective as the density interfaces did not always have the same number of horizontal nodal lines and these were not vertically aligned (e.g., V1H5, Figure 5.4g). 5.2.4 Modal Decomposition To identify the modes which dominate the baroclinic response to wind, the ELCOM-simulated response was decomposed in terms of the modes obtained from the ThVDC model. The decomposition followed a procedure simplified from Shimizu et al. (2007) for non-rotating modes. The vertical eigenvalue problem (Münnich et al., 1992) for the KN temperature profile of day 225 gives phase speeds of 0.26 and 0.12 m s-1 for the V1 and V2 modal structures, respectively. The corresponding internal Rossby radii are 2.3 km and 1.0 km, respectively. The average width of curtain C2 at the level of the upper interface is 1.7 km which suggests that the effect of rotation on V1 modes is secondary while the effect on V2 modes may be more significant (Hutter et al., 2011). Since the basins are long and narrow and since we are primarily concerned with the relative contribution of modes to the Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 113 baroclinic response to wind, we neglect rotation in ELCOM and in modal decomposition of its results. Compared to the rotating modal analysis of Lemmin et al. (2005) and Shimizu et al. (2007), the eigenvalue problem solved by the non-rotating ThVDC (equation 5.1) is 1/6 in size since the two components of flow are eliminated from the problem and only real coefficients are considered compared to complex coefficients in the rotating analysis. With the same computational resources, this allows for significant improvement in the grid resolution2 and thus better representation of constricted reaches such as the Narrows. Assuming no rotation, the modal decomposition of Shimizu et al. (2007) can be sim- plified as follows. Interface deflections and layer flows associated with mode r are given by, ηi,r = arZi,re−1r (5.4) (qi,r, pi,r) = br (Qi,r,Pi,r)e−1r (5.5) where ηi,r (x,y, t) is the displacement of interface i associated with mode r (ηi,r > 0 indicates upward deflection); qi,r (x,y, t) and pi,r (x,y, t) are the mode-r depth-integrated flows in layer i along the x and y directions, respectively; t is time; ar (t) = 2 ¨ ( Σ3i=1ρiciiZi,rηi ) dxdy; (5.6) br (t) = 2 ¨ ( Σ3i=1ρi (Qi,rqi+Pi,rpi)h −1 i ) dxdy; (5.7) and er = ¨ Σ3i=1ρi ( ciiZ2i,r+ ( Q2i,r+P 2 i,r ) h−1i ) dxdy (5.8) in which ηi (x,y, t) is the simulated displacement of interface i and qi (x,y, t) and pi (x,y, t) are the simulated depth-integrated flows in layer i along the x and y directions, respectively. The total energy Et,r contained in mode r is given by the summation of the potential Ep,r and kinetic Ek,r energies, Ep,r = 0.25a2re −1 r Ek,r = 0.25b2re −1 r , (5.9) and the work done by wind on mode r is given by, Wr (t) = ¨ (q1,rτx+ p1,rτy)h−11 dxdy (5.10) 2See Appendix C for the resolution of numerical grids used in previous studies. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 114 where τx (x,y, t) and τy (x,y, t) are the components of surface wind stress along the x and y directions, respectively. For mode r, the difference between the work done by wind Wr (t) and the rate of change of the modal energy Et,r (t) gives the rate of modal dissipation Dr (t). For linear damping, DrE−1k,r = 4κr (Arneborg and Liljebladh, 2001) in which κr is the decay rate of the response of mode r to wind. The dimensionless damping ratio ζr = κrω−1r characterizes the strength of damping as underdamped, 0 < ζr < 1, or overdamped, ζr > 1. To decompose the baroclinic response to wind in Knewstubb and Natalkuz Lakes (i.e., determine ar and br), the simulated thermal and velocity fields were output from ELCOM at 3 h intervals. The depths of the 7°C and 16°C isotherms were determined by linear interpolation in each column of the ELCOM grid. The spatially-averaged depths of the 7°C and 16°C isotherms were calculated at each time step (dashed lines in Figure 5.5g,h) and were subtracted from the respective isotherm depths in each grid column to obtain the interface deflections η2 and η3. The Layer flows qi and pi were determined by vertical integration of the x and y velocity components, respectively, from the bottom to the 7°C isotherm, between the 7°C and 16°C isotherms, and from the 16°C isotherm to the surface. Isotherm deflections and layer flows were interpolated to the ThVDC grid using bilinear interpolation. To determine ar (t) and br (t) (equation 5.6 and 5.7), the spatial structures Zi,r for all modes were normalized to er = 1 GJ (equation 5.8). With this common normalization, modal contributions to interface deflections ηi,r (x,y, t) can be compared by examining a single snapshot of ηi,r (at a given time) and the time series of ar as will be explained later. Using the calculated ar and br, total modal energies Et,r = Ep,r+Ek,r were obtained from equation 5.9. Only modes with periods Tr > 0.8 d were considered as these modes contained most of the total energy of the simulated response. Finally, to characterize the typical features of the baroclinic response of Knewstubb and Natalkuz Lakes, we examined the metalimnetic deflections 0.5(η2,r+η3,r), changes in metalimnetic thickness η2,r−η3,r, and damping ratios ζr for the modes which dominated the response (i.e. with the largest Et,r). In estimating decay rates κr to obtain ζr, we excluded segments of the simulation with energy dissipation rates Dr < 5% of its maximum value. This eliminated unreliable estimates of ζr when energy dissipation was negligible (∼ 7 to 10 d of the two-week simulation). Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 115 5.3 Results 5.3.1 Observed Baroclinic Response Response to summer storm We focus on the baroclinic response to wind during two weeks in August 1994 (days 225.5- 240) when the strongest storm in the summer of 1994 occurred. During the storm (days 233-236.5), wind at raft KN and NA blew predominantly from the WSW (Figure 5.5a,b) which is the prevailing wind direction in the summer (Chapter 4). Prior to the storm (days 225.5-233), there was relatively low wind activity (Figure 5.5c). On day 233.2, wind speed increased and 14 h later reached a peak (P1) of 8 m s-1 at raft KN after which wind speed decreased to 1 m s-1. Wind speed increased again on day 234.3, quickly reaching 9 m s-1 by day 234.5, and peaking at 9.3 m s-1 on day 234.8 (peak P2). Wind speed decreased to 3 m s-1 by day 235.2 but increased once again to 6 m s-1 by day 235.4, peaking at 7 m s-1 on day 235.6 (peak P3), and then weakening slowly over the next 24 h. Following the storm (days 236.5-240), wind activity was low averaging 0.6±1.3 m s-1. At raft NA, wind revealed a similar pattern to KN wind but with magnitudes about 20% higher and a slight delay in peak P3 (day 235.8) (Figure 5.5d). Corresponding to wind forcing during days 225.5-240, variations in the metalimnion depth and thickness are observed at mooring KN and NA (Figure 5.5e,f). At mooring KN, the 7°C to 16°C isotherms appear to deflect in-phase which indicates the dominance of V1 modes in the metalimnetic deflections at KN. The prominent feature of the V1 modes is the ∼ 7 m deepening of metalimnetic isotherms between day 233.2 and 235.4 and the subsequent gradual rebound. Superimposed on this long time-scale (∼7 d) deflection, smaller upward deflections (1-3 m) of shorter time scales (∼ 1 d) are observed at or following the wind peaks P1-P3. The influence of V2 modes at KN is evident from the separation of the 7°C and 16°C isotherms which varied between 6 m and 18 m with a standard deviation of 2.5 m about the average thickness. At mooring NA, the upper metalimnetic isotherms (12-16°C) during days 233-234.8 showed deepening similar to the metalimnetic isotherms at mooring KN albeit with a more rapid and abrupt rebound and revealing no upward deflections corresponding to the wind peaks P1 and P2. In contrast to mooring KN, successive contractions and expansions in the metalimnion can be visually discerned at NA. The separation between the 7°C and 16°C isotherms varied between 2 m and 17 m with a standard deviation of 3.6 m. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 116 W in d Sp ee d (m s− 1 ) (c) P1 P2 P3 0 5 10 (d) W in d D ire ct io n (a) N E S W N (b) D ep th (m ) (e) 0 10 20 30 (f) Julian Day in 1994 D ep th (m ) (g) S1 S2 S3 226 228 230 232 234 236 238 240 0 10 20 30 Julian Day in 1994 (h) 226 228 230 232 234 236 238 240 Figure 5.5: Wind forcing and thermocline response between day 225.5 and 240 in 1994 at mooring KN (left panels) and NA (right panels): (a,b) observed wind direction, (c,d) observed wind speed, (e,f) observed isotherm depths, and (g,h) simulated isotherm depths. Isotherms in (e-h) are displayed at 1°C intervals with solid black lines denoting the 7°C and 16°C isotherms. In panels (c-h), the dotted vertical lines, P1, P2, and P3, indicate features which are described in the text. Note that peak P3 in (c) is slightly shifted to the left. Arrows in (g,h) denote the times of snapshots S1-S3 at which the simulated metalimnetic response is investigated. Dashed lines indicate the domain-wide average depths of the simulated 7°C and 16°C isotherms. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 117 The larger range and standard deviation of variations in the metalimnetic thickness at NA indicates a stronger V2 signal at mooring NA than at KN. Spatial structure We now look at temperature profiles from the June 1994 CTD casts. Because wind direction during CTD profiling was similar to that during the storm of days 233-236.5 (Figure 5.6a), it is expected that the CTD temperature profiles reflect a spatial structure similar to the baroclinic response during the storm of days 233-236.5. The wind record at station SK started on the same day as the profiling (day 159) and shows westerly wind of ∼6 m s-1 (Figure 5.6a,b). For the next three days (days 160-162), the SK record shows diurnal westerly wind and, given this strong diurnal periodicity3, the condition prior to profiling on day 159 could be similar. Temperature profiles from the downcasts show differences in the metalimnion depth and thickness depending on the cast location (Figure 5.6c). Within Knewstubb Arm, Knewstubb mid-reach, and the eastern part of Natalkuz Lake (casts 1 to 3), the top of the metalimnion was deeper by 3-18 m than in the western part of Natalkuz Lake (casts 4 and 5). Also, at the center of Knewstubb mid-reach (cast 2), the metalimnion was only 3 m thick compared to a thickness of ∼ 20 m in Knewstubb and Chelaslie Arms (casts 1 and 5). As will be shown, these features that are revealed by the June 1994 CTD casts are also seen in the ELCOM response to the storm of days 233-236.5. In the next section, we validate the ELCOM response and then use this response to understand the observed dynamics of the metalimnion. 5.3.2 Simulated Baroclinic Response Model validation At mooring KN, ELCOM qualitatively reproduced the observed variations in metalimnetic depth and thickness (Figure 5.5g). For the 16°C isotherm, there is good agreement between the simulated and observed depth as indicated by a correlation coefficient of 0.97 and a root- mean square (RMS) error of 0.8 m. This RMS error amounts to 8% of the observed depth range and to 0.3 of the standard deviation of observed depth variations. Similarly, for the 3Also, see Figure 4.5a in Chapter 4. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 118 W in d Sp ee d (m s− 1 ) (b) 0 3 6 Julian Day in 1994 W in d D ire ct io n (a) 15 9 16 0 16 1 16 2 16 3 E S W N E 4 5 6 7 8 9 10 110 10 20 30 40 50 60 D ep th (m ) Temperature (°C) (c) 1 2 3 4 5 Figure 5.6: (a) Wind direction and (b) speed during CTD profiling on 8 June 1994 (day 159). (c) Temperature profiles from the CTD downcasts. Locations of the casts in Knewstubb and Natalkuz Lakes are marked on Figure 5.1. Note the colder surface temperature and shallower metalimnion in the west side of Natalkuz Lake (casts 4 and 5). In (a,b), grey hatching indicates the duration of profiling. Wind conditions are from station SK. metalimnetic thickness (the separation between the 7°C and 16°C isotherms), there is good agreement between the observed and simulated thickness at mooring KN. The correlation coefficient is 0.9 and the RMS error is 1.5 m which equates to 13% of observed range and 0.6 of the standard deviation of observed variations. Overall, the ratio of RMS error to range at mooring KN (8-13%) is comparable to that obtained by Okely et al. (2010) using ELCOM for the Valle de Bravo Reservoir. Further, the small ratios of RMS errors to standard deviations (0.3 and 0.6) indicate that errors from the model are smaller than the observed variations in metalimnetic depth and thickness (Jin et al., 2000). At mooring NA, the agreement between the observed and simulated metalimnetic vari- ations in depth and thickness is slightly worse than at KN (Figure 5.5e-h). The RMS error in the simulated depth of the 16°C isotherm is 1.4 m and the RMS error in metalimnetic thickness is 2.4 m. Nevertheless, the simulation reasonably reproduces the trends in the NA response as indicated by correlation coefficients of 0.94 and 0.85 for the depth of the 16°C isotherm and the metalimnetic thickness, respectively. The good agreement in trend between the simulated metalimnion depth and thickness and the observed series indicates that the dominant processes are captured in ELCOM. The Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 119 simulated response to wind adequately characterizes the basin-wide baroclinic response of Knewstubb and Natalkuz Lakes and will be used to identify the dominant baroclinic modes through modal decomposition. Simulated spatial structure Snapshots of the simulated isotherms were generated along curtains C1 and C2. The times of snapshots S1, S2, and S3 correspond to the wind peak P2, between the P2 and P3 peaks, and following the storm, just after the metalimnion at NA expanded, respectively (Figure 5.5). In Natalkuz Lake, curtain C2 shows that, during the storm of days 233-236.5, the 16°C isotherm was consistently sloped downward from Euchu Arm towards the Narrows (38-62 km in S1 and S2, Figure 5.7a,c). The 7°C isotherm was less tilted than the overlying isotherms which led to compression of the metalimnion at the downslope end (38 km) and expansion at the upslope end (62 km). East of the Narrows, within Knewstubb mid-reach, a similar pattern of metalimnetic response is observed though with smaller amplitudes (9-28 km in S1 and S2, Figure 5.7a,c). At the most constricted section of the Narrows (28-30 km in S1 and S2, Figure 5.7a,c) isotherm levels changed abruptly because of the upward metalimnetic deflection and expan- sion east of the Narrows and downward deflection and compression west of the Narrows. In Knewstubb and Big Bend arms, curtain C1 shows opposite tilts of the metalimnion between snapshots S1 and S2 (Figure 5.7b,d). During the wind peak P2 (corresponding to S1), the metalimnion was titled downward towards Big Bend Arm (Figure 5.7b) and, when the wind relaxed between peaks P2 and P3 (corresponding to S2), the metalimnion reversed its tilt (Figure 5.7d). After the wind storm (day 239.8 corresponding to S3), all isotherms relaxed towards horizontal equilibrium levels (Figure 5.7e,f) which were deeper than the pre-storm levels (dashed lines in Figure 5.5g,h). This deepening was the result of mixing driven by wind stirring and density instabilities in upwelling fronts of epilimnetic and upper metalimnetic isotherms (not shown). The compression of the metalimnion at the eastern side of Natalkuz Lake (38 km) in S1 and S2 was reversed and the metalimnion expanded (Figure 5.7e) as observed at mooring NA (Figure 5.5f). As will be explained below, this expansion was caused by flow within the metalimnion from the previously dilated metalimnion in the western part of Natalkuz Lake. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 120 D ep th (m ) (a) 40 20 0 D ep th (m ) (c) 40 20 0 Distance (km) D ep th (m ) (e) 102030405060 60 40 20 0 (b) (d) Distance (km) (f) 0 10 Figure 5.7: Snapshots of simulated isotherms along curtains (a,c,e) C1 and (b,d,f) C2. Times of snapshots (a,b) S1, (c,d) S2, and (e,f) S3 are marked on Figure 5.5g,h. Solid black lines denote the 7°C and 16°C isotherms at the time of the snapshot. Dotted lines indicate the domain-wide average depths of the 7°C and 16°C isotherms. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 121 The ELCOM-simulated spatial structure described above reveals that, during the storm, the component of the WSW wind acting along Knewstubb and Big Bend Arms pushed warm surface water downwind towards Big Bend Arm. As a result, the metalimnion deflected downward in Big Bend Arm and upward in Knewstubb Arm (Figure 5.7b). Similarly, in Knewstubb mid-reach, the WSW wind drove surface water northeastward raising the metalimnion east of the Narrows and deepening it in the northeast part of the mid-reach (Figure 5.7a,c). A similar response is observed in Natalkuz Lake where the WSW wind pushed surface water out of the lake towards the Narrows, raised the metalimnion in the west part of Natalkuz Lake, and deepened the metalimnion in the east part of the lake. As the metalimnion deflected vertically, it changed in thickness. These changes were most pronounced in Natalkuz Lake, where under the action of the WSW wind, the metal- imnion compressed west of the Narrows and expanded at the opposite side, in Euchu and Chelaslie Arms (Figure 5.7a,c). As observed at mooring NA, the metalimnetic thickness decreased to ∼ 4 m between day 235.0 and 238.5 (Figure 5.5f). In laboratory experi- ments by Kranenburg (1985) in a flat-bottom flume, the compressed metalimnetic thick- ness at the downwind endwall was found to be about 0.2 of the upper layer thickness h1. For Natalkuz Lake, the value of this compressed metalimnetic thickness equates to 2.4 m (h1 = 12 m) which is in good agreement with the observed thickness (∼ 4 m), given that the wind stress over Natalkuz Lake was not steady and, more importantly, that the Narrows permitted through flow as opposed to the no through-flow at the downwind endwall in the experiments of Kranenburg (1985). After the storm (day 239), the metalimnion at mooring NA expanded (Figure 5.5f). The ELCOM simulation shows that this expansion was part of re-adjustment of isotherms (Figure 5.7e) whereby the horizontal density gradients associated with metalimnetic tilt and upwelled fronts relaxed driving west to east metalimnetic flow in Natalkuz Lake. The W to E advection of metalimnetic water masses caused a so-called metalimnetic intrusion (Appt et al., 2004), a relatively rapid divergence of metalimnetic isotherms progressing along the thalweg as registered at mooring NA. Similar observations of adjustment of metalimnetic thickness after removal of wind forcing were made in laboratory flumes (Monismith, 1986), the sub-basin of Upper Lake Constance (Appt et al., 2004), and the single-basin Mundaring Weir Reservoir (Okely and Imberger, 2007). Observational support for the simulated response can be seen in the temperature profiles from the June 1994 CTD casts. The similarity in features should, however, be viewed as Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 122 suggestive in light of the weaker stratification in June and the uncertainty in determining wind conditions prior to profiling (Figure 5.6a,b). Comparison of the ELCOM response during the storm of days 233-236.5 to the temperature profiles from the CTD casts shows that the simulated compression of the metalimnion (9-15 km in S1 and S2, Figure 5.7a,c) is in agreement with the sharp metalimnion observed at the center of Knewstubb mid-reach (cast 2, Figure 5.6c). Further, the simulated metalimnetic expansion (0-9 km in S1, Figure 5.7b) agrees with the diffuse metalimnion observed in Knewstubb Arm (cast 1, Figure 5.6c) and the simulated metalimnetic upwelling at the west end of Natalkuz Lake (>58 km in S2, Figure 5.7c) agrees with the observed broad metalimnion in Chelaslie Arm which extended to the surface and led to cooler surface temperatures (cast 6, Figure 5.6c). In this section, we have described the observed and simulated baroclinic response to wind. Next, we decompose the simulated response as a proxy of the observed response to identify the dominant baroclinic modes which caused the observed dynamics. Decomposition of response From the model decomposition of the simulated baroclinic response, seven modes were identified as energetically dominant (Table 5.1). During one week from the start of the storm (days 233-240), these dominant modes had average modal energies Et,r that exceeded 0.25 GJ and reached maximum values of Et,r > 1 GJ. The V1H1 mode contained the highest average energy over days 233-240, Et,r = 6.6 GJ, and the highest maximum energy of 11.4 GJ on day 236.50. Next in order, the V2H1 and V1H2 modes had comparable average energies, Et,r ' 1.5 GJ. However, V1H2 reached a higher maximum energy (Et,r = 4.8 GJ on day 236.12) than V2H1 (Et,r = 1.9 GJ on day 237.0). The V2H3, V2H5, V2H6, and V1H5 modes followed with average Et,r between 0.3 and 0.8 GJ. The maximum energy of these modes reached relatively high values between 1.1 and 1.9 GJ. Other baroclinic modes with Tr > 0.8 d (not listed in Table 5.1) had average Et,r < 0.25 GJ and maximum Et,r < 0.5 GJ. We examined the adequacy of the seven dominant modes in resolving isotherm deflec- tions at mooring KN and NA (Figure 5.8). For the upper density interface, superposition of contributions from the dominant modes, Σrη2,r, at the moorings captures most of the deflections of the 16°C isotherm simulated using ELCOM (Figure 5.8a,b) which, in turn, is representative of the observed deflections. For the lower density interface, the superposition of η3,r for the dominant modes does not capture all the simulated deflections of the 7°C Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 123 Table 5.1: Characteristics of baroclinic modes which dominated the simulated baroclinic response. The listed modes have average and maximum modal energies > 0.25 GJ and > 0.5 GJ, respectively. Mode Period Tr (d) Total Modal Energy, Et,r (GJ) Average Maximum V2H1 11.9 1.4 1.9 V1H1 6.8 6.6 11.4 V2H3 3.9 0.7 1.5 V1H2 3.3 1.6 4.8 V2H5 3.1 0.8 1.9 V2H6 2.3 0.4 1.1 V1H5 1.7 0.3 1.6 isotherm at the moorings (Figure 5.8c,d). At mooring NA, the superposition of η3,r is in agreement with the low frequency deflections of the 7°C isotherm but higher modes are needed to resolve the high frequency deflections (Figure 5.8d). Overall, because the lower interface has smaller deflections than the upper interface, the dominant modes can explain the observed large amplitude variations in the depth and thickness of the metalimnion. To characterize the typical features of the baroclinic response of Knewstubb and Natalkuz Lakes, we examined the individual contributions of dominant modes to metalimnetic deflec- tions, 0.5(η2+η3), and changes in the metalimnion thickness, η2−η3. For snapshot S2 (day 235.4), metalimnetic deflections are mainly attributed to the V1 modes; V1H1, V1H2, and V1H5 (Figure 5.9). Contributions from V2 modes are much smaller. In Knewstubb Lake, the V1H1 mode deepened the metalimnion while the V1H2 mode raised the metal- imnion but with much smaller amplitudes. The V1H5 mode had a nodal line at the north end of Knewstubb mid-reach. South of this nodal line, the metalimnion deflected upwards while, north of the nodal line, the metalimnion was depressed down with an upward tilt from Knewstubb Arm to Big Bend Arm as observed in Figure 5.7d. In the east part of Natalkuz Lake, the V1H1 displaced the metalimnion upwards which was counteracted by deepening due to the V1H2 and V1H5 modes. In the west part of Natalkuz Lake, both V1H1 and V1H2 caused upward deflection of the metalimnion as observed in Figure 5.7a,c. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 124 (a) η 2 (m ) −9 −6 −3 0 3 (b) (c) η 3 (m ) 22 7 23 0 23 3 23 6 23 9 Julian Day in 1994 −9 −6 −3 0 3 (d) 22 7 23 0 23 3 23 6 23 9 Julian Day in 1994 Figure 5.8: Deflections of the upper (η2) and lower (η3) interfaces at mooring KN (a,c) and NA (b,d). The superposition of deflections from the dominant modes listed in Table 5.1 (black lines) are compared to deflections simulated with ELCOM (grey lines). Figure 5.9: Snapshot S2 (day 235.4) of metalimnetic deflections 0.5(η2,r+η3,r) of dominant modes. Black contours indicate zero deflections (i.e. nodal lines). Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 125 226 228 230 232 234 236 238 −1 −0.5 0 0.5 1 Julian Day in 1994 a ⋆ r (a) V1H1 V2H1 226 228 230 232 234 236 238 Julian Day in 1994 (c) V1H2 V1H5 226 228 230 232 234 236 238 Julian Day in 1994 (b) V2H3 V2H5 V2H6 Figure 5.10: Normalized modal amplitudes of interface deflection between day 225.5 and 240 for modes which dominate the baroclinic response of Knewstubb and Natalkuz Lakes. Vertical dotted lines indicate the time of snapshot S2. Modal amplitudes are normalized by the amplitude at S2. The above features correspond to snapshot S2. The metalimnetic deflections 0.5(η2,r+η3,r) at any time have the same structure as for snapshot S2 with a magnitude proportional to the ratio of ar at the investigated time and ar at time S2. We denote this ratio as a∗r and, if a∗r < 0, the structure of the metalimnetic deflections is the inverse of that at S2; i.e. has the inverse sign of the structures shown in Figure 5.9. To track the evolution of modal contributions to metalimnetic deflections, time-series of a∗r for each of the dominant modes was examined (Figure 5.10). The largest magnitude of a∗r was about 1.2 which indicates that ηi,r for snapshot S2 nearly corresponds to the maximum deflections associated with the dominant modes. For V1H1, a∗r > 0 during the entire simulation (days 225.5-240, Figure 5.10a) which reflects a persistent metalimnetic deepening in Knewstubb Lake and lifting in Natalkuz Lake due to V1H1. For V1H2 and V1H5, a∗r shows irregular fluctuations with reversals in sign (Figure 5.10c). During the storm (days 233.5-236.5), a∗r > 0 which reveals persistence of the V1H2 and V1H5 metalimnetic deflections observed at S2. Noteworthy, the relative phase and magnitude of the deflections caused by the V1H1 and V1H2 modes at the KN and NA moorings are in agreement with the results of a previous application of a two- layer variable cross-section (TVC) model to a longitudinal axis running from Knewstubb Arm to Chelaslie Arm (Chapter 4). We now focus on changes in the metalimnetic thickness η2,r−η3,r (Figure 5.11). Both V1 and V2 modes contributed to changes in thickness. While the upper and lower interfaces for V1 modes mostly deflect in-phase, the deflections at the same location have different Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 126 amplitudes which result in changes in the thickness of the intermediate density layer. In Knewstubb and Big Bend Arms, the metalimnion thickness was most affected by V2H1 and V2H5 which caused metalimnetic compression and by V2H3 and V1H5 which caused metalimnetic expansion. In Knewstubb mid-reach, the metalimnion was compressed due to V2H1 and V2H3 and expanded due to V2H5. The V1H5 caused expansion in the south part of the mid-reach and compression in the north part as observed in Figure 5.7a,c. In the east part of Natalkuz Lake, the V1H1 and V2H5 modes caused metalimnetic compression which was counteracted by the V1H2 mode. In the west part of Natalkuz Lake, the V1H1 and V2H1 modes were associated by metalimnetic expansion which was augmented in Euchu Arm by modes V2H5 and V2H6. The a∗r series for V2H1 is very similar to that for V1H1 (Figure 5.10a) indicating that these two modes are responsible for the long time-scale asymptotic motions observed at mooring KN. In contrast, mode V2H6 was only energetic during the storm as indicated by a∗r ∼ 0 prior to day 234 and following day 238. For the V2H3 and V2H5 modes, the a∗r series show similar and nearly regular fluctuations (Figure 5.10b). On day 238, a∗r for V2H5 reversed sign reflecting a structure opposite to that in Figure 5.11d which is associated with metalimnetic expansion in the east part of Natalkuz Lake. Examination of the magnitudes of a∗r for the dominant modes between day 238 and 239 with the corresponding structures in Figure 5.11 indicates that the V2H5 mode was the main contributor to the metalimnetic intrusion which caused the observed expansion at mooring NA around day 239. Most of the energy input by wind to the modes was dissipated and this dissipation oc- curred rapidly. This is shown in Figure 5.12 by the small difference between the cumulative energy input by wind and the cumulative dissipated energy for each mode. Damping ratios ζr calculated for each mode showed relatively high variability which is reflected in Figure 5.13 by the 10th, 50th and 90th percentiles. The main feature conveyed by the percentiles is the general decreasing trend in ζr from low to high modes. Assuming a constant decay rate κr for all modes, damping ratios should be linearly proportional to modal periods Tr as given by ζr = κrTr (2pi)−1. A linear best fit to the median ζr with κr = 3.0 d-1 satisfactorily reproduces the decreasing trend in ζr (solid line in Figure 5.13). The 10th and 90th percentiles have a similar trend to the median ζr. This reflects a fixed relationship between damping of different modes which supports the assumptions made in Chapter 2 about relative damping of horizontal modes and in Chapter 3 about relative damping of vertical modes. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 127 Figure 5.11: Snapshot S2 (day 235.4) of variations in metalimnetic thickness η2,r−η3,r for dominant modes. Black lines indicate no changes in thickness. The modal structure in (g) was classified as H5 according to the number of nodal lines in the metalimnetic deflection (Figure 5.9g). Julian Day in 1994 E n er g y (G J ) 226 228 230 232 234 236 238 240 0 10 20 30 40 50 60 70 80 90 V1H1 V1H2 V2H1 V2H5 V2H3 V2H6 V1H5 Figure 5.12: Cumulative energy input by wind (solid) and dissipated energy (dashed) for the dominant modes. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 128 Tr (d) ζ r V 1 H 1 V 1 H 2 V 2 H 1 V 2 H 5 V 2 H 3 V 2 H 6 V 1 H 5 24681012 0.1 1.0 10 10th 50th 90th Figure 5.13: Damping ratios of the dominant baroclinic modes corresponding to the 10th, 50th (median), and 90th percentiles. For each mode, damping ratios estimated when energy dissipation rate Dr was less than 10% of its maximum value were excluded from calculating the percentiles. The solid line represents a linear best fit to the median ζr as explained in the text. The damping ratios in Figure 5.13 also shows that the V1H1 mode was overdamped (ζr > 1). This is qualitatively in agreement with the previous application of the TVC model (calibrated ζr = 1.3) (Chapter 4). Further, Figure 5.13 shows that the V1H2 mode has a median ζr = 0.8 which is in good agreement with ζr = 0.7 predicted by the TVC (Chapter 4). 5.4 Discussion In single basin lakes, wind aligned along the lake thalweg is expected to induce a monotonic tilt of the metalimnion from the upwind endwall towards the downwind endwall. In lakes with multiple basins and arms, this expected metalimnetic response is modified by bathymet- ric features including geometric constrictions and changes in thalweg orientations relative to the wind field. In addition to the domain-wide tilt, these bathymetric features act to divide the domain into sub-basins in which the wind induces local tilts. The resulting metalimnetic response consists of these sub-basin scale tilts superimposed on the domain-wide tilt. In Knewstubb and Natalkuz Lakes, the geometric constriction in the Narrows leads to sub-basin tilts in each of the lakes, along with a metalimnetic tilt that extends over the two Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 129 lakes (Figure 5.14). Additionally, the response of Knewstubb Lake is modified by the sharp bend in the thalweg between Knewstubb Arm and Knewstubb mid-reach. Because the wind field over Nechako Reservoir is nearly uniform, wind projected on the thalweg does not follow the same direction from the Narrows to the west end of Knewstubb Arm; the wind blows from SW to NE along Knewstubb mid-reach and from W to E along Knewstubb Arm 4. This leads to a sub-basin tilt along Knewstubb Arm that is not a continuation of the metalimnetic tilt along Knewstubb mid-reach (Figure 5.14). The response would be different if the northeast wind along Knewstubb mid-reach turned towards the west at the junction with Knewstubb Arm. In this case, we expect the metalimnetic tilt along Knewstubb mid-reach to extend into Knewstubb Arm. In summary, the metalimnetic response in Knewstubb and Natalkuz Lakes shows that sub-basin tilts are caused by the following bathymetric features: geometric constrictions (the Narrows) and changes in thalweg orientation relative to wind direction (bend between Knewstubb Arm and mid-reach). The effect of geometric constrictions can be isolated by examining the response of Amisk Lake which is composed of two basins connected by a narrow elevated channel and has a straight thalweg (Chapter 3). The V1H1 mode of Amisk lake is responsible for the interbasin exchange and lake-wide metalimnetic tilt (see Figure 3.9b,c in Chapter 3). The higher horizontal modes V1H2 and V1H3 which resemble V1H1 modes of the decoupled basins are responsible for metalimnetic oscillations in the individual sub-basins. Similar to Amisk Lake, modal decomposition of the simulated response in Knewstubb and Natalkuz Lakes revealed that the domain-wide tilt is caused by the V1H1 mode while the sub-basin tilts are associated with higher horizontal modes, V1H2 and V1H5 (Figure 5.14). In particular, V1H2 significantly contributes to the local metalimnetic tilt along Natalkuz Lake. The excitation of this mode can be explained by decoupling Natalkuz Lake from Knewstubb Lake. The effect of the WSW wind is to induce a metalimnetic tilt within the decoupled sub-basin that would be attributed to its V1H1 mode. Examination of Figure 5.4f shows that the spatial structure of V1H2 of Knewstubb and Natalkuz Lakes is similar to the expected V1H1 modal structure of the latter lake uncoupled from Knewstubb Lake. The metalimnetic tilts in the sub-basins of the domain are generally associated with compression of the metalimnion at the downwind end of the sub-basins and expansion at the upwind end. These variations in thickness are caused by a combination of V1 and V2 modes. In particular, considerable thickness variations are attributable to the V1H1 4See Figure 4.2d in Chapter 4. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 130 Figure 5.14: Schematic of typical features of the baroclinic response to wind and modal contributions. Black lines indicate the deflected upper and lower bounds of the metalimnion (solid) and their initial levels (dashed). Modes which to contribute to the observed metalimnetic tilts along the sub-basins are denoted in ellipses while modes which contribute to the observed metalimnetic compression or expansion at sub-basin ends are denoted in rectangles. Arrows indicate along-thalweg wind direction. mode (Figure 5.11a) which is not usually credited for inducing changes in the metalimnion thickness (Münnich et al., 1992). The changes in thickness associated with V1H1 can be explained by two features of the spatial structure of the interface deflections Zi,r for this mode (Figure 5.4a). 1) The upper and lower interfaces mostly deflect in-phase but with different magnitudes. In Knewstubb Lake, deflections of the upper interface are larger than the lower interface. This relative deflection is reversed in the west part of Natalkuz Lake. The difference in deflection magnitude causes changes in metalimnetic thickness. 2) The nodal lines of Z2,r and Z3,r for V1H1 are not vertically aligned. Between the two nodal lines (in Intata Reach and its junction with the main basin of Natalkuz Lake), the upper and lower interfaces deflect in anti-phase and hence induce changes in metalimnetic thickness. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 131 While several modes are excited along with V1H1, modal decomposition shows that the V1H1 mode is the most energetic mode in Knewstubb and Natalkuz Lakes (Table 5.1). This is similar to the observation that, in many lakes of simpler geometry, V1H1 dominates the baroclinic response to wind (Mortimer, 1953; Hutter et al., 1983; Lemmin, 1987; Wiegand and Carmack, 1986). Mortimer (1952) attributed the dominance of V1H1 to a smaller damping rate than higher horizontal modes but also acknowledged that the spatial wind distribution can play a role in suppressing the V1H1 mode and preferentially exciting higher modes (Mortimer, 1952, 1953). In the north basin of Lake Lugano (which is baroclinically decoupled from the south basin), Hutter et al. (1983) attributed the dominance of V1H1 to resonance with wind. Similarly, in Kootenay Lake, Wiegand and Carmack (1986) noted the dominance of V1H1 and described the amplification and damping of this mode due to seasonal tuning in and out with energy-containing wind frequencies. In the main basin of Lake Zurich (which is baroclinically decoupled from the sub-basin), Horn et al. (1986) observed the dominance of V1H1 over higher gravity modes and focused on the role of the timing of wind pulses relative to existing V1H1 oscillations in strengthening or damping these oscillations. In Lake Baldegg, Lemmin (1987) hypothesized that the dominance of V1H1 resulted from an energy cascade where the V1H1 is directly energized by wind. Following cessation of wind, gradual leakage of the V1H1 energy feeds energy into higher modes. As the receiver of wind energy and the energy source for higher modes, V1H1 emerges as the most dominant mode (Lemmin, 1987). Along with amplification due to resonance and possible weaker damping, we argue that the prominence of the V1H1 in these lakes is correlated to the relative spatial uniformity of the wind field. In the the north basin of Lake Lugano, wind data shown by Hutter et al. (1983) reveals reasonable along-thalweg uniformity (their Figure 8). In Kootenay Lake, Stevens et al. (1996) argues that, for powerful synoptic-scale storms, wind is relatively uniform along the 100 km thalweg of Kootenay Lake. Similarly, Lemmin et al. (2005) expects a uniform wind field over Lake Baldegg due to its small size. In Lake Zurich, although Horn et al. (1986) notes the spatial non-uniformity of the wind field over the lake, this is not reflected by wind measurements that they present. For the 29 km long main basin of Lake Zurich, the along-thalweg wind stress at two stations ∼ 14 km apart shows good agreement (see Horn et al., 1986, their Figure 3) suggesting relative along-thalweg uniformity of wind forcing. In all these lakes where the V1H1 mode is prominently excited, there is evidence that the wind distribution is characterized by some uniformity. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 132 In lakes of large lateral extent, the dominance of rotation-affected V1H1 mode (ex- hibiting Kelvin-like behavior) has been frequently observed. Examples include Lake Biwa (Kanari, 1974; Shimizu and Imberger, 2008), Lake Tiberias (Antenucci et al., 2000; Gomez- Giraldo et al., 2006; Shimizu and Imberger, 2008), Upper Lake Constance (Appt et al., 2004), and Lake Geneva (Lemmin et al., 2005). Though the wind field over these lakes is spatially non-uniform, the response of the V1H1 kelvin-like mode is well explained by models forced by the domain-average wind speed and direction; i.e. by a uniform wind distribution (Laval et al., 2003; Umlauf and Lemmin, 2005). The evidence compounded above indicates that the dominance of V1H1 mode is related to the spatial uniformity of the wind field or to large domain-average wind relative to spatial variability. The observations cited above are from lakes of simple geometries. The dominance of the V1H1 in Knewstubb and Natalkuz Lakes extends these observations to lakes with multiple basins and arms that are oriented at different angles relative to the wind field. In disagreement with the hypothesis of Mortimer (1952, 1953), the V1H1 mode in Knewstubb and Natalkuz Lakes was found to be more damped than higher horizontal modes even though the V1H1 is the most energetic. This is attributable to the large energy input by wind to this mode (Figure 5.12) which, in turn, is related to the considerable spatial uniformity of the wind field indicated by wind measurements at KN, NA, and SK (Chapter 4). As shown by Shimizu et al. (2007), amplitudes (ar and br) of baroclinic modes are proportional to how strong the modal velocity structures in the upper layer ( Q1,rh−11 ,P1,rh −1 1 ) are represented in the spatial distribution of wind stress (τx,τy). Decomposition of a uniform wind distribution in terms of the structures Q1,rh−11 and P1,rh −1 1 will generally result in the lowest mode being the strongest component. This is supported by decomposition of uniform wind along the thalweg of a flat-bottom two-layer narrow basin where the magnitude of the wind components cor- responding to Q1,rh−11 and P1,rh −1 1 of the V1Hn modes is largest for V1H1 and diminishes with increasing the rank of the modes (i.e. n) in proportion to n−1 (Heaps and Ramsbottom, 1966). In summary, the V1H1 mode is expected to dominate the baroclinic response to wind even in lakes of complex geometry as Knewstubb and Natalkuz Lakes provided that the wind is spatially uniform or has a significant spatial-average magnitude. Chapter 5. The Baroclinic Response to Wind in a Multi-arm Multi-basin Reservoir 133 5.5 Conclusions Field observations, three-dimensional hydrodynamic numerical modeling, and modal de- composition were used to characterize the typical features of the baroclinic response in Knewstubb and Natalkuz Lakes to prevailing WSW winds and to identify the baroclinic modes which dominate the response. The study methodology is suited for elongated multi- arm lakes where rotation of the earth has a secondary effect on the response. The methodol- ogy is also useful in inferring modal damping particularly if field observations are insufficient for reliable calibration of damping rates. By using three-dimensional numerical modeling to expand the field data, we were able to infer the dominant modes and their damping rates at a level of confidence close to that in the model representation of the actual response to wind. For Knewstubb and Natalkuz Lakes, the model revealed features of the baroclinic re- sponse that have not been well studied in the context of multi-arm, multi-basin water bod- ies, namely, the establishment of sub-basin metalimnetic tilts that are superimposed on a domain-wide tilt and the creation of along-thalweg variations in metalimnetic thickness associated with the sub-basin tilts. These features are attributable to geometric constric- tions and changes in thalweg orientation relative to the wind field. The model also shows the occurrence of a metalimnetic intrusion which explains the observed expansion of the metalimnion subsequent to a wind storm. For the complex bathymetry of Knewstubb and Natalkuz Lakes, the results of modal decomposition indicated the dominance of the V1H1 mode which is similar to the case in lakes of less complicated geometries. The V1H1 mode was found to be responsible for the domain-wide tilt while higher horizontal modes caused the sub-basin tilts. The modal decomposition also indicated that low modes are more strongly damped than higher modes. To improve the characterization of the baroclinic response, the next step is to identify the effect of rotation and dissipation on baroclinic modal structures and how they influence the response to wind. 134 Chapter 6 Conclusions 6.1 Summary This study examined the baroclinic response to wind in narrow, elongated lakes. The exam- ined lake bathymetries include idealized single basins of uniform width and variable depth (Chapter 2), a two-basin lake of non-uniform cross-section and a nearly straight thalweg (Amisk Lake, Chapter 3), and two connected lakes with multiple arms (Knewstubb and Natalkuz Lakes in Nechako Reservoir, Chapter 4 and 5). For these lakes, the modal composi- tion of the baroclinic response was analyzed and linked to lake bathymetry and stratification, temporal and spatial characteristics of wind forcing, and damping. To achieve the goals of the study, field observations, modal analysis models (MAMs), modal-based forced models (MFMs), and the three-dimensional hydrodynamic numerical model ELCOM were used. Neglecting the thickness of the metalimnion, the two-layer variable cross-section MFM model (TVC) was developed to simulate the baroclinic response to wind by superposition of horizontal baroclinic modes. The TVC was used in different forms; with uniform width and step forcing in Chapter 2, with Fourier forcing components in Chapter 3, and with impulsive forcing and excluding side arms in Chapter 4. The thickness of the metalimnion was accounted for by using a three-layer box MFM in Chapter 3, and using the three-layer variable-depth complete MAM (ThVDC) in Chapters 4 and 5. Because field observations in Knewstubb and Natalkuz Lakes were not sufficient for full characterization of the energetically dominant baroclinic modes, ELCOM was used to simulate the evolution of the three-dimensional thermal and velocity fields which were then spatially decomposed into constituent baroclinic modes (Chapter 5). In Chapter 2, the effect of along-thalweg depth variability on the response was examined for idealized curved-bottom and wedge-shaped basins. For Amisk Lake (Chapter 3), the ex- citation and modulation of horizontal baroclinic modes were linked to characteristics of wind forcing. Also, the evolution of the vertical structure for the exchange between the basins of Chapter 6. Conclusions 135 Amisk Lake were linked to the interaction of vertical modes excited due to the thickness of the metalimnion. In Nechako Reservoir (Chapter 4), the longitudinal modes which dominate the baroclinic response to wind in Knewstubb and Natalkuz Lakes were identified and the strong modal damping was demonstrated. For the strongest summer storm over the lakes, the baroclinic response was modeled using ELCOM and spatially decomposed into constituent horizontal and vertical modes (Chapter 5). The analysis revealed the energetic modes which contributed to metalimnetic deflections and changes in metalimnetic thickness on domain- wide and sub-basin scales. 6.2 Contributions This study contributed the following to the understanding of the effects of bathymetry, stratification, wind forcing characteristics, and damping on the baroclinic response to wind: • Bathymetry and spatial wind distribution – Similar to lakes of simple geometries, the V1H1 mode dominates the baroclinic response in lakes with baroclinically-connected multiple basins. This V1H1 dominance is driven by uniform along-thalweg wind or wind with a significant spatial-average magnitude that preferentially couples with the spatial distribution of layer flow for H1 since it has the most uniform distribution among horizontal modes. In Amisk Lake, because of the absence of topographic features that can cause sheltering, wind along the thalweg of the lake is expected to be uniform. Based on the TVC with uniform wind stress, the exchange flow between the two basins of Amisk Lake was shown to be dominated by V1H1. Thermistor mooring data also show the dominance of V1H1 (Chapter 3). For Nechako Reservoir, records from wind stations indicated the near-uniformity of the wind field over Knewstubb and Natalkuz Lakes. The resulting along-thalweg wind was directed from west to east with a significant along-thalweg average magnitude. In response, the V1H1 mode was found to be the most energetic mode (Chapter 4 and 5). – In lakes and reservoirs, bathymetric features such as geometric constrictions (e.g., sills and contractions) or changes in thalweg orientation relative to the wind Chapter 6. Conclusions 136 field act to create sub-basins in which wind forcing causes local metalimnetic tilts. These sub-basin scale tilts deepen (raise) the metalimnion at the downwind (upwind) ends of the sub-basins and are superimposed on the lake-wide V1H1 tilt. The metalimnetic tilts in the sub-basins are attributable to higher horizontal modes that are equivalent to V1H1 modes of the sub-basins. Because of the com- plex bathymetry in Nechako Reservoir, strong summer storms from the WSW cause local metalimnetic tilts in the sub-basins (Knewstubb and Big Bend Arms, Knewstubb mid-reach, and Natalkuz Lake) on top of the domain-wide V1H1 tilt from Natalkuz to Knewstubb Lake (Chapter 5). The most energetic mode after V1H1 was V1H2 which can be viewed as the V1H1 mode of Natalkuz Lake (if decoupled from Knewstubb Lake) (Chapter 5). Also, in Amisk Lake, thermistor mooring data and the TVC showed the excitation of V1H2 which is equivalent to V1H1 of the north basin (Chapter 3). – Along the thalweg, reservoirs commonly have a wedge-shape formed by a ver- tical dam at one basin end and a bottom shoaling towards the other end. Under uniform wind forcing, asymmetry in the wedge-shaped bottom profile leads to asymmetry in transient metalimnetic deflections between basin ends; the deflec- tions at the shallow end are larger than at the deep end. The response associated with the wedge profile is also characterized by high localized velocity shear at the shallow end. The potentially large metalimnetic deflections and velocity shear may have implications on entrainment of streams inflows that typically enter reservoirs from the shallow end. This interaction will depend on the along- thalweg direction of the wind as to whether it deflects the metalimnion upward or downward at the shallow end (Chapter 2). – The baroclinic response in symmetric curved-bottom basins is associated with more uniform along-thalweg distributions of velocity shear compared to flat- bottom basins. This more uniform velocity shear reduces variability in shear- induced entrainment and raises the confidence in the results of vertical one- dimensional mixed layer models which assume horizontal uniformity of entrain- ment processes (Chapter 2). Chapter 6. Conclusions 137 • Bathymetry and stratification – If the thickness of the metalimnion in two-basin lakes is significant, the verti- cal structure of interbasin exchange flow is determined by interaction between vertical modes. At the beginning of a wind storm, the inter-basin exchange is initially dominated by the first vertical mode, V1, which causes a two-layer exchange flow. The second and higher vertical modes which have longer periods than V1 gradually shift the exchange to three and more layers (Chapter 3). By the time the exchange shifts to more than two layers, the exchange flows are significantly smaller than during two-layer exchange because of the stronger damping of higher vertical modes (e.g., V2) than V1 (Chapters 3 and 5). – When the thickness of the metalimnion is significant in lakes with geometric constrictions and changes in thalweg orientation relative to the wind field, the sub-basin metalimnetic tilts mentioned above are associated with variations in the metalimnetic thickness along the sub-basins; the metalimnion is compressed at the downwind end and expanded at the other end of each sub-basin (Chapter 5). • Bathymetry and damping – In basins connected by elevated and laterally-constricted reaches, modes with flow anti-nodes located in the geometrically constricted reach can be strongly damped. Evidence from Knewstubb and Natalkuz Lakes and Amisk Lake reveal the overdamping of the V1H1 which has its flow anti-node within the connecting reach (Chapter 3 and 4). • Temporal wind pattern and damping – If wind forcing is aperiodic, the strength of modal damping can be qualitatively inferred from spectral and wavelet coherence between the wind and the baro- clinic response. For periodic forcing, coherence is expected to be independent of damping and hence inference of the strength of damping from coherence is difficult (Chapter 4). – Modal response to periodic wind of different frequencies depends on the strength of modal damping. If damping is strong, a baroclinic mode resonates with a Chapter 6. Conclusions 138 wider range of forcing frequencies around the modal period than if damping was weak. The implication is that a strongly damped mode can maintain resonance with periodic wind even as the modal period changes during the stratified season (Chapter 3). This study also contributed to the methods used in the analysis of the baroclinic response to wind by the development of the two-layer variable cross-section (TVC) model for nar- row, elongated lakes with thalwegs that are not too sinuous. The TVC is a simple modal- based forced model which simulates the baroclinic response from the decoupled evolution of baroclinic modes. The usefulness of the TVC is in linking the excitation and modulation of baroclinic modes to lake bathymetry and the along-thalweg distribution and temporal pattern of wind forcing. 6.3 Future Work The modal-based models used in this study idealized the stratification as layers of homoge- neous density. To some extent, placing the interfaces between density layers is subjective and represents a source of error in modal analysis and decomposition. In modal analysis of lakes, accounting for continuous stratification (CS) has been limited to one horizontal dimension (longitudinal) with no consideration of the lateral dimension (Münnich, 1996; Fricker and Nepf, 2000; Vidal et al., 2005, 2007). Extension to two dimensions is hampered by the need for intensive computational resources (Fricker and Nepf, 2000; Shimizu, 2008), robust numerical schemes to avoid artificial modes (Maas and Lam, 1995), and greater effort in processing and interpreting the three-dimensional output (Shimizu, 2008). Until two- dimensional CS modal analysis becomes practical, analysis of baroclinic modes of narrow, elongated lakes could be improved by adding width variability to longitudinal CS models. Moreover, the CS variable-width model (CSVW) could be used in modal decomposition of spatial temperature and velocity fields or in simulating evolution of individual baroclinic modes as outlined by Shimizu et al. (2007) and Shimizu and Imberger (2008) for layered stratification. The necessary condition for this is the orthogonality of baroclinic modes in longitudinal variable width domains. Fulfillment of this condition is expected given the modal orthogonality for continuous stratification in two horizontal dimensions (Shimizu, 2008). Chapter 6. Conclusions 139 Proper understanding of flow and transport dynamics in a specific lake requires extensive studies. For example, several studies and field campaigns were needed to advance the understanding of the baroclinic response of Lake Tiberias including Antenucci et al. (2000); Hodges et al. (2000); Laval et al. (2003); Antenucci and Imberger (2003); Gomez-Giraldo et al. (2006); Marti and Imberger (2006); Shimizu and Imberger (2008, 2010). This large body of research on Lake Tiberias shows how much work is needed to fully resolve the baroclinic response in any lake. In particular, field observations of high spatial and temporal resolution are needed. In Knewstubb and Natalkuz Lakes, collection of more field data could help validate simulated features including compression and deepening of the metalimnion in Big Bend Arm, metalimnetic upwelling in the west part of Natalkuz Lake, and upward metalimnetic deflection in the south part of Knewstubb mid-reach. In Amisk Lake, the exci- tation of higher horizontal modes and the second vertical mode could be further investigated with thermistor moorings deployed at multiple locations along the thalweg and having faster sampling rates and smaller vertical spacings. With respect to modeling, a natural extension to this work is to consider the effect of rotation and dissipation on the periods and structures of the baroclinic modes of Knewstubb and Natalkuz Lakes and Amisk Lake and to use the rotating dissipative modes to decompose the baroclinic response of the lakes to wind. For Amisk Lake, analysis of the baroclinic response and interbasin exchange can be improved by applying a continuous-stratification variable-width (CSVW) model. Baroclinic exchange between connected basins has important implications for water quality (Lawrence et al., 1997; Umlauf and Lemmin, 2005; Laval et al., 2008) and affects the performance of lake oxygenation systems (Rueda et al., 2010). In Nechako Reservoir, more field observations are needed to properly characterize the exchange between Knewstubb and Natalkuz Lakes through the Narrows and between Natalkuz Lake and the north basins of Nechako Reservoir through the Intata Reach. To quantify the exchange flow, velocity mea- surements in the entire water column are needed. Additionally, within the Narrows, ELCOM reveals abrupt changes in the levels of metalimnetic isotherms. To validate this simulated feature, closely-placed thermistor moorings are needed within the Narrows. Finally, the Narrows represents an opportunity for field studies of dissipation and mixing associated with relatively rapid flow around single trees and tree arrays particularly those penetrating the metalimnion. In Amisk Lake, deeper understanding of the damping of the fundamental mode which controls the interbasin exchange requires quantifying the dissipation in the bottom boundary layer within the connecting channel. An important aspect of such study Chapter 6. Conclusions 140 is to observe the changes in the interbasin exchange and bottom boundary layer dissipation as the thermocline deepens relative to the sill of the connecting channel. Metalimnetic upwelling at upwind basin endwalls has been frequently observed and discussed in lakes (Mortimer, 1952; Laval et al., 2008) and in laboratory experiments (Moni- smith, 1986; Stevens and Imberger, 1996). However, upwelling at the connection between basins in multi-basin lakes, the so-called split layer (Appt et al., 2004), has rarely been observed and analyzed. In response to strong wind forcing, Hollan and Simons (1978) and Appt et al. (2004) observed upwelling of deep water between the main and side basins of Upper Lake Constance. In explaining the upwelling, Appt et al. (2004) suggested that for a split-layer to occur, the removal of epilimnetic water from the upwind basin has to be slower than upwelling in the downwind basin; i.e. the following condition has to be satisfied (Appt et al., 2004) vu < 2ci Lu Ld (6.1) where vu is the average epilimnetic velocity in the upwind basin and Lu and Ld are the lengths of the upwind and downwind basins, respectively. In Upper Lake Constance, a non-uniform wind field due to topography shadowing over the side basin (which is the upwind basin) led to small vu which satisfied the condition in equation 6.1. For Knewstubb and Natalkuz lakes, equation 6.1 necessitates that the average epilimnetic velocity in Natalkuz Lake (which is the upwind basin) be less than 0.8 m s-1 for a split-layer to occur at the south end of Knewstubb mid-reach. ELCOM simulated velocities show that this condition was satisfied; however, a split-layer was not observed in the simulations. The inadequacy of the condition in equation 6.1 indicates the need for more work to be done on split-layers particularly on collecting field observations around basin connections in multi-basin lakes. Perhaps a more comprehensive condition for the occurrence of split layers than equation 6.1 isWd < 1