Offline Implementation and Evaluation of the Canadian Small Lake Model with the Canadian Land Surface Scheme over Western Canada

<jats:title>Abstract</jats:title>                <jats:p>The Canadian Small Lake Model (CSLM), version 2, was run with the Canadian Land Surface Scheme (CLASS), version 3.6.1, in an offline regional test over western Canada. Forcing data were derived from ERA-Interim and downscaled using the fifth-generation Canadian Regional Climate Model (CRCM5). The forcing precipitation field was adjusted using monthly data from the Canadian Gridded Temperature and Precipitation Anomalies (CANGRD) observation-based dataset. The modeled surface air temperature was evaluated against CANGRD data, the modeled albedo against MODIS data, and the modeled snow water equivalent against Canadian Meteorological Centre (CMC) and Global Snow Monitoring for Climate Research (GlobSnow) data. The lake simulation itself was evaluated using the Along Track Scanning Radiometer (ATSR) Reprocessing for Climate: Lake Surface Water Temperature and Ice Cover (ARC-Lake) dataset. Summer surface lake temperatures and the lake ice formation and breakup periods were well simulated, except for slight warm/cold summer/fall surface temperature biases, early ice breakup, and early ice formation, consistent with warm/cold biases in the climate simulation. Tests were carried out to investigate the sensitivity of the CSLM simulation to the default values assigned to the shortwave extinction coefficient and the average lake depth, and changing the former from 0.5 to 2.0 m−1 and the latter from 10.0 to 50.0 or 5.0 m had minimal effects on the simulation. Comparisons of the average annual variations of the simulated net shortwave radiation, turbulent fluxes, snowpack, and maximum and minimum daily surface temperatures between the land and the lake fractions for tundra, boreal, and southern regions showed patterns consistent with those expected. Finally, a test of the CSLM over the large resolved lakes in the model domain demonstrated a performance comparable to that for subgrid lakes.</jats:p>

section_content: Inland water bodies are regionally important components of the climate system. Their large thermal mass leads to local damping of the diurnal and seasonal temperature cycles (e. g., Martynov et al. 2012). The development of a warm surface layer over the course of the summer leads to large surface evaporation fluxes in the fall when the surface of the lake cools and overturning of the water mass occurs. 

Snow-covered lakes have a much higher albedo than adjacent forested areas, leading to marked differences in the surface energy balance (e. g., Verseghy et al. 2017). The effects of large lakes on local weather and climate have long been recognized (e. g., Notaro et al. 2013; Balsamo et al. 2012; Steenburgh et al. 2000; Lofgren 1997 ). 

As a result, efforts have been underway for a number of years to incorporate parameterizations for inland lakes into regional and global climate models. To avoid undue computational complexity, these efforts have generally made use of 1D (single column) lake models. 

Two such models have primarily been used for this purpose: the Freshwater Lake model (FLake ; Mironov 2008 ) and the Hostetler model (Hostetler et al. 1993). Bonan (1995) analyzed the effect of coupling the Hostetler model to the National Center for Atmospheric Research Community Climate Model, version 2 (CCM2). Samuelsson et al. (2010) studied the effect on the simulated European regional climate of adding FLake to the Rossby Centre Regional Atmospheric Model, version 3. 5 (RCA3. 5). Martynov et al. (2012) coupled both FLake and the Hostetler model to the fifthgeneration Canadian Regional Climate Model (CRCM5), demonstrating strengths and weaknesses in the lake simulations of each over North America. Bennington et al. (2014) evaluated and improved the simulation of large lakes by the Hostetler model in the International Centre for Theoretical Physics (ICTP) Regional Climate Model, version 4 (RegCM4). Le Moigne et al. (2016) coupled FLake to the global climate model of the Centre National de Recherches Météorologiques (CNRM) and investigated the effect of doing so on the climate simulation. 

In the above studies, attention has tended to focus on the simulation of large lakes, and the lake models used have generally only been evaluated using local-scale or in situ observational data. Le Moigne et al. ( 2016 ) evaluated FLake at a global scale using the satellitederived Along Track Scanning Radiometer (ATSR) Reprocessing for Climate: Lake Surface Water Temperature and Ice Cover (ARC-Lake) dataset (MacCallum and Merchant 2013), but their analysis was done on a per-lake basis and again focused on comparatively large lakes (. 500 km 2 ). More recently, attention has begun to turn to water bodies that are subgrid in scale, that is, those that occupy less than 50% of any given climate model grid cell. Even at regional climate model scales (typically tens of kilometers), in certain parts of the world such as Canada and northern Europe there are thousands of small lakes that remain unresolved, and in some areas the fractional coverage of such lakes can exceed 40%. There has therefore been growing interest in investigating the effects of such lakes in climate models. 

