A General Lake Model (GLM 3.0) for linking with high-frequency sensor data from the Global Lake Ecological Observatory Network (GLEON)

The General Lake Model (GLM) is a onedimensional open-source code designed to simulate the hydrodynamics of lakes, reservoirs, and wetlands. GLM was developed to support the science needs of the Global Lake Ecological Observatory Network (GLEON), a network of researchers using sensors to understand lake functioning and address questions about how lakes around the world respond to climate and land use change. The scale and diversity of lake types, locations, and sizes, and the expanding observational datasets created the need for a robust community model of lake dynamics with sufficient flexibility to accommodate a range of scientific and management questions relevant to the GLEON community. This paper summarizes the scientific basis and numerical implementation of the model algorithms, including details of sub-models that simulate surface heat exchange and ice cover dynamics, vertical mixing, and inflow–outflow dynamics. We demonstrate the suitability of the model for different lake types that vary substantially in their morphology, hydrology, and climatic conditions. GLM supports a dynamic coupling with biogeochemical and ecological modelling libraries for integrated simulations of water quality and ecosystem health, and options for integration with other environmental models are outlined. Finally, we discuss utilities for the analysis of model outputs and uncertainty assessments, model operation within a distributed cloud-computing environment, and as a tool to support the learning of network participants.

Given the diversity of lakes among continents, regionspecific pressures, and local management approaches, the Global Lake Ecological Observatory Network (GLEON; http://gleon.org, last access: 14 January 2019) was initiated in 2004 as a grass-roots science community with a vision to observe, understand, and predict freshwater systems at a global scale . In doing so, GLEON has been a successful example of collaborative research within the hydrological and ecological science disciplines. GLEON aims to bring together environmental sensor networks, numerical models, and information technology to explore ecosystem dynamics across a vast range of scalesfrom individual lakes or reservoirs  to regional (Read et al., 2014;Klug et al., 2012) and global extents O'Reilly et al., 2015). Ultimately, it is the aim of the network to facilitate discovery and synthesis and to provide an improved scientific basis for sustainable freshwater resource management.
Environmental modelling forms a critical component of observing systems as a way to make sense of the "data deluge" (Porter et al., 2012), allowing users to build virtual domains to support knowledge discovery at the system scale (Ticehurst et al., 2007;Hipsey et al., 2015). In lake ecosystems, the tight coupling between physical processes and water quality and ecological dynamics has long been recognized. Modellers have capitalized on a comprehensive understanding of physical processes (e.g. Imberger and Patterson, 1990;Imboden and Wüest, 1995) to use hydrodynamic models as an underpinning basis for coupling to ecological models. Such models have contributed to our understanding of lake dynamics, including applications associated with climate change , eutrophication dynamics (Matzinger et al., 2007), harmful algal bloom dynamics (Chung et al., 2014), and fisheries (Makler-Pick et al., 2011).
In recent decades, a range of one-, two-, and threedimensional hydrodynamic models has emerged for lake simulation. Depending on the dimensionality, the horizontal resolution of these models may vary from metres to tens of kilometres with vertical resolutions from sub-metre to several metres. As in all modelling disciplines, identifying the most parsimonious model structure and degree of complexity and resolution is challenging, and users in the lake modelling community often tend to rely on heuristic rules or practical reasons for model choice (Mooij et al., 2010). Highresolution models are suited to studying events that occur at the timescale of flow dynamics, but are not always desirable for ecological studies over longer timescales due to their computational demands and level of over-parameterization. On the other hand, simple models may be more agile for a particular application and more suited to parameter identification and scenario-testing workflows. However, it has been the case within GLEON that simple models are often less applicable across a wide variety of domains, making them less generalizable, which is a key requirement of synthesis studies across many waterbodies. Despite the fact that there is a relatively large diversity of models and approaches for aquatic ecosystem simulation (Janssen et al., 2015), it is generally agreed that to improve scientific collaboration within the limnological modelling community, there is an increasing need for flexible, open-source community models (Trolle et al., 2012). Whilst acknowledging that there is no single model suitable for all applications, a range of open-source community models and tools can enhance scientific capabilities and foster scientific collaboration and combined efforts (Read et al., 2016). There are examples of such initiatives being successful in the oceanography, hydrology, and climate modelling communities.
With this in mind, the General Lake Model (GLM), a onedimensional (1-D) hydrodynamic model for enclosed aquatic ecosystems, was developed. The lake modelling community has often relied on 1-D models, which originated to capture lake water balance and thermal stratification dynamics (e.g. Imberger and Patterson, 1981;Saloranta and Andersen, 2007;Perroud et al., 2009;Stepanenko et al., 2013). The use of 1-D structure is justified across a diverse range of lake sizes given the dominant role of seasonal changes in vertical stratification in lake dynamics, including oxygen, nutrient, and metal cycling and plankton dynamics (Hamilton and Schladow, 1997;Gal et al., 2009). Despite advances in computing power and more readily available 3-D hydrodynamic drivers, 1-D models continue to remain attractive as they are easily linked with biogeochemical and ecological modelling libraries for complex ecosystem simulations. This allows 1-D models to be used to capture the long-term trajectory and resilience of lakes and reservoirs to climate change, hydrologic change, and land use change. For example, such models have been used to study long-term changes to oxygen, nutrient cycles, and the changing risk of algal blooms (e.g. Peeters et al., 2007;Hu et al., 2016;Snortheim et al., 2017). Furthermore, the low computational requirements of this approach relative to 3-D models is more suited to parameter identification and uncertainty analysis, making it an attractive balance between process complexity and computational intensity.
GLM emerged as a new open-source code in 2012, with the design goal of balancing the complexity of dimensional representation, applicability to a wide range of standing waters, and availability to a broad community (e.g. GLEON has > 700 members from around 50 countries). The scope and capability of the model have developed rapidly with application to numerous lakes and lake types within the GLEON network and beyond (e.g. Read et al., 2014;Bueche et al., 2017;Snortheim et al., 2017;Weber et al., 2017;Menció et al., 2017;Bruce et al., 2018). It is unique in that its suitability now ranges from ephemeral wetlands and ponds to deep lakes, from natural systems to heavily managed man-made reservoirs, and across climatic regions. Given that individual applications rarely engage the full array of features or describe the full details of the model structure, the aim of this paper is to present a complete description of GLM, including the scientific background (Sect. 2) and model code organiza-tion (Sect. 3). The approach to coupling with biogeochemical models is also discussed (Sect. 4), since a main objective of GLM's development is to link its hydrodynamic simulation with water quality models to explore the effects of stratification and vertical mixing on biogeochemical cycles and lake ecology. Finally, an overview of the use of the model within the context of GLEON-specific requirements for model analysis, integration, and education (Sects. 5-6) is described. In order to better define the typical level of model performance across these diverse lake types, a companion paper by Bruce et al. (2018) has undertaken a systematic assessment of the model's error structure against 31 lakes.
2 Model overview 2.1 Background and layer structure The 1-D approach adopted by GLM resolves a vertical series of layers that capture the variation in water column properties. Users may configure any number of inflows and outflows, and more advanced options exist for simulating aspects of the water and heat balance (Fig. 1). Depending on the context of the simulation, either daily or hourly meteorological time series data for surface forcing are required, and daily time series of volumetric inflow and outflow rates can also be supplied. The model is suitable for operation in a wide range of climate conditions and is able to simulate ice formation, as well as accommodating a range of atmospheric forcing conditions.
Although GLM is a new model code written in the C programming language, the core layer structure and mixing algorithms are founded on principles and experience from model platforms including the DYnamic REservoir Simulation Model (DYRESM; Imberger and Patterson, 1981;Hamilton and Schladow, 1997) and the Dynamic Lake Model (DLM; Chung et al., 2008). Other variations have been introduced to extend this underlying approach through applications to a variety of lake and reservoir environments (e.g. Hocking and Patterson, 1991;McCord and Schladow, 1998;Gal et al., 2003;Yeates and Imberger, 2003). The layer structure is numbered from the lake bottom to the surface and adopts the flexible Lagrangian layer scheme first introduced by Imberger et al. (1978) and Imberger and Patterson (1981). The approach defines each layer, i, as a "control volume" (Fig. 1) that can change thickness by contracting and expanding in response to inflows, outflows, mixing with adjacent layers, and surface mass fluxes. As the model simulation progresses, density changes due to surface heating, vertical mixing, and inflows and outflows lead to dynamic changes in the layer structure associated with layers amalgamating, expanding, contracting, or splitting. Notation used throughout the model description is provided in Table 1.
As layers change, their volumes change based on the sitespecific hypsographic curve, whereby the overall lake vol-ume, V max , is defined as H max H 0 A [H ] dH , with the elevation (H ) and area (A) relationship provided as a series of points based on bathymetric data. This computation requires the user to provide a number, N BSN , of elevations with corresponding areas. The cumulative volume at any lake elevation is first estimated as where 2 ≤ b ≤ N BSN . Using these raw hypsographic data, a refined height-area-volume relationship is then internally computed using finer height increments (e.g. H mi ∼ 0.1 m), giving N MORPH levels that are used for subsequent calculations. The area and volume at the height of each increment, H mi , are interpolated from the supplied information as where V mi and A mi are the volume and area at each of the elevations of the interpolated depth vector, and V b and A b refer to the nearest b level below H mi such that H b < H mi .
The interpolation coefficients are computed as Within this lake domain, the model solves the water balance by including several user-configurable water fluxes that change the layer structure. Initially, the layers are assumed to be of equal thickness, and the initial number of layers, N LEV [t = 0], is computed based on the initial water depth. Water fluxes include surface mass fluxes (evaporation, rainfall, and snowfall), inflows (surface inflows, submerged inflows, and local run-off from the surrounding exposed lakebed area), and outflows (withdrawals, overflow, and seepage). Surface mass fluxes operate on a sub-daily time step, t, by impacting the surface layer thickness (described in Sect. 2.2), whereby the dynamics of inflows and outflows modify the overall lake water balance and layer structure on a daily time step, t d , by adding, merging, or removing layers (described in Sect. 2.7). Depending on whether a surface (areal) mass flux or volumetric mass flux is being applied, the layer volumes are updated by interpolating changes in layer heights, whereby V i = f [h i ] and i is the layer number, or layer heights are updated by interpolating changes in layer volumes, whereby Each layer also contains heat, salt (S), and other constituents (C), which are generically referred to as scalars. These are subject to mass conservation as layers change thickness or are merged or split. The specific number of other constituents depends on the configuration of the associated water quality model, but typically includes attributes such as oxygen, nutrients, and phytoplankton. Layer density is computed from the local salinity and temperature according to TEOS-10 (http://teos-10.org, last access: 16 December 2018), whereby ρ i = ρ [T i , S i ]. When density instabilities occur between adjacent layers or when sufficient turbulent kinetic energy becomes available to overcome stable density gradients, then layers merge, thereby accounting for the process of mixing (Sect. 2.6). For deeper systems, a stable vertical density gradient forms seasonally in response to periods of high solar radiation creating warm, buoyant water overlying cooler, denser water, separated by a metalimnion region which includes the thermocline. Layer volumes change due to depth-specific changes in mixing, inflows, and outflows. Thickness limits, z min and z max , are enforced to adequately resolve the vertical density gradient, generally with fine resolution occurring in the metalimnion and thicker cells where gradients are weak. The number of layers, N LEV [t], is adjusted throughout the simulation to maintain homogenous properties within a layer. It has been reported that numerical diffusion at the thermocline can be restricted using this layer structure and mixing algorithm (depending on the minimum and maximum layer thickness limits set by the user), making it particularly suited to longterm investigations and ideally requiring limited site-specific calibration (Patterson et al., 1984;Hamilton and Schladow, 1997;Bruce et al., 2018).
Because this approach assumes layer properties are laterally averaged, the model is suitable for investigations in which resolving the horizontal variability is not a requirement of the study. This is often the case for ecologists and biogeochemists studying central basins of natural lakes (e.g. Gal et al., 2009), managers simulating drinking water reservoirs (e.g. Weber et al., 2017), mining pit lakes (e.g. Salmon et al., 2017), or for analyses exploring the coupling between lakes and regional climate (e.g. Stepanenko et al., 2013). Further, whilst the model is able to resolve vertical stratification, the approach is also able to be used to simulate shallow lakes, wetlands, wastewater ponds, and other small waterbodies that experience well-mixed conditions. In this case, the layer resolution, with upper and lower layer bounds specified by the user, will automatically be reduced, and the mass of water, constituents, and energy will continue to be conserved. The remainder of this section outlines the model components and provides example outputs for five waterbodies that experience a diverse hydrology.

Water balance
The general nature of the model to accommodate a wide diversity of lake types has necessitated flexibility in the configuration of water inputs and outputs (schematically depicted in Fig. 1). The net water flux over the entire lake is summarized as where V S is the total lake volume, t is time, and A S is the lake surface area; the changes due to fluxes at the water sur-face, h S , are expanded upon below, and the remaining inflow and outflow terms are described in detail in Sect. 7. For practical reasons the equation is numerically solved in two stages with different times steps for the surface flux change and all other fluxes. Furthermore, in any given application, not all the inputs and outputs are relevant and users may customize the water balance components accordingly; examples demonstrating lake hydrology from wetlands to reservoirs to deep lakes are presented in Fig. 2. Note that Eq. (4) accounts for the liquid water balance, and in cold climates the model will also track the amount of water allocated into an overlying ice layer (Sect. 2.4), which interacts with the surface water balance as indicated next. The mass balance of the surface layer is computed at each model time step ( t; usually hourly) by modifying the surface layer height, h S , according to where E is the evaporation mass flux computed from the latent heat flux φ E , described below (E = φ E /λ v ρ s ; m s −1 ), R F is rainfall, and S F is snowfall (m s −1 ). Depending on the meteorological conditions, precipitation will either be added to the water volume or to the surface of the ice cover (see Sect. 2.4), and R F and S F therefore influence the water surface height depending on the presence of ice cover according to and Here, f R and f S are user-defined scaling factors that may be applied to adjust the input data values R x and S x , respectively. The surface height of the water column is also impacted by ice formation or melting of the ice layer sitting on the lake surface according to d z ice /dt, as described in Sect. 2.4. Q R is an optional term to account for run-off to the lake from the exposed riparian banks, which may be important in reservoirs with a large drawdown range or wetlands where periodic drying of the lake may occur. The run-off volume generated is averaged across the area that the active lake surface area (A s ) is not occupying, and the amount is calculated using a simple model based on exceedance of a rainfall intensity threshold, R L (m day −1 ), and run-off coefficient: where f ro is the run-off coefficient, defined as the fraction of rainfall that is converted to run-off at the lake's edge, and A max is the maximum possible area of inundation of the lake (the area provided by the user as the N BSN value).
Note that mixing dynamics (i.e. the merging or splitting of layers to enforce the layer thickness limits) will impact the thickness of the surface mixed layer, z SML , but not change the overall lake height. However, in addition to the terms in Eq. (5), h S is modified due to volume changes associated with river inflows, withdrawals, seepage, or overflows, which are described in subsequent sections.

Surface energy balance
A balance of shortwave and longwave radiation fluxes and sensible and evaporative heat fluxes (all W m −2 ) determines the net cooling and heating across the surface. The general heat budget equation for the uppermost layer is described as where c w is the specific heat capacity of water, T s is the surface temperature, and z s and ρ s are the depth and density of the surface layer (i = N LEV ), respectively. The right-hand side (RHS) heat flux terms are numerically computed at each time step and include several options for customizing the individual surface heat flux components, which are expanded upon below.

Solar heating and light penetration
Solar radiation is the key driver of lake thermodynamics and may be input based on daily or hourly measurements from a nearby pyranometer. If data are not available then users may choose to have GLM compute surface irradiance from a theoretical approximation based on the Bird Clear Sky Model (BCSM) (Bird, 1984) modified for cloud cover and latitude. The options for input are summarized as Option 2: sub-daily insolation data provided Option 3: insolation computed from the BCSM where φ SW 0 is the solar radiation flux entering the surface layer, φ SW x is the incoming shortwave radiation flux supplied by the user, f SW is a scaling factor that may be applied and adjusted as part of the calibration process (for example, to capture the effects of shading), and α SW is the albedo for shortwave radiation. If daily data are supplied (option 1), the model continues to run at a sub-daily time step, but applies the algorithm outlined in Hamilton and Schladow (1997) to distribute the daily solar energy flux over a diurnal cycle based on the day of the year, d, and time of day, t − t . For   Figure 2. A 2-year times series of the simulated daily water balance for five example lakes (a-e) that range in size and hydrology. The water balance components summarized are depicted schematically in the inset and partitioned into inputs and outputs. The daily net water flux is computed from Eq. (4). For more information about each lake, the simulation configuration, and input data, refer to the "Data availability" section. option 3 the BCSM is used (Bird, 1984;Luo et al., 2010): where the total irradiance,φ SW , is computed from direct beamφ DB and atmospheric scatteringφ AS components (refer to Appendix A for a detailed outline of the BCSM equations and parameters). In GLM, the clear-sky value is then reduced according to the cloud cover data provided by the user, C x , according to which is based on a polynomial regression of cloud data from Perth Airport, Australia, compared against nearby sen-sor data (R 2 = 0.952; see also a similar relationship by Luo et al., 2010).
The albedo, α SW , is the reflected fraction of the incoming radiation and depends on surface conditions, including the presence of ice, waves, and the angle of incident radiation. For open water conditions, users may configure the following.
Option 1 : daily approximation from Hamilton and Schladow (1997) Here, zen is the solar zenith angle (radians) as outlined in Appendix A, RH x is the relative humidity, ς is the percentage of atmospheric diffuse radiation, d is the day of year, and U 10 is wind speed. The second (oceanic) and third (lacustrine) options are included to allow for diel and seasonal variation of albedo from approximately 0.01 to 0.4 depending on the sun angle (Fig. 3). Option 4 may be better for higher-latitude sites; also note that albedo is calculated separately during ice cover conditions using a customized algorithm, as outlined below in Sect. 2.4. The depth of penetration of shortwave radiation into the lake is wavelength specific and depends on the water clarity via the light extinction coefficient, K w (m −1 ). Two approaches are supported in GLM. The first option assumes that the photosynthetically active radiation (PAR) fraction of the incoming light is the most penetrative and follows the Beer-Lambert law: where z is the depth of any layer from the surface. K w may be set by the user as constant, read in from a time series file, or linked with the water quality model (e.g. FABM Figure 3. Variation of albedo (α SW ) with solar zenith angle (SZA = zen 180/π , degrees) for options 2 and 3 (Eq. 13). For option 3, settings of RH x = 80 % and U 10 = 6 m s −1 were assumed. or AED2; see Sect. 4); in the latter case the extinction coefficient will change as a function of depth and time according to the concentration of dissolved and particulate constituents. For this option Beer's law is only applied for the photosynthetically active fraction, f PAR , which is set as 45 % of the incident light. The amount of radiation heating the surface layer, φ SW S , is therefore the photosynthetically active fraction that is attenuated across z s , plus the en- , which implicitly assumes that the near-infrared and ultraviolet bandwidths of the incident shortwave radiation have significantly higher attenuation coefficients (Kirk, 1994). The second option adopts a more complete light absorption algorithm that integrates the attenuated light intensity across the bandwidth spectrum: where l is the bandwidth index and φ SW i l [z i ] is the radiation flux at the top of the ith layer for the lth bandwidth fraction. For this option, the model by Cengel and Ozisk (1984) is adopted to compute the penetration of individual bandwidth fractions, which more comprehensively resolves the incident and diffuse radiation components of the light climate, taking into account the angle of incident light, transmission across the light surface (based on the Fresnel equations), and reflection off the bottom. These processes are wavelength specific and the user specifies the number of simulated bandwidths, N SW , their respective absorption coefficients, K w l , and reflectivity of light at the sediment, α sed . The light reaching the benthos is relevant in some applications as an indicator of benthic productivity or as a proxy for the type of benthic habitat that might emerge. In addition to the light profiles, GLM therefore predicts the benthic area of the lake in which light intensity exceeds a user-defined fraction of the surface irradiance, f BEN crit , (Fig. 4): where h BEN = h S −z BEN , z BEN is calculated from Beer's law, and the daily average benthic area above the threshold is then reported as a percentage (100 × A BEN /A s ).

Longwave radiation
Longwave radiation can be provided as a net flux, an incoming flux, or, if there are no radiation data from which longwave radiation can be computed, then it may be calculated by the model internally from the cloud cover fraction and air temperature. Net longwave radiation is described as where σ is the Stefan-Boltzmann constant and ε w the emissivity of the water surface, assumed to be 0.985. If the net or incoming longwave flux is not provided, the model will compute the incoming flux from where α LW is the longwave albedo (0.03). The emissivity of the atmosphere can be computed considering emissivity for cloud-free conditions (ε a ) based on air temperature (T a ) and vapour pressure and extended to account for reflection from clouds such that ε * a = f [T a , e a , C x ] (see Henderson-Sellers, 1986;Flerchinger et al., 2009). Options adapted from a range of authors include the following: where C x is the cloud cover fraction (0-1) and e a the air vapour pressure calculated from relative humidity. Note that cloud cover is typically reported in octals (0-8), and thus a value of 1 would correspond to a fraction of 0.125. Some data may also include cloud type and their respective heights. If this is the case, good correspondence has been reported by averaging the octal values for all cloud types to get an average cloud cover. If longwave radiation data do not exist and cloud data are also not available, but solar irradiance is measured, then GLM rad_mode setting 3 will instruct the model to compare the measured and theoretical clear-sky solar irradiance (estimated by the BCSM; Eq. 11) to approximate the cloud cover fraction by assuming that φ SW x /φ SW = f [C x ]. Note that if neither shortwave or longwave radiation is provided, then the model will use the BCSM to compute incoming solar irradiance, and cloud cover will be assumed to be 0 (noting that this is likely to overestimate downwelling shortwave radiation).

Sensible and latent heat transfer
The model accounts for the surface fluxes of sensible heat and latent heat using commonly adopted bulk aerodynamic formulae. For sensible heat, where c a is the specific heat capacity of air, C H is the bulk aerodynamic coefficient for sensible heat transfer, T a the air temperature, and T s the temperature of the water surface layer. The air density (kg m −3 ) is computed from ρ a = 0.348 (1 + r)/(1 + 1.61r) p/T a , where p is air pressure (hPa) and r is the water vapour mixing ratio, which is used to compute the gas constant. For latent heat, where C E is the bulk aerodynamic coefficient for latent heat transfer, e a the air vapour pressure, e s the saturation vapour pressure (hPa) at the surface layer temperature ( • C), ω the ratio of the molecular mass of water to the molecular mass of dry air (= 0.622), and λ v the latent heat of vaporization. The vapour pressure is calculated by the linear formula from Tabata (1973): The net heat fluxes for the example lakes are shown in Fig. 5. . Example light data outputs from a GLM application to Woods Lake, Australia, showing (a) the ratio of benthic to surface light, 100φ PAR BEN /φ SW 0 (%), overlain on the lake map based on the bathymetry, with the area where f BEN crit < 0.2 (i.e. less than 20 % of surface irradiance) depicted in grey, (b) a time series of the depth variation in light (W m −2 ), and (c) a time series of A BEN /A s (as %) for various f BEN crit . Note that the 2-D projection of the 1-D lake model in (a) can assist in managing the lake condition but assumes uniformity of K w .

Correction for non-neutral atmospheric stability
For long time integrations (e.g. seasonal), the bulk-transfer coefficients for momentum, C D , sensible heat, C H , and latent heat, C E , can be assumed approximately constant because of the negative feedback between surface forcing and the temperature response of the waterbody (e.g. Strub and Powell, 1987). At finer timescales (hours to weeks), the thermal inertia of the waterbody is too great, so the transfer coefficients should be specified as a function of the degree of atmospheric stratification experienced in the internal boundary layer that develops over the water (Woolway et al., 2017). Monin and Obukhov (1954) parameterized the stratification in the air column using the now well-known stability parameter, z/L, which is used to define corrections to the bulk aerodynamic coefficients C H and C E using the numerical scheme presented in Appendix B. The corrections may be optionally applied within a simulation, and if enabled, the transfer coefficients used above are automatically updated. To ensure that the data provided are from within the internal boundary layer over the lake surface, users should preferably provide wind speed, air temperature, and relative humidity data that have been collected over the lake surface (at a height of 2-10 m, depending on lake size), supplied at approximately hourly resolution.

Wind sheltering
Wind sheltering may be important depending on the lake size and shoreline complexity and is parameterized according to several methods based on the context of the simulation and data available. For example, Hipsey and Sivapalan (2003) presented a simple adjustment to the bulk-transfer equation to account for the effect of wind sheltering in small reser-voirs using a shelter index to account for the length scale associated with the vertical obstacle relative to the horizontal length scale associated with the waterbody itself. Markfort et al. (2010) estimate the effect of a similar sheltering length scale on the overall lake area calculated based on surrounding topography and canopy heights relative to the water surface. Therefore, within GLM, users may specify the degree of sheltering or fetch limitation using either constant or direction-specific options for computing an "effective" area.
Option 3: user-defined shelter index Here, A WS is a user-defined critical lake area for wind sheltering to dominate, x WS is a user-defined sheltering distance, and L D the lake diameter (L D = 0.5(L crest +W crest )). For option 1, the sheltering factor is held constant for the simulation based on the size of the lake, whereas options 2 and 3 require users to additionally input wind direction data and a direction function, f WS [ wind ], to allow for a variable sheltering effect over time. In the case of option 2, this function scales the sheltering distance, x WS , as a function of wind direction, , 1)), whereas in the case of option 3 the function reads in an effective area scaling fraction directly based on a precalculated shelter index. The ratio of the effective area to the total area of the lake, A E /A S , is then used to scale the wind speed data input by the user, U x , as a means of capturing the average wind speed over the entire lake surface such that U 10 = f U U x A E /A S , where f U is a wind speed adjustment factor that can be used to assist calibration or to correct the raw wind speed data to the reference height of 10 m.

