Lakes and other lentic (standing) waters support extensive ecosystem services such as water supply, flood mitigation, hydropower, aesthetic and cultural benefits, as well as fisheries and biodiversity (Mueller et al., 2016). Lakes are often considered to be "sentinels of change", providing a window into the sustainability of activities in their associated river basins (Williamson et al., 2009). They are also particularly susceptible to impacts from invasive species and land use development, which often lead to water quality deterioration and loss of ecosystem integrity. Recent estimates have demonstrated their significance in the earth system, contributing to heterogeneity in land surface properties and feedbacks to regional and global climate through energy, water and biogeochemical transfers (Martynov et al., 2012, Cole et al., 2007). For example, Tranvik et al. (2009) suggested carbon burial in lakes and reservoirs is substantial on the global scale, on the order of 0.6 Pg yr -1 , or four times the oceanic burial rate.
Given the diversity of lakes among continents, region-specific pressures and local management approaches, the Global Lake Ecological Observatory Network (GLEON: gleon.org) was initiated in 2004 as a grass-roots science community with a vision to observe, understand and predict freshwater systems at a global scale (Hanson et al., 2016). In doing so, GLEON has been a leading 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 scales -from an individual lake or reservoir (Hamilton et al., 2015) to regional (Read et al., 2014;Klug et al., 2012), and even global trends (Rigosi et al., 2015;O´Reilly et al., 2015). Ultimately, it is the aim of the network to facilitate primary discovery and synthesis 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 recognised, and models have capitalized on 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 that have contributed to our understanding of lake dynamics, including aspects such as mixing regimes, eutrophication dynamics (Matzinger et al., 2007), harmful algal bloom dynamics (Chung et al., 2014), and fisheries (Makler-Pick et al., 2009).
In recent decades a range of 1, 2, and 3-dimensional hydrodynamic models has emerged for lake simulation across a diverse range of time scales. Depending on the dimensionality, the horizontal resolution of these models may vary from metres to tens of kilometres, and the spatial resolution from sub-metre to several metres. As in all modelling disciplines, identifying Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2017-257 Manuscript under review for journal Geosci. Model Dev. Discussion started: 20 November 2017 c Author(s) 2017. CC BY 4.0 License. application to numerous lakes 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., 2017). GLM has been designed to be an opensource community model suited to modelling studies across a broad spectrum of lakes, reservoirs and wetlands. It balances complexity of dimensional representation, applicability to a wide range of standing waters, and availability to a broad community (e.g., GLEON). Given that individual applications of the model are not able to describe the full array of features and details of the model structure, the aim of this paper is to present a complete description of GLM, including the scientific background (Section 2), model code organization (Section 3), approach to coupling with biogeochemical models (Section 4), and to overview use of the model within the context of GLEON specific requirements for model analysis, integration and education (Section 5-6).
GLM adopts a 1D approach for simulating lake mixing processes by resolving a vertical series of layers that describe 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. Depending on the context of the simulation, either daily or hourly meteorological time-series data for surface forcing is 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 is 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, to which the reader is also referred (e.g., Hocking & Patterson, 1991;McCord & 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 & Patterson (1981). The approach defines each layer as a 'control volume' that can change thickness by contracting and expanding in response to inflows, outflows, mixing with adjacent layers, and surface mass fluxes. Layer thickness limits are enforced to adequately resolve the vertical density gradient with fine resolution occurring in the metalimnion and thicker cells where mixing is occurring, as depicted schematically in Figure 1. Unlike fixed-grid (Eulerian) design of most 1D lake and ocean models, where mixing algorithms are typically based on resolving vertical velocities, it has been reported that numerical diffusion at the thermocline can be restricted by this approach (depending on the user-defined minimum (ℎ ) and maximum (ℎ
layer thickness limits set by the user), making it particularly suited to long-term investigations, and requiring limited site- specific calibration (Patterson et al., 1984;Hamilton & Schladow, 1997;Bruce et al., 2017). Layers each have a unique density computed based on the local salinity and temperature, and when sufficient energy becomes available to overcome density instabilities between adjacent layers, they will merge, thereby accounting for the process of mixing. For deeper systems, a stable vertical density gradient will form in response to periods of high solar radiation creating warm, less-dense conditions near the surface with cooler conditions deeper in the water, separated by a metalimnion region which includes the thermocline. The number of layers, 𝑁 (𝑡), is adjusted throughout the simulation to maintain homogenous properties within a layer. Initially, the layers are assumed to be of equal thickness, and the initial number of layers, 𝑁 (𝑡 = 0). 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. whereby 𝐴 = 𝑓(ℎ ), and 𝑖 is the layer number. This computation requires the user provides a number 𝑁 of depths with corresponding areas, and the volumes are estimated as:
where 1 < 𝑏 ≤ 𝑁 . Using the raw hyposgraphic data, a refined depth-area-volume relationship is internally computed using finer depth increments (e.g., ~ 0.1 m), giving 𝑁 levels that are used for subsequent calculations. The area and volume at the depth of each increment, ℎ is interpolated from the supplied information as:
where 𝑉 and 𝐴 are the volume and area at each of the elevations of the refined depth vector, and 𝑉 refers to the nearest b level below ℎ such that ℎ < ℎ . The interpolation coefficients are computed as:
The density in each layer i is computed based on the temperature, 𝑇, and salinity, 𝑆, at any given time according to the UNESCO (1981) equation of state whereby 𝜌 = 𝜌 𝑇 , 𝑆 . Density calculations can also be customised as required.
Because this approach assumes layer properties are laterally averaged, the model is suitable for investigations where resolving the horizontal variability is not a requirement of the study. This is often the case for ecologists and biogeochemists studying natural lakes (e.g., Gal et al., 2009), managers simulating drinking water reservoirs (e.g., Weber et al., 2017), or 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, it may also 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 simplify, and mass and energy will continue to be conserved.
The model solves the water balance of the lake domain by including several user-configurable water fluxes. The components include surface mass fluxes (evaporation, rainfall and snowfall), inflows (surface inflows, submerged inflows and local runoff from the surrounding exposed lake bed area) and outflows (withdrawals, overflow and seepage). The dynamics of inflows and outflows modify the overall lake water balance on a daily time-step, and may impact upon the layer structure by adding, merging or removing layers (described in Sect. 2.6). In addition, the mass balance the surface layer is computed at each model time step (usually hourly), by modifying the surface layer height according to: where h S is the top height of the surface layer (m), t is the time (days), E is the evaporation mass flux computed from the heat flux 𝜙 (W m -2 ) described below, R is rainfall and S is snowfall (m day -1 ), and 𝑓 is a user-definable scaling factor that may be applied to increase or reduce the rainfall data (default = 1). 𝑄 is an optional term to account for runoff to the lake from the exposed banks, which may be important in reservoirs with a large drawdown range, or wetlands where periodic drying of the lake may occur. The runoff volume generated is averaged across the current lake area (𝐴 ), and the amount is calculated using a simple model based on exceedance of a threshold rainfall intensity, R L (m day -1 ), and runoff coefficient:
where 𝑓 is the runoff coefficient, defined as the fraction of rainfall that is converted to runoff at the lake's edge, and 𝐴 is the maximum possible area of inundation of the lake (as defined by the area provided by the user at 𝑁 area value).
Note 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. 4, h S will also be modified as a result of ice formation/melt, and river inflows, withdrawals, seepage or overflows impacting upon the surface layer, which are described in subsequent sections; these are in addition to the above described terms.
A balance of shortwave and longwave radiation fluxes, and sensible and evaporative heat fluxes determine the net cooling and heating for GLM. The general heat budget equation is described as:
where c p is the specific heat capacity of water (4186 J kg -1 °C-1 ), 𝑇 is the surface temperature of the surface mixed layer and 𝑧 is the depth of the surface mixed layer. The RHS heat flux terms, including several options for customizing the individual surface heat flux components, are expanded upon individually below.
Solar radiation is the key driver of the lake thermodynamics, however, data may not always be available from a nearby pyranometer. Users may choose to either have GLM compute surface irradiance from a theoretical approximation based on the Bird Clear Sky insolation model (BCSM) (Bird, 1984), modified for cloud cover and latitude, or alternatively, hourly or daily solar radiation intensity data may be specified directly. If the BCSM is used, then 𝜙 is calculated from (Bird, 1984;Luo et al., 2010): where the model computes total irradiance, 𝜙 (W m -2 ) from direct beam 𝜙 , and atmospheric scattering 𝜙 components (refer to Appendix A for a detailed outline of the BCSM equations and parameters). In GLM, the clear sky value is reduced according to the cloud cover, C, according to:
𝑓 𝐶 = 0.66182 𝐶 -1.5236 𝐶 + 0.98475 (8)
which is based on a polynomial regression of cloud data from Perth Airport, Australia, compared against nearby sensor data (R 2 = 0.952; see also Luo et al., 2010). 5 where z is the depth of the layer from the surface, 𝑓 is a scaling factor that may be applied and adjusted as part of the calibration process, and K w is the light extinction coefficient (m -1 ). K w may be set by the user as constant or linked to the water quality model (e.g. FABM or AED2, see Sect 4) in which case the extinction coefficient will change as a function of depth and time according to the concentration of dissolved and particulate constituents. Beer's Law is only applied for the photosynthetically active fraction (PAR) component, 𝑓 , which is set as 45% of the incident light. The amount of light heating the surface layer, 𝜙 , is therefore the above photosynthetically average fraction that is attenuated across z SML , plus the remaining 1 -𝑓 fraction which accounts for near infra-red and ultraviolet bandwidths of the incident shortwave radiation with significantly higher attenuation coefficients (Kirk, 1994).
In some applications, the extent to which the benthos has a suitable light climate is a good indicator of benthic productivity, and a proxy for the type of benthic habitat that might emerge. In addition to the light profiles, GLM predicts the benthic area of the lake where light intensity exceeds a user defined value (Figure 4), 𝜙 .
where ℎ = ℎ -𝑧 , and 𝑧 is calculated from Beer's law:
The daily average benthic area above the threshold is reported in the lake.csv summary file as a percentage (𝐴 /𝐴 ).
Longwave radiation can either be specified as net flux, incoming flux or, if there is no radiation data from which longwave radiation can be computed, then it may be calculated by the model internally based on the cloud cover fraction and air temperature. Net longwave radiation is described as:
where
and s is the Stefan-Boltzman constant and e 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 𝛼 is the longwave albedo (0.03), and the emissivity of the atmosphere is computed considering emissivity of cloudfree conditions (𝜀 ), based on air temperature (𝑇 ) and vapour pressure, extended to account for reflection from clouds, such that 𝜀 * = 𝑓 𝑇 , 𝐶 calculated from (Henderson-Sellers, 1986):
where, C is the cloud cover fraction (0-1), e a the air vapour pressure calculated from relative humidity, and options 1-4 are chosen via the cloudmode variable. Note that cloud cover is typically reported in octals (1-8) with each value depicting a fraction of 8. So 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 results have been reported by averaging the octal values for all cloud types to get an average value of cloud cover If longwave radiation data does not exist and cloud data is also not available, but solar irradiance is measured, then it is possible to get GLM to compare the measured and theoretical (BCSM) solar irradiance to approximate the cloud cover fraction. This option utilises the above relation in Eq. 7 to compute 𝜙 , and clouds are approximated by assuming that 𝜙 𝜙 = 𝑓(𝐶). Note that if neither shortwave or longwave radiation are provided, then the model will use the BCSM to compute incoming solar irradiance and cloud cover will be assumed to be 0.
The model accounts for the surface fluxes of sensible heat and latent heat using commonly adopted bulk aerodynamic formulae. For sensible heat:
where c p is the specific heat capacity of air (1005 J kg -1 °C-1 ), C H is the bulk aerodynamic coefficient for sensible hear transfer (~1.3x10 -3 ), T a the air temperature (°C) and T s the temperature of the surface layer (°C). The air density is in kg m -3 and computed from 𝜌 = 0.348 (1 + 𝑟)/(1 + 1.61𝑟) 𝑝/𝑇 , where 𝑝 is air pressure (hPa) and r is the 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), 𝜅 is the ratio of molecular weight of water to molecular weight of air ( = 0.622) and 𝜆 is the latent heat of vaporisation. The vapour pressure can be calculated by the following formulae:
Correction for non-neutral atmospheric stability : For long-time integrations (i.e., 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 water body (e.g. Strub and Powell, 1987). At finer timescales (hours to weeks), the thermal inertia of the water body is too great and so the transfer coefficients must 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) parameterised 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 applied as options in the model, and requires that the measurement of wind speed, air temperature and relative humidity within the internal boundary layer over the lake surface are supplied at an hourly resolution.
Still-air limit : The above formulations only 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. 16-17 are no longer 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 have small fetches that limit the build-up of wind speed. These water bodies are often sheltered from the wind and will develop surface temperatures warmer than the atmosphere for considerable periods.
Therefore, we optionally augment Eqs. 16-17 with calculations under low wind-speed conditions by calculating the evaporative and sensible heat flux values for both the given U x and for an assumed U x = 0. The chosen value is found by taking the maximum value of the two calculations:
where 𝜙 is the zero-wind flux, given below, and applies for both evaporative and sensible heat fluxes, and 𝜙 , is calculated from Eqs. 16-17. The two zero-wind speed heat flux equations are taken from TVA ( 1972), but modified slightly to return energy flux in SI units (W m -2 ):
where 𝐶 = 𝜅 𝑒 𝑝, with the appropriate vapour pressure values, e, for both surface and ambient atmospheric values. Here, v is the molecular heat conductivity of air (0.1 kJ m -1 h -1 K -1 ), n is the kinematic viscosity of the air (0.0548 m 2 h -1 ], r o is the density of the saturated air at the water surface temperature, r s is the density of the surface water, x is a roughness correction coefficient for the lake surface (0.5) and a is the molecular heat diffusivity of air (0.077 m 2 h -1 ). Note that the impact of low wind speeds on the drag coefficient is captured by the modified Charnock relation (Eq. A24), which includes an additional term for the smooth flow transition (see also Figure A1). where A C is the critical area. In GLM, the ratio of the effective area to the total area of the lake 𝐴 𝐴 is then used to scale 𝑈 as a means of capturing the average wind speed over the entire lake surface.
The algorithms for GLM ice and snow dynamics are based on previous ice modelling studies (Patterson and Hamblin, 1988;Gu and Stefan, 1993;Rogers et al., 1995;Vavrus et al., 1996;Launiainen and Cheng, 1998;Magee et al., 2016). To solve the heat transfer equation, the ice model uses a quasi-steady assumption that the time scale for heat conduction through the ice is short relative to the time scale of meteorological forcing (Patterson and Hamblin, 1988;Rogers et al., 1995). The steady-state conduction equations, which allocate shortwave radiation into two components, a visible (A1=70%) and an infra-red (A2=30%) spectral band, are used with a three-component ice model that includes blue ice (or black ice), white ice (or snow ice) and snow (see Eq. 1 and Fig. 5 of Rogers et al., 1995). White ice is generated in response to flooding, when the mass of snow that can be supported by the buoyancy of the ice cover is exceeded (see Eq. 13 of Rogers et al., 1995). By assigning appropriate boundary conditions to the interfaces and solving the quasi-steady state of heat transfer numerically, the model computes the upward conductive heat flux from the ice or snow cover to the atmosphere, 𝜙 . The estimation of 𝜙 involves the application of an empirical equation (Ashton, 1986) to estimate snow conductivity (K s ) from its density (Figure 6).
At the ice (or snow) surface, a heat flux balance is employed to provide the condition for surface melting:
where L is the latent heat of fusion (see physical constants, Table 1), h i is the height of the upper snow or ice layer, t is time, r is the density of the snow or ice, determined from the surface medium properties, T 0 is the temperature at the solid surface, T m is the melt-water temperature (0 o C) and f net (T 0 ) is the net incoming heat flux, at the solid surface:
where f LWin and f LWout are incoming and outgoing longwave radiation, f H and f E are sensible and evaporative heat fluxes between the solid boundary and the atmosphere, and f R is the heat flux due to rainfall. These heat fluxes are calculated as above with modification for determination of vapor pressure over ice or snow (Gill, 1982) and the addition of the rainfall heat flux (Rogers et al., 1995). T 0 is determined using a bilinear iteration until surface heat fluxes are balanced (i.e. f 0 (T 0 ) = -f net (T 0 )) and T 0 is stable (± 0.001 o C). In the presence of ice (or snow) cover, surface temperature T 0 > T m indicates that Accretion or ablation of ice is determined through the heat flux at the ice-water interface, q f , Solving for heat conduction through ice yields:
where 𝜙 is the shortwave radiation penetrating the surface, K refers to the light attenuation coefficient of the ice and snow components designated with subscripts s, w and e for snow, blue ice and snow ice respectively, and h refers to the thickness of snow, white ice (snow ice) and blue ice. Q white is a volumetric heat flux for formation of snow ice, which is given in Eq. 14 of Rogers et al. (1995). Ice and snow light attenuation coefficients in GLM are fixed to the same values as those given by Rogers et al. (1995). Shortwave albedo for the ice or snow surface is a function of surface medium (snow or ice), surface temperature and ice or snow thickness (see Table 2, Vavrus et al., 1996). Values of albedo derived from these functions vary from 0.08 to 0.6 for ice and from 0.08 to 0.7 for snow.
The imbalance between q f and the heat flux from the water to the ice, q w , gives the rate of change of ice thickness at the interface with water:
where r blue is the density of blue ice and q w is given by a finite difference approximation of the conductive heat flux from water to ice:
where K m is molecular conductivity and DT is the temperature difference between the surface water and the bottom of the ice, which occurs across an assigned depth Dz. A value for D z of 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).
Figure 6 summarizes the algorithm to update ice cover, snow cover and water depth. The ice cover equations are applied when water temperature first drops below 0 °C. 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 After the change in ice thickness due to heat exchange is calculated, the effects of snowfall, rainfall, and compaction of snow 5 are calculated through appropriate choice of one of several options, depending on the air temperature and whether ice or snow is the upper solid boundary (Figure 6). 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: r s,max = 300 kg m -3 , r s,min = 50 kg m -3 ) truncated to the appropriate limit. The snow compaction model is based on the exponential decay formula of McKay (1968), with selection of snow compaction parameters based on air temperature (Rogers et al., 10 1995) as well as on rainfall or snowfall. The approach of snow compaction used by Rogers et al. (1995) is to set the residual snow density to its maximum value when there is fresh snowfall. This method is found to produce increases in snow density that are too rapid when there is only light snowfall. As a result, GLM uses a gradual approach where the new snowfall and where 𝜉 is the K-H billow length scale (described below), 𝑢 is the shear velocity at the interface of the mixed layer, and 𝐶 , 𝐶 , and 𝐶 are mixing efficiency constants. For mixing to occur, the energy must be sufficient to lift up water at the bottom of the mixed layer, denoted here as the layer 𝑘 -1, with thickness ∆ℎ , and accelerate it to the mixed layer velocity. This also accounts for energy consumption associated with K-H production and expressed as, E PE :
where z is the thickness of the surface mixed layer. To numerically resolve the above equations we sequentially compute the different components of the above expressions in light of the layer structure. Here GLM follows the algorithm outlined in Imberger and Patterson (1981) whereby cooling is computed so that layers are combined due to convection, then stirring, and then shear and K-H mixing are computed.
To compute mixing due to convective cooling, the value for 𝑤 * is calculated, which is the turbulent velocity scale associated with convection. The model adopts the algorithm used in Imberger and Patterson (1981), whereby the potential energy that is released by mixed layer deepening is computed by looking at the moments of the different layers in the surface mixed layer (from layers K to N LEV ):
where 𝜌 is the mean density of the mixed layer including the combined layer, 𝜌 is the density of the k th layer, Δz is the height difference between two consecutive layers within the loop (Δz = h -h ), h is the mean height of layers to be mixed ( h = 0.5[ h + h ] ), and h is the epilimnion (surface mixed layer) mid height, calculated from: h = 0.5 ℎ + h .
The velocity scale 𝑢 * is associated with wind stress and calculated according to the wind strength:
where 𝐶 is the drag coefficient for momentum. and 𝑔 = is the reduced gravity between the mixed layer and k-1 layer. If the condition is not met the energy is stored for the next time-step.
Once stirring is completed, mixing due to velocity shear is applied. Velocity shear at the interface is approximated from:
where t is characteristic time scale over which the shear has been operating, considered relative to t shear , which is the time beyond which there is no shear production (i.e., 𝑢 = 0 if the time exceeds t shear ). This cut-off time assumes use of only the energy produced by shear at the interface during a period equivalent to half the basin-scale seiche duration, T i , and modified to account for damping:
where 𝑇 is the time-scale of damping (see Spigel, 1978). The wave period is approximated based on the stratification as
where 𝐿 is the length of the basin at the thermocline and c is the internal wave speed. Once the velocity is computed, 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: where the K-H length scale is 𝜉 = 𝐶 𝑢 𝑔 and ∆𝜉 = 2 𝐶 𝑢 ∆𝑢 𝑔 ; in this case the reduced gravity is computed from the difference between the epilimnion and hypolimnion, and 𝐶 is a measure of the billow mixing efficiency.
Once shear mixing is done, the model checks the resultant density interface to see if it is unstable to shear (i.e., K-H billows would be expected to form). This occurs if the gradient is less than the K-H length scale. If K-H mixing is required, layers are further split and subject to mixing using an algorithm similar to above.
Mixing below the SML in lakes, in the deeper stratified regions of the water column, is modelled using a characteristic vertical diffusivity, 𝐾 = 𝐾 + 𝐾 , where 𝐾 is the fixed molecular diffusivity of scalars. Three hypolimnetic mixing options are possible in GLM including: (i) no diffusivity (ii) a constant vertical diffusivity 𝐾 over the water depth below the thermocline or (iii) a derivation by Weinstock (1981), 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 𝛼 is interpreted as a vertical diffusivity (m²/s). For the Weinstock (1981) model, the diffusivity is computed according to:
where 𝛼 is the mixing efficiency of hypolimnetic TKE (~0.8 in Weinstock, 1981) and 𝑘 is the turbulence wavenumber:
and 𝑢 * = 1.612 ×10 𝑈 . The term 𝑁 is the Brunt-Väisälä (buoyancy) frequency defined as:
Estimating the turbulent dissipation rate can be complex and GLM adopts a simple approach as described in Fischer et al.
(1980) where a "net dissipation" is approximated by assuming dissipation is in equilibrium with energy inputs from external drivers:
which is expanded and calculated per unit volume as: The diffusivity is calculated according to Eq. 42, but since the dissipation is assumed to concentrate close to the level of strongest stratification, the "mean" diffusivity is modified to decay exponentially with distance from the thermocline:
where 𝜎 is the variance the N 2 distribution below ℎ and scales with the depth over which mixing decays.
Once the diffusivity is approximated (for either model 1 or 2 in Eq. 43), the diffusion of any scalar, 𝐶, between two layers is numerically accounted for by the following mass transfer expressions:
where 𝐶 is the weighted mean concentration of 𝐶 for the two layers, and ∆𝐶 is the concentration difference between them.
𝑓 is related to the diffusivity according to: 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 to both the upper and lower layers.
Inflows can be specified as local runoff from the surrounding (dry) lake domain (Q R described above, Eq. 5), rivers entering at the surface of the lake that will be buoyant or plunge depending on their momentum and density (Sect 2.6.1), or submerged inflows including groundwater (Sect 2.6.2). Any number of inflows to the lake body can be specified and these are applied daily. Four options for outflows are included in GLM, including withdrawals from a specified depth (Sect 2.6.3), adaptive offtake (Sect 2.6.4), vertical groundwater seepage (Sect 2.6.5).
For river inflows, depending on the density of the river water, the inflow will form a positively or negatively buoyant intrusion that will enter the lake and insert at a depth of neutral buoyancy. As the inflow inserts it will entrains water depending on the rate of mixing created by th inflowing water. In GLM, as the inflow crosses layers it will entrain water from each, until it reaches a level of neutral buoyancy and undergoes insertion. Therefore, when it reaches its point of neutral buoyancy a new layer of thickness dependent on the inflow volume at that time (including additions from entrainment) is created. Following insertion, the inflow layer may then amalgamate with adjacent layers depending on numerical criteria within the model for combining or splitting layers. The rate of entrainment of the intrusion, 𝐸, can be calculated in a number of ways. For simplicity, in GLM the rate has been adapted from the first approximation given in Fischer et al. (1979):
where 𝐶 is the user specified drag coefficient for the inflow. The Richardson's number is adapted from Fischer et al.
(1979) as:
where 𝛼 is the stream half angle and 𝜙 is the slope of the inflow at the point where it meets the water body (Figure 10). As the inflow parcel travels through the layers, the increase in inflow thickness due to entrainment is estimated as:
where ℎ is the inflow thickness, 𝐸 is the entrainment rate and 𝑑𝑥 is the distance travelled by the inflowing water, calculated from the flow rate and inflow thickness. The initial estimate of the intrusion height is computed from Imberger and Patterson (1981) where 𝑄 is the inflow discharge provided as a boundary condition and 𝑔′ is the reduced gravity of the inflow given as:
where 𝜌 is the density of the inflow and 𝜌 the density of the surface layer. The distance travelled by the inflow aliquot, 𝑑𝑥, is estimated as the distance travelled in the vertical and the slope of the inflow river, 𝜙 and given by:
where 𝑑𝑧 is the distance travelled in the vertical. The velocity of the inflow aliquot for that day is then calculated as:
Following conservation of mass, the flow is estimated to increase according to Fischer et al. (1979), as in Antenucci et al.
(2005):
The above entrainment and insertion algorithm is repeated for each inflow. Aside from importing mass into the lake, river inflows also contribute turbulent kinetic energy to the hypolimnion, as discussed in the Sect 2.5.2 (e.g., see Eq. 42), and contribute to the scalar transport in the water column (Figure 11a).
Submerged inflows are inserted at the user-specified depth with zero entrainment. The submerged inflow volume is added as a new layer which may then be mixed with adjacent layers (above or below) depending on the density difference, until neutral buoyancy is attained (Figure 11b). This option can be used across one or more layers to account for groundwater inputs, or for capturing a piped inflow, for example.
Outflows from a specific depth can be accommodated including outlets from a dam wall offtake, or other piped withdrawal, or for removing water that may be lost due to groundwater discharge. For a stratified water column, the water will be removed from the layer corresponding to the specified withdrawal depth, as well as layers above or below depending on the strength of discharge and stability of the water column. Accordingly, the model assumes a line-sink algorithm where the thickness of the withdrawal layer is dependent on the internal Froude (𝐹𝑟) and Grashof (𝐺𝑟) numbers, and the parameter, R (see Fischer et al., 1979;Imberger and Patterson, 1981):
where 𝑊 , 𝐿 and 𝐴 are the width, length and area of the lake at the outlet elevation, and 𝑣 is the vertical diffusion coefficient averaged over the withdrawal layer. The Brunt-Väisälä frequency averaged over the thickness of the withdrawal layer, 𝑁 , is calculated as:
where 𝑑𝑧 is the thickness of the withdrawal layer, 𝜌 is the density of the lake at the height of withdrawal and 𝜌 is the density of the lake at the edge of the withdrawal layer. The thickness of the withdrawal layer is then calculated depending on the value of R (Fischer et al. 1978), such that:
The proportion of fluid withdrawn from each layer either above or below the layer of the outlet elevation is determined using a curve that fits the region of fluid drawn in a given time. To calculate the width and length of the lake at the height of 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, 𝐿 and the lake width, 𝑊 are given by:
Seepage of water from the bottom layer is also configurable 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. Seepage is configured to leave the lake at a constant rate:
where h B is the depth of the bottom-most layer at any time, and G is the seepage rate (m day -1 ). 𝐺 is constrained within the model to ensure no more than 50% of the layer can be reduced in any one time-step. Note that in shallow-lake or wetland simulations, the layer structure may simplify to a single layer, in which case the surface and bottom layer are the same, and
Eqs. 4 and 60 are effectively combined.
Wind-induced resuspension of sediment from the bed of shallow lakes is sporadic and occurs as the waves created at the water surface create oscillatory currents that propagate down to the lake-bed. GLM does not predict resuspension and sediment concentration directly, but computes the bottom shear stress for later use by sediment and water quality modules.
Nonetheless, even without this explicit formulation, the model can identify the areal extent and potential for bed-sediment resuspension by computing the area of the lake over which the bed shear stress exceeds some critical value required for resuspension to occur. To compute the stress at the lake bottom we estimate the surface wave conditions using a simple, fetch-based, steady state wave model (Laenen and LeTourneau, 1996;Ji 2008). The wave geometry (wave period, significant wave height and wave length), is predicted based on the wind speed and fetch over which the waves develop (Figure 13), calculated as:
Using this model, the wave period, T, is calculated from fetch as: 5
where:
and ℎ is the average lake depth. Wave length is then estimated from: The total shear stress at the lake bed is calculated as: 5
where 𝑈 is the mean velocity of the layer, computed during the mixing calculations (Eq. 33). The friction factors use D (a typical particle diameter). For the current stress we compute 𝑓 = 0.24 log 12𝑑 2.5𝐷 and for waves:
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 code structure outlined in Figure 15. The model is written in C, with a Fortran-based interface module to link with Fortran-based water quality modelling libraries in Sect. 4. The model compiles with gcc, and gfortran, and commercial compilers, with support for Windows, OS X and Linux.
To facilitate the use of the model in teaching environments and for users with limited technical support, the model 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 (Figure 16). The outputs from predictions are stored into a structured netcdf file, and this can be visualised in real-time through the simple inbuilt plotting library (libplot), or may be opened for post-processing in MATLAB or R (see Sect. 5.1). Parameters and configuration details are input through the main glm.nml text file (Figure 16) and default parameters and their associated description are outlined in Table 1.
Beyond modelling the vertical temperature distribution, the water, ice and heat balance, as well as 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., 2013) and the Framework for Aquatic Biogeochemical Models (FABM; Bruggeman and Bolding, 2014). Through connection with these libraries, GLM can simulate the seasonal changes in vertical profiles 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 are highlighted here relevant to managing physical-ecological interactions.
Firstly, the model is designed to allow a user-defined number of sediment zones that span the depth of the lake. Using this approach, the current setup allows for depth-dependent sediment properties, both for physical properties such as roughness or sediment heat flux, and also biogeochemical properties such as sediment nutrient fluxes and benthic ecological
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" (https://github.com/GLEON/glmtools) and the MATLAB library is called "GLMm"
(https://github.com/AquaticEcoDynamics/GLMm). Both tools have utilities for model 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 results figures above), comparisons to field observations, and various evaluations of model performance.
In order to compare the performance of the model for varied types of lakes, numerous different metrics of model performance are relevant. These include simple measures like surface or bottom temperature, or ice thickness, however, it is also possible to compare the model's performance in capturing higher-order metrics relevant to lake dynamics, including Schmidt Stability, thermocline depth, ice on/off dates (see also Bruce et al., 2017, 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 high frequency sensor data for calculation of key lake stability metrics (Read et al., 2011), then continuous wavelet transform comparisons are also possible (Kara et al., 2012), allowing assessment of the time-scales over which the model is able to capture the observed variability within the data.
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 attempts at physically based descriptions where there is relatively little variation (Bruce et al., 2017), thereby reducing the number of parameters that remain uncertain, however, for others their variation reflects imperfect formulation of some processes that are not fully considered. 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 (Figure 18). Wrappers and examples for use of GLM within the openDA framework and PEST are also being tested, giving users access to a wide range of uncertainty assessment and data assimilation algorithms. The PEST framework allows for calibration of complex model using highly-parameterised regularisation with pilot-points (Doherty, 2015), and sensitivity matrices derived from the calibration process can also be utilised in linear and non-linear uncertainty analysis. 5
Questions relevant to land use and climate change are driving scientists to develop scenarios for how lake ecosystem services 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 10 consequences from the emergent properties of the simulation, such as changes in water clarity, extent of anoxia, mixing regime, or habitability to fishes (Hipsey et al., 2015). Often, it takes years or even decades for lakes to respond fully to changes in exogenous drivers, requiring simulations to recreate lake behavior over extended periods. While most desktop computers can run a decade-long, low-resolution simulation in less than one minute, 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), 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 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 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, as well as 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 setup.
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.
Similarly, scripts exist for coupling GLM with the Weather Research Forecasting (WRF) model, or similar climate models, for specification of the meteorological input file from weather prediction simulations.
The above coupling approaches require the models to be run in sequence, however, for the simulation of lake-wetlandgroundwater systems, 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 can be simulated using GLM coupled with the 3D groundwater flow model, FEFLOW (https://www.mikepoweredbydhi.com/products/feflow). For this case the GLM code is compiled as a Dynamic Link Library (DLL) 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 a significant error, however, error can be assessed and reduced using smaller time step lengths. FEFLOW models can be simulated for flow-only, or including heat and/or solute transport. Depending on the
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 development (Stewart et al., 2005;Zohar and Dori 2011). For example, modelling forces students to analyze patterns in data, create evidence-based hypotheses for those patterns and make their hypotheses explicit, and develop predictions 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 analyze hundreds of model simulations through distributed computing overlay networks embedded via the GRAPLEr interface, described above. Hence, students participating in the module learn computing and quantitative skills in addition to improving their understanding of how climate change is affecting 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; all skills that are transferrable for many other applications. Thus, following previous studies (Schwarz andWhite 2005, Schwarz et al. 2009), we predict that increased experience with GLM modelling will not only build students' understanding of lake ecosystem concepts but also their ability to interpret natural phenomena and generate new scientific knowledge.
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 for their individual lake simulation requirements, in addition to being suitable for application in a distributed way across tens to thousands of lakes for regional and global scale assessments. Recent examples include an application of the model for assessing how the diversity of >2000 lakes in lake-rich landscape in Wisconsin respond to meteorological conditions and projected warming (Read et al., 2014), and given its computationally efficient nature it is envisioned to be made available as a library for use within in land-surface models (e.g., the Community Land Model, CLM), allowing 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 the simulation of lake stratification, air-water interaction of momentum and heat, and also biogeochemically relevant variables associated with contemporary questions about greenhouse gases 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 their application of the model is suited to this 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 longerterm 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 on lake ecosystem dynamics. However, the model can also be used for shallow lakes, ponds and wetland environments where the water column is relatively well mixed. In order to better define the typical level of model performance across these diverse lake types, a companion paper by Bruce et al. (2017), has undertaken a systematic assessment of the model's error structure against 31 lakes from across GLEON, to which readers can refer to for further guidance. In cases where the assumption of one-dimensionality is not met for a particular lake application, then it is possible for users to apply two or three dimensional models; potentially these can be temporally nested within a longer term GLM simulation. This paper has focused on description of the hydrodynamic model, but we highlight 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 evapoconcentration. 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 its interpretation. The 1D nature of the model makes the run-times modest and therefore the model is suitable for application within more intensive parameter identification and uncertainty assessment procedures. This is particularly relevant as the needs for network participants to expand model configurations to further include biogeochemical and ecological state variables. It is envisioned that continued application of the model to lakes of GLEON 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 cyberinfrastructure, 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 cross-disciplinary collaboration (Weathers et al., 2017). Initial use of GLM within the classroom has found that teaching modules integrating GLM into classes improves students' understanding of how climate change is altering lake ecosystems.