Technical challenges and solutions in representing lakes when using WRF in downscaling applications

The Weather Research and Forecasting (WRF) model is commonly used to make high-resolution future projections of regional climate by downscaling global climate model (GCM) outputs. Because the GCM fields are typically at a much coarser spatial resolution than the target regional downscaled fields, lakes are often poorly resolved in the driving global fields, if they are resolved at all. In such an application, using WRF’s default interpolation methods can result in unrealistic lake temperatures and ice cover at inland water points. Prior studies have shown that lake temperatures and ice cover impact the simulation of other surface variables, such as air temperatures and precipitation, two fields that are often used in regional climate applications to understand the impacts of climate change on human health and the environment. Here, alternative methods for setting lake surface variables in WRF for downscaling simulations are presented and contrasted.


Introduction
When using global climate model (GCM) fields to drive finer-scale regional climate model (RCM) runs, typically the RCM does not have an oceanic or lake physics component and relies on the GCM output to provide all water surface temperatures and ice cover.Within a downscaling simulation, by design, the GCM is at a coarser spatial resolution than the RCM, so inland water bodies in the region being simulated are either poorly resolved or not resolved by the GCM.Prior to 2013, the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008) required exogenously prescribed water surface temperatures, as there was no capability to prognosticate water temperatures.WRF has included an optional coupled ocean component since version 3.5 was released in April 2013 (WRF User's Guide, 2014).Other RCMs have been coupled to ocean models in order to simulate regions around the Arctic, Mediterranean Sea, and Indian Ocean (e.g., Rinke et al., 2003;Ratnam et al., 2009;Artale et al., 2010;Gualdi et al., 2013).However, when using WRF's default configuration, the sea surface temperature (SST) fields used during the simulation are calculated from the driving data during the preprocessing steps performed before WRF runs the simulation; during the model run, these prescribed water temperatures are input at a userspecified frequency which is usually daily or sub-daily.Similarly, lake surface temperatures (LSTs) and lake ice cover are prescribed by spatial interpolation from the SST and sea ice fields in the driving data.In this study, we examine the use of the Advanced Research WRF (Skamarock and Klemp, 2008) model applied as an RCM in regions where the driving larger-scale data have a poor representation of lakes.
When the WRF Preprocessing System (WPS) interpolates skin temperatures from the coarser global data set (where both land and water temperatures are included in a single field), masks are applied such that water temperatures from the GCM are used to set water temperatures on the finer, target grid.Using the standard methods in WPS, interpolation is first attempted using 16 surrounding grid cells in the coarser grid; if this method fails due to a lack of the requisite 16 valid data points, WPS attempts other interpolation Figure 1.The ocean mask from the 1 • CESM data (which is used by WPS to determine the locations of land and water points from CESM), as shown in the area corresponding to a WRF 36 km continental US domain (left), and the 36 km WRF grid's land-water mask (right).Labels are placed to indicate the locations of Lake Superior ("S"), Lake Michigan ("M"), Lake Huron ("H"), Lake Erie ("E") and Lake Ontario ("O"), as well as Hudson Bay ("HB").
Figure 2. The skin temperature (K) processed from CESM to the 36 km WRF grid using WPS and valid at 00:00 UTC, 1 December 1994.White circles indicate the locations of Pyramid Lake, Great Salt Lake, and Lake Sakakawea, from west to east, respectively.
techniques using as many as four grid cells and as few as one.While a full description of all WPS interpolation techniques is beyond the scope of this study, more information is available in the WRF User's Guide (2014, pp. 3-56 to 3-59).When all other methods fail due to the lack of nearby water grid cells, WPS defaults to the "search" approach, in which the nearest water point is used to set LSTs.When employing the search option, water cells in the driving data are often distant from and unrepresentative of the target cell in the WRF domain.The search option in WPS performs no interpola-tion or averaging, sometimes resulting in abrupt, nonphysical temperature discontinuities.
Here we show the result of using this default methodology to downscale 1 • Community Earth System Model (CESM) fields to a 36 km WRF domain (199×127) covering the continental US, and subsequently similar examples in other downscaling studies are discussed.However, it should be noted that the use of CESM as an example is arbitrary because similar results have been obtained with other global data sets as well.The CESM ocean mask, used to interpolate the GCM's SST fields to the WRF grid, has no water grid cells over the North American interior (Fig. 1).As a result, water temperatures in Hudson Bay are used to set temperatures over the larger westernmost areas of the Laurentian Great Lakes, while LSTs in the southeastern areas of the Great Lakes are set by Atlantic SSTs (Fig. 2).At the time shown in Fig. 2, the LSTs interpolated from CESM onto the 36 km WRF grid contain discontinuities of approximately 17 K between adjacent grid cells in Lake Michigan and Lake Huron, while a smaller discontinuity of approximately 3 K is created in Lake Superior.It should be noted that various interpolation options are available in WPS and can be specified by the user.The description in the paragraph above is representative of the interpolation process as defined by WPS's default settings.Even though this process could be changed by the model user, the key issue remains that when lakes are poorly represented or completely absent, the problem of how to specify the lake state is not amenable to any interpolation method.
The problems of using larger-scale data to define LSTs with the default options in WPS are not limited to the Great Lakes.None of the lakes resolved by WRF at 36 km have valid LSTs in the CESM ocean mask (Fig. 1).Using the search option in WPS results in setting the LSTs to unrealistic values throughout the domain.Temperatures in Pyramid Lake, Great Salt Lake, as well as several smaller lakes east of the Rocky Mountains in both Canada and the US are assigned from the Pacific Ocean (Fig. 2), while lake temperatures in the southeastern and central US are set from SSTs in the Gulf of Mexico and Atlantic Ocean.Two adjacent grid cells representing Lake Sakakawea in North Dakota are assigned LSTs differing by approximately 10 K because the western cell is set from the Pacific while the eastern cell is prescribed from Hudson Bay (Fig. 2).Using any interpolation method to assign LSTs when no suitable data are available will adversely affect the accuracy of downscaled simulations that are based on forcing from those LSTs.Mallard et al. (2014;hereafter M14) also discuss problems that arise when downscaling coarse global data to a 12 km grid covering the eastern US.In M14, the National Centers for Environmental Prediction (NCEP) -Department of Energy Atmospheric Model Intercomparison Project (AMIP-II) reanalysis (hereafter R2; Kanamitsu et al., 2002) is used to drive historical simulations as a proxy or stand-in for a similarly coarse GCM.In contrast to the CESM example discussed above, R2 has at least a partial representation of the western Great Lakes, but nevertheless has only three inland water points to represent all five of the Great Lakes (Fig. 1 of M14).Therefore, using the standard interpolation methods with R2 results in unrealistically large, abrupt, and nonphysical LST discontinuities in eastern Lake Erie and Lake Ontario, where water temperatures are set using Atlantic SSTs, while the LSTs in western Lake Erie and in the three western Great Lakes are interpolated from the three lake cells in R2 (M14).
In WRF, ice cover can either be interpolated from the driving data and assigned to cover some fraction of a grid cell, or it can be treated as a binary field that is set to 100 % at grid cells where the water surface temperature drops below a specified threshold.The default threshold value was 271 K (slightly below the freshwater freezing temperature of approximately 273 K), but it was changed to 100 K as of version 3.5.1 to avoid the unintended creation of ice by this method when using WRF's default settings (Table 1).When fractional ice values are prescribed from the driving data set, the WPS methods applied to interpolate sea and lake ice differ from those used for SSTs and LSTs.If there are no surrounding water grid cells in the driving data set, an ice cover value of zero is assigned rather than employing the search method.When M14 downscaled ice cover from R2, it was shown that ice concentrations of zero were applied to points through Lake Huron, Lake Erie and Lake Ontario throughout a 2-year simulation (Fig. 3 of M14), even though partial ice coverage was observed on all three lakes during that historical period.Moreover, almost complete ice coverage of Lake Superior and Lake Michigan occurred in a single day (M14).Wang et al. (2012b) conducted a climatology of ice cover in the Great Lakes over the period 1973 to 2010 and showed that, in the average seasonal cycle of ice cover, the maximum fractional coverage of Lake Superior was approximately 50 % (their Fig. 3).Although Wang et al. noted that the standard deviation of ice cover is quite large (exceeding the mean values in some of the Great Lakes), the seasonal cy-cles in their study showed the accumulation of ice coverage over months, not the abrupt appearance of lake-wide ice over daily periods.Ultimately, M14 improved the representation of the Great Lakes in their downscaled simulations by applying a coupled lake model, which will be discussed further in a subsequent section.Whereas M14 showed the results of using a single lake model, the current work presents a broader range of approaches, recognizing that the most preferable method to represent lake fields may vary between different RCM applications.
Prior studies downscaling other global data sets and GCMs have also noted findings similar to the example shown here (Fig. 2) and the results of M14.Using WRF as an RCM over eastern Africa, Argent (2014) showed that the use of WPS's default interpolation methods resulted in oceanic temperatures from a global SST data set applied to set LSTs throughout Lake Victoria.Discontinuities in LSTs with WRF were noted in the Great Lakes basin by Bullock et al. (2014), who downscaled R2 to 12 km, and by Gao et al. (2012), who downscaled CESM to a 4 km grid.Within the downscaled simulations produced for the North American Regional Climate Change Assessment Program (NARCCAP; Mearns et al., 2012), problems with producing realistic LSTs and ice cover for the Great Lakes region are documented using several approaches with various RCMs, including WRF (NAR-CCAP, 2014).For some NARCCAP model configurations, caution is recommended when using surface variables in the region surrounding the Great Lakes.Previous work examining the value of dynamical downscaling has noted that downscaled simulations have the most potential to add value relative to GCM simulations in areas of complex topography and along coastlines because of increased resolution in regional models (e.g., Feser et al., 2011).Although RCMs better resolve the coastlines (and therefore the presence of lakes) than the driving GCMs, using erroneous LSTs and lake ice cover could impair the simulation of interactions between lakes and overlying air masses.The potential benefits gained by downscaling to a grid spacing that better resolves land-water interfaces may not be realized if the lake state (defined here by LSTs and ice) is unrealistically represented.Even as additional computing resources allow GCMs to increase in resolution and better represent lakes, RCMs will also be run at finer scales; therefore, it can be expected that smaller lakes with important effects on mesoscale and microscale climatology will continue to be unresolved by the driving data sets.
The purposes of this paper are to describe various techniques that can be used to set LSTs and lake ice cover in the WRF model for downscaling, and to discuss the benefits and possible shortcomings of each approach.The effects of these techniques on simulated lake-atmosphere interactions, both in the present climate and in future climate states, are discussed in context with relevant previous literature.

WRF version Released
Updates of interest 3.3 April 2011 "Alternative initialization of lake SSTs" option included in WPS so users can set LSTs from temporally averaged 2 m temperatures.

3.5
April 2013 CLM available as an LSM within WRF, but with its lake model disabled.

3.5.1
September 2013 Default surface water temperature at which WRF prescribes ice ("seaice_threshold") is lowered from 271 to 100 K.

3.6
April 2014 CLM lake model available with any choice of LSM.Lake depths can be prescribed as a constant or as a spatially varying 2-D field.

Comparison of methods
As will be shown below, choice of the appropriate methodology for representing a lake in a downscaling configuration is dependent on what interactions must be simulated between the atmospheric fields and the lake state and how the lake state is expected to be impacted by climate change when downscaling future GCM projections.In regional climate simulations conducted over the continental US, the Laurentian Great Lakes are a prominent feature, as Lake Superior is the largest freshwater lake in the world (by surface area) at over 82 000 km 2 .Several studies have concluded that the Great Lakes strongly influence the surrounding regional climate, moderating extremes in near-surface temperatures, and affecting precipitation and passing cyclones and anticyclones on an annual cycle (e.g., Wilson, 1977;Bates et al., 1993;Scott and Huff, 1996;Notaro et al., 2013).Climatologically, the greater heat capacity of the lakes serves to enhance precipitation and convection during September to March, when warmer surface water (relative to low-level atmospheric temperatures) reduces atmospheric stability (e.g., Notaro et al., 2013).Conversely, the slower warming of the lakes in boreal spring results in the opposite effect during the April-August period, when the relatively cool lakes enhance atmospheric stability and reduce precipitation and convection.These periods are referred to as the lake unstable and lake stable seasons, respectively.Lake-effect precipitation has also been documented outside the Great Lakes as well, such as in Lake Champlain (Tardy, 2000;Laird et al., 2009), Lake Tahoe (Cairns et al., 2001), and the Great Salt Lake (Carpenter, 1993;Steenburgh and Onton, 2001).A review by Schultz et al. (2004) states that lake-effect snowfall has been observed to occur over lakes with fetches of only 30 to 50 km, citing prior studies over Bull Shoals Lake of Arkansas (Wilken, 1997) as well as Lake Tahoe and Pyramid Lake in Nevada (Cairns et al., 2001;Huggins et al., 2001).Interactions between the lakes and surrounding regions are also strong in tropical environments as well.For example, the immediate region surrounding Lake Victoria in Africa has the highest recorded frequency of thunderstorms in the world with approximately 300 storm days per year (Asnani, 1993).Overall, while a comprehensive review of the impact of each lake on regional climate is beyond the scope of this study, prior work indicates that even lakes that are smaller than the Great Lakes can be anticipated to have substantial effects on regional climate.
Prior studies have also illustrated that even relatively small errors in prescribed LSTs in a downscaling configuration can adversely affect simulated precipitation in regions surrounding lakes.The sensitivity study of Wright et al. (2013) showed significant changes in lake-effect snowfall over the Great Lakes in idealized simulations where LSTs were uniformly warmed by 3 • C. Anyah and Semazzi (2004) simulated changes in the spatial patterns and intensity of precipitation, as well as the amount of evaporation, over Lake Victoria in a modeling study where LSTs were uniformly changed by only 1.5 • C.
Interactions between the lakes and overlying air masses are also governed by the amount of lake ice in climates that permit lakes to freeze.Previous studies have found the presence of ice suppresses turbulent latent and sensible heat fluxes from the lake to the air mass (e.g., Zulauf and Krueger, 2003;Gerbush et al., 2008).As shown in the lake-effect snow case studies simulated by Wright et al. (2013), the presence of ice coverage over the lake's surface inhibits downstream precipitation.As a result, lake-effect snowfall decreases in some areas surrounding the Great Lakes during the later portion of the lake unstable season, as the water's surface freezes during the winter and early spring months.Overall, past studies indicate that if LSTs and ice are not properly prescribed, inaccurate values of precipitation and temperature in the lee of lakes result from a downscaled simulation.

WRF's alternative lake setting
Since the release of WRF version 3.3 in April 2011, an "alternative initialization of lake SSTs" option is provided in WPS to set LSTs (WRF User's Guide, 2014; Table 1).When employing this method, LSTs can be set using temporally averaged 2 m air temperatures from the driving data set, with the averaging period set by the user.Bullock et al. (2014), when downscaling a proxy GCM (R2) over a 12 km grid covering the Great Lakes, attempted to use the alternative lake setting to account for the greater thermal inertia of the Great Lakes by incorporating seasonal temperature changes after a 1-month time lag.Following the procedure of Bullock et al., if a user were to perform a simulation over the month of May, a single LST field would first be generated by temporally averaging air temperatures during the previous month of April; subsequently this static LST field would be used to set inland water temperatures throughout the month of May.Because Bullock et al. (2014) preprocessed the driving data in monthly segments, the LST field was prescribed to vary with time on a monthly basis.Using this method may imitate the seasonal changes observed over the Great Lakes, producing a lake stable and unstable season during the appropriate months.A drawback to this methodology is that the same lag time is used throughout the model grid, regardless of lake depth.Therefore, in this approach, large, deep lakes are implied to heat and cool on the same timescale as small, shallow lakes.Meanwhile, it is expected that observed seasonal temperature changes over smaller and shallower lakes would more closely follow atmospheric temperature changes than in large, deep lakes.If employed for simulations outside the Great Lakes, the procedure used by Bullock et al. (2014) should be modified to imitate the observed relationship between changing air temperatures and LSTs.
In its default configuration used prior to the release of version 3.5.1,WRF prescribes ice cover at grid cells where LST is less than 271 K (Table 1).This value is applied at all water points regardless of salinity.As winter 2 m air temperatures are frequently below freezing in the Great Lakes area, Bullock et al. (2014) found that unrealistically large spatial coverage of ice occurred when using the alternative lake setting in WRF version 3.4.1,with all five Great Lakes completely frozen for most of the winter.Such erroneous ice cover would be expected to negatively impact the simulation of precipitation, 2 m temperatures, and other variables influenced by sensible and latent heat fluxes supplied by the Great Lakes.Therefore, the use of the alternative lake setting in WRF may not be appropriate in some regions where sub-freezing air temperatures would result in unrealistic temporal and spatial coverage of sub-freezing LSTs and ice.
However, this is not a concern for tropical lakes where air temperatures would not be sufficiently low enough to result in frozen lakes.Argent (2014, Sect.3) demonstrated the utility of the alternative lake setting in WRF simulations over Lake Victoria in eastern Africa, finding that it improved the accuracy of simulated rainfall relative to the use of the default interpolation in which oceanic SSTs were used to set Lake Victoria's LSTs.

Climatological LSTs and ice
Another approach for setting LSTs and lake ice coverage when downscaling with WRF is to prescribe these variables from higher-resolution data sets of climatologically averaged quantities.This can be viewed as assuming stationarity for the lake state as is frequently done for other input variables in an RCM, such as land use and vegetation.Even for retrospective climate simulations, using this approach could be detrimental because the interannual variability of LSTs and ice -and its effects on the prediction of extreme eventswould not be captured using this method.When making future projections, it must be considered that prior studies have shown that LSTs cannot be assumed to be stationary in future warmer climates; in fact, some studies conclude that nonlinear feedbacks exist between regional climate change and LSTs and ice for some lakes.An observational study by Austin and Colman (2007) found that the multi-decadal warming trend in the Great Lakes region was amplified in the lake temperatures, relative to surrounding inland temperatures, because of the earlier breakup of ice and earlier springtime warming of surface water.In the downscaling simulations of Gula and Peltier (2012), increased snowfall was simulated in the lee of the Great Lakes in a warmer, mid-century climate because lake ice forms later in the winter.Gula and Peltier conclude that the impact of having the lakes remain free of ice is that increased latent and sensible heat fluxes are present for a longer time period during the lake unstable season, lessening the stability of the overlying air mass and enhancing precipitation.Magnuson et al. (2000) concluded that observed ice coverage is decreasing in lakes and rivers throughout the Northern Hemisphere.Such a decrease in ice coverage has been linked by observational studies to increases in lake-effect precipitation in the Great Lakes region (Assel and Robertson, 1995;Burnett et al., 2003;Kunkel et al., 2009).Because ice suppresses fluxes of latent and sensible heat (e.g., Zulauf and Krueger, 2003;Gerbush et al., 2008), decreasing ice cover in a warmer climate allows larger fluxes of latent and sensible heat to modify the overlying air mass, increasing downstream precipitation during the lake unstable season.None of the impacts on the lake state reviewed here (the warming of LSTs and more open water from which to produce fluxes) would be considered in the WRF model using LSTs and ice based on present-day climatology, and the effects of changing lake conditions on atmospheric stability, humidity, precipitation and convection would not be simulated.This approach could be improved by adding a linear increase to observed LSTs over time, which may be a valid approximation for the effect of climate change on some lakes.However, such an approach would not capture the nonlinear impacts of climate change (as described by Austin and Colman, 2007) on the Great Lakes.Overall, the efficacy of using a climatologically based approach is dependent on the amount of interannual variability, as well as the impacts of climate change on the lake state and whether those effects can be accounted for by the inclusion of a linear LST anomaly.

Land mask modification
To avoid the issues with LSTs discussed in Sect. 1 and illustrated in Fig. 2, Gao et al. (2012) modified the GCM land mask in the Great Lakes area so that skin temperatures from land points in the GCM were used to set LSTs on the WRF grid in their downscaled simulations.This treatment successfully eliminated the abrupt temperature discontinuities (such as those in Fig. 2) produced by interpolating a coarse data set.However, the effects of the lakes themselves are lost if GCM land temperatures are used to prescribe RCM water temperatures and the lake-land temperature contrasts, with their associated mesoscale phenomena such as lake breezes and lake-effect precipitation, are eliminated.Notaro et al. (2013) conducted an idealized modeling experiment in which the Great Lakes were replaced with forest and field land cover types.They found that the presence of the lakes affected precipitation, 2 m air temperatures and their variability, water vapor, cloud cover, incoming shortwave radiation, the hydrological budget, and the intensity of passing cyclones and anticyclones.The approach used by Gao et al. (2012), in which land surface temperatures from the GCM are used to specify water temperatures, partially accounts for some lake effects (such as changes in surface friction and albedo) because WRF would recognize the presence of a water surface.However, all processes related to the LST (e.g., ice formation, latent and sensible heat flux, 2 m temperature and moisture values, outgoing longwave radiation from the surface) would be negatively impacted by this treatment.Additionally, some impacts of climate change on the future lake state could be lost.For example, the amplification of Great Lakes LSTs, relative to over-land temperatures, observed by Austin and Colman (2007) will not be captured if land temperatures are used to set LSTs.

Use of simulated lake fields from GCM
A more sophisticated class of approaches for better representing the lake state in a downscaling configuration involves the use of a lake model.This can be done either by using outputs from the GCM's lake model (if available), driving a stand-alone lake model offline with GCM fields to simulate LSTs and ice, or by coupling a lake model to the RCM when downscaling.The CESM has a lake model embedded within its land surface model (LSM), version 4 of the Community Land Model (CLM4).CLM4 accounts for the presence of subgrid-scale lakes using the one-dimensional lake model described in Oleson et al. (2010).It is a column model partially based on the Hostetler lake model (e.g., Hostetler and Bartlein, 1990;Hostetler et al., 1993Hostetler et al., , 1994)), and it simulates 10 water layers through the depth of the lake, as well as additional layers for thermally active soil underneath and snow and ice above.However, when producing the downscaled simulation shown in Fig. 2, output from CLM's lake model was not easily accessible with other CESM outputs from the same simulation within archiving systems such as the Earth System Grid Federation.Lake temperatures and ice from CESM, and other GCMs with embedded lake models, could be leveraged by RCMs such as WRF to account for the impact of climate change on the lake state.In areas where lakes are at least partially resolved by the GCM, this approach would be effective at driving the RCM with simulated changes in LSTs and ice cover consistent with future projections and at keeping the RCM solution in the regions affected by lakes consistent with the GCM simulation.However, some small lakes may remain unrepresented by GCM data.

Use of a stand-alone lake model
If lake model outputs from the GCM are unavailable, one alternative is to use a stand-alone lake model driven by GCM fields to downscale the lake state in a manner which is consistent with the GCM's atmospheric fields.In the downscaling experiments performed by Gula and Peltier (2012) over the period 2050-2060, the Freshwater Lake (FLake) model was utilized to provide simulated LSTs and lake ice to WRF in the Great Lakes basin.GCM fields from the Community Climate System Model, with a spectral resolution of T85 (∼ 1.4 • grid spacing), were used to drive a FLake simulation on a 10 km regional grid, and the LSTs and ice cover simulated by FLake were subsequently used to drive the downscaled WRF simulation.In this one-way WRF-FLake model configuration, changes in LSTs and ice respond to changes in atmospheric variables in the driving GCM, but the lake model output is produced on the higher-resolution regional WRF grid.FLake is a 1-D column model which is highly reliant on empirical relationships and has been used in several studies with other RCMs (e.g., Mironov, 2008;Kourzeneva et al., 2008;Martynov et al., 2008;Mironov et al., 2010;Samuelsson et al., 2010).FLake requires a 2-D field of lake depths and the 1-D column model is called at each point.Therefore, the simulated LSTs are sensitive to lake depth, as well as the driving GCM fields.

Use of a coupled lake model within an RCM
In WRF version 3.6 a CLM-based lake model can be utilized with other non-CLM land surface models (WRF User's Guide, 2014; Table 1).This lake model is taken from CLM version 4.5 (Subin et al., 2012;Oleson et al., 2013) with some modifications by Gu et al. (2015) as discussed further below.Although a version of CLM4 was available as an LSM option within WRF version 3.5, the lake model in CLM4 was disabled in WRF (Table 1).In WRF version 3.6, CLM's Hostetler-based lake model can be applied by using horizontally varying lake depths (which are available in WPS version 3.6) or a uniform lake depth can be assigned to all lakes at runtime.Gu et al. (2015)  of this model configuration (WRF 3.2 and CLM 3.5) to simulate a 16-month period from 2001 to 2002 at 10 km grid spacing.It was shown that the lake model simulated LSTs well in Lake Erie but generated large biases in LSTs when compared to buoy observations in Lake Superior.However, the LST bias was reduced by reformulating the eddy diffusivity parameter in the CLM lake model, and it was concluded that the updated lake model within WRF-CLM was reasonably able to reproduce observed LSTs.However, no ice was observed during the period and the ability of WRF-CLM to accurately simulate ice cover was not examined in Gu et al. (2015).
In an alternative coupled approach, the prior work of Gula and Peltier (2012) has been updated with the option of using WRF-FLake as a two-way coupled model, where atmospheric variables simulated by WRF are used by FLake at each time step in the WRF model, and simulated LSTs and ice thicknesses are provided back to WRF by FLake.M14 concluded that the use of WRF-FLake resulted in a more accurate representation of LSTs and lake ice, relative to interpolation from the R2.Substantial improvements were shown in the simulation of the temporal and spatial variability of ice cover, and errors in LSTs were reduced by the use of the coupled model.Similar to Martynov et al. (2010), M14 found that FLake performed worst in the largest and deepest lake (Lake Superior) and best for the smallest and shallowest (Lake Erie).
When using an embedded lake model within an RCM, it can be anticipated that the period of time needed for spin-up could be larger than it is when all water conditions are simply prescribed.To spin up the WRF-FLake model in M14, the stand-alone version of the FLake model was driven with atmospheric conditions from the proxy GCM in a spin-up procedure recommended by Mironov et al. (2010) when using FLake.In this methodology, the initial year of the simulation is "looped" over 10 annual cycles with meteorological variables from the initial year repeatedly used to force the lake model, and the lake state at the end of each year used to initialize FLake for the start of the next year, ensuring that the simulated lake state converges to equilibrium with these atmospheric conditions by the end of the 10-cycle simulation.Output from the first year of this offline simulation is shown in Fig. 3 illustrating the adverse effects of using FLake output without adequate spin-up time.A time series taken from a representative point in Lake Superior shows unrealistically cool LSTs (below 200 K) occurring during the initial months of the simulation.Also during this period, unrealistically large ice coverage formed, freezing over all five Great Lakes.The observed ice cover plotted in Fig. 3 is much more limited in its spatial extent.Observed ice cover is plotted from National Ice Center (NIC) ice charts, which are processed and provided by the Great Lakes Environmental Research Laboratory (GLERL; Wang et al., 2012a).The FLake model results obtained after the spin-up period showed realistic values of LSTs and ice cover (M14).
To examine how WRF-CLM reacts during the initial months of a simulation, without any spin-up time, output from a 12 km WRF-CLM simulation (version 3.6) is shown in Fig. 4. In this simulation, the same methods as in M14 are followed but with the following changes: the model version is updated from 3.4.1 to 3.6, the CLM lake model is used in place of FLake, and no spin-up procedure is employed for initialization of the lake model (initial LSTs are interpolated from R2).As in M14, the Noah LSM (Chen and Dudhia, 2001) is used.Similar to the example shown in Fig. 3, significant overestimation of ice coverage occurs during the first year (Fig. 4).Although some adverse effects in this simulation are introduced due to the use of LSTs interpolated from the coarse R2 data to provide an initial state, the similarity of these results to FLake's fields in Fig. 3 suggests that the lack of spin-up time is a common problem to both model runs.It is also implied by the methodologies of other CLM-based studies, which do use spin-up or initialization procedures.Previous work by Subin et al. (2012) with the lake model in CLM4 used a 110-year period for the spin-up of their reference simulation.In their experiments with WRF-CLM, Gu et al. (2015) used an observed LST field for initialization.The nine sub-surface layers in their model were initialized based on the shape of an observed profile of lake temperatures, valid during that period of the year and taken from Lake Superior.Using this initialization methodology for a future downscaled simulation is not possible due to lack of observations, but simulated future lake profiles could possibly be utilized for initialization of downscaled runs.Overall, when using an embedded lake model in a downscaling application, users should consider how the lake model is being initialized or spun up in order to achieve results with accuracy similar to the prior studies discussed above.If the lake state is initially poorly prescribed from the GCM (with results similar to those shown in Fig. 2), a protracted spin-up could be required to reach equilibrium with the driving fields in the RCM and obtain more realistic results.
It has been noted previously that both WRF-FLake and WRF-CLM, as well as other 1-D lake models, tend to exhibit Table 2.A summary of the pros and cons of each method of treating lake surface temperatures and ice coverage described in the text.All approaches were found to eliminate unrealistic temperature discontinuities resulting from WRF's default interpolation methods as shown in Fig. 2.

Methodology
Positives Potential drawbacks WRF's alternative lake setting Effective at representing LSTs when lake temperatures are closely coupled with atmospheric temperatures.
Unrealistic ice formation possible when 2 m temperatures are below freezing.Cannot account for varying lake depths and differing timescales of warming and cooling throughout lakes.

Climatological
Observed LSTs and ice taken from high-resolution analyses.
For long-term simulations, user must include temperature trend or LSTs will not be in equilibrium with future climate state.Does not represent interannual variability of lake state.
Land mask modification Future LSTs can be taken from projected GCM temperatures.
Eliminates land-lake temperature contrasts.
Lake model component Models have ability to simulate future changes in LSTs and ice.
Additional preprocessing needed to provide lake model spin-up for RCM run or to use lake fields simulated by GCM.
difficulty in simulating deep lakes (e.g., Martynov et al., 2010;Stepanenko et al., 2010;Gu et al., 2015;M14).Some model error can be attributed to the fact that one-dimensional column models cannot represent 2-and 3-D processes (e.g., currents, drifting ice, and formation of a thermal bar).While more sophisticated lake models could be coupled with WRF, using computationally efficient 1-D models is advantageous in downscaling applications, where computational resources are taxed by the use of finer resolution.Additionally, Martynov et al. (2010) noted that more complex 3-D lake models are generally run with much finer grid spacing (∼ 2 km) than typical RCMs.Martynov et al. (2010) also compared the simulated water temperatures and ice coverage from the Hostetler and FLake models, finding that FLake generally performed better, but that the Hostetler model provides more opportunity to improve model performance because it utilizes more vertical layers and is less reliant on parameterization.A comparison of 1-D lake models by Thiery et al. (2014) showed favorable results for both FLake and Hostetler-based models (including the lake model found in CLM4) and noted their computational efficiency.When making regional climate projections with these models it should be noted that both WRF-FLake and WRF-CLM assume that lake depths are constant in time, which could be a poor assumption depending on the lake being modeled and the future period.Also, more complex lake models may be appropriate for higher-resolution (∼ 2 km grid spacing) RCM simulations focused on regions where lake dynamics are not adequately captured by the column lake models discussed here.

Conclusion
It has been shown in the present study and in previous work (e.g., Gao et al., 2012;Bullock et al., 2014; M14) that downscaling typically coarse GCM data, using WRF's default interpolation methods, to finer-resolution WRF grids results in LST discontinuities and spurious ice formation in the Great Lakes (Fig. 2).Although the default interpolation methods in WRF can easily be modified to alter the interpolation scheme or to eliminate the search option, none of these simple changes will overcome the challenges of setting the LSTs for inland water bodies that are not resolved by driving data when WRF is used as an RCM.Various alternate methods have been presented, and a summary of the positives and potential drawbacks to each approach is shown in Table 2. Using WRF's "alternative" lake setting instead of the default interpolation method in WPS eliminates unrealistically large and abrupt spatial discontinuities in temperature, but it causes large, deep lakes (such as Lake Superior) to erroneously freeze when ice is set based on an air-temperature threshold.All the other approaches discussed above can simulate more realistic ice cover than the default interpolation.However, the simulation of ice cover is obviously not a factor in downscaling studies where the environment does not become sufficiently cold to produce lake ice, such as those focusing on tropical regions.For example, the alternative lake setting has been used to improve rainfall results (relative to the use of WRF's default interpolation techniques) over Lake Victoria in eastern Africa by Argent (2014).Using climatological values in a future warmer climate will adversely affect results because LSTs cannot be assumed to be stationary over time.A warming trend could be applied to observed LST fields in order to improve this approach; however, a realistic trend may be complex to derive for some lakes as Austin and Colman (2007) have shown an observed nonlinear amplification of warming LSTs relative to inland temperatures in the Great Lakes region.The land mask alteration method of Gao et al. (2012) is effective at preventing discontinuities in surface temperatures, but the use of temperatures from land grid cells in the GCM to set LSTs in the RCM eliminates the presence of land-lake temperature contrasts which impact precipitation, winds (i.e., land-sea breeze), and other near-surface fields.The use of a lake model (either coupling a lake model to the RCM or using outputs from the GCM's lake model to drive the RCM) can improve the representation of the lakes in retrospective simulations and has the ability to simulate nonlinear impacts of climate change on LSTs and ice cover (e.g., Gula and Peltier, 2012, M14).For downscaling applications using WRF, we recommend setting LSTs and ice cover from either an RCM-or GCMdriven lake model, especially when simulating mid-latitude regions.In their studies focused on the Great Lakes, Notaro et al. (2013) and Wright et al. (2013) state that accurate predictions of changes in LSTs and ice cover from lake models are needed when simulating changes in regional climate.Zhao et al. (2012) also recommended the use of a lake model for simulating changes in regional precipitation in the Great Lakes basin.Including prognostic changes in the lake state is also possible if GCM data sets include predicted lake surface temperatures and ice within their publicly available outputs.For regional climate modeling efforts in which the RCM data is being archived for various end-user applications, we recommend the use of GCM-or RCM-driven lake modeling approaches.If such an approach is not used, the potential adverse effects of setting LSTs and ice cover using interpolation from the GCM should be documented, as is currently done in NARCCAP (2014).
The accuracy of the various approaches presented here is sensitive to the characteristics of the lakes to which they are being applied.Approaches which set LSTs as a function of over-land temperatures (such as the land mask modification approach or WRF's alternative lake setting) may perform adequately when applied to smaller, shallower lakes where LST changes are more closely coupled to air temperature changes.Investigators performing RCM experiments should consider both the present-day interactions between the lake and overlying air masses as well as the potential climate change impacts on the lakes within their model domain when choosing an approach.

Figure 3 .
Figure 3. Surface temperature from the initial year of a 10-year FLake spin-up simulation, taken from a point near the north shore of Lake Superior (48.47 • N, 87.54 • W) and shown hourly from 1 January to 31 December 2005 (top).LSTs at all lake cells are initialized with a default value of 274.15 K, and the time series shows either ice or water surface temperatures depending on whether ice is present.Simulated ice thickness (m) taken from day 30 of the same FLake simulation, valid 30 January 2005 (bottom left).Fractional ice values observed on this date plotted from the NIC ice analysis (bottom right).

Figure 4 .
Figure 4. Simulated ice cover (%) taken from a WRF simulation (valid 2 March 2006, after ∼ 4 months of simulation time) with the same model configuration as described in M14, but simulated with WRF version 3.6 and the use of the CLM lake model in place of FLake (left).A 2-D field of lake depths (instead of a single default value) was used from WPS to set the lake depth in this simulation.Ice coverage observed on this date is plotted from the NIC ice analysis (right).

Table 1 .
List of WRF versions discussed in the text, ordered chronologically by the date of release and with relevant model updates summarized.