Still-air limit
The above formulations apply when sufficient wind exists to create a defined boundary layer over the surface of the water. As the wind tends to zero (the "still-air limit"), Eqs. (22)-(23) become less appropriate as they do not account for free convection directly from the water surface. This is a relatively important phenomenon for small lakes, cooling ponds, and wetlands since they tend to have small fetches that limit the energy input from wind. These waterbodies may also have large areas sheltered from the wind and will develop surface temperatures warmer than the atmosphere for considerable periods. Therefore, users can optionally augment Eqs. (22)-(23) with calculations for low wind speed conditions by calculating the evaporative and sensible heat flux values for both the given U 10 and for an assumed U 10 = 0. The chosen value for the surface energy balance (as applied in Eq. 9) is found by taking the maximum value of the two calculations: Option 1: no-sheltering area max φ X , φ X 0 Option 2: still-air sheltered area where φ X 0 is the zero-wind flux for either the evaporative or sensible heat flux (φ E 0 and φ H 0 , respectively) and φ X is calculated from Eqs. (22)-(23). The two zero-wind-speed heat flux equations are from TVA (1972), but modified to return energy flux in SI units (W m −2 ).
Here, ν = κ e/p, with the appropriate vapour pressure values, e, for both surface and ambient atmospheric values. K air is the molecular heat conductivity of air (J m −1 s −1 C −1 ), ν a is the kinematic viscosity of the air (m 2 s −1 ), ρ o is the density of the saturated air at the water surface temperature, ρ s is the density of the surface water, f 0 is a dimensionless roughness correction coefficient for the lake surface, and D a is the molecular heat diffusivity of air (m 2 s −1 ). Note that the impact of low wind speeds on the drag coefficient is captured by the modified Charnock relation (Eqs. B2-B3), which includes an additional term for the smooth flow transition (see also Fig. A1).