In this paper, we present a study evaluating the performance of the recently developed Canadian Small Lake Model (CSLM), version 2 (MacKay et al. 2016, manuscript submitted to J. Hydrometeor. ), for the first time on a regional scale, and analyzing the effect of adding it in a regional offline simulation of the land surface climate over western Canada using the Canadian Land Surface Scheme (CLASS), version 3. 6. 1 (Verseghy 2017). Sensitivity tests of background parameters are also presented. This study extends the work presented by MacKay et al. (2016, manuscript submitted to J. Hydrometeor. ), which reports testing of the CSLM at the local scale using field data from the Experimental Lakes Area in northern Ontario. A major difference between version 1 of the CSLM (MacKay 2012) and version 2 is the addition of the CLASS snow model above the lake ice, so particular attention is paid to the wintertime simulation. Like Le Moigne et al. (2016), we make use of the ARC-Lake data in our evaluation of the CSLM, but focus on subgrid-scale lakes, which the CSLM was designed to model.
section_name: Introduction

section_name: Models

section_content: The CSLM is a one-dimensional thermodynamic lake model with a surface mixed layer scheme based on the integrated turbulent kinetic energy (TKE) approach of Rayner (1980) and Imberger (1985). In this scheme the epilimnion is mixed via well-known processes (surface stirring, buoyancy production, shear production, etc. ) with conventional parameterizations, but the closure assumption of Rayner (1980) retains the nonsteady term in the TKE budget so that a truly time-varying diurnal mixed layer is modeled. In addition, a fully nonlinear energy balance is computed in a thin (currently set to 5 cm) ''skin'' layer facilitating a faster response to sudden changes in atmospheric forcing. Complete details can be found in MacKay (2012). 

These features address some weaknesses in the FLake and Hostetler models, especially regarding the simulation of small lakes in global and regional climate and numerical weather prediction applications. The Hostetler model does not represent surface mixing processes per se but rather parameterizes an eddy diffusivity based on a current profile deduced from Ekman theory (Henderson-Sellers 1985). Such an approach has been found to produce inadequate mixing, partly due to the scheme's inability to simulate individual storm mixing events (Bennington et al. 2014; Martynov et al. 2012 ). In addition, Ekman theory cannot be relevant for the millions of small lakes whose length scales must be smaller than the internal Rossby radius (e. g., Downing et al. 2006; Verpoorter et al. 2014). For such lakes other processes (e. g., seiching and inflows) must govern the current profile. 

FLake takes a similar bulk layer approach to the surface mixed layer as in CSLM, but it does not retain the detailed structure in the metalimnion that ensues after the various daily excursions of the diurnal mixed layer. Instead, the metalimnion below the surface mixed layer is represented by an empirical fourth-order polynomial, and the hypolimnion is absent. As such, and as noted by several authors, the model is relatively sensitive to the specified lake depth and is really only appropriate for lakes shallower than about 50 m (Mironov et al. 2010; Samuelsson et al. 2010). The advantage of such a two-layer approach is computational efficiency, but applications (e. g., ecological) that require more detailed thermal structure beneath the surface layer might not be well served. 

For lakes that freeze seasonally, snow and ice cover will be important aspects for lake simulation. Neither FLake nor the Hostetler model currently includes a very sophisticated treatment of snow processes, though both models do produce ice. Recently, the complete snowpack physics of CLASS was added to the CSLM along with a new scheme for the production of snowice, that is, ice that is formed when the weight of snow depresses the lower portion of the snowpack below the water line, causing the inflow and freezing of lake water (MacKay et al. 2016, manuscript submitted to J. Hydrometeor. ). That study demonstrated that the new version of the model simulated realistic development of both the mixed layer and wintertime evolution of snow and ice.
section_name: a. The Canadian Small Lake Model

section_content: CLASS was originally developed in the late 1980s for the Canadian Global Climate Model and has been under continuous development since then (Verseghy 2017). It simulates the energy and water balances of the land surface vegetation, soil, and snow (and also includes a parameterization for land ice). The vegetation canopy is treated as a single layer with phenological characteristics simulated separately for the major vegetation classes of needleleaf trees, broadleaf trees, crops, and grass. Background parameters such as maximum and minimum plant area index, albedo, and rooting depth are assigned according to the vegetation types present on each grid cell. Radiation transmission and precipitation interception are dependent on leaf area index, and transpiration is regulated by temperature, solar radiation, vapor pressure, and water stress. Intercepted snow is unloaded at an exponentially decreasing rate. The soil column in its standard configuration is discretized as three layers of thicknesses 0. 10, 0. 25, and 3. 75 m (although since the release of version 3. 0, CLASS has had the flexibility to handle any number of soil layers to whatever depth the user specifies). Heat transmission, water infiltration and redistribution, and soil freezing are modeled using physically based algorithms. Snow is modeled as a single thermal layer with a bulk temperature and a surface skin temperature. The snow albedo decreases exponentially and the density increases exponentially with time. The snow thermal conductivity is a function of density. 

Meltwater produced at the surface infiltrates into the snow and refreezes until the snowpack reaches isothermal conditions at 08C, at which point meltwater can reach the soil surface and infiltrate. Grid cells are divided into subareas of vegetation over soil, vegetation over snow, snow over soil, and bare soil; vertical fluxes are calculated separately for each subarea.
section_name: b. The Canadian Land Surface Scheme

section_content: The regional domain selected for this study was centered on western Canada (Fig. 1 ). This region represents a globally significant area of high subgridscale lake coverage. The selected model grid spacing was 0. 258 (about 25 km), a standard spacing for regional climate models. At this spacing, many grid cells have fractional subgrid lake coverages of 0. 3 or more (Fig. 2 ). The domain stretches from the high Arctic tundra in the north, through the wide band of boreal forest running from the northwest to the southeast, to croplands and grasslands in the south, and is bounded on the west by the Rocky Mountains, the Coast Mountains, the Mackenzie Mountains, and neighboring ranges. 

The atmospheric forcings required to run CLASS and the CSLM over this domain were obtained by dynamically downscaling data from the European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (Dee et al. 2011). The model used for the downscaling was CRCM5 (Zadra et al. 2008) ; the technique is described in Martynov et al. (2013). By this means, forcing data for our offline run were obtained at a temporal interval of 15 min (compared with 6 h for the ERA-Interim data) and a grid spacing of 0. 258 (compared with 0. 758 for the ERA-Interim data). The saved fields were those required to run CLASS and the CSLM: incoming shortwave and longwave radiation; atmospheric temperature, humidity, and wind speed at the lowest model level; the height of the lowest model level; fractional cloud cover; rainfall and snowfall; and surface pressure. A modeling period from June 1990 to May 2011 was selected, with the first year treated as a spinup interval. 

Gridded land cover and soil information were generated over the model grid as background information for CLASS. The land-cover data were obtained from the 1-km Global Land Cover 2000 (GLC2000) dataset (Bartholome and Belward 2005). Over Canada, soil texture and depth were obtained from the ''Soil Landscapes of Canada'' dataset published by the Canadian Soil Information Service of Agriculture and Agri-Food Canada, and over the United States from the USGS dataset (Miller and White 1998). The initial soil temperature was set to the average forcing air temperature for each grid cell over the modeling period, and initial soil moistures were set to field capacity. 

For the CSLM, the fractional lake coverage (shown in Fig. 2 ) was provided by the GLC2000 dataset. Since we are primarily interested in subgrid-scale lakes, for the initial part of this model validation exercise we neglected large lakes, that is, grid cells where the lake coverage is 100%. Thus, over each grid cell the CSLM was run for the lake fractional area and CLASS for the land fractional area. The initial lake temperature was set to 148C everywhere, since offline tests have shown that the modeled temperature profile took less than a year to spin up. The CSLM requires assigned gridded values for three additional parameters: the average lake horizontal length scale, the extinction coefficient for shortwave radiation, and the average lake depth. The length scale for each distinct lake at the original GLC2000 grid spacing of 1 km was obtained as the square root of the lake area, from which the average length scale for each cell on the model grid was obtained by averaging over the lakes present within the grid cell. For the extinction coefficient, measured values are not available globally, so we chose a set value of 0. 5 m 21, consistent with the value used by Bennington et al. (2014) for the Hostetler model and by Le Moigne et al. (2016) for FLake. For the lake depth we made use of the global dataset of Kourzeneva (2010), but since over our domain it only provides depth data for the largest lakes, we had to select a default value. We chose 10. 0 m, the same as that adopted in the studies of Samuelsson et al. (2010), Martynov et al. (2012), and Le Moigne et al. (2016).
section_name: Modeling domain and input data

section_content: Table 1 gives an overview of the datasets used for validation. To evaluate the precipitation and temperature fields, we made use of the Canadian Gridded Temperature and Precipitation Anomalies (CANGRD) dataset (Milewska et al. 2005) 3a. Corrections have been applied to account for systematic errors such as gauge undercatch and to address the effects of observation system changes, station relocation, and spatial inhomogeneities. A statistical optimal interpolation procedure was used to generate monthly grids of temperature and precipitation anomalies (departures from 1961-91 normals) on a 50-km polar stereographic grid. Actual temperature and precipitation fields were then generated by combining the anomalies with the 1960-91 normals, gridded using the ANUSPLIN model (based on thin plate smoothing splines, developed at the Australian National University). For ease of comparison, the CANGRD data were interpolated to our model grid by means of bilinear interpolation. 

To assess the modeled albedo, we made use of the whitesky albedo product derived from MODIS data, available from the year 2000 at a global grid spacing of 0. 058 (Schaaf et al. 2002). The data were transferred to the model grid using box averaging and stored as monthly averages. To minimize the problem of cloud cover contamination of the MODIS images, only high-quality (grade 2 or better) retrievals were included in the averaging. 

To evaluate the modeled snow water equivalent (SWE), we made use of two datasets. The first was derived from the daily global snow depth analysis that has been produced by the Canadian Meteorological Centre (CMC) since 1998 and archived on a 24-km grid (Brown and Brasnett 2010). It is based on a simple snow accumulation and melt model forced using reanalyzed temperatures and forecast precipitation from the Canadian Global Environmental Multiscale (GEM) weather prediction model. 

Observed daily snow depths from surface stations are incorporated into the product using optimal interpolation. The locations of the surface stations within our model domain are shown in Fig. 3b. Red stars indicate stations reporting snow depth on 15 April in any year between 2000 and 2009; as an arbitrary example, blue circles indicate stations reporting on that date in 2007. The number of stations reporting on any given date in a year can vary between 30 and 50, depending on snow conditions. SWE is derived from the snow depth using a climatologically derived lookup table based on density (Brown and Mote 2009). The SWE fields were interpolated to our model grid using the nearest-neighbor approach and stored as monthly averages. 

The second SWE dataset that we used is the European Space Agency (ESA) Global Snow Monitoring for Climate Research (GlobSnow) product, version 2
section_name: Validation data a. Validation of climate simulation

section_content: To evaluate the CSLM simulation over our domain, we made use of the recently developed ARC-Lake database (MacCallum and Merchant 2012, 2013). ATSR is a remote sensing instrument installed on board ERS-2 and Envisat. ARC-Lake is a project funded by the European Space Agency to derive lake temperature and ice cover globally; version 3. 0 of the product covers the period 1995-2012. The data are available either on a FIG. 5. As in Fig. 4, but for seasonally averaged near-surface temperature T (8C). 0. 058 grid, or on a per-lake basis, from daily daytime and nighttime satellite passes. Since substantial spatial and temporal gaps exist in the data, we first computed monthly daytime and nighttime averages over the 0. 058 grid, then averaged the two together and regridded the data to our model grid using box averaging.
section_name: b. Validation of lake simulation

section_name: Results and discussion

section_content: Prior to evaluating the performance of the lake model, we examined the overall climate simulation, as biases in the background forcing fields will naturally affect the surface fields. The precipitation field in particular is of vital importance for the winter evolution of the snowpack; without realistic snowfall forcing it is impossible to obtain a realistic snow simulation. We therefore adopted the method of Verseghy et al. (2017) in applying an adjustment factor to the precipitation forcing field prior to performing the model run, on the basis of observation-based data. For each month k of the simulation period, an adjustment factor F adj was calculated as the ratio of the CANGRD precipitation P CANGRD divided by the model-downscaled precipitation P mod for that month:

This factor was then used to scale the precipitation at each model time step within the month, so that for each month over the simulation, the monthly average modeled precipitation was exactly equal to the CANGRD value. Verseghy et al. (2017) made use of the CRU Time Series, version 3. 21 (CRU TS3. 21), dataset from the Climate Research Unit of the University of East Anglia (Harris et al. 2014) to perform this adjustment, but they found that it tended to overcorrect the precipitation field, resulting in a distinct underestimation of the snowpack in some areas. We also found, in preliminary trial runs over our domain, that precipitation adjustment using the CRU data resulted in a serious underestimation of the snowpack over areas in Alaska such as the Brooks Range. For this reason, we selected the CANGRD precipitation data to calculate the monthly adjustment factor. Although it is only available over Canada, for that part of the domain it should provide more realistic precipitation forcing fields, as it incorporates corrections for snow gauge undercatch that are not included in the CRU data (see section 4a). This is important since, as will be shown below, most of the lake validation data are in fact located in Canada. Figure 4 shows the seasonal average unadjusted precipitation compared with the CANGRD values over the model domain for the 20-yr modeling period. It can be seen that the precipitation adjustment results in large corrections over the Coastal, Rocky, and Mackenzie Mountains, together with more modest downward adjustments in the fall and spring in the high Arctic and the southern plains area, respectively, and an upward adjustment in the summer in southern Saskatchewan and western Ontario. 

Turning to the outputs of the model run, Fig. 5 shows the seasonal average modeled surface air temperature compared with CANGRD over the 20-yr modeling period. It can be seen that over the continental interior, the biases are in general on the order of 1. 08C, with a larger cold bias in the Northwest Territories and Nunavut in the fall, and larger warm biases in the high Arctic in the winter and spring, in western Ontario in the winter, and in the southern plains in the summer. All year-round, but especially in the winter and spring, large apparent cold biases sometimes exceeding 5. 08C are found over the mountainous regions. Complex topography is difficult to model well and observations in mountain areas are sparse and strongly biased toward valleys, so for the purposes of this study we excluded from the analysis areas of elevation exceeding 800 m. To examine the effects of lakes on an ecosystem basis, we further defined three major ecozones, shown in Fig. 6 : the tundra ecozone, consisting of areas north of 608N with vegetation coverage less than 50%; the southern ecozone, consisting of areas south of 558N with vegetation coverage dominated by crops and grasses; and the intervening FIG. 6. Locations of tundra, boreal, and southern ecozones used for averaging of surface fields. Mountainous regions are excluded from the averaging, as is Alaska, because of the unknown magnitude of biases in the simulated climate in these areas due to the sparse nature of the surface observations. boreal forest ecozone, dominated by evergreen needleleaf forest. We included in the southern and boreal ecozones areas south of the Canada-U. S. border, since according to Fig. 4 the precipitation adjustment in these regions would probably have been small, but we excluded Alaska from the boreal and tundra ecozones since it was uncertain what magnitude of precipitation adjustment might have been required there (and the fractional lake coverage, as seen in Fig. 2, is in any case small). 

Figure 7 shows plots of the modeled average surface albedo compared with MODIS for two selected months: April, as representative of the spring peak snow period, and August, as representative of midsummer conditions, over the 10-yr period of overlap between the modeling period and the MODIS data (2000-10). Figure 7 (top) shows the modeled values, Fig. 7 (middle) shows the differences between the modeled and MODIS values, and Fig. 7 (bottom) shows the differences between a second model run neglecting the presence of lakes and the MODIS values. Figure 7 (bottom) is included to provide a comparison with the results reported in Verseghy et al. (2017). For a modeling domain centered on eastern Canada, they found a general negative albedo bias over boreal forest areas during the winter snow cover period, which was halved when the lake areas were treated as bare soil with a bright snow cover. In our case, when we included lakes, which are similarly bright in the winter, our April albedo bias in these areas was likewise halved. (Generally, when climate models are run without including subgrid lakes, the vegetation fractions are normalized by 1 minus the lake coverage, thus in effect extrapolating the dark boreal forest over lake areas. ) In August, conversely, including lakes has little effect on the albedo of boreal forest areas (since both forest and lakes are dark), but causes a slight improvement in the albedo of the tundra area. In this region we have a substantial positive albedo bias, exceeding 0. 2 in places. Verseghy et al. (2017) reported the same high albedo bias over snow-free tundra areas, and we trace this to the high values assigned in the model lookup table to dry sandy soil, dominant in this area, which in retrospect are actually more characteristic of desert regions. Since the vegetation coverage over much of the tundra area is less than 25%, the bare soil albedo will have a large effect. 

Finally, we compared the modeled evolution of SWE with the observation-based CMC and GlobSnow data for the 12-yr period of overlap between the modeling period and the two datasets (1998-2010). Verseghy et al. (2017) found the CLASS-modeled SWE to be consistently higher than the CMC values over their eastern Canada modeling domain, but a comparison with GlobSnow was not possible because of the large wintertime snow accumulation in the Quebec interior coupled with the fact that the GlobSnow product begins to underestimate SWE at values greater than 100-150 mm due to signal saturation (Dong et al. 2005). Over our western Canada domain, with its lower snow accumulations (excluding the mountains), we were able to compare the modeled SWE with both the CMC and the GlobSnow data, and on that basis we found that the model performed very well. Figure 8 (right) shows the modeled average annual variation of SWE for the tundra, boreal, and southern zones compared with the CMC and the GlobSnow data, and it can be seen that in all cases the modeled values lie between the two observationbased datasets except for March and April in the southern zone, when the modeled SWE is overestimated by a few millimeters. The total SWE for the model run is slightly less than the SWE averaged over the land only, because of the ongoing conversion of snow to snow-ice over lakes, as described in section 2a. 

Both lines are included in the plot because the CMC product, being based on a simple accumulation-melt model coupled with station observations, should presumably be more comparable to the CLASS land-only SWE, while the GlobSnow SWE, being a remote sensing product, should be more comparable to the total land-plus-lake SWE. On the other hand, it may be that the GlobSnow passive microwave retrieval algorithms are not able to distinguish very well between snow and snow-ice, or for that matter lake ice, which complicates the analysis. It must further be noted that as shown in Fig. 8 (left), the spatial differences between the CMC and the GlobSnow SWE fields are quite substantial. The March simulated SWE, for example, shows a local maximum around Reindeer Lake, which is present in the GlobSnow SWE but not in the CMC SWE. Conversely, the simulated SWE shows some lake-effect snow accumulation to the east of Great Bear Lake, which occurs (to a lesser extent) in the CMC SWE but is absent from the GlobSnow SWE. Bearing in mind these uncertainties, however, on the whole the SWE simulation seems satisfactory.
section_name: a. Climate simulation

section_content: We turn now to a direct evaluation of the performance of the CSLM. To this end, we made use of the ARC-Lake dataset described in section 4b above. We compared the modeled subgrid lake surface temperature during the comparatively ice-free months of July-September over our domain, and the evolution of the lake ice cover during the spring breakup period of May-June and the fall freeze-up period of October-November, for the 16-yr interval of overlap between the modeling period and the dataset (1995) (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) (2006) (2007) (2008) (2009) (2010) (2011). Figure 9 shows the monthly average subgrid lake temperatures for August and scatterplots and statistics for July-September. It can be seen that the model simulates the seasonal progression of lake-surface warming very well, with RMSE values on the order of 28C. The scatterplots show some areas of temperature overestimations in July and August, changing to an overall mean bias error (MBE) of 21. 88C in September, consistent with the warm summer surface air temperature bias of 18-28C across the model domain and the onset of negative temperature biases in the fall (Fig. 5 ). Values of the coefficient of determination r 2 are very high, around 0. 9, in August and September; they are lower in July, but this reflects the smaller number of data points for that month (given the fact that lakes north of 608N are still largely ice covered). Figure 10 shows the monthly average lake ice cover for June and October, and Table 2 shows the statistics for May, June, October, and November. (Scatterplots are not very informative in this case, as there is such a strong 0-1 dipole signal in the ice cover. ) Consistent with the overall summer high surface air temperature bias (Fig. 5 ), the modeled ice cover in June shows a tendency to disappear slightly too early. However, the statistics for May and June indicate excellent performance on the part of the model, with MBE values of magnitude 0. 11 or less, RMSE values of around 0. 2, and r 2 values of around 0. 8. In the fall, conversely, the modeled ice cover tends to form somewhat too early, as shown in Fig. 10. As with the lake surface temperatures, this is consistent with the fall cold bias over much of the domain (Fig. 5 ). The r 2 values are lower (dropping to 0. 26 in November), but the MBE and RMSE values are still comparatively small. We may conclude that, given the known biases in the modeled climate, the CSLM is performing very well in its simulation of surface subgrid lake temperature and ice cover.
section_name: b. Subgrid lake simulation

section_content: As noted in section 3, two of the background parameters for the CSLM, the shortwave extinction coefficient and the average lake depth, were assigned default values of 0. 5 m 21 and 10. 0 m respectively, since there are no global datasets of the former available, and since our FIG. 11. Average annual variation over the period 1991-2011 of land (solid line) vs lake (dashed line) variables over tundra, boreal, and southern areas: net shortwave radiation, snow water equivalent, sensible heat flux, latent heat flux, daily max surface temperature, and daily min surface temperature. lake depth database includes only the larger lakes in the model domain. We therefore conducted a few tests to investigate the sensitivity of the lake surface temperature and ice cover simulation to varying these two parameters. Since the value of 0. 5 m 21 that we used for the extinction coefficient is characteristic of a fairly clear lake, we reran the simulation with a default value of 2. 0 m 21 instead, representing a more turbid lake. The statistics corresponding to Figs. 9 and 10 and Table 2 are listed in Table 3. It can be seen that for the July-September temperatures, the magnitude of the MBE varies little, remaining at around 0. 18, 0. 28, and 2. 08C; the RMSE varies by at most 0. 38C, and the difference in the r 2 values is negligible. For the May-June and October-November ice coverages, the MBE and RMSE vary by at most 0. 05, and the r 2 by at most 0. 08. This confirms the results of earlier single-column sensitivity tests, which showed that the effect of varying the extinction coefficient is seen mainly in the deeper lake temperatures. 

For the average lake depth, we tried a test increasing the default value from 10. 0 to 50. 0 m, which is close to the value often used as the limiting depth for FLake (e. g., Martynov et al. 2012; Le Moigne et al. 2016). We also tried a test decreasing the default value to 5. 0 m, to see whether this would enhance the annual amplitude of the lake surface temperature as found with other single-column models, even though doing so tends to violate the assumption of zero heat flux at the lake bottom (made in the absence of modeled heat conduction from the bottom sediments): with a lake depth of 10. 0 m, less than 0. 7% of the incoming shortwave radiation at the surface of the lake reaches the bottom, whereas for a lake depth of 5. 0 m, the percentage is greater than 8%. The statistics corresponding to those in Table 3 are shown in Tables 4 and 5. It can be seen that for the 50. 0-m lake depth test, the lake temperature MBE varies from the base run by less than 0. 068C, the RMSE by 0. 078C, and the r 2 imperceptibly, and the ice cover statistics vary by at most 0. 02 throughout. For the 5. 0-m test, the greatest differences are found in the September temperature statistics, with MBE and RMSE both varying by 0. 68C; for most of the other statistics, the differences are 0. 05 or less. It is therefore evident that the simulation is not particularly sensitive to variations in the lake depth beyond 5. 0 m. As with the variation in the extinction coefficient, this had been confirmed earlier by offline single-column tests.
section_name: c. Sensitivity of simulation to lake parameters

section_content: Previous studies (e. g., Bonan 1995) have investigated the differing annual regimes of various climate parameters for land and lake areas, demonstrating how changes in surface fluxes and temperatures can be expected for the overall land areas as the lake fraction increases. By way of comparison, we show in Fig. 11 the annual variation in surface net shortwave flux, sensible heat flux, latent heat flux, SWE, and daily maximum and minimum temperatures for the land tile and the lake tile separately from our model run, averaged over the tundra, boreal, and southern parts of the domain. Note that the effect of the lake tile is muted in offline runs such as the present one, with prescribed atmospheric forcings; when run in coupled mode with a climate model, when surface feedbacks are included, the effect of subgrid lakes can be expected to be considerably larger, as shown, for example, by Le Moigne et al. (2016). The net shortwave radiation plot shows summertime values for lakes that are consistently slightly higher than for land, reflecting the generally lower albedo of lakes; the difference is smallest for the boreal zone, where the albedo difference is the least. In the winter, the net shortwave radiation over lakes is lower than over land, but the differences are small for the southern and tundra areas, where the vegetation is partly or completely buried by snow. For the dark boreal forest region the differences are much larger (as discussed in section 5a), by up to 63% in March. For SWE, the land values are considerably larger than the lake values over the whole year, by up to more than a factor of two for the boreal and tundra areas, because of the ongoing conversion of snow to snow-ice described in section 2a. The latent heat flux for land peaks in the early summer for all three regions, when vegetation is actively growing and water availability stresses are generally low. Conversely, for lakes the latent heat flux peak occurs in the late summer or early fall, when cooler weather in conjunction with the large thermal inertia of the lakes leads to a large pulse in the surface turbulent fluxes. The sensible heat fluxes show similar patterns, with peak land values occurring in the early summer; the boreal forest fluxes are the earliest to ramp up because of the large surface roughness. For lakes, the sensible heat fluxes, like the latent heat fluxes, peak in the fall, with a slight lag as the lakes cool. 

Overall, the Bowen ratios are higher over the land, reflecting the vegetation stomatal resistance to transpiration that is absent over water bodies. Finally, the daily average maximum and minimum temperature patterns demonstrate the lake thermal inertia effect during the summer period, with land maxima/minima consistently higher/lower than those for lakes, by as much as 108C. During the period of winter snow cover the land-lake differences are much smaller, with lake maximum and minimum temperatures consistently lower than land values, with a maximum difference of 48C for boreal regions where the land-lake contrast is the greatest. 

e. Applying the CSLM to large lakes

The CSLM was designed specifically to model small lakes, and as a single-column model it is naturally not equipped to properly account for the effects of horizontal advection and lake-wide circulations that occur in large lakes. Previous studies (e. g., Martynov et al. 2012; Bennington et al. 2014) have suggested that such large lakes probably need to be represented using 3D models [such as the Great Lakes Ice-Circulation Model (GLIM), described in Wang et al. (2010) ], so that such processes can be included. However, because other single-column lake models such as FLake and the Hostetler model have nevertheless been applied and evaluated over large lakes, we decided to redo our model run, this time including the large inland lakes that had previously been excluded (shown as masked out in Fig. 2 ), in order to investigate whether the model performance was in fact markedly different. A list of the largest lakes in our We may conclude that for this modeling domain at least, the realism of the simulation of the large lakes is quite comparable to that of the subgrid lakes.
section_name: d. Differences between land and lake surface fluxes and variables

section_content: We have presented the first regional offline test of the Canadian Small Lake Model, in conjunction with the Canadian Land Surface Scheme. The modeling domain chosen was centered over western Canada, a globally significant area of high subgrid lake coverage. Atmospheric forcing fields for the CSLM and CLASS were obtained by dynamically downscaling ERA-Interim data using CRCM5. The rainfall and snowfall fields were adjusted to match monthly totals from the CANGRD observation-based precipitation and temperature database. Background vegetation and lake coverage fields were obtained from the GLC2000 database, and the lake depths were assigned on the basis of the Kourzeneva (2010) dataset, where available. The lake shortwave radiation extinction coefficient and the default average lake depth were assigned values of 0. 5 m 21 and 10. 0 m, respectively. 

Evaluation of the simulated surface air temperature field against CANGRD data showed moderate biases over the continental interior, but large negative biases over mountainous regions; these were therefore excluded from the subsequent analysis. The modeled spring and summer albedo fields compared with MODIS data were shown to be consistent with those reported by Verseghy et al. (2017) over an eastern Canada domain, with the addition of lakes improving the springtime simulation over the boreal forest region and the summertime simulation over the tundra region. The snow water equivalent was well simulated, with the seasonal evolution of the snowpack lying between observation-based CMC and GlobSnow values, again reinforcing the findings of Verseghy et al. (2017). 

The surface temperature and ice cover of subgrid lakes over the domain, as simulated by the CSLM, were compared with data from the ARC-Lake dataset. The model showed a slight tendency to overestimate early summer and to underestimate early fall lake surface temperatures, and to simulate ice breakup and formation too early, but as these were consistent with temperature biases in the climate simulation, overall the model was deemed to have performed well. Tests showed little sensitivity to the values assigned as defaults to the shortwave light extinction coefficient and the average lake depth. Overall land-lake differences showed expected patterns: higher net shortwave radiation over lakes in summer overall, and lower in winter for boreal forest areas; lower SWE values over lakes because of the conversion from snow to snow-ice; maximum turbulent fluxes from land in the summer and from lakes in the fall, owing to the delayed release of heat storage from lakes; and damped surface maxima and minima over lakes compared with land, likewise owing to the greater thermal inertia of lakes. A final test was carried out to ascertain whether the CSLM could do a reasonable job of simulating large lakes, despite having been developed primarily for subgrid lakes, and its performance turned out to be quite comparable to that for small lakes. Our findings thus demonstrate the robustness of the CSLM and its applicability for future operational implementation in regional climate and numerical weather prediction models. 

Acknowledgments. Thanks are due to Katja Winger at l'Université du Québec à Montréal, who performed the CRCM5 downscaling runs; to Daniel Robitaille, who interpolated the observation-based data onto the model grid, converted them to a standard format, and performed the initial single-column sensitivity tests of the lake model; and to Ed Chan, who provided muchappreciated support with software applications and with data manipulation and visualization. Richard Harvey reviewed an initial draft of this paper and supplied helpful comments.
section_name: Summary and conclusions

section_name: FIG.

section_content: 12. Monthly average large lake surface temperature (8C): (top) ARC-Lake August, (middle) modeled August, and (bottom) scatterplots of modeled vs ARC-Lake for July, August, and September.
section_name: FIG.