Patterns and drivers of deep chlorophyll maxima structure in 100 lakes: The relative importance of light and thermal stratification

The vertical distribution of chlorophyll in stratified lakes and reservoirs frequently exhibits a maximum peak deep in the water column, referred to as the deep chlorophyll maximum (DCM). DCMs are ecologically important hot spots of primary production and nutrient cycling, and their location can determine vertical habitat gradients for primary consumers. Consequently, the drivers of DCM structure regulate many characteristics of aquatic food webs and biogeochemistry. Previous studies have identified light and thermal stratification as important drivers of summer DCM depth, but their relative importance across a broad range of lakes is not well resolved. We analyzed profiles of chlorophyll fluorescence, temperature, and light during summer stratification from 100 lakes in the Global Lake Ecological Observatory Network (GLEON) and quantified two characteristics of DCM structure: depth and thickness. While DCMs do form in oligotrophic lakes, we found that they can also form in eutrophic to dystrophic lakes. Using a random forest algorithm, we assessed the relative importance of variables associated with light attenuation vs. thermal stratification for predicting DCM structure in lakes that spanned broad gradients of morphometry and transparency. Our analyses revealed that light attenuation was a more important predictor of DCM depth than thermal stratification and that DCMs deepen with increasing lake clarity. DCM thickness was best predicted by lake size with larger lakes having thicker DCMs. Additionally, our analysis demonstrates that the relative importance of light and thermal stratification on DCM structure is not uniform across a diversity of lake types.

lakes and oceans (Moll et al. 1984;Weston et al. 2005;Giling et al. 2017) and influence nutrient cycling (Jamart et al. 1977;Letelier et al. 2004). DCMs create a vertical resource gradient for primary consumers and thus influence zooplankton vertical distributions and diel vertical migration (Williamson et al. 1996a;Winder et al. 2003), as well as predator aggregations (Tiselius et al. 1993). Additionally, DCMs are often associated with elevated densities of mixotrophic (Bird and Kalff 1989) and heterotrophic (Dolan and Marras e 1995) protozoans, and bacterial biomass that can be 10X greater than elsewhere in the water column (Auer and Powell 2004). Consequently, the presence and characteristics of DCMs, such as their depth or size, can have important implications for predator-prey interactions, biogeochemistry, and energy flow in aquatic ecosystems.
Many characteristics contribute to regulating the depth of DCM formation in lakes (Durham and Stocker 2012;Cullen 2015). Light availability, the vertical distribution of nutrients, and the location of thermal gradients in the water column are frequently identified as the main abiotic drivers of the depth of DCM formation (e.g., Fee 1976;Abbott et al. 1984). DCMs form because phytoplankton must balance opposing gradients of light from above with the availability of nutrients, which often increase with depth under stratified conditions (Abbott et al. 1984;Carney et al. 1988;Clegg et al. 2007;Descy et al. 2010). Other biological factors can contribute to DCM formation, including high in situ growth at depth (Coon et al. 1987;Lande et al. 1989), or behavioral aggregations to nutrient-replete depths via physiologicalinfluenced swimming in motile species (Durham and Stocker 2012), or buoyancy-controlled movements in non-motile species (e.g., Oliver 1994;Villareal et al. 1996). Horizontal movement of water masses such as intrusions or those that cause shear may also contribute to DCM formation by laterally transporting deep, nutrient-or phytoplankton-replete waters (Durham and Stocker 2012). Under oligotrophic conditions, DCMs may also form independently of phytoplankton biomass peaks due to increased chlorophyll (Chl) per unit biomass because in these highly transparent lakes sufficient light may penetrate to depths with relatively higher nutrient conditions, a process known as photoacclimation or photoadaptation (Fennel and Boss 2003). These mechanisms are not mutually exclusive and may act in concert to contribute to the formation of a DCM (Mignot et al. 2014).
Thermal stratification determines DCM depth because it can regulate nutrients mixing from deep waters into the euphotic zone (Abbott et al. 1984;Varela et al. 1994;Mellard et al. 2011;White and Matsumoto 2012) and is an important factor in determining the depth of strong vertical gradients in nutrients (i.e., the nutricline). Non-motile phytoplankton may also settle at depths of neutral density along a vertical temperature gradient (Alldredge et al. 2002;Durham and Stocker 2012). Accordingly, DCMs are frequently associated with both the thermocline and areas of low light, often reported as 1-3% of incident surface photosynthetically active radiation (PAR; Fee 1976; Banse 1987Banse , 2004Perez et al. 2002;Morel et al. 2007), or as daily integrated PAR values of 0.1-1.2 mol quanta m 22 d 21 (Letelier et al. 2004;Mignot et al. 2014).
Other characteristics of DCM structure, aside from DCM depth, have received far less attention. Though analytical definitions have varied, DCM thickness is generally defined as the depth range over which a Chl peak occurs (Platt et al. 1988;Beckmann and Hense 2007;Hanson et al. 2007;Jobin and Beisner 2014). DCM thickness can be quantified as the depth range where Chl values are within a certain percentage of the maximum Chl concentration (e.g., Jobin and Beisner 2014), or where a certain percentage of the integrated Chl occurs (e.g., Platt et al. 1988). Although, when specifically considering deep biomass layers, the depth range of net primary production can also be used (Beckmann and Hense 2007).
The few studies that have examined factors controlling DCM thickness have concluded that DCMs are thinner and more sharply peaked with stronger stratification at the depth at which they develop (i.e., where the relative availability of light vs. nutrients is optimized; Varela et al. 1994;Mellard et al. 2011). In areas of strong stratification, non-motile cells are hypothesized to aggregate in thin layers by sinking or rising to the depths of neutral buoyancy (Alldredge et al. 2002;Durham and Stocker 2012). Once these cells reach neutral buoyancy, vertical dispersion is inhibited by the strong density gradient and in situ growth (Lande et al. 1989) and photoacclimation (Fennel and Boss 2003) may further contribute to high Chl concentrations over a small vertical depth range. Theoretical models predict that the exponential decay of light in the water column, combined with a gradient of nutrient concentrations below the thermocline, means that optimal conditions for high phytoplankton growth can exist over an increasingly broad depth range with increasing lake transparency, which may also influence DCM thickness (Beckmann and Hense 2007). However, because our understanding of mechanistic controls of DCM thickness primarily stems from theoretical modeling (e.g., Varela et al. 1994;Beckmann and Hense 2007;Mellard et al. 2011), not empirical studies (but see Mignot et al. 2011), the way in which these mechanisms interact across broad gradients of lake morphometry and transparency remains an open question.
DCM features (i.e., depth and thickness) are likely controlled by simultaneously operating characteristics that covary with one another. Specifically, the strong co-variation of light and temperature with depth (Houser 2006;Read and Rose 2013) complicates our ability to differentiate the relative importance of each in controlling DCM structure. Both light attenuation (Read and Rose 2013) and convective vs. wind-driven mixing (Read et al. 2012) control the distribution of incident thermal energy throughout the water column. As a result, both transparency and surface area may interact to determine the strength and location of water column stability, as well as the relative position of light and nutrient gradients. Moreover, lake surface area and depth interact to influence transport processes such as mixing, sedimentation, resuspension, and diffusion (Håkanson 2005), and thus may also modulate the influence of light and thermal stratification on DCM structure.
A fundamental challenge in predicting DCM characteristics and identifying the drivers of DCM structure in lakes has been a lack of data across broad gradients of lake transparency and morphometry. Most previous work, by necessity, has focused on DCM depth in one lake (Abbott et al. 1984), a small set of lakes through time (e.g., Hamilton et al. 2010), or in a small set of lakes in a constrained geographical region (Fee 1976; but see Longhi and Beisner 2009). These studies have provided insight into the roles of lake trophy or clarity on DCM depth. For example, DCMs deepen with increasing water transparency (Fee 1976;Longhi and Beisner 2009;Hamilton et al. 2010). However, the relatively small number of lakes limits our understanding of generalized patterns in DCM structure across broad gradients of lake size or transparency. Fortunately, recent advances in fluorometric sensors and their widespread use have improved our ability to collect information on the vertical distribution of Chl fluorescence across a large number of lake ecosystems .
We collated a 100-lake dataset of light, temperature, and Chl fluorescence profiles collected during the summer stratified period from the Global Lake Ecological Observatory Network (GLEON; Weathers et al. 2013;Rose et al. 2016) to assess the relative importance of variables associated with light availability vs. those associated with the depth and strength of thermal stratification to predict DCM depth and thickness. Additionally, we assessed how light-and stratification-associated variables scaled across broad gradients of lake morphometry and transparency. To deal with the analytical challenges of correlated explanatory variables and non-linear relationships that often occur in large-scale ecological studies, we used a random forest algorithm (Breiman 2001) to contrast the relative importance of predictors for the two focal DCM characteristics and to examine how their importance varied across lakes.