Snow and ice dynamics
The extent of ice and snow cover can significantly impact the lake water balance and mixing regime depending on the prevailing environmental conditions. The algorithms for GLM ice and snow dynamics are based on previous ice modelling studies that adopt a three-layer scheme for resolving ice and snow split into blue ice (or black ice), white ice (or snow ice), and snow layers (Patterson and Hamblin, 1988;Gu and Stefan, 1993;Rogers et al., 1995;Vavrus et al., 1996;Launiainen and Cheng, 1998;Magee et al., 2016). Blue ice is formed through direct freezing of lake water into ice, whereas white ice is generated in response to seeping of lake water onto the ice surface when the mass of snow that can be supported by the buoyancy of the ice cover is exceeded (see below; Rogers et al., 1995). The snow layer is subject to compaction and melting based on surface meteorological conditions and the ice layers are affected by the lake water temperature at the lower boundary. Blue ice initially forms when the water at the lake surface goes below 0 • C. Once fresh snow deposits on the surface it is subject to densification, which depends on the air temperature and amount of rainfall (Fig. 6); the density of fresh snowfall is determined as the ratio of measured snowfall height to water-equivalent height, with any values exceeding the assigned maximum or minimum snow density (defaults: ρ s,max = 300 kg m −3 , ρ s,min = 50 kg m −3 ) truncated to the appropriate limit. The snow compaction equation is based on the exponential decay formula of McKay (1968), with the selection of snow compaction parameters based on air temperature and depending on whether rainfall or snowfall is being added. When the weight of snow exceeds the buoyancy of the ice layer, the ice will be forced downward and lake water will seep into the snow layer, leading to formation of white ice. This downward movement and white ice formation is limited to the snow amount matching the buoyancy deficit of the ice layer, and the lake height is reduced accordingly. Woods Lake To capture the changing thickness of the ice and snow layers due to melting or freezing, the model employs a quasisteady-state assumption to solve the heat transfer equation through the layers by assuming that the timescale for heat conduction is short relative to the timescale of changes in meteorological forcing (Patterson and Hamblin, 1988;Rogers et al., 1995). By assigning appropriate boundary conditions at the ice-atmosphere and ice-water interfaces, the model computes the upward conductive heat flux through the ice and snow cover to the atmosphere, termed φ 0 .
At the upper surface (which could be ice or snow), a heat flux balance is employed to provide the condition for surface melting: where λ f is the latent heat of fusion, z ice,snow is the height of either the upper snow or ice layer, ρ ice,snow is the density of the relevant snow or ice layer determined from the surface medium properties, T 0 is the temperature at the solid surface, and T m is the meltwater temperature (0 • C). φ net [T 0 ] is the net incoming heat flux for non-penetrative radiation at the solid surface: where the heat fluxes between the solid boundary and the atmosphere are calculated as outlined previously, but with modification for the determination of vapour pressure over ice or snow (e s ice [T 0 ] = e s [T 0 ] 1 + 9.72 × 10 −3 T 0 + 4.2 × 10 −5 T 2 0 ; Jeong, 2009) Figure 6. (a) Decision tree describing updates to the snow cover each time step according to the amount of incident rainfall ( * R F ) and snowfall ( * S F ), air temperature (T a ), and snow compaction rules. (b) Schematic of ice and snow layers and heat fluxes. Refer to the text and Table 1 for definitions of other variables. Here, * R F = f R R x /c secday and * S F = f S S x /c secday if ice cover is present; otherwise they are set to 0 and the model reverts to Eqs. (6)-(7). Rogers et al., 1995). To determine the flow of heat through the layers, Rogers et al. (1995) derived the following.
Here, = z snow K snow + z white K white + z blue K blue , φ SW 0 is the shortwave radiation penetrating the ice-snow surface, K refers to the light attenuation coefficient of the ice and snow components designated with subscripts s, w, and b for snow, white ice, and blue ice, respectively, and the z terms refer to the thickness of snow, white ice, and blue ice. This is rearranged and solved for T 0 and φ 0 by using a bilinear iteration until surface heat fluxes are balanced (i.e. φ 0 [T 0 ] = − φ net [T 0 ]) and T 0 is stable (±0.001 • C). In the presence of ice (or snow) cover, a surface temperature T 0 > T m indicates that energy is available for melting. The amount of energy for melting is calculated by setting T 0 = T m to determine the reduced thickness of snow or ice (as shown in Eq. 32). The estimation of φ 0 applies an empirical equation to estimate snow conductivity, K snow , from its density (Ashton, 1986): The heat flux in the ice at the ice-water interface is where φ si is a volumetric heat flux for the formation of white ice, which is given in Eq. (14) of Rogers et al. (1995), and ice and snow light attenuation coefficients in GLM are also fixed to the same values as those given by Rogers et al. (1995). Shortwave albedo for the ice or snow surface (required for Eq. 10) is a function of surface medium (see Table 1 of Vavrus et al., 1996) with values varying from 0.08 to 0.6 for ice and from 0.08 to 0.7 for snow, depending on the surface temperature and the layer thicknesses; an additional scaling factor for the snow albedo, f α , is also implemented to aid calibration. Accretion or ablation of blue ice occurs at the ice-water boundary based on the conductive heat flux from water into the ice, φ w , as given by the finite-difference approximation: where K water is the molecular conductivity of water (assuming the water is stagnant under the ice), and T is the temperature difference between the surface water of the lake and the bottom of the blue ice layer, T m − T s . This occurs across an assigned length scale, δ wi , for which a value of 0.1-0.5 m is usual based on the reasoning given in Rogers et al. (1995) and the typical vertical water layer resolution of a model simulation (0.125-1.5 m). Note that a wide variation in techniques and values is used to determine the basal heat flux immediately beneath the ice pack (e.g. Harvey, 1990), which suggests that this may need careful consideration during calibration.
The imbalance between φ f moving through the blue ice layer and the heat flux from the water into the ice, φ w , gives the rate of change of ice thickness at the interface with water: The ice thickness is set to its minimum value of 0.05 m, which is suggested by Patterson and Hamblin (1988) and Vavrus et al. (1996). The need for a minimum ice thickness relates primarily to horizontal variability of the ice cover during formation and closure (ice-on) periods. The ice cover equations are discontinued and open water conditions are restored in the model when the thermodynamic balance first produces ice thickness < 0.05 m. Example outputs are shown in Fig. 7; see also Yao et al. (2014) for a previous application.

Sediment heating
The water column thermal budget may also be affected by heating or cooling from the soil-sediment below. For each layer, the rate of temperature change depends on the temperature gradient and the relative area of the layer volume in contact with bottom sediment: where K soil is the soil-sediment thermal conductivity and δz soil is the length scale associated with the heat flux. The temperature of the bottom sediment varies seasonally and also depending on its depth below the water surface such that where z is the soil-sediment zone that the ith layer overlays (see Sect. 4 for details), T z is the temperature of this zone, T z mean is the annual mean sediment zone temperature, δT z is the seasonal amplitude of the sediment temperature variation, and d T z is the day of the year when the sediment temperature peaks. By defining different sediment zones, the model can therefore allow for a different mean and amplitude of littoral waters compared to deeper waters. A dynamic sediment temperature diffusion model is also under development, which will be suitable when empirical data for the above parameters in Eq. (40) are not available.

Stratification and vertical mixing
Mixing processes in lakes are varied and depend upon the degree of meteorological and hydrological forcing, the lake morphometry, and the nature of thermal stratification experienced by the lake at the time of forcing. Numerous models adopt an eddy-diffusivity approach whereby mixing is captured using the advection-dispersion equation (e.g. Riley and Stefan, 1988). GLM adopts an energy balance approach as used in DYRESM whereby the mixing dynamics are based on estimating the amount of turbulent energy available, which is separately computed for the surface mixed layer (surface mixing) and for mixing below the thermocline (deep mixing).

Surface mixed layer
To compute mixing of layers, GLM works on the premise that the balance between the available energy, E TKE , and the energy required for mixing to occur, E PE , provides the surface mixed layer (SML) deepening rate dz SML /dt, where z SML is the depth from the surface to the bottom of the surface mixed layer. For an overview of the dynamics, readers are referred to early works on bulk mixed layer depth models by Kraus and Turner (1967) and Kim (1976), which were subsequently extended by Imberger and Patterson (1981) and Spigel et al. (1986) as a basis for hydrodynamic model design. Using this approach, the available kinetic energy is calculated due to contributions from wind stirring, convective overturn, shear production between layers, and Kelvin-Helmholtz (K-H) billowing. Overall, the turbulent energy generated for mixing is summarized as (Hamilton and Schladow, 1997) where δ KH is the K-H billow length scale (described below), u b is the shear velocity at the interface of the mixed layer, and C K , C W , and C S are mixing efficiency constants. For mixing to occur, the energy must be sufficient to lift up water in the layer below the bottom of the mixed layer, denoted here as the layer k − 1 with thickness z k−1 , and accelerate it to the mixed layer velocity, u * . This must also account for energy consumption associated with K-H billowing. In total, the energy required to entrain a layer into the mixed layer is expressed as E PE : where C T is a mixing efficiency constant to account for unsteady turbulence. To numerically resolve Eqs. (41) and (42) the model sequentially computes the different components of the above expressions with respect to the layer structure, thereby checking the available energy relative to the required amount (depicted schematically in Fig. 8). GLM follows the sequence of the algorithm presented in detail in Imberger and Patterson (1981), whereby layers are combined due to convection and wind stirring first, and then the resultant mixed layer properties are used when subsequently computing the extent of shear mixing and the effect of K-H instabilities.
Plots indicating the role of mixing in shaping the thermal structure of the example lakes are shown in Fig. 9.
To compute the mixing energy available due to convection, in the first step, the value for w * is calculated, which is the turbulent velocity scale associated with convection brought about by cooling at the air-water interface. The model adopts the algorithm used in Imberger and Patterson (1981), whereby the potential energy that would be released by mixed layer deepening is computed from the first moment of layer masses in the epilimnion (surface mixed layer) about the lake bottom relative to the well-mixed condition. This is numerically computed by summing from the bottom-most layer of the epilimnion, k, up to N LEV : where ρ SML is the mean density of the mixed layer including the combined layer, ρ i is the density of the ith layer, z i is the height difference between two consecutive layers within the loop ( The velocity scale u * of the surface layer is associated with wind stress and calculated according to the wind strength: where C D is the drag coefficient for momentum. The model first checks to see if the energy available (Eqs. 43 and 44) can overcome the energy required to mix the k − 1 layer into the surface mixed layer (Fig. 8e); i.e. mixing of k − 1 occurs if where g k = g ρ ρ o is the reduced gravity between the mixed layer and the k − 1 layer calculated as g (ρ SML − ρ k−1 ) /(0.5 (ρ SML + ρ k−1 )). If the mixing condition is met, the layers are combined, the energy required to combine the layer is removed from the available energy, k is adjusted, and the loop continues to the next layer. When the mixing energy is substantial and the mixing reaches the bottom layer, the mixing routine ends. If the condition in Eq. (45) is not met, then any residual energy is stored for the next time step, and the mixing algorithm continues as outlined below.
Once stirring is completed, mixing generated due to velocity shear is then accounted for. Parameterizing the shear velocity, denoted u b , in a one-dimensional model can be problematic; however, the approximation used in Imberger and Patterson (1981) is applied: where u b old is from the previous time step and zeroed between shear (wind) events. Therefore, this model yields a simple linear increase in the shear velocity over time for a constant wind stress. This is considered relative to δt shear , which is the cut-off time beyond which it is assumed that no further shear-induced mixing occurs for that event. This cut-off time assumes the use of only the energy produced by shear at the interface during a period equivalent to half the basin-scale seiche duration, δt iw , which can be modified to account for damping (Spigel, 1978): where δt damp is the timescale of damping. The wave period is approximated based on the stratification as δt iw = L META /2c, where L META is the length of the basin at the thermocline calculated from √ A k−1 (4/π) (L crest /W crest ), whereby an ellipse shape is assumed and c is the internal wave speed, where δ epi and δ hyp are characteristic vertical length scales associated with the epilimnion and hypolimnion: where V epi and V k−1 are the associated volumes.
The time for damping of internal waves in a two-layer system can be parameterized by estimating the length scale of the oscillating boundary layer, through which the wave energy dissipates, and the period of the internal standing wave (see Spigel and Imberger, 1980): Once the velocity is computed from Eq. (46), the energy for mixing from velocity shear is compared to that required for lifting and accelerating the next layer down, and layers are combined if there is sufficient energy (Fig. 6f) where the billow length scale is δ KH = C KH u 2 b /g EH and δ KH = 2 C KH u b u b /g EH ; in this case the reduced gravity is computed from the difference between the bulk epilimnion and hypolimnion waters (Eq. 49), and C KH is a measure of the billow mixing efficiency.
Once energy from shear mixing is exhausted, the model checks the resultant density interface to see if it remains unstable to shear such that K-H billows would be expected to form, i.e. if the metalimnion thickness is less than the K-H length scale, δ KH . If this condition is met, a six-layer set is created about the thermocline to relax the stratification, set to have a linear density profile over δ KH (Fig. 8g), and the surface layer properties are updated.

Deep mixing
Mixing below the thermocline in lakes, in the deeper hypolimnion, is modelled using a characteristic vertical diffusivity, D Z = D ε +D m , where D m is a constant molecular diffusivity for scalars and D ε is the turbulent diffusivity. Three hypolimnetic mixing options are possible in GLM including (1) no diffusivity, D Z = 0, (2) a constant vertical diffusivity D Z over the water depth below the surface mixed layer, or (3) a derivation by Weinstock (1981) used in DYRESM, which is described as being suitable for regions displaying weak or strong stratification, whereby diffusivity increases with dissipation and decreases with heightened stratification. For the constant vertical diffusivity option, the coefficient C HYP is interpreted as the vertical diffusivity (m 2 s −1 ), i.e. D z = C HYP , and applied uniformly below the surface mixed layer. For the Weinstock (1981) model, the diffusivity varies depending on the strength of stratification and the rate of turbulent dissipation according to where C HYP in this case is the mixing efficiency of hypolimnetic TKE (∼ 0.8 in Weinstock, 1981), k TKE is the turbulence wavenumber defined below, and u * is defined as above. The stratification strength is computed using the Brunt-Väisälä (buoyancy) frequency, N 2 , defined for a given layer i as where ρ ref is the average of the layer densities. This is computed from layer three upwards, averaging over the span of five layers until the vertical density gradient exceeds a set tolerance. N 2 varies following an approximate normal distribution with height, centred at the height at which the centre of buoyancy is located and computed each time step from the first moment of the vertical N 2 distribution. Additionally, GLM estimates the vertical length scale associated with 1 standard deviation about the centre of the N 2 distribution, denoted δz σ . The diffusivity increases in line with the turbulent dissipation rate. This can be complex to estimate in stratified lakes; however, GLM adopts a simple approach as described in Fischer et al. (1979) in which a "net dissipation" is approximated by assuming that dissipation is in equilibrium with energy inputs from external forcing: which is expanded and calculated per unit mass as where ρ = 0.5 ρ 1 + ρ N LEV is the mean density of the water column. The work done by inflows is computed based on the flow rate and considers the depth to which the inflow plunges and the difference in density between the inflow water and layer into which it inserts, summed over all configured inflows (refer Sect. 2.7). These sources are normalized over the mass of water contained above the area of mixing. This is estimated as V N 2 , the fractional volume of the lake that is contained above the height that corresponds to 1 standard deviation below the centre of buoyancy and is therefore the volume of the lake over which 85 % of the N 2 variance is captured. The turbulence wavenumber, k TKE , is then estimated from where c wn is a coefficient. Since the dissipation is assumed to concentrate close to the level of strongest stratification, the "mean" diffusivity suggested by Eq. (52) is modified to decay exponentially within the layers as they increase their distance from the thermocline: where δz σ is used to scale the depth over which the mixing is assumed to decay below the bottom of the mixed layer, Once the diffusivity is approximated (either using a constant value or Eq. 57), the diffusion of any scalar, C (including temperature, salinity, and any water quality attributes), between two layers is numerically accounted for by the following mass transfer expressions:  Figure 8. Schematic depiction of layer changes during stratification and mixing. Consecutive panels show changes from (a) the initial layer and thermal profile, to (b) heating due to solar radiation, to (c) evaporative cooling, which creates (d) convective mixing followed by (e) a wind event causing stirring and (f) shear mixing across the thermocline. If the metalimnion remains unstable to shear it may be subjected to mixing from K-H billowing, which opens up the thermocline as depicted in panel (g).
where C is the weighted mean concentration of C for the two layers, and C is the concentration difference between them.
The smoothing function, f dif , is related to the diffusivity according to and the above diffusion algorithm is run once up the water column and once down the water column as a simple explicit method for capturing diffusion of mass to both the upper and lower layers. An example of the effect of hypolimnetic mixing on a hypothetical scalar concentration released from the sediment to the water column layers and accumulating in the hypolimnion is shown in Fig. 10.

Inflows and outflows
Aside from the surface fluxes of water described above, the water balance of a lake is controlled by inflows and outflows. Inflows can be specified as local run-off from the surrounding (dry) lake domain (Q R described separately above; Eq. 8), rivers entering at the surface of the lake that will be buoyant or plunge depending on their momentum and density (Sect. 2.7.1), or submerged inflows (including groundwater) that enter at depth (Sect. 2.7.2). Four options for outflows are included in GLM. These include withdrawals from a specified depth (Sect. 2.7.3), adaptive offtake (Sect. 2.7.4), vertical groundwater seepage (Sect. 2.7.5), and river outflowoverflow from the surface of the lake (Sect. 2.7.6). Any number of lake inflows and outflows can be specified, and, except for the local run-off term, all are applied at a daily time step. Depending on the specific settings of each, these water fluxes can impact the volume of the individual layers, V i , and the overall lake volume (Eq. 4). Inflows have a prescribed composition (temperature, salinity, and scalars), except local run-off, which is assumed to be at air temperature with zero salinity.

River inflows
As water from an inflowing river connects with a lake or reservoir environment, it will form a positively or negatively buoyant intrusion depending on the density of the incoming river water in the context of the water column stratification. As the inflow progresses towards insertion, it will entrain water at a rate depending on the turbulence created by the inflowing water mass (Fischer et al., 1979). For each configured inflow the entrainment coefficient, E inf , is computed based on the bottom drag being experienced by the inflowing water, C D inf , and the water stability using the approximation given in Imberger and Patterson (1981) as written in Ayala et al. (2014): where the inflow Richardson number, Ri inf , characterizes the stability of the water in the context of the inflow. Imberger and Patterson (1981) derived a simple estimate of Ri inf based on the drag coefficient by assuming the velocity (and Froude number) is typically small and considering the channel geometry, which is adapted in GLM as where α inf is the stream half-angle assuming an approximate triangular cross section, and inf is the angle of the slope of the inflow thalweg relative to horizontal in the region where it meets the waterbody (Fig. 11). Therefore, using Eqs. (60) and (61) bottom roughness can be used to parameterize the characteristic rate of entrainment as it enters the waterbody. On entry, the inflow algorithm captures two phases: first, the inflowing water crosses the layers of the lake until it reaches a level of neutral buoyancy, and second, it then undergoes insertion. In the first part of the algorithm, the daily inflow parcel is tracked down the lakebed and its mixing with layers is updated until it is deemed ready for insertion. The initial estimate of the intrusion thickness, z inf 0 , is computed as in Antenucci et al. (2005) and Ayala et al. (2014): where Q inf 0 = f inf Q inf x is the inflow discharge entering the domain based on the data provided as a boundary condition, Q inf x , and g inf is the reduced gravity of the inflow as it enters. Here, ρ inf is the density of the inflow computed from the supplied inflow properties of temperature and salinity (T inf x , S inf x ), and ρ s is the density of the surface layer.
If the inflowing water is deemed to be positively buoyant (ρ inf < ρ s ) or the model only has one layer (N LEV = 1), then the inflow water over the daily time step is added to the surface layer volume ( V N LEV = Q inf 0 t d ), and h S is updated accordingly. Otherwise, this inflow volume is treated as a parcel which travels down through the lake layers, and its properties are subsequently incremented over each time step, j (currently daily), until it inserts. The thickness of an inflow parcel increases over each increment due to entrainment, assuming where z inf j is the inflow thickness and x inf j is the distance travelled by the inflowing water parcel over the j th time step. The distance travelled is estimated based on the change in the vertical height of the inflow, δz inf , and the angle of the inflow river, φ inf , as given by The vertical excursion for the step is approximated as the difference between its starting height and the bottom of the nearest layer that it sits above, h i j −1 , such that δz inf j = h S − z inf j −1 − h i j −1 , where z inf j −1 is the depth of the inflow from the surface at the start of the time step, and this is subsequently updated from z inf j = z inf j −1 + x inf j sin inf . The average velocity of the inflow parcel is updated based on the incoming flow rate from where the denominator links the relationship between inflow height and channel width in order to define the crosssectional area of the flow. This velocity is used to estimate the timescale of transport of the parcel (δt d = x inf j u inf j ). Following the conservation of mass, the flow is estimated to increase according to Fischer et al. (1979) (see also Antenucci et al., 2005): whereby Q inf j is removed from the volume of the corresponding layer, i j , and added to the previous time step inflow Q inf j −1 to capture the entrainment effect on the inflow for the next increment. The properties associated with Q inf j are assumed to match those of the i j layer and mixed into the inflow parcel to update temperature, salinity, and density, ρ inf j . The inflow travel algorithm (Eqs. 62-67) increments through j until the density of the inflow first reaches its depth of neutral buoyancy: ρ inf j ≤ ρ i j . Once this condition is met, the insertion depth is defined as z inf ins I , its density as ρ inf ins I , and the second part of the algorithm then creates a new layer of thickness dependent on the inflow's volume at that time, Q inf ins I c secday , which includes the successive layer additions from entrainment; Eq. (67). Since a new inflow parcel is created each day and the user may configure multiple inflows, N INF , a complex set of parcels being tracked via Eqs. (60)-(67), and a queue of new layers to be inserted are created. Following the creation of a new layer for an inflow parcel, N LEV is incremented and all layer heights above the new layer are updated, paying attention to the lake hypsography. The new inflow layer is then subject to the thickness limit criteria within the layer-limitchecking routine and may amalgamate with adjacent layers or be divided into thinner layers.
Aside from importing mass into the lake, river inflows also contribute turbulent kinetic energy that may dissipate in the hypolimnion, as discussed in Sect. 2.6.2 (e.g. see Eq. 55), and they contribute to the scalar transport in the water column Inflow inserted as new layer at depth of neutral bouyancy when E n tr a in m e n t, Figure 11. Schematic showing inflow insertion depth, entrainment, E inf , slope, inf , and bottom slope angle, α inf , of an inflowing river, I , entering at a flow rate of Q inf 0 and an estimated starting thickness of z inf 0 .
by adding mass contained within the inserted inflow parcels and contributing to mixing of properties via entrainment as described above (Fig. 12a); see also Fenocchi et al. (2017).

Submerged inflows
Submerged inflows are inserted at the user-specified depth, h inf , with zero entrainment by utilizing the second part of the algorithm described in Sect. 2.7.1. Once the submerged inflow volume is added as a new layer it may then be mixed with adjacent layers (above or below) depending on the density difference and layer thickness criteria (Fig. 12b). This option can be used across one or more inflow elevations to account for groundwater input to a lake or for capturing a piped inflow, for example.

Withdrawals
Outflows from a specific depth can include outlets from a dam wall offtake, other piped withdrawal, or removal of water that may be lost due to groundwater recharge. For a stratified water column, the water will be removed from the layer corresponding to the specified withdrawal height, h outf , and layers above or below depending on the strength of discharge and stability of the water column. Accordingly, the model assumes an algorithm in which the thickness of the withdrawal envelope is dependent on the internal Froude (F r) and Grashof (Gr) numbers and the parameter R (see Fischer et al., 1979;Imberger and Patterson, 1981): where W outf , L outf , and A outf are the width, length, and area of the lake at the outlet elevation, and D 2 outf is the vertical diffusivity averaged over the layers corresponding to the withdrawal thickness, δ outf (described below). To calculate the width and length of the lake at the height of the outflow, it is assumed, firstly, that the lake shape can be approximated as an ellipse and, secondly, that the ratio of length to width at the height of the outflow is the same as that at the lake crest. The length of the lake at the outflow height, L outf , and the lake width, W outf , are given by where A outf is the area of the lake at the outflow height, L crest is the length, and W crest the width of the lake at the crest height.
The thickness of the withdrawal layer is calculated depending on the value of R (Fischer et al., 1979) such that If stratification is apparent near h outf , either above or below this elevation, then the thickness computed in Eq. (73) may not be symmetric about the offtake level (Imberger and Patterson, 1981); therefore the algorithm separately computes the thickness of the withdrawal layer above and below, denoted δ outf top and δ outf bot , respectively. The Brunt-Väisälä frequency is averaged over the relevant thickness, N 2 outf , and calculated as where ρ outf is the density of the layer corresponding to the height of the withdrawal, i outf , and ρ i is the density of the water column at the edge of the withdrawal layer, as determined below. The flow of water taken from each layer influenced by the withdrawal, Q outf i , either above or below the layer of the outlet elevation requires identification of the uppermost and lowermost layer indices influenced by the outflow, denoted i top and i bot . Once the layer range is defined, Q outf i is computed for the layers between i outf and i top and between i outf and i bot by partitioning the total outflow using a function to calculate the proportion of water withdrawn from any layer that fits the region of water drawn in a given time (Q outf i = f f outf Q outf x , h i , h i−1 , h outf , δ outf bot , δ outf top ; see Imberger and Patterson, 1981, Eq. 65). Given that users configure any height for a withdrawal outlet and flow rates of variable strength, the upper (h outf + δ outf top ) and lower (h outf − δ outf bot ) elevation limits computed by the algorithm are limited to the lake surface layer or bottom layer. Once computed, the volumes are removed from the identified layer set, and their height and volumes updated accordingly. Q outf i is constrained within the model to ensure no more than 90 % of a layer can be removed in a single time step. Depending on the fractional contribution from each of the layers from which the water is withdrawn, the water taken will have the associated weighted average of the relevant scalar concentrations (T outf , S outf , C outf ), which are reported in the outlet file for the particular withdrawal. This routine is repeated for each withdrawal considered, denoted O, and the model optionally produces a summary file of the combined outflow water and its properties.

Adaptive offtake dynamics
For reservoir applications, a special outflow option has been implemented that extends the dynamics in Sect. 2.7.3 to simulate an adaptive offtake or selective withdrawal. This ap-proach is used for accommodating flexible reservoir withdrawal regimes and their effects on the thermal structure within a reservoir. For this option, a target temperature is specified by the user and GLM identifies the corresponding withdrawal height within a predefined (facility) range to meet this target temperature during the runtime of the simulation; i.e. the withdrawal height adaptively follows the thermal stratification in the reservoir. The target temperature can be defined as a constant temperature or a time series (via a *.csv file), such as a measured water temperature from an upstream river that could be used to plan environmental releases from the reservoir to the downstream river. The selected height of the adaptive offtake is printed out in a *.txt file to assist reservoir operation. In addition to the basic adaptive offtake function, GLM can also simulate withdrawal mixing; i.e. water from the adaptive offtake is mixed with water from another predefined height (e.g. the bottom outlet). For this option, the discharges at both locations need to be predefined by the user (via the standard outflow *.csv files) and GLM chooses the adaptive withdrawal from a height at which the water temperature is such that the resulting mixing temperature meets the target temperature. This withdrawal mixing is a common strategy in reservoir operation in which deep water withdrawal and temperature control are required simultaneously to prevent deleterious downstream impacts. An example of the adaptive offtake function with and without withdrawal mixing, assuming a constant water temperature of 14 • C for the outflow water, shows that GLM is able to deliver a constant outflow temperature of 14 • C during the stratified period (Fig. 13). In winter when the water column is cooler than 14 • C, the model withdraws surface water. The adaptive offtake functionality can be used in a standalone mode or also linked to the dissolved oxygen concentration (when operated with the coupled water quality model AED2; see Sect. 4). In the latter case, the effect of the withdrawal regime on the oxygen dynamics in the hypolimnion can be simulated (see Weber et al., 2017). In this setting, the simulated hypolimnetic dissolved oxygen concentration at a specified height is checked against a user-defined critical threshold. If the hypolimnetic oxygen falls below the critical threshold, the height of the adaptive offtake will be automatically switched to a defined height (usually deep outlets in order to remove the oxygen-depleted water) to withdraw water from this layer until the oxygen concentrations have recovered.

Seepage
Seepage of water from the lake can also be configured within the model, for example, as might be required in a wetland simulation or for small reservoirs perched above the water table that experience leakage to the soil below. The seepage rate, Q seepage , can be assumed constant or dependent on the overlying lake head:

Option 2: Darcy flux based on water height
where G is the seepage rate (m day −1 ), K seep is the sediment hydraulic conductivity (m day −1 ), and δz soil is an assumed sediment thickness over which the seepage is assumed to occur. The water leaving the lake is treated as a "vertical withdrawal" whereby the water exits via the bottommost layer(s), and the amount V G = Q seepage t d is generally all taken from the bottom-most layer (i = 1); however, it is constrained within the model to ensure no more than 90 % of the layer can be reduced in any one time step; where V G > 0.9V i=1 , the routine sequentially loops up through the above layers until enough lake volume has been identified to cover the seepage demand. Once the individual layer volumes are incremented due to the seepage flux, V G i , the heights of all layers (h 1 : h s ) are recomputed based on the hypsographic curve using h i = f [V i ]. Where seepage reduces the lake below 0.05 m, the lake becomes dry and will continue to have zero volume until there are new inputs from rain or inflows (e.g. Fig. 9a).

Overflows
Once the lake volume exceeds the maximum volume, the excess water is assumed to leave the domain as an overflow. The flow rate, Q ovfl , is computed based on the interim volume, V * S , prior to the end of the daily time step, where V * and h S is the cumulative change in the daily water level over the day. Users can optionally also specify a crest elevation which lies below the elevation of maximum lake volume and support a rating curve linking the height of water above the crest level with the overflow volume: where h * S is the interim update to the water surface height prior to the overflow computation, C D weir is a coefficient related to the drag of the weir, W weir is the width of the weir crest, and h crest is the height of the crest level. The overflow rate is then computed as the sum of the flow over the weir crest and the volume of water exceeding the volume of the domain: If no crest is configured below the maximum lake height, then Eq. (77) assumes Q weir = 0.

Wave height and bottom stress
Resuspension of sediment from the bed of lakes depends on the stresses created by water movement across the lake bottom. Wind-induced resuspension, in particular, is sporadic and occurs as the waves at the water surface create oscillatory currents that propagate down to the lakebed and exceed a critical threshold. The wave climate that exists on a lake can be complex and depend on the fetch over which the wind has blown, the time period over which the wind has blown, and complicating factors such as wind sheltering and variations in bottom topography. The horizontally averaged nature of GLM means that only a single set of wave characteristics across the entire lake surface can be computed for a given time step and these are assumed to be at steady state. Note that GLM does not predict resuspension and sediment concentration directly, but computes the bottom shear stress for later use in coupled sediment and water quality modules. Since each layer has a component that is considered to overlay sediment (Sect. 4), the stress experienced at the sediment-water interface is able to approximated as a function of depth in relation to the surface wave climate. The model can therefore identify the depth range and areal extent to which there is potential for bed-sediment resuspension to occur, i.e. by computing the area of the lake over which the bed shear stress exceeds some critical value required for resuspension.
The model estimates surface wave conditions using a simple, fetch-based, steady-state wave model (Laenen and Le-Tourneau, 1996;Ji, 2008). The average wave geometry (wave period, significant wave height, and wavelength) is predicted Figure 13. Adaptive offtake reservoir simulation; water temperatures of the adaptive offtake model assuming a constant target temperature of 14 • C (a, b) without and (c, d) with mixing with the bottom outlet withdrawal. The black dashed line (a, c) represents the height range of the variable withdrawal facility (AOF) and the magenta lines the adaptive offtake and second withdrawal height (here: bottom outlet). In the scenario with the second withdrawal activated (c), the bottom outlet was periodically opened during flooding conditions. Panels (b) and (d) indicate where the actual withdrawal temperature (DrawTemp, T outf ) was able to meet the target (TargetTemp). based on the wind speed and fetch over which the waves develop (Fig. 14), whereby the average fetch is approximated in the one-dimensional model formulation from the lake area, Based on these properties the orbital wave velocity at the surface can be translated down the depth of the water column such that in the ith layer it is calculated as (Sheng and Lick, 1979) For each layer, the total shear stress experienced at the lakebed portion of that layer (equivalent in area to A i −A i−1 ) is calculated from where U m is the mean layer velocity, which for simplicity is assumed based on the velocity estimate made during the mixing calculations (Eq. 44) in the surface mixed layer such that The friction factors depend upon the characteristic particle diameter of the lake bottom sediments, δ ss , and the fluid velocity. For the current-induced stress, we compute f c i = 0.24/ log 12z avg /2.5δ ss z i and for waves (Kleinhans and Grasmeijer, 2006) where δ ss z i is specific for each layer i, depending on which sediment zone it overlays (see Sect. 4). The current-and wave-induced stresses at the lake bottom manifest differently within the lake, as demonstrated in Fig. 15 for a shallow lake.

Code organization and model operation
Aside from the core water balance and mixing functionality, the model features numerous options and extensions in order to make it a fast and easy-to-use package suitable for a wide range of contemporary applications. Accommodating these requirements has led to the modular code structure outlined in Fig. 16. The model is written in C, with a Fortran-based interface module to link with the Fortran-based water quality modelling libraries described in Sect. 4. The model compiles with GCC, GFortran, and commercial compilers, with support for Windows, OS X, and Linux. The model may also be compiled as a library, termed libGLM, that can be called as a plug-in to other models (e.g. see Sect. 5.4). Whilst the model is not object oriented, users may easily customize the specific modules described in Sect. 2 by adding or extending options for alternate schemes or functions.
To facilitate the use of the model in teaching environments and for users with limited technical support, it may be operated without any third-party software, as the input files consist of "namelist" (nml) text files for configuration and csv files for meteorological and flow time series data (Fig. 17). The outputs from predictions are stored into a structured NetCDF file, which can be visualized in real time through the simple inbuilt plotting library (libplot) or may be opened for post-processing in MATLAB, R, or any other tool supporting the open NetCDF format (see Sect. 5.1). Parameters and configuration details are input through the main glm.nml text file (Fig. 17) and default parameters and their associated descriptions are outlined in Table 1.

Dynamic coupling with biogeochemical and ecological model libraries
Beyond modelling the vertical temperature distribution, the water, ice, and heat balance, and the transport and mixing in a lake, the model has been designed to couple with biogeochemical and ecological model libraries. Currently the model is distributed pre-linked with the AED2 simulation library (Hipsey et al., 2019c) and the Framework for Aquatic Biogeochemical Models (FABM; Bruggeman and Bolding, 2014). Through connection with these libraries, GLM creates a set C of scalar variables, where C ∈ C, which resolve the vertical profiles and mass balance of turbidity, oxygen, nutrients, phytoplankton, zooplankton, pathogens, and other water quality variables of interest. Documentation of these models is beyond the scope of the present paper; however, two features associated with the coupling are highlighted here as relevant to managing physical-ecological interactions. Firstly, the model is designed to allow for a user-defined number of sediment zones that span the depth of the lake. Using this approach, the current set-up allows for depthdependent sediment properties, both for physical properties such as roughness or sediment heat flux (as outlined in previous sections) and also biogeochemical properties such as sediment nutrient fluxes and benthic ecological interactions. Since the GLM layer structure is flexible over time (i.e. layer heights are not fixed), any interactions between the water and sediment-benthos must be managed at each time step. The model supports the disaggregation and/or aggregation of layer properties for mapping individual water layers to one or more sediment zones (Fig. 18). The weightings provided by each layer to the sediment are based on the relative depth overlap of a layer with the depth range of the sediment zone, with the heights of zone boundaries denoted h z . This approach makes the model suitable for long-term assessments of wetland, lake, and reservoir biogeochemical budgets, as is required for C, N, and other attribute balances (Stepanenko et al., 2016).
Secondly, the water quality modules feed back to GLM properties related to the water and/or heat balance. Feedback options are included for external libraries to provide water density additions and to modify bottom friction, f w , the light attenuation coefficient, K w , solar shading, f SW , and rainfall interception, f R .

Workflow tools for integrating GLM with sensor data and supporting models
The GLM has been designed to support the integration of large volumes of data coming from instrumented lakes, including many GLEON sites. These data consist of highfrequency and discrete time series observations of hydrologic fluxes, meteorology, temperature, and water quality (e.g. Hamilton et al., 2015). To facilitate research that re- volume of the lake at the top of the surface layer interim calculation of the volume of the lake at the top of the surface layer m 3 variable used to estimate lake volume prior to overflow calculation V N 2 a fractional volume of the lake that contains 85 % of N 2 variance m 3 variable computed as the volume of the lake above the height which is 1 standard deviation (δz σ ) below the height at the centre of buoyancy volume of the layer below the surface mixed layer-epilimnion depth that an inflow parcel associated with inflow I inserts m from water surface variable depth from the surface at which an inflow reaches its level of neutral buoyancy z i thickness of the ith layer m variable z k−1 thickness of the layer below the epilimnion m variable z min minimum layer thickness m 0.5 should be estimated relative to lake depth; set in &glm_setup  waveband 2, white ice light extinction m −1 20.0 K b1 waveband 1, blue ice light extinction m −1 1.5 K b2 waveband 2, blue ice light extinction m −1 20.0 K s1 waveband 1, snow light extinction m −1 6 K s2 waveband 2, snow light extinction m −1 20 K snow molecular heat conductivity of snow J m −1 s −1 • C −1 computed dependent on snow density according to Eq. (35) K white molecular heat conductivity of white ice J m −1 s −1 • C −1 2.3 K blue molecular heat conductivity of blue ice J m −1 s −1 • C −1 2.0 K water molecular heat conductivity of water J m −1 s −1 • C −1 0.57 K air molecular heat conductivity of air J m −1 s −1 • C −1 2.8 × 10 −3 reported as 0.1 kJ m −1 h −1 K −1

TVA (1972)
K soil heat conductivity of soil-sediment J m −1 s −1 • C −1 1.2 varies from 0.25 for organic soil to 2.9 for inorganic particles K seep hydraulic conductivity of soil below the lake m day −1 configurable set in &outflows  rate of a single water inflow provided by the user as input to the model m 3 ss −1 time series input based on the non_avg flag in &glm_setup; the supplied value at the current time step or the average of the current and past time step is used, depending on whether the daily data are referenced from midday or midnight Q outf x rate of a single water outflow provided by the user as input to the model m 3 s −1 time series input Q inf 0 rate of a single water inflow prior to the inflow entering the lake Q outf rate of a single water outflow exiting the lake flow rate of water being extracted from the ith layer R x rainfall rate supplied in the input file m day −1 time series input user-supplied rainfall rate R L rainfall intensity threshold before run-in occurs m day −1 0.04 depends on land slope and soil type R snow critical rainfall rate incident on snow that controls densification m day −1 configurable set in &snow ice R F rainfall rate entering the water column m s −1 computed Eq. (6) * R F rainfall rate incident on the ice-snow layer m s −1 computed  quires running the model using these data sources, we have created GLM interfaces in the R and MATLAB analysis environments. These tools support user-friendly access to the model and include routines that streamline the process of calibrating models or running various scenarios. In addition, for assessment of lake dynamics in response to catchment or climatic forcing, it is desirable to be able to connect GLM with other model platforms associated with surface and groundwater simulation and weather prediction (Read et al., 2016).

R and MATLAB libraries for model set-up and post-processing
The R and MATLAB scientific languages are commonly used in aquatic research, often as part of automated modelling and analysis workflows. GLM has a client library for both, and these tools are shared freely online. The R package is called "glmtools" and the MATLAB library is called "GLMm"(available via the GLM website). Both tools have utilities for model output pre-and post-processing. The preprocessing components can be used to format and modify data inputs and configuration files and define options for how GLM executes. Post-processing tools include visualizations of simulation results (as shown in the figures above), comparisons to field observations, and various evaluations of model performance.

Utilities for assessing model performance, parameter identification, and uncertainty analysis
In order to compare the performance of the model for various types of lakes, numerous metrics of model performance are relevant. These include simple measures like surface or bottom temperature and ice thickness. It is also possible to assess the model's performance in capturing higher-order metrics relevant to lake dynamics, including Schmidt stability, thermocline depth, and ice on-off dates (see also Bruce et al., 2018, for a detailed assessment of the model's accuracy across a wide diversity of lakes across the globe). With particular interest in the model's ability to interface with highfrequency sensor data for the calculation of key lake stability metrics (Read et al., 2011), continuous wavelet transform comparisons are also possible (Kara et al., 2012), allowing for the assessment of the timescales over which the model is able to capture the observed variability within the data.  Figure 18. Schematic of a lake model layer structure (indicated by layers i = 1 : 7), in conjunction with five sediment "zones" (Z1-Z5) activated when benthic_mode=2. The dynamically varying layer structure is remapped to the fixed sediment zone locations at each time step in order for the sediment zone to receive the average overlying water properties and for the water to receive the appropriate information from benthic-sediment variables.
As part of the modelling process, it is common to adjust parameters to get the best fit with available field data and, as such, the use of a Bayesian hierarchical framework in the aquatic ecosystem modelling community has become increasingly useful (e.g. Zhang and Arhonditsis, 2009;Romarheim et al., 2015). Many parameters described throughout Sect. 2 are physically based descriptions in which there is relatively little variation (Bruce et al., 2018), thereby reducing the number of parameters that remain uncertain. For others, however, their variation reflects the imperfect formulation of some processes that are not completely described numerically. Therefore, within MATLAB, support scripts for GLM to work with the Markov chain Monte Carlo (MCMC) code outlined in Haario et al. (2006) can be used to provide improved parameter estimates and uncertainty assessment ( Fig. 19; see also Huang et al., 2017). Example setups for use of GLM within the PEST (Parameter Estimation Tool) have also been developed, giving users access to a wide range of assessment methodologies. The PEST framework allows for the calibration of complex models using highly parameterized regularization with pilot points (Doherty, 2015). Sensitivity matrices derived from the calibration process can also be utilized in linear and non-linear uncertainty analysis.

Operation in the cloud: GRAPLEr
Questions relevant to land use and climate changes are driving scientists to develop numerous scenarios for how lake ecosystems might respond to changing exogenous drivers.
An important approach to addressing these questions is to simulate lake or reservoir physical-biological interactions in response to changing hydrology, nutrient loads, or meteorology and then infer consequences from the emergent properties of the simulation, such as changes in water clarity, extent of anoxia, mixing regime, or habitability to fishes . Often, it takes years or even decades for lakes to respond fully to changes in exogenous drivers, requiring simulations to recreate lake behaviour over extended periods. While most desktop computers can run a decade-long, low-resolution simulation in less than 1 min, high-resolution simulations of the same extent may require minutes to hours of processor time. When questions demand hundreds, thousands, or even millions of simulations, the desktop approach is no longer suitable.
Through access to distributed computing resources, modellers can run thousands of GLM simulations in the time it takes to run a few simulations on a desktop computer. Collaborations between computer scientists in the Pacific Rim Applications and Grid Middleware Assembly (PRAGMA) and GLEON have led to the development of GRAPLEr (GLEON Research and PRAGMA Lake Expedition in R), a software written in R that enables modellers to distribute batches of GLM simulations to pools of computers (Subratie et al., 2017). Modellers use GRAPLEr in two ways: by submitting a single simulation to the GRAPLEr Web service, along with instructions for running that simulation under different climate scenarios, or by configuring many simulations Figure 19. Depiction of parameter uncertainty for a GLM simulation of Lake Kinneret, Israel, following calibration against observations (green circles) via MCMC for (a) epilimnion temperature, (b) hypolimnion temperature, (c) thermocline depth, and (d) Schmidt number. The black line indicates the 50th percentile likelihood of the prediction, and the grey bands depict the 40th, 60th, and 80th percentile. on the user's desktop computer and then submitting them as a batch to the Web service. The first approach provides a high degree of automation that is well suited to training and instruction, and the second approach has the full flexibility often needed for research projects. In all approaches, GRAPLEr converts the submitted job to a script that is used by the scheduling programme HTCondor (Thain et al., 2005) to distribute and manage jobs among the computer pool and ensure that all simulations run and return results. An iPOP overlay network (Ganguly et al., 2006) allows the compute services to include resources from multiple institutions and cloud-computing services.
GRAPLEr's Web service front end shields the modeller from the compute environment, greatly reducing the need for modellers to understand distributed computing; they therefore only need to install the R package, know the URL of the GRAPLEr Web service, and decide how the simulations should be set up.

Integration with catchment and climate models
GLM simulations may be coupled with catchment models, such as the Soil Water Assessment Tool (SWAT) or similar catchment models, simply by converting the catchment model output into the inflow file format via conversion scripts (e.g. Bucak et al., 2018). Similarly, scripts exist for coupling GLM with the Weather Research Forecasting (WRF) model or similar climate models for the specification of the meteorological input file from weather prediction simulations (e.g. Hansen et al., 2017).
The above coupling approaches require the models to be run in sequence. For the simulation of lake-wetlandgroundwater systems, however, two-way coupling is required to account for the flow of water into and out of the lake throughout the simulation. For these applications, the interaction has been simulated using GLM coupled with the 3-D groundwater flow model FEFLOW (https://www. mikepoweredbydhi.com/products/feflow, last access: 16 December 2018). For this case, the GLM code is compiled as a dynamic link library (DLL), termed libGLM, and loaded into FEFLOW as a plug-in module. The coupling between GLM and FEFLOW is implemented using a one-step lag between the respective solutions of the groundwater and lake models. This approach in most simulations does not introduce significant error; however, error can be assessed and reduced using smaller time step lengths. The GLM module was designed to accommodate situations of variable lake geometry by using a dry-lake-wet-lake approach, whereby dry-lake areas are defined as those above the current lake level and wet-lake areas as below the current lake level. Different boundary types in FEFLOW are assigned to dry-lake and wet-lake areas (Fig. 20). The calibration of such coupled models is often complex given the large number of parameters and sensitivities when different sources of information are utilized (for example, flow and water level measurements). The FEFLOW-GLM coupling structure allows for a relatively straightforward integration with PEST based on existing FEFLOW workflows.

GLM as a tool for teaching environmental science and ecology
Environmental modelling is integral for understanding complex ecosystem responses to anthropogenic and natural drivers and also provides a valuable tool for engaging students learning environmental science (Carey and Gougis, 2017). Previous pedagogical studies have demonstrated that engaging students in modelling provides cognitive benefits, enabling them to build new scientific knowledge and conceptual understanding (Stewart et al., 2005;Schwarz et al., 2009). For example, modelling forces students to analyse patterns in data, create evidence-based hypotheses for those patterns, make their hypotheses explicit, and develop predic- tions of future conditions (Stewart et al., 2005). As a result, the U.S. National Research Council has recently integrated modelling into the Next Generation Science Standards, which provide recommendations for primary and secondary school science pedagogy in the United States (NRC, 2013). However, it remains rare for undergraduate and graduate science courses to include the computer-based modelling that environmental scientists need to manage natural ecosystems.
A teaching module for the use of GLM within undergraduate and graduate classrooms has been developed to explore lake responses to climate change (Carey and Gougis, 2017). The GLM module, called the "Climate Change Effects on Lake Temperatures", teaches students how to set up a simulation for a model lake within R. After they are able to successfully run their lake simulations, they force the simulation with climate scenarios of their own design to examine how lakes may change in the future. To improve computational efficiency, students also learn how to submit, retrieve, and analyse hundreds of model simulations through distributed computing overlay networks embedded via the GRAPLEr interface (Sect. 5.3). Hence, students participating in the module learn computing and quantitative skills in addition to improving their understanding of how climate change affects lake ecosystems.
Initial experiences teaching GLM as well as pre-and post-assessments indicate that participation in the module improves students' understanding of lake responses to climate change (Carey and Gougis, 2017). By modifying GLM boundary condition data and exploring model output, students are able to better understand the processes that control lake responses to altered climate and improve their predictions of future lake change. Moreover, the module exposes students to computing and modelling tools not commonly experienced in most university classrooms, building competence with manipulating data files, scripting, creating figures and other visualizations, and statistical and time series analysis; these are all skills that are transferrable to many other applications.

Conclusions
As part of GLEON activities, the emergence of complex questions about how different lake types across the world are responding to climate change and land use change has created the need for a robust, accessible community code suitable for a diverse range of lake types and simulation contexts. Here, GLM is presented as a tool that meets many of the needs of network participants with suitability for a wide array of lake types and sizes, whilst also meeting the need for a distributed simulation across tens to thousands of lakes as is required for regional-and global-scale assessments (e.g. Kirillin et al., 2011). Recent examples have included an application of the model for assessing how the diversity of > 2000 lakes in a lake-rich landscape in Wisconsin respond to climate, including projected warming (Read et al., 2014;Winslow et al., 2017). Given its computationally efficient nature, it is envisioned that GLM can be made available as a library for use within in land surface models (e.g. the Community Land Model, CLM), allowing for an improved representation of lake dynamics in regional hydrological or climate assessments. With further advances in the degree of resolution and scope of Earth system models, we further envisage GLM as an option suitable to be embedded within these models to better allow for the simulation of lake stratification, air-water interaction of momentum and heat, and also biogeochemically relevant variables associated with contemporary questions about greenhouse gas emissions such as CO 2 , CH 4 , and N 2 O.
Since the model is one-dimensional, it assumes no horizontal variability in the simulated water layers and users must therefore ensure that their application of the model is suited to this simplifying assumption. For stratified systems, the parameterization of mixing due to internal wave and inflow intrusion dynamics is relatively simple, making the model ideally suited to longer-term investigations ranging from weeks to decades (depending on the domain size) and for coupling with biogeochemical models to explore the role that stratification and vertical mixing play in lake ecosystem dynamics. However, the model can also be used for shallow lakes, ponds, and wetland environments in which the water column is relatively well mixed. In cases in which the assumption of one-dimensionality is not met for a particular lake application, a two-or three-dimensional model may be preferred. This paper has focused on a description of the hydrodynamic model, but we highlight the fact that the model is a platform for coupling with advanced biogeochemical and ecological simulation libraries for water quality prediction and integrated ecosystem assessments. As with most coupled hydrodynamic-ecological modelling platforms, GLM handles the boundary conditions and transport of variables simulated within these libraries, including the effects of inflows, vertical mixing, and evapo-concentration. Whilst the interface to these libraries is straightforward, the Lagrangian approach adopted within GLM for simulation of the water column necessitates the adoption of sediment zones on a static grid that is independent from the water column numerical grid.
More advanced workflows for operation of the model within distributed computing environments and with data assimilation algorithms is an important application when used within GLEON capabilities related to high-frequency data and their interpretation. The 1-D nature of the model makes the runtimes modest and therefore the model suitable for application within more intensive parameter identification and uncertainty assessment procedures. This is particularly relevant to the needs of network participants to expand model configurations to further include biogeochemical and ecological state variables. It is envisioned that continued application of the model will allow us to improve parameter estimates and ranges, and this will ultimately support other users of the model in identifying parameter values and assigning parameter prior distributions. Since many of the users the model is intended for may not have access to the necessary cyber-infrastructure, the use of GLM with the open-source GRAPLEr software in the R environment provides access to otherwise unavailable distributed computing resources. This has the potential to allow non-expert modellers within the science community to apply good modelling practices by automating boundary condition and parameter sensitivity assessments, with technical aspects of simulation management abstracted from the user.
Finally, the role of models in informing and educating members of the network and the next generation of hydrologic and ecosystem modellers has been identified as a critical element of synthesis activities and supporting crossdisciplinary collaboration . Initial use of GLM within the classroom has shown that teaching modules integrating GLM into classes improves students' understanding of lake ecosystems. Data availability. The five example lakes used to demonstrate the model operation are described along with model input files (and the associated hydrologic and meteorological forcing data) within the GitHub repository: https://github. com/AquaticEcoDynamics/GLM_Examples (last access: 14 January 2019) (https://doi.org/10.5281/zenodo.2538489; Hipsey et al., 2019b).
The direct beam (DB) radiation on a horizontal surface at ground level on a clear day is given bŷ φ DB = 0.9662φ ETR T rayleigh T ozone T mix T watvap T aerosol cos [ zen ] , (A19) and the atmospheric scattering (AS) component iŝ φ AS = 0.79φ ETR T ozone T mix T watvap T aa cos [ zen ] . (A20) The total irradiance hitting the surface (W m −2 ) is thereforê .
The albedo is computed for the sky as Appendix B: Non-neutral bulk-transfer coefficients The iterative procedure used in this analysis to update-correct the bulk-transfer coefficients based on atmospheric conditions is conceptually similar to the methodology discussed in detail in Launiainen and Vihma (1990). The first estimate for the neutral drag coefficient, C DN , is specified as a function of wind speed as it is commonly observed to increase with U 10 . This is modelled by first estimating the value referenced to 10 m of height above the water from Option 1 : Francey and Garratt (1978), Hicks (1972) 0.001 U 10 ≤ 5 0.001 (1 + 0.07 (U 10 − 5)) U 10 > 5 Option 2 : Babanin and Makin (2008) 1.92 × 10 −7 U 3 10 + 0.00096, and then computing the Charnock formula with the smooth flow transition (e.g. Vickers et al., 2013): where a is the Charnock constant, and here u * is the approximated friction velocity of the atmosphere near the surface ( C DN−10 U 2 10 ) initially estimated using Eq. (B1). The drag is recomputed using where κ is the von Karman constant (Fig. B1). Note that the neutral humidity-temperature coefficient, C HWN−10 , is held constant at the user-defined C H value and is assumed not to vary with wind speed.  During stable stratification (L > 0) they are assumed to take a form modified from Hicks (1976).
where the N subscript denotes the neutral value and X signifies either D, H , or E for the transfer coefficient and o, θ , or q for the roughness length scale. Inclusion of the stability functions into the substitution and some manipulation (Imberger and Patterson, 1990;Launianen and Vihma, 1990) yields the transfer coefficients relative to these neutral values. Hicks (1975) and Launianen and Vihma (1990) suggested an iterative procedure to solve for the stability-corrected transfer coefficient using Eq. (B10) based on some initial estimate of the neutral values (as input by the user). The surface flux is subsequently estimated according to Eqs. (22)-(23) and used to provide an initial estimate for L (Eq. B4). The partially corrected transfer coefficient is then recalculated and so the cycle goes. Strub and Powell (1987) and Launiainen (1995) presented an alternative based on estimation of the bulk Richardson number, Ri B , defined as www.geosci-model-dev.net/12/473/2019/ and related as a function of the stability parameter, z/L, according to where it is specified that C HN = C EN = C HEN . Figure B2 illustrates the relationship between the degree of atmospheric stratification (as described by both the bulk Richardson number and the Monin-Obukhov stability parameter) and the transfer coefficients scaled by their neutral value.