Study sites
In situ Chl fluorescence (ChlF) and temperature profiles were collated from lakes representing 13 different ecological regions in 8 countries (Abell et al. 2008, Fig. 1). The lakes spanned a broad range of lake size, maximum depth, transparency, dissolved organic carbon concentrations, and trophic state (Table 1; Supporting Information Table S1). We focused on thermally stratified temperate lakes during summer because the water column must be stratified for a DCM to form (Durham and Stocker 2012 and references therein;White and Matsumoto 2012). All data were collected during mid-summer when surface water temperatures were the warmest in the water column and greater than 48C, which are the requirements for thermal stratification (Boehrer and Schultze 2008). While it is important to recognize that peaks in ChlF may not necessarily represent a biomass peak (Fennel and Boss 2003), ChlF is a widely used proxy for phytoplankton biomass. One ChlF profile collected during the mid-summer stratification was used for each lake. Although this is a potential limitation due to temporal variability in DCMs that have been observed in some lakes (Hamilton et al. 2010;Brentrup et al. 2016), we chose to focus here on cross-lake patterns in DCM characteristics, rather than seasonal dynamics, to maximize the number of available lakes that had at least one ChlF vertical profile.

Chlorophyll fluorescence data processing
We classified lakes as having a DCM if a depth could be identified below the top 5% of the water column depth, where ChlF was at least 1.5 times the average fluorescence in the top 5% of the water column. Our DCM filter was followed by visual assessment of each profile to verify DCM presence and to distinguish whether a lake had a single or multiple peaks, resulting in a dataset of 100 lakes that exhibited a DCM. Although we could not compare lakes with and without a DCM, these criteria allowed us to focus specifically on the drivers of variation in DCM structure within a population of lakes that exhibited DCMs. Note that while nonphotochemical quenching (NPQ) can reduce ChlF in high light areas near the surface of lakes (Serra et al. 2009;Huot and Babin 2010), we did not correct for NPQ because we are not aware of appropriate methods, given the diverse data and sensors used in our dataset. While NPQ may have caused us to identify a lake as having a DCM when none was present, this is unlikely because only three lakes fell close to our 1.5x fluorescence-depth threshold (1.5-1.6) for DCM classification. The median ratio of average ChlF in the top 5% of the water column to the ChlF at the DCM was 6 (mean 5 10.9, 1 st -3 rd quartile 5 3.5-12.0, range 5 1.5-76.2). Further, Mignot et al. (2011) examined the effects of NPQ on DCM depth and thickness in marine systems by comparing estimated Chl derived from in situ fluorescence vs. high performance liquid chromatography. The authors found that NPQ caused an underestimation of DCM thickness in 24% of their dataset and that identification of the DCM depth was generally robust to NPQ.
Fluorometer values were not comparable across lakes due to differences in the algorithm and calibration procedures among different sensor models. We therefore normalized all profiles to a maximum ChlF value of 100 within each lake to facilitate comparisons among lakes. While normalization meant that we could not compare the magnitude of the ChlF peak across lakes differences in ChlF values and ranges between fluorometer manufactures would have made that unfeasible even before the normalization. The depth intervals of the ChlF data within each lake were standardized to 0.1 m using a piecewise cubic hermite polynomial function (R Package pracma; Borchers 2015), which prevents local overshooting of data values, and thus preserved the shape of the ChlF profiles (Moler 2004). Prior to interpolation, the median vertical resolution of the ChlF profiles across the  Table S1 in the Supporting Information for more detailed information on each lake. Points in North American and Europe are semi-transparent so more darkly colored points represent overlapping points. We calculated depth and thickness of the ChlF peak for each profile and used these as derived response variables in all analyses (Fig. 2). The depth of the DCM corresponded to the maximum ChlF observation in each lake. The thickness of the DCM was the depth range over which a ChlF peak occurred (Fig. 2). We fit a curve to each ChlF profile to characterize DCM thickness. Our approach built on previous studies that have characterized Chl profiles in a similar manner. Lewis et al. (1983) used a Gaussian distribution to represent vertical profiles of Chl a, and, Platt et al. (1988) and Morel and Berthon (1989) superimposed a Gaussian distribution onto a constant background Chl concentration. Recent studies have improved the curve fitting by superimposing the Gaussian profile onto a linearly (Uitz et al. 2006) or exponentially decreasing background Chl concentration (Mignot et al. 2011).
We fit an exponential power distribution (Subbotin 1923), also known as a p-generalized normal distribution (Nadarajah 2006), with a constant background fluorescence to our ChlF profiles. The exponential power distribution is more flexible than a Gaussian distribution due to an additional shape parameter, which improved our characterization of the ChlF profiles. Our varied dataset prevented quantification of profile amplitude or magnitude however. For example, expressing the peak ChlF value as a quotient of the average or median ChlF in the profile was problematic because some of the profiles in the dataset did not cover the entire depth of the water column. Restricting the calculation of median or average ChlF values to within the euphotic zone (above the 1% PAR depth) did not provide an accurate or biologically relevant baseline ChlF estimate as the DCM extended below the 1% PAR depth in some lakes.
We used an optimization procedure (R Developement Core Team 2015; function optim) to fit Eq. 1 to the normalized and 0.1 m resolution ChlF profile in each lake.
is an additive parameter that represents the background ChlF as implemented by Platt et al. (1988) and Morel and Berthon (1989). The multiplicative parameter, b (unitless), adjusts the amplitude of the DCM above the background ChlF. The remainder of the right-hand term in Eq. 1 is the exponential power distribution (R package pgnorm; Nadarajah 2005; Kalke and Richter 2013;Kalke 2015) where l (m) is the mean of the distribution. DCM thickness was defined as the depth range that encompassed 1r (m) on either side of l. The unitless shape parameter, p, controls the tail region and thus defines the shape of the distribution (Kalke and Richter 2013). When p 5 2, the curve is a Gaussian distribution, p < 2 indicate sharply peaked distributions and when p > 2 peaks are more broadly distributed (Fig. 2b, Nadarajah 2005). C represents that gamma function.
Optimization was based on minimizing the total sum of squares between the estimated ChlF values from the curve fit and the normalized ChlF values for each profile. To provide more ecologically meaningful descriptive metrics of the DCM, constraints were placed on the exponential power cence (ChlF) as defined by Eq. 1 where the background ChlF is represented by a, the mean of the exponential power distribution is represented by l, and the thickness is defined by 2r (i.e., 1r on either side of l). (b) The shape parameter, p, controls the tail weight and thus the shape of the profile (see "Methods" section for details), and profiles with several values of p and (r 5 3) are shown to demonstrate the variety of shapes that could be characterized with our curve-fitting procedure. (c-f) Example ChlF profiles (green) normalized to 100 and overlaid with the exponential power distribution (black). The solid point represents the depth of the ChlF maximum (Z DCM ), which is the depth of the maximum ChlF value and the open circle represents, l. The thickness of the deep chlorophyll maximum (DCM) is shown by the dashed gray lines and represents 1r on either side of l. Profile (c) is from Emerald Lake (Wyoming, U.S.A.), Z DCM 5 14.1 m, l 5 14.9, r 5 3.29, p 5 3.300, profile (d) is from Heart Lake (Wyoming, U.S.A.) with Z DCM 5 12.8 m, l 5 13.0, r 5 3.20, p 5 0.686, (e) is from Pincher Lake (Ontario, Canada), Z DCM 5 5.8 m, l 5 6.3, r 5 1.86, p 5 2.029 and (f) is from Lake Truite (Quebec, Canada), Z DCM 5 4.3 m, l 5 4.4, r 5 0.99, p 5 0.661. The y-axes in (c-f) are consistent to improve comparisons but note that the depths of the lakes differ. distribution parameters including: 0 < l < maximum lake depth, r < one half of the maximum profile depth and a > -0.1. In cases where l-r was less than 0 (i.e., the fitted DCM extended above the lake surface), we defined the top of the DCM as the lake surface. This scenario occurred only six times in our dataset and decreased the thickness of the DCM an average of 3.3 m (median 5 1.9 m). When l1r extended below the maximum depth of the lake, the DCM bottom and thus the DCM thickness was not defined (five lakes were excluded because of this criterion).

Limnological data
Based on mechanistic understanding of the factors that control DCM structure (e.g., Klausmeier and Litchman 2001;Mellard et al. 2011;Durham and Stocker 2012;Cullen 2015), we selected a suite of predictors of DCM depth and thickness associated with three categories: light attenuation, thermal stratification, and lake morphometry. Light attenuation variables included the 1% PAR depth and epilimnetic dissolved organic carbon (DOC) concentration because DOC is a major regulator of light attenuation in aquatic ecosystems (Williamson et al. 1996b;Bukaveckas and Robbins-Forbes 2000). While DOC and light attenuation are closely related to one another, we included DOC as a separate predictor in our model because recent ecological studies have shown that the effects of DOC on algal biomass can be positive at low DOC concentrations but negative at high DOC concentrations (Seekell et al. 2015). Thus the role of DOC in predicting DCM characteristics likely varies across lakes. Given the importance of stratification for DCM formation (Cullen 1982(Cullen , 2015, we included several metrics related to the location and strength of thermal density gradients in the water column including the thermocline depth, thickness of the metalimnion, and the stratification strength (i.e., buoyancy frequency) at the thermocline and within the DCM (as defined by the depths of l6r). These metrics may predict the relative positions of light and nutrient gradients in the water column and thus the DCM structure. The strength of stratification at such depths may also be an important predictor of DCM thickness by influencing nutrient availability or the ability of cells to aggregate in the water column (Klausmeier and Litchman 2001;Durham and Stocker 2012;Cullen 2015). Lake surface area and maximum depth were also included in our analysis because they can influence transport processes in the water column (e.g., mixing, or circulation-aided diffusion; Håkanson 2005), and modulate the influence of the light and thermal metrics on DCM structure. Light, temperature, and DOC data for each lake were collected on the same day as the ChlF profiles except for Lake Champlain, where the Secchi measurement was taken 4 d later.
We determined the light attenuation coefficient of PAR, K d (m 21 ), for each lake in the dataset using two methods. For 49 of the lakes, radiometer data were available and we estimated K d , as the slope of the linear regression of the natural logarithm of downwelling planar irradiance (E d ) vs. depth. The regression was performed over the entire portion of the light profile that showed a visibly log-linear relationship between the ln(E d (z)) vs. depth (Kirk 1994). In some lakes, the depth range of the regression was altered to avoid boat shadows or surface waves that cause the light measurements in the near surface to be particularly noisy. In the remaining 51 lakes, where radiometer data were unavailable, K d was estimated by dividing 1.7 by the Secchi depth (Poole and Atkins 1929;Idso and Gilbert 1974). The depth at which 1% of the surface PAR remained was estimated as: While light levels at the DCM are often reported as a percent of surface irradiance, daily-integrated light level (mols photons m 22 d 21 ) at a given depth is a better quantitative descriptor of the amount of light at the DCM (Letelier et al. 2004;Mignot et al. 2014;Cullen 2015). However, in lakes where light attenuation (K d ) was estimated using a Secchi disk, irradiance measurements were not available. Because almost all study lakes were located in temperate latitudes where surface irradiance during summer is similar, our estimated 1% PAR depths are likely a reasonable proxy for comparing the amount of light available in the water column among lakes. We estimated multiple metrics of thermal structure for each lake in the dataset, including the depth of the thermocline, the thickness of the metalimnion, buoyancy frequency at the thermocline, and the average buoyancy frequency in the DCM peak or within the metalimnion. Raw temperature profile data were converted to 0.1 m resolution either by averaging or linear interpolation, depending on the original data resolution. The only exception was for Lake Stechlin, Germany, where we interpolated the temperature data to 0.5 m due to a low spatial sampling resolution (approximately 1 m). From these temperature profiles, we estimated the depth of the seasonal thermocline following procedures outlined in Read et al. (2011) using the rLakeAnalyzer package (Winslow et al. 2016). To examine the effects of stratification strength, we also used the temperature profiles to estimate profiles of buoyancy frequency (N) in each lake, which was defined as: where g is the gravitational constant, p 0 is the density at each depth, and @p=@z is the density gradient (Kalff 2002). Density gradients were determined by applying finite center differencing to the density profiles, which were calculated using the 0.1 m resolution temperature data and thermodynamic equations specific to freshwater conditions (Chen and Millero 1986). The top and bottom depths of the metalimnion were defined as the depths above and below the thermocline where N reached at least 35% of the value of N at the thermocline; these depths were found by moving away from the thermocline depth along the profile of N (Pernica et al. 2014) and were verified with visual inspection. Finally, from the N profiles, we derived three additional metrics of thermal stratification: (1) the value of N at the thermocline; (2) the mean buoyancy frequency within the DCM (determined by averaging all N values from the density profiles that fell between the top and bottom depths of the DCM peak, i.e., l6r, see above); and (3) the mean buoyancy frequency within the metalimnion (determined by averaging all N values from the density profiles that fell within the metalimnion). Last, we collated several other limnological variables for each lake that may influence light attenuation or thermal structure including, lake surface area and maximum depth (both transformed using a logarithm base 10). Depth profiles and surface data of soluble or total nutrients were not available for approximately one third of the lakes in our dataset. Thus, we did not include nutrients as predictors in our random forest (RF) analyses because the lakes with nutrient data did not represent a random subsample of the entire dataset and were biased toward the low elevation and low transparency lakes. Results from a separate RF analysis of just the lakes with nutrient data are included in the supplemental information (Supporting Information Fig. S4).

Random forest analysis
We used a random forest (RF) algorithm to assess the relative importance of predictors for each of the two DCM characteristics. RF analyses are robust to overfitting and insensitive to outliers (Breiman 2001) and are able to handle datasets with relatively large numbers of candidate predictors. Importantly, RF analyses are able to detect non-linear relationships without the need to specify them in advance (Jones and Linder 2015), accurately identify important predictors of the response when there is high collinearity among predictors (Archer and Kimes 2008), and provide accurate predictions of the response (Prasad et al. 2006;Cutler et al. 2007).
A RF algorithm is a machine learning technique based on classification, or in our case, regression tree analyses (Breiman 2001;Cutler et al. 2007). The RF algorithm builds many regression trees using a unique bootstrapped sample of the data to construct each individual tree (Liaw and Wiener 2002). Results of each tree are combined to produce a model-averaged fit (Liaw and Wiener 2002). The method differs from traditional regression trees in that only a small number of randomly selected predictor variables are available for selection at each node, which decreases the correlation between trees and therefore the error rate of the entire forest (Archer and Kimes 2008).
A measure of importance was calculated for each predictor by cross-validating each tree with data that were not used when the tree was constructed, referred to as the outof-bag data (Breiman 2001). For a given predictor variable within a single tree, the out-of-bag values were randomly permuted and the mean square error (MSE) of the tree was compared between the original out-of-bag and permuted out-of-bag datasets. A large percent increase in MSE from the original out-of-bag data to the permuted out-of-bag data indicated that the predictor variable had high explanatory power for the response. Increase in percent MSE cannot be compared across RF with different response variables, so we calculated a relative increase in percent MSE (Rel. IMSE) for each RF analysis by dividing the percent MSE for each variable by the highest observed percent MSE for a given response variable, following Read et al. (2015).
We constructed a RF of 1500 trees for each of the two response variables using 1% PAR depth (m), DOC concentration (mg L 21 ), thermocline depth (m), metalimnion thickness (m), buoyancy frequency at the thermocline (s 21 ), lake surface area (log 10 (km 2 )), and maximum depth (log 10 (m)) as predictors included in each analysis. For the RF analyses with DCM thickness as the response variables, only profiles with a fitted R 2 value 0.8 were used and buoyancy frequency within the DCM (s 21 ) was included as an additional predictor variable. We also excluded six additional lakes because the bottom depth of the DCM, as defined by the curve fitting, extended below the bottom of those lakes. Because DCM thickness was derived from the curve fitting procedures, our criteria ensured that only profiles where the thickness was reasonably well characterized were used in the analysis (n 5 75 lakes total). All 100 lakes were included when DCM depth was the response variable because it did not use outputs from the curve fitting. For lakes where predictor variables such as maximum depth or DOC concentration were not available, the lake was excluded from that individual tree, but not the entire analysis, following recommendations from Liaw and Wiener (2002). Our data exclusion method affected 2.4% of the complete data matrix (700 observations) for the RF when DCM depth was the response variable and 2.2% for the RF when DCM thickness was the response (587 observations).
To visualize the direction and nature of the relationship between the most important predictors and each response, and thus assess how the influence of each predictor scales across lakes, we constructed partial dependency plots (PDP) for the most important predictors from each RF model. Partial dependence is computed by predicting the response (e.g., DCM depth or thickness) from the RF over a range of values for the variable of interest, while holding all other variables in the dataset constant (Friedman 2001; also see Auret and Aldrich 2012 for an excellent description). In essence, a PDP represents the relationship between a single variable and the response, after accounting for the average effects or interactions of the other predictors in the model (Carlisle et al. 2009).
All analyses were completed in R version 3.2.2 (R Developement Core Team 2015) and the RF analyses were done with an implementation written by Liaw and Wiener (2002; R package randomForest).

Cross-lake patterns in DCM characteristics
We found a strong positive relationship between DCM depth and thickness, and the 1% PAR depth (Fig. 3a,b), thermocline depth (Fig. 3c,d), lake surface area (Fig. 3e,f) and maximum lake depth (Supporting Information Fig. S1). The depth of the DCM ranged from 1.0 m to 53.9 m (mean 5 9.2 6 0.8 (1 SE), median 5 6.5; Table 2). In the 75 lakes where the ChlF distributions were well described by the curve-fitting procedure, DCM thickness ranged from 0.7 m to 29.3 m (mean 5 6.8 6 0.8 (1 SE), median 5 4.3; Table 2). Eleven of the 100 lakes exhibited distinct double ChlF peaks in the profile, although no profile showed more than two distinct peaks. Measured limnological variables did not differ between lakes with a single vs. double peak. DCM depth and thickness were also positively related to each other (Supporting Information Fig. S2). The depth of the DCM was located at, or shallower than, the 1% PAR depth in 87 of the 100 lakes (Fig. 4a). The median percent of incoming PAR that was available at the depth of the DCM across the entire dataset was 4.8%. With one exception, the percent of PAR at the DCM was less than 23% in all lakes (interquartile range 5 1.6-7.05%, min.-max. 5 0.0001-52.9%). In the majority of the lakes, the depth of the DCM was at or deeper than the top of the metalimnion (85 lakes). Of that majority, 33 lakes had a DCM depth that was located entirely below the metalimnetic layer (Fig. 4b).
The 1% PAR depth was located below the top of the metalimnion in 99% of the lakes. In the one exception, Lake Atitlan (Guatemala), the 1% PAR depth was located 2 m above the estimated top of the metalimnion, although the DCM in Lake Atitlan was located below the 1% PAR depth. The 1% PAR depth was located below the thermocline in 96 of the 100 lakes and below the bottom of the metalimnion in 71 of the lakes (Fig. 4c).
For a small and geographically biased subset of 37 lakes (31 located in Quebec and Ontario, Canada, one in Virginia, U.S.A., five in western Europe), a spectrofluorometer (Fluoroprobe: bbe Moldaenke GmbH, Schwentinental, Germany) was used to collect ChlF. The spectrofluorometer parses the ChlF of four spectral groups of phytoplankton based upon differences in accessory pigments and emission spectra (Beutler et al. 2002). These spectral groups represent the following broad taxonomic groupings: Browns (diatoms, dinoflagelletes, and chrysophytes), Greens (chlorophytes), Cyanobacteria (phycocyanin containing cyanobacteria), and Cryptophytes (cryptophytes and cyanobacteria containing phycoerythrin). Browns were the dominant spectral group at the depth of the DCM in a majority of the lakes (26 of 37 lakes). Cryptophytes largely dominated the spectral composition of the DCM in the remaining lakes (7 of 37), although in one lake the DCM was evenly comprised of Browns and Cryptophytes. Cyanobacteria dominated the DCM in only three lakes and Greens were never the dominant spectral group in the DCM (Supporting Information Fig. S3).

Random forest analyses
Across all lakes, variables associated with light attenuation (i.e., 1% PAR depth and DOC concentration) were more important than thermal stratification metrics to explain DCM depth (Fig. 5a). Together, all predictors in the RF model explained 63% of the variation in DCM depth. The 1% PAR depth was by far the most important variable explaining DCM depth, followed by DOC concentration, and the depth of the thermocline (Fig. 5a). The partial dependency plots (PDP) showed that the predicted DCM depth increased rapidly with increasing 1% PAR depth until the 1% PAR depth reached approximately 30 m, after which the response of the DCM depth to further increases in the 1% PAR depth was weaker (Fig. 6a). There was a weak negative relationship between DOC concentration and the DCM depth at very low concentrations of DOC, and a stronger negative relationship at intermediate DOC levels. At DOC concentrations above 7-8 mg L 21 , further increases in DOC did not predict further decreases in DCM depth (Fig. 6b). DCM depth increased with increasing depth of the thermocline (Fig. 6c).
Lake surface area and maximum lake depth were the most important variables explaining DCM thickness, followed closely by the depth of the thermocline (Fig. 5b). All predictors in the RF model explained 70% of the variation in DCM thickness. Predicted DCM thickness changed very little with log-transformed lake surface area for lakes smaller than approximately 1 km 2 , but showed a substantial increase in thickness in lakes larger than approximately 10 km 2 (Fig.  7a). Similarly, predicted DCM thickness showed only moderate increases when maximum lake depth was shallower than 30 m, but a steeper positive relationship when maximum lake depth was greater than 30 m (Fig. 7b). The DCM also thickened with increasing thermocline depth, and this relationship was strongest when the thermocline depth ranged from approximately 5-15 m (Fig. 7c,d).

Discussion
Across a diverse and globally distributed set of lakes, we identified strong but markedly different patterns in the relative importance of light vs. thermal stratification in explaining DCM depth and thickness. Light and thermal stratification are important predictors of DCMs (e.g., Fee 1976;Cullen 1982Cullen , 2015. However, their strong covariation had previously precluded a sufficient evaluation of their relative importance and influence over the different elements of DCM structure within and across lake types. We found that light attenuation and associated variables (e.g., DOC concentration) were more important predictors of DCM depth than were predictors associated with thermal stratification. In contrast, DCM thickness was predicted by lake size, maximum depth, and other limnological variables related to thermal stratification. Importantly, our analysis also suggests that the relative importance of light and thermal stratification on DCM structure is not uniform across broad gradients of lake morphometry and transparency, but overall these abiotic characteristics explain a majority of the variation in DCM depth or thickness among lakes.

DCM depth
The positive relationship that we observed between DCM and 1% PAR depths across a diverse set of lake types (Fig. 3a) is consistent with results from previous studies, which encompassed a narrower range of lake types (e.g., Fee 1976). The positive relationship has largely been attributed to phytoplankton balancing the need for light from above and nutrients from hypolimnetic waters below to maximize growth (Abbott et al. 1984;Durham and Stocker 2012), although photoacclimation also plays a role (Fennel and Boss 2003). Our results, from a wider range of lake types, also show that the relationship between 1% PAR and DCM depth is not uniform across a gradient of lake transparency (Fig. 6a). Predicted DCM depths from the RF analysis were more sensitive to increases in 1% PAR depth between depths of 10-30 m in the water column (K d range: 0.46-0.19), but were less sensitive when 1% PAR was deeper than 30 m (K d < 0.19).
Differences in the relative position of the 1% PAR depth and the nutricline may be responsible for the variability in the relationship between DCM and 1% PAR depths across lakes observed in our dataset (Abbott et al. 1984;Barbiero and McNair 1996). The relative positions of the 1% PAR and nutricline depths will alter the depth at which light availability from above and nutrient availability from below are balanced and thus interact to control the DCM depth. For example, Barbiero and McNair (1996) observed that the DCM in a small oligotrophic lake deepened over time as the depth of the nutricline increased relative to the 1% PAR depth. Similarly, in a modeling study, Mellard et al. (2011) demonstrated that when the depth of the nutricline was held constant, DCM depths were sensitive to small decreases in light attenuation and thus to the relative position of the 1% PAR and nutricline depths. Nutrient concentrations generally increase with increasing depth in a stratified water column (Holm-Hansen et al. 1976;Cullen 1982;Moll et al. 1984) and likely play an important role in regulating DCM depths. While depth profiles of nutrient data were not available for all of the lakes in our dataset, these previous studies strongly suggest the need to examine the role of vertical distribution of nutrient availability in future empirical studies. However, we note that our RF model explained a substantial amount of variation in DCM depth among lakes in the absence of nutrient predictors, indicating that physical characteristics alone can provide the ability to explain a majority (63%) of variation in DCM depth across lakes with diverse characteristics.
DOC is a major regulator of water clarity and light attenuation in aquatic ecosystems due to its light-absorbing properties (Williamson et al. 1996b;Bukaveckas and Robbins-Forbes 2000). Therefore, we were not surprised that DOC emerged as an important explanatory variable of DCM depth in our analysis. As with 1% PAR depth, the relationship between predicted DCM depth and DOC concentration was not uniform across lakes (Fig. 6b), suggesting that the DOC concentration within a lake will alter the role that DOC plays in influencing DCM depth. Because the relationship between DOC concentration and the 1% PAR depth is characterized by a negative power function, a unit increase in DOC reduces the 1% PAR depth much more in low DOC lakes compared to high DOC lakes (Morris et al. 1995). Supporting this "diminishing returns" concept, our results show that the predicted DCM depth became shallower with increasing DOC until approximately 7-8 mg L 21 , after which continued increases in DOC had little influence on the predicted DCM depth. Williamson et al. (1996b) identified a similar threshold in DOC concentration ( 8-9 mg L 21 ), above which the influence of DOC in regulating the 1% PAR depth decreased substantially.
DCMs are often considered characteristic of highly transparent, oligotrophic lakes (Fee 1976;Moll and Stoermer  1982). Our results, while supporting previous findings that DCMs may be more common in oligotrophic lakes, show that they are not exclusive to these systems (Williamson et al. 1996b). We found that DCMs can also be found in meso-to eutrophic and dystrophic lakes, as was also noted by two recent studies (Simmonds et al. 2015;Brentrup et al. 2016; Supporting Information Table S1). Previous studies also indicated that light must penetrate below the surface mixed layer for a DCM to form (e.g., Hamilton et al. 2010;Simmonds et al. 2015;Brentrup et al. 2016). Consistent with this, in all 100 lakes in our dataset, including the most eutrophic lakes (Trophic State Index, TSI > 50: Supporting Information Table S1), we observed that the 1% PAR depths extended below the thermocline.

DCM thickness
Our RF analysis found that larger and deeper lakes have thicker DCMs. Modeling studies predict that weak thermal stratification around the depth of the DCM will result in thicker Chl peaks (Mellard et al. 2011). The importance of lake size in predicting DCM thickness in our study is likely due to differences in stratification strength between large and small lakes. Lake surface area controls several aspects of thermal stratification and the location of density gradients in the water column (MacIntyre and Melack 2009). The controls of lake stratification reflect the balance of energy between solar radiation inputs that enhance stratification and wind and convective mixing, which serve to weaken stratification (MacIntyre and Melack 2009;Read et al. 2012). Larger lakes are less sheltered from the wind (Markfort et al. 2010) and have a longer fetch (Mazumder and Taylor 1994). Thus, larger lakes are subject to higher levels of wind-mixing and weaker stratification than smaller lakes (Rueda-Valdivia and Schladow 2009), creating conditions for a thicker DCM. Interestingly, our analysis also showed that the effect of lake surface area on predicted DCM thickness is stronger in lakes > 1 km 2 , which is a similar size cutoff identified in previous studies that examined the importance of lake size on physical processes (1-10 km 2 ; e.g., Hondzo and Stefan 1993;Read et al. 2012).
Weaker stratification with increasing lake size may explain why lake size emerged as an important predictor of DCM thickness in our model. However, the mechanism remains unclear because buoyancy frequency within the DCM, a metric of thermal stratification strength, was not a particularly important predictor of DCM thickness in our RF analysis. Previous studies have found that the strength of stratification around the depth of the DCM was negatively related to DCM thickness (Varela et al. 1992;Klausmeier and Litchman 2001;Jobin and Beisner 2014). While our dataset showed exclusively narrow DCMs when buoyancy frequency within the DCM was high (i.e., when stratification was strong), our data showed highly variable DCM thickness in weakly stratified conditions (i.e., at low buoyancy frequency). The high variability in DCM thickness when buoyancy frequency was low is likely the reason that buoyancy frequency within the DCM does not appear as an important predictor of DCM thickness in our RF analysis. This high variability in DCM thickness under weakly stratified conditions Fig. 4. Frequency polygons of the ratios of (a) DCM depth to 1% PAR depth, (b) DCM depth, and (c) 1% PAR to metalimnion top, bottom, and thermocline depth for each lake in the dataset, respectively. Note that the x-axis scale differs among panels but that the bin size is 0.3 for all panels. In (b, c), a single lake is not shown that had DCM depth: thermocline and DCM depth: metalimnion bottom depth ratio of 44 and a PAR: thermocline and PAR: metalimnion top value of 40.
suggests that other factors unaccounted for in our study were likely affecting DCM thickness (e.g., behavioral aggregation, (Alldredge et al. 2002), zooplankton grazing, (Pilati and Wurtsbaugh 2003) or phytoplankton competition (Carney et al. 1988)) while thin DCMs observed under strongly stratified conditions suggests that stratification strength may become a more important controller of DCM thickness with increasing buoyancy frequency at the depth of the DCM (Fig. 8).
Maximum lake depth was also an important explanatory variable of DCM thickness. It is likely that lake depth itself is not a direct causal driver of DCM thickness. Rather, lake , only 75 lakes where the ChlF curve fitting procedure yielded R 2 0. 8 and where the DCM bottom depths were defined above the bottom of the lake were included. The DCM thickness was defined by 2r from the DCM depth (see "Methods" section for full details). Buoyancy frequency in the DCM represents the average buoyancy frequency within that depth range. Metalimnion is abbreviated as meta., buoyancy frequency as buoy. freq, and thermocline as thermo. The predictors in the random forest model (a) explained 63% of the variation in DCM depth and 64% of the variation in DCM thickness (b). depth and basin shape can influence and are correlated with many other physical, chemical, and biological characteristics of lakes (Håkanson 2005;Johansson et al. 2007). For example, maximum lake depth is negatively correlated with lake productivity (Carpenter 1983), and in a study of over 1000 lakes in the United States, maximum lake depth was the most important explanatory variable of epilimnetic levels of total phosphorus, total nitrogen, turbidity, and DOC (Read et al. 2015). The exact mechanism behind maximum lake depth as an important predictor of DCM thickness in our study, and its role in other previous empirical studies, remains unclear. The influence of lake depth on transport processes such as mixing, sedimentation, resuspension, or diffusion (Håkanson 2005), which influence vertical gradients in nutrients and light, and their relative vertical positions, may be important. In a coupled bio-physical process model, Håkanson (2005) showed that when lake surface area is held constant, Secchi disk depths were deeper in a deep vs. a shallow lake. This occurred because the proportion of bottom sediment area that is susceptible to turbulent resuspension is lower in deeper lakes compared to shallower lakes, resulting in a greater flux of nutrients from the sediments to the water column in a shallower lake. The increased nutrient concentration due to this resuspension fueled an increase in epilimnetic primary production, and thus a decrease in Secchi disk depths in the shallower lake (Håkanson 2005). Differences in maximum lake depth may also contribute to variation in the vertical distribution of thermal energy in the water columns of lakes, with important implications for DCM structure.

Phytoplankton community composition of the DCM
The composition of the DCM from the spectrofluorometer data are consistent with previous lakes studies that have shown that DCMs are typically dominated by Browns (diatoms, dinoflagelletes, and chrysophytes) (Fee 1976;Abbott et al. 1984;Barbiero and Tuchman 2004;Saros et al. 2005;Hamilton et al. 2010;Simmonds et al. 2015), cryptophytes (Barbiero and Tuchman 2004;Camacho 2006), and less frequently cyanobacteria (Padis ak et al. 2003;Camacho 2006). Browns, cryptophytes and cyanobacteria are hypothesized to dominate DCMs because they can maintain their vertical position in the water column either through active in the random forest model where DCM depth was the response variable; (a-c) are in descending order of importance. The partial dependency is the dependency of the predicted mean DCM depth, (the mean response) on each explanatory variable after averaging out the effects of all other predictors in the model. Tick marks along the x-axis represent the 10 th -90 th percentiles of each predictor, in 10 percentile increments, from the empirical dataset. Inset in (b) shows the full range of DOC concentrations in the dataset, while the main plot shows a subset of these data where the majority of empirical values exist for improved visualization of the relationship. Horizontal relationships indicate values of the explanatory variables where changes in the variable have little influence on the response after accounting for all other variables in the analysis.
swimming or buoyancy regulation, show decreased settling rates under nutrient-replete and low light conditions (i.e., diatoms, Richardson and Cullen 1995), or because they show high growth rates under low light conditions due to the presence of phycoerythrin (Gervais et al. 1997;Camacho et al. 2001). Our results, however, should be interpreted with caution as the subset of lakes where spectrofluorometer data were available represented a geographically constrained subset of our entire dataset (31 of the 37 lakes were located in southeastern Canada). Future studies that examine species composition of the DCM across a broad range of lakes may provide more insights into the generality of the patterns observed in our data and the relative importance of different mechanism of DCM formation (e.g., behavioral aggregations vs. in situ growth).

Additional considerations
While care was taken to collect high-quality, representative data from across our study systems, some fundamental data collection challenges warrant consideration. Approximately half of the light attenuation data for this study was derived from Secchi depth observations. Differences in attenuation controlled by particulate scattering vs. DOC absorption may lead to variation in 1% PAR depth not explained by Secchi depth (Koenings and Edmundson 1991). We were not able to comprehensively compare differences in the estimates of 1% attenuation depths between radiometer and Secchi readings because only six lakes in our dataset had both radiometer and Secchi readings. Within those six lakes, we found no consistent pattern in over-or underestimation of 1% PAR depths based on values derived from Secchi vs. radiometer-based measurements. The difference in 1% PAR depth estimates between the Secchi vs. radiometer-based methods within each lake varied by an average of 20% (6 10% SD). Despite the added variability in 1% PAR depth estimates that was introduced in our data, 1% PAR depth was by far the most important explanatory variable for DCM depth, suggesting that the methods of light measurements were sufficient to observe this important pattern. Broad-scale use of optical profiling techniques over Secchi disks will improve our understanding of the relationship between DCM depth and the underwater light environment by providing estimates of irradiance or daily radiant exposure throughout the water column, spectral composition at the depth of the DCM, and more accurate estimates of light attenuation.
To our knowledge we are the first to incorporate an exponential power distribution into a model to characterize ChlF profiles (Eq. 1). We found that the additional parameter improved (based on higher R 2 ) the characterization of the ChlF profiles compared with a curve fit using a Gaussian distribution. Further, visual inspection did not indicate the model overfit profiles, rather, the model was more versatile and better fit the range of DCM shapes observed (Fig. 2c-f). Using the exponential power distribution allowed us to expand the breadth and diversity of lakes that we were able to include in our study. For the subset of lakes where R 2 >5 0.8 from the ChlF curve fitting procedure, the median value of the shape parameter was 1.7 (mean 5 8.0 (3.3 standard error), 1 st -3 rd quartile 5 1.0-2.6; note that the mean is driven by a small number of large values) indicating that the DCMs in our dataset were similar in shape to a Gaussian distribution but tended to have heavier tails. Characterizations of DCMs in future studies may be improved by using the exponential power distribution with a shape parameter that is less than two.

Conclusions
Our results show that the depth of light attenuation, lake size, and maximum lake depth explain a substantial amount of variation in the DCM depth and thickness. Light attenuation was a more important predictor of DCM depth than was thermal stratification, while DCM thickness was best predicted by lake size and maximum depth. Long-term increasing or decreasing trends in water clarity (Lottig et al. 2014;Williamson et al. 2015) will alter the depths of light attenuation (Williamson et al. 1996b) and thermal stratification (Read and Rose 2013) with important consequences for DCM structure. In an important extension of previous studies, our analysis also suggests that changing environmental conditions (e.g., eutrophication, oligotrophication, brownification, changing thermal structure) may affect the depth and thickness of DCMs differently across broad gradients of lake size, maximum depth, and transparency. For example, small decreases in the 1% PAR depth may cause a more dramatic reduction in the DCM depth in high transparency lakes but have a much smaller influence in low transparency lakes. The diversity of lakes in our dataset demonstrates that DCMs occur across many different types of lakes, not just clear oligotrophic lakes. Future ecosystem studies may need to consider physical, chemical, and biological processes that occur in deeper parts of the water column, as well as in the surface mixed layer of many lake types. Our study demonstrates that recent advances in high-resolution profiling, in combination with collaborative, data-sharing networks such as GLEON, provide excellent opportunities to improve the understanding of important ecological processes by allowing for studies that encompass a large number of diverse lakes.

Author contribution statement
All authors contributed to the conception of the work as part of the Global Lake Ecological Observatory Network (GLEON) meeting in Orford, Qu ebec, Canada (2014) and/or subsequent conference calls with discussions led by BEB, CCC, and THL. After the lead author, authors are listed in two tiers according to contribution. THL, BEB, CCC, PP, KCR, and YH contributed substantially to discussions regarding the analytical approaches and conceptual framing of this project. The authors after these initial six are listed alphabetically and contributed to interpretation of the results and the discussion. THL led the writing of the manuscript, and conducted all statistical analyses. PP estimated all metrics associated with buoyancy frequency. BEB, CCC, ID, HPG, THL, KCR, JAR, JDS, DS, and PV contributed data. All authors critically revised the manuscript for intellectual content and approved of the final version of the manuscript for publication. Fig. 8. Average buoyancy frequency (s 21 ) in the depth range of the DCM vs. the thickness of the DCM. Note the high variability in DCM thickness when buoyancy frequency is low, but that the thickness is less variable and thinner at high buoyancy frequency values.