Urban nocturnal cooling mediated by bluespace

The spatiotemporal characteristics of air temperature and humidity mediated by urban bluespace are investigated using a combination of dense network of climatological observations in a medium-sized US city, computational fluid dynamics, and analytical modeling approaches. Both numerical simulation and observational results show that the rate of change of hourly averaged air temperature and humidity at 3.5 m over urban areas peaks 2 h after sunset, while it decreases with time monotonically over greenspace, indicating different impacts due to presence of urban lakes. The apparent temperature decreases with distance to lakes in urban area due to higher near-shore humidity. This highlights that urban lakes located near city center can deteriorate the nighttime cooling effects due to elevated humidity. Finally, two analytical models are presented to explain the connection between the surface and air temperature as well as the spatial variation of air temperature and humidity adjacent to the urban lakes. These simplified models with parameters being inferred from the network of measurements have reasonably good performance compared to the observations. Compared to other sophisticated numerical simulations, these analytical models offer an alternative means that is easily accessible for evaluating the efficacy of bluespace on urban nocturnal cooling.


L v
Latent heat of vaporization of water Q Available energy Q l Available energy of lake surface Q u Available energy of urban surface q Mass fraction of water vapor q l Humidity above lakes q ls Humidity of lake surface q ref Reference humidity q us Humidity of urban surface R n Net all wave radiation RH Relative humidity T Ambient temperature T i Initial surface temperature T l Air temperature above lakes T ls Temperature of lake surface T ref Reference temperature T s Surface temperature T s_average Nighttime averaged surface temperature T us Temperature of urban surface T 3.5 Air temperature at 3.5 m above surface t Time after sunset t sunrise∕sunset Sunrise/sunset time u Velocity u x Streamwise velocity u y Vertical velocity u 10 Velocity at 10 m above surface

Introduction
Nearly 70% of the world population will be living in urban areas by 2050 (United Nations 2018), imposing multifaceted challenges to urban sustainability (UN DESA 2018). The urban heat island (UHI) effect since its first discovery in the nineteenth century has been an increasingly significant concern for the urban thermal environment under climate change (Cheval and Dumitrescu 2009;Lee and Baik 2010;Cui and Shi 2012;Howard 2012). There have been increasing efforts in cities worldwide in adopting the "nature-based solutions" for UHI mitigation (Kleerekoper et al. 2012;Santamouris et al. 2017;Yu et al. 2020), such as making use of the greenspace (i.e., urban vegetation and forestry) and bluespace (natural or man-made water bodies). Extensive studies by means of field observation, numerical simulation, and theoretical scaling have been conducted to characterize the mechanisms and efficacy of urban cooling by greenspace and bluespace (Robitu et al. 2004;Chatzidimitriou et al. 2013;Theeuwes et al. 2013;Zoras et al. 2015;Gu et al. 2016;Gunawardena et al. 2017;Fan et al. 2019;Hu and Li 2020). It is found that the size, shape, and configuration of the bluespace and greenspace as well as the mean wind speed in the surrounding area can largely influence the cooling effect. Despite the consensus that evapotranspiration from the urban greenspace is one of the dominant pathways of reducing UHI (Taha 1997;Zhang et al. 2010;Kleerekoper et al. 2012;Santamouris et al. 2017;Yu et al. 2020), research on the cooling effect of urban water bodies shows a diverse range of conclusions (Völker et al. 2013;Hu and Li 2020). Compared to the vast literature on urban greenspace, more research and observational evidence on the impact of bluespace is needed. Due to the large heat capacity of water, urban water bodies may even impede cooling in urban areas especially during nighttime (Theeuwes et al. 2013;Hu and Li 2020). High humidity around the urban bluespace is found to decrease evaporation (Theeuwes et al. 2013), hence reducing its effectiveness in UHI mitigation and negatively impacting the human comfort level due to synergistic effects of heat and humidity (Hu and Li 2020).
Despite that complex factors simultaneously control the cooling effects of urban bluespace, they modulate the nocturnal urban thermal environment via two key mechanisms. Firstly, over distinct types of the urban land and open water surfaces, the surface temperature T s decreases from an initial value T i with a square root dependence of time (e.g., T s − T i ∼ √ t ) at night when net longwave radiation is in balance with the ground heat flux (Wang and Li 2017). However, as surface temperature decreases in time, the surface humidity (assuming no surface ice) can simultaneously vary significantly over the urban bluespace, while remaining fully saturated according to the Clausius-Clapeyron relation. Thus, the synergistic effect of heat and humidity of urban bluespace on the thermal environment, in particular, their co-evolution with time, needs to be better quantified.
Secondly, if considerable spatial gradients of air temperature exist over different types of surfaces, air circulations will be generated, such as the well-known rural-urban circulation. Within cities, micro-scale circulations (Richiardone and Brusasca 1989;Fan et al. 2016;Wang and Li 2016;Omidvar et al. 2020) between built-up surfaces, greenspace, and bluespace can also contribute to the urban surface cooling. To quantify the contribution of such convective cooling on the urban thermal environment at night, previous studies have proposed general analytical models for city-scale circulations. For example, Wang and Li (2017) scaled the UHI intensity and uncovered that it is related to horizontal flow velocity as well as vertical temperature gradient (Wang and Li 2017). Previous research (Lee et al. 2012;Wang and Li 2017) also applied a theoretical scaling analysis to the air temperature and found that the nocturnal UHI decreases exponentially, in reasonably good agreement with observations. Compared to more sophisticated and costly mesoscale weather models, the analytical modeling approach is found to be reasonably accurate, thus being valuable as an accessible tool for quick evaluation and assessment of urban climatology. However, extending analytical approach to characterize the impact of urban bluespace on the intra-urban variability of temperature and humidity remains relatively underexplored.
Given the increasing availability of urban meteorological measurements (Muller et al. 2013;Warren et al. 2016;Konstantinov et al. 2018), these data not only will be useful in validating numerical models in urban climate research but may also provide an alternative means to constrain the parameters in physically based analytical models. In addition, numerical simulations can compensate the deficiency of field observations, where physical quantities are usually measured at a limited number of locations. Thus, in this study, we analyze observational data from a spatially dense meteorological network from the US Long-Term Ecological Research Network (LTER) of North Temperate Lakes over the metropolitan area of Madison, Wisconsin (Schatz et al. 2019;Hu and Li 2020). Combining the field observations and numerical simulations using computational fluid dynamics and analytical approaches, we aim to address the following three main research questions: 1. How do the local meteorological quantities (wind, temperature, and humidity) change with time, given distinct temporal evolution of surface temperatures over different land surface types? 2. How does the relative location of bluespace with respect to the urban land surface affect the spatial characteristics of the nocturnal thermal environment? 3. How can we formulate analytical models to assess the spatial-temporal characteristics of air temperature and humidity over different land cover types?
The paper is organized as follows: Observational data are first analyzed (Section 2.1 and Section 3.1). Numerical simulation using ANSYS Fluent 19.4.0 commercial software is conducted to reveal the structures of micro-scale air circulations (Section 2.2 and Section 3.2). Finally, two analytical models, which explain the connection between the time-varying surface and air temperature, and the development of urban thermal boundary layer, respectively, are proposed to represent the spatiotemporal characteristics of the urban climatological thermal environment (Section 4). Concluding remarks are presented in Section 5.

Observational dataset
For the US LTER of North Temperature Lakes observations (https:// lter. limno logy. wisc. edu/), 152 HOBO U23 Pro V2 sensors are attached on streetlight poles at a height of about 3.5 m over various landscapes. The network provides the half-hourly record of the ambient temperature ( T( • C) ) and relative humidity ( RH(%) ) from March 2012 to December 2017 (Schatz et al. 2019;Hu and Li 2020). After eliminating the sites with incomplete data availability, the data collected from 137 sites is used in the current research. Only data during clear sky and calm weather conditions (i.e., wind speed at 10 m being less than 3.13 m/s) from Jun 1st to Aug. 31st are chosen for analysis to focus on the convective effect arising from different types of surfaces and minimize any mean advective effect. In addition, understanding the potential cooling effect of urban bluespace is also the most pertinent under such meteorological condition. Figure 1a illustrates the satellite image of the observational domain. The center of the urban area in Madison locates between two lakes and both the urban area and lakes are surrounded by the rural area. Figure 1b shows some of the sites that are in the urban area between two lakes, which are selected for more detailed analysis in Section 3.1.  ANSYS Fluent 19.4.0 commercial software is used in this research. A semi-idealized setting, analogous to the spatial configuration of land use types in Madison, is set up. The dimensions of the computational domain are set to be 14 km in streamwise (x) and spanwise (z) direction and 1 km in vertical (y) direction, respectively, as illustrated in Fig. 2a.

Setup of computational fluid dynamics modeling using Fluent
The bottom of the domain is divided into three different types of surfaces, namely, the impervious area, lake, and greenspace; each is characterized by distinct surface properties (more details will be discussed). The dimensions of these areas are shown in Fig. 2b. The impervious area represents the built-up urban land area in the observation. The distance between the two lakes is 1 km. The mesh contains about 2,000,000 cells. A total of 140 grids with each size of 100 m are imposed in horizontal direction and 100 unevenly distributed grids with the bias factor of 50 being imposed in vertical direction. The sizes of these grids are between 1 m near the bottom and 40 m near the top of the domain. The Renormalization Group (RNG)-k-ε model (Ansys 2019) with standard wall function is used as the turbulent model. The governing equations for air and water vapor (both as ideal gases) are solved in the simulation, using a secondorder upwind scheme. For pressure-velocity coupling, the standard SIMPLE algorithm is adopted (Ansys 2019). Boundary conditions for the domain are set as follows. The top boundary is set to be symmetric implying that the flows are parallel to that face. Lateral boundary conditions are periodic, which indicates that the influence of the lakes and impervious areas to those boundaries is small. For momentum boundary condition, no slip boundary condition is imposed to all three areas in the bottom of the domain. The roughness height K S and roughness constant C S of the three types of surfaces are displayed in Table 1, where K S is the roughness height representing the average height of surface roughness elements. C S indicates the degree of randomness of the distribution of surface roughness elements and it is in the range of [0.5, 1] . Here, the value of K S and C S for water surface follows that of Tominaga et al. (2015). For impervious area and greenspace, the surface is considered uniform; thus, C S = 0.5 is adopted (Ansys 2019). Because grass and crops are the dominant vegetation of greenspace, K S of greenspace is set to 0.1 m which is close to their estimated heights. K S of impervious area is set to be 0.2 m, which is larger than that of greenspace but still much smaller than half of the vertical grid size such that the urban area is treated as a rough surface (Ansys 2019).
For the temperature boundary conditions, the timedependent surface temperature T s of the three types of urban surfaces can be obtained from surface energy balance and Green's function approach (Wang and Li 2017). Assuming that the net long wave radiation balances the conductive ground heat flux, the time evolution of nighttime T s can be solved as (Wang and Li 2017): The air temperature and humidity in Section 3.2 are averaged in the spanwise direction within the two dash lines  where T i is the initial surface temperature, L net is the net longwave radiation, and and are molecular thermal conductivity and diffusivity, respectively. The values of these parameters are shown in Table 2.
The choice of T i and L net in Table 2 is informed by the following reasoning based on previous field observations. First, L net roughly remains constant during nighttime (Offerle et al. 2003;Wang et al. 2013;Li et al. 2018). Second, the total change of air temperature over the night can be approximated by the change of surface temperature because a negative feedback between the surface and air temperature exists in regulating both temperatures (more details will be given in Section 4.1 in the proposed analytical model). Thus, from the observation data of air temperature drop for different dominant land surface types near the measuring stations (i.e., lake, urban area, and greenspace) (Hu and Li 2020), the surface temperature drop in these three areas can be estimated. Then L net can be calculated from Eq. (1). In addition, the absolute value of L net in rural area can be 1 − 9W∕m 2 lower than that in urban area during summer (Oke and Fuggle 1972), which is consistent with the lower magnitude of L net for greenspace estimated from this method. To impose the initial surface temperature T i for these three types of areas, first from Eq. (1), the nighttime averaged surface temperature of three areas T s_average is derived as a function of T i : where t sunrise(sunset) are the local time (LT) of sunrise (sunset). To estimate T s_average , the remotely sensed surface temperature obtained from the ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS) is used as a surrogate. Surface temperatures at three locations in urban, greenspace, and lake are averaged over summer nights as T s_average in three areas, respectively. Albeit a rather crude estimation of T s_average , as being limited by the lack of continuous temporal observations of the surface temperature from satellite data, this method gives a firstorder estimate and is straightforward to implement. With L net computed from the first step and the known local time of sunrise and sunset, T i as the only unknown can be obtained by solving Eq. (2). Figure 3 shows the temporal change of T s for the three types of surfaces computed using Eq. (1) for the respective values of L net and T i . Boundary conditions for mass fraction of water vapor, q , are set as follows. Very close to the lake surface, the air is assumed to be fully saturated (Li et al. 2018). Thus, humidity on the lake surface varies with time as a result of the time-varying lake surface temperature, according to the Clausius-Clapeyron relation. Observation from previous research suggests that the latent heat flux in urban and rural areas during summer night is very close to zero (Yap and Oke 1974;Cleugh and Oke 1986). Thus, zero diffusive mass flux is imposed on computational domain corresponding to impervious and greenspace, implying that there is zero evaporation and no water vapor can go through these boundaries.
The background wind across the domain is set to be 0 m/s because we are interested in the clear and calm weather conditions, where the mean flow is dominated by circulations driven by contrasts in surface temperature. To spin up the simulation, the model is first run under steady state for 240,000 iterations, until the velocity and temperature field do not change by more than 10% every 40,000 iterations. In this steady state, the surface temperatures of the three surface types are set as constant, at 1.5 • C higher than T i of those areas, respectively. Such a setting is necessary to achieve the desired values of air temperature at the end of the spin up period and it will not impact the subsequent results for further analysis. The initial value of q is set to be 0.0109 for the entire domain, which is estimated from the observed humidity in Madison averaged over the 137 observational sites (c.f. Sec. 2.1). Note that the observations are only for near-surface humidity at 3.5 m and extrapolating to higher altitudes may incur inaccuracy. However, the nearsurface humidity is assumed to significantly impacted by the surface source from the lakes, thus justifying the setting of initial q . After the spin-up period, velocity and temperature fields in the entire domain are fully developed. The timedependent boundary conditions for T s discussed above are then imposed. The time step of the unsteady simulation is set to be Δt = 4s . The simulation is then run for 18,000 s.

Validation of Fluent simulations
A simulation of air flow over open water is performed and compared with results in Tominaga et al. (Tominaga et al. 2015) to validate that Fluent is able to reproduce the development of boundary layer in a wind tunnel experiment. For the validation simulation, the dimensions of the computational domain are set to 3 m in streamwise direction, 1 m in spanwise direction, and 0.98 m in vertical direction, respectively. The temperature T and humidity q of the inflow are 20 ℃ and 0.008, respectively, while that of the water surface are 16 ℃ and 0.01143, which is the saturated humidity at T = 16 • C . Therefore, a boundary layer will develop above the water surface, where T and q vary with x and y. Other setup of the simulation follows that of Tominaga et al. (2015). Figure 4 shows the comparison of normalized temperature and humidity distribution along streamwise direction at 3 different heights in our results and those presented by Tominaga et al. (2015). The reference temperature and mass fraction of water vapor are T ref = 20 • C and q ref = 0.008 , respectively. The present results agree well with those in Tominaga et al. (2015). The validation of grid convergence is presented in the supplementary material.

Observational analysis
The air temperature at 3.5 m, T 3.5 , changing with time and location in urban area between the two lakes is shown in Fig. 5a and b, respectively. T 3.5 decreases almost linearly with time as shown in Fig. 5a. Figure 5b illustrates that T 3.5 at locations far from the lakes (distance to west lake = 500 m, for example) is higher than locations close to the lakes (distance to west lake = 0 or 1000 m, for example). For the same measurement locations, the temporal and spatial variation of humidity at 3.5 m, q 3.5 , is illustrated in Fig. 5c and d. The value of q 3.5 is higher near the lakes and decreases with distance to the lake shore. In addition, for all locations between the two lakes, q 3.5 increases dramatically in the first 2 h, and after 20:00 LT, q 3.5 at locations far from the lakes still increases with time while for near-shore sites, a slight decrease with time is observed.
The 137 sites around Madison can be divided into three categories, namely, the "lake," "urban," and "greenspace" surface types, respectively, based on the criteria in Hu and : the "lake" sites have zero distance to open water surface; land surfaces with green area fraction less than 60% are categorized as "urban"; and those of greenspace have green area fractions greater than 60%. Figure 6 displays T 3.5 of these three types of surfaces changing with time after sunset. Value of T 3.5 averaged over all "lake" sites decreases the slowest while that of the greenspace surface type decreases the fastest, which is likely due to different cooling rates of the surfaces, resulting from differences in thermal properties (heat capacity, conductivity, and emissivity). The synergistic effect of temperature, humidity, and wind can be represented by apparent temperature (AT), which is given by: where T is in ℃, e in hPa is water vapor pressure, and u 10 is the mean wind speed at the reference level of 10 m (Hu and Li 2020). The observational results of hourly apparent temperature at 3.5 m, AT 3.5 , over urban, greenspace, and lake show similar trend as those of T 3.5 in Fig. 6, respectively (results not shown, u 10 is represented by the (3) AT = T + 0.33e − 0.7u 10 − 4 wind speed obtained from the local airport, rather than the measuring sites).
To gain a qualitative understanding of the effects of both surface types and distance to water to the cooling rate of air temperature, we first consider the average cooling rate. We binned the observational data according to ranges of distance from the lake. The cooling rate at 3.5 m above the surface at 20:00 LT shows distinct trends with respect to distance to water (see Fig. 7). Since there are only two measuring sites over lake, their results are not shown in Fig. 7. Surfaces over greenspace cool down faster whereas lower cooling rate is generally observed for sits over urban area close to the water (i.e., 0-500 m). The cooling rates at 20:00 LT averaged over Fig. 5 Observational results of (a) air temperature at 3.5 m T 3.5 changing with local time and (b) T 3.5 distribution in the urban area between two lakes (sites shown in Fig. 1b). The solid and dash lines in (a) and (b) show locations closer to west and east lake, respectively. (c) and (d) are the same as (a) and (b), respectively, except for humidity Fig. 6 Observation of the hourly T 3.5 averaged over sites of "lake," "urban," and "greenspace" surface types after sunset for the chosen days in summer months Fig. 7 Observational results of 3.5 m air cooling rate at 20:00 LT changes with distance to water. The line in the box is the median, the top and bottom of the box are 75th percentile and 25th percentile, respectively, and the black lines out of the box are maximum and minimum surface types (i.e., greenspace, urban, and lake) are 1.82 , 1.46 , and 1.19K∕h , respectively, implying that the greenspace not only has the lowest air temperature, but the temperature also decreases the most rapidly.
The dense network of observational sites records only air temperature and humidity, but no anemometers are installed for measurement of wind speed and direction. Numerical simulation with computational fluid dynamics modeling is able to overcome this observational limitation and reveal how microscale circulations develop and possibly impact the observed temperature and humidity patterns as being analyzed in Section 3.2.

Spanwise and time averaged numerical analysis
Due to symmetry of the setup of different surface types, we consider the cross-section of a vertical plane along the streamwise direction at z = 0m (x-y plane) (see dotted line in Fig. 2b), where both early and late periods after the spin-up time are shown in Fig. 8. During the first 1000 s, the average horizontal flows are from greenspace to impervious area near the surface and from impervious area to greenspace near the top (Fig. 8a) while a strong updraft is found in the center of the domain (Fig. 8b). Hence, there are two large circulations in the domain. These circulations result from the temperature difference between different surfaces: air temperature near the impervious surface is higher, causing updrafts, then the air near the greenspace will move towards impervious area to fill the low pressure area, thus resulting in the large circulations in the whole domain, which are ubiquitous and were found in many previous studies (Richiardone and Brusasca 1989;Fan et al. 2016;Wang and Li 2016;Omidvar et al. 2020). However, at about 8000 s, the large circulations gradually weaken and finally disappear (results not shown). After that, two other circulations, which are in the reverse direction of the previous ones during earlier time, form and occupy most of the domain, as can be seen from Fig. 8c and d. The initially strong updraft developed near the impervious area is suppressed as air temperature decreases (Fig. 8d) and magnitudes of both u x and u y are lower than earlier period. Emergence of the weakened and reversed circulations may be related to the transition of atmospheric boundary layer from unstable to stable as the surface temperature gradually decreases below the air temperature (i.e., the surfaces cool down faster than the air temperature far above the surface). Even though flow close to the surface (i.e., y ≤ 20m ) is still from the greenspace to the impervious area due to its higher surface temperature, but the vertical extent of this microscale circulation is limited to the near surface.
Note that semi-idealized setup was implemented in the modeling, where urban building morphology, details of the urban greenspace, and effects of varying meteorological conditions are neglected. Thus, simulation results are useful to understand the generic wind circulation phenomenon due to contrasting surface temperatures between different surface types, but they should not be interpreted as simulation corresponding to any particular day in the chosen study area. Next, we examine the simulated air temperature and humidity. The magnitude of T 3.5 over impervious area increases with distance to water while that of greenspace decreases with distance to water, as can be seen from Fig. 9a and c, which are both consistent with the observation. Also, the influence distance of the lake on the greenspace is about 500 m, agreeing well with the finding in Hu and Li (2020) using the same observational data set. However, the temperature variation between different places in impervious area is smaller than that of the observation, which is likely caused by an oversimplification of the surface morphology and actual landscape types in the current semi-idealized simulation. Figure 9e shows that the greenspace influences the lake up to a distance about 1000 m. The magnitude of q 3.5 over greenspace is close to the background magnitude at 0.0109 (Fig. 9d), due to the existing pattern of circulation that features a consistent near-surface wind direction from greenspace to the lake. Thus, the advective transport of q from the source (i.e., the lake) upwind is ineffective and there is no additional evaporation by construction of our simulation setup. Nevertheless, if nighttime evaporation from local sources of water vapor in the greenspace is non-zero, it is likely to dominate over the observed humidity and thus impact of the lakes is negligible with the given circulation pattern. For q 3.5 over majority of the lake surfaces (500 m < x < 2750 m) (Fig. 9f), at first, it increases with time and reaches the maximum value at 2-3 h, which is close to its saturated value given the local air temperature T 3.5 . After that, q 3.5 starts to decrease with time as a result of the decreasing lake surface temperature, dictated by the Clausius-Clapeyron relation. This also explains the trend of q 3.5 over the impervious surface far away from the lake: for instance, for distance to water being 500 m, the air remains above saturation and water vapor from the lake is being Fig. 9 Left column: spanwise and hourly averaged T 3.5 for (a) impervious, (c) greenspace, and (e) lake surfaces simulated by Fluent. Right column: same as left column, respectively, except for humidity. City center in (e) and (f) is at z = 0m plane transported by advection, such that q 3.5 always increases with time. However, impervious region close to the lake (distance to water < 100 m, for example) is strongly influenced by humidity variation over the lake surface. Therefore, the trend is very similar to that of the lake:q 3.5 first increases with time and reaches the maximum around 2-3 h before starting to decrease. This is also consistent with our findings in observational results (see Fig. 5c and d).
In summary, both observational and simulation results support that T 3.5 for the impervious area (i.e., urban area) increases with distance to water. The results of simulation further indicate that T 3.5 of greenspace decreases with distance to water. For q 3.5 , both observational and simulation results indicate that humidity at locations with high humidity (i.e., above or close to the lakes) first increases with time until nearly saturated and then decreases with time but remains nearly saturated, while for locations with low humidity (i.e., far from the lakes), q 3.5 always increases with time.
The rate of change of T 3.5 computed from hourly averaged values in both simulation and observation are compared and shown in Fig. 10. Simulation results capture the magnitude and trend of the observation reasonably well. The non-monotonic trend in cooling rate over impervious surface and its maximum at 20:00 LT are observed in both simulation and observational results (Fig. 10a). Such an increasing cooling rate during the first 2-h magnitude of T 3.5 can be explained by considering the temporal change of the difference between T s and T 3.5 . During the first hour of the simulation, T s , which is initially higher than T 3.5 in impervious area, gradually decreases and becomes lower than the corresponding T 3.5 due to the large surface cooling rate. However, during this period, T 3.5 − T s has undergone a sign change and thus only a small positive difference between T 3.5 and T s exists, giving rise to a slow cooling rate of the air. As surface cooling proceeds to the second hour, the larger magnitude in T 3.5 − T s leads to an increased cooling rate that subsequently reduces their difference and thus slower cooling rates in later hours. These observed temporal change in cooling rate of air, dT 3.5 ∕dt , hints the existence of a negative feedback between dT 3.5 ∕dt and the magnitude of T 3.5 − T s , which will be explored further in Section 4.1. In Fig. 10b and c, both simulation and observation indicate that the cooling rates of the greenspace and lake decrease with time. The relative error of cooling rate between simulation and observation of impervious area is about 10%; while that of greenspace is about 20%, which may be caused by lack of considering evaporation effect in the simulation. The relative error of lake reaches 25%, because the measuring locations in the observation are near the lake shore, which are also influenced by the land properties, therefore, it may not be fully representative of the temporal change of T 3.5 above the lakes.
The spanwise and hourly averaged AT 3.5 simulated by Fluent in three areas are shown in Fig. 11. Because the simulated u 10 is about 0.2 m/s, the impact of u 10 is no more than 0.2 ℃; thus, it is negligible compared with that of temperature and humidity. Although T 3.5 of impervious area slightly increases with distance to water (Fig. 9a), AT 3.5 close to lakes can be 0.8 ℃ higher than that far from the lake because Fig. 10 Air cooling rate at 3.5 m above the surfaces for both simulation and observation in (a) impervious, (b) greenspace, and (c) lake of high humidity near the lakes (Fig. 11a). The temporal change of the cooling rate of AT 3.5 in impervious area has a similar trend to that of T 3.5 , which first increases with time during the first 2 h and then decreases. For AT 3.5 of greenspace as shown in Fig. 11b, due to 3.5-m air humidity being close to the background value of 0.0109 throughout the simulation time, the spatial-temporal behavior of AT 3.5 is very close to that of T 3.5 (Fig. 9c). AT 3.5 over majority of the lake surfaces (500 m < x < 2750 m) shown in Fig. 11c is constant with x. AT 3.5 between 2750 m < x < 3500 m decreases with x because of the same trend of both T 3.5 and q 3.5 ( Fig. 9e and f). It should be noted that because of the initial condition for q 3.5 in all areas is set to be the background humidity 0.0109, which is likely underestimated over the lakes, thus, AT 3.5 during the first hour may be underestimated, leading to an underestimated cooling rate of AT 3.5 in the first hour. From the second to the fifth hour, the cooling rate of AT 3.5 monotonically decreases with time.

Analytical models
In this section, analytically tractable models are developed for two different purposes. Model A (Section 4.1) aims to quantitatively connect the temporal evolution of the surface temperature to that of air temperature. Model B (Section 4.2) is developed to understand how urban open water surfaces impact the nearshore temperature and humidity.

Model A: Connection between seasonally averaged time-varying surface and air temperature
Under calm weather conditions at night, where the surface energy budget over a generic urban land surface is assumed to be balanced by the net longwave radiation and ground heat flux, the linkage between cooling rate of near surface air (e.g., at 3.5 m) and surface cooling rate is mediated by the difference between the air and surface temperatures. It is also assumed that advection only affects limited areas near the boundaries of the surface. The following relation can be proposed for the cooling rate of air: where the surface temperature T s = T i − a √ t , k < 0 is a constant that indicates air cooling rate, which should depend on the height at which T is measured (assuming near surface, similar to the observational setup at 3.5 m). The value of a is 2L net √ ( ∕ ) (see Eq. (1)), representing the linkages of T(t) with surface cooling rate, which is controlled by the surface material properties (κ and λ) and the net longwave radiation. Equation (4) represents a negative feedback between surface and air temperatures: a decrease in T s leads to increased dT dt that intends to restore the difference between T s and T . It is important to emphasize that the proposed relation is by no means a predictive one . 11 Simulation results of spatial and temporal distribution of AT 3.5 in (a) impervious area, (b) greenspace, and (c) lake for any particular summer night. Instead, it is meant to represent the climatological nocturnal cooling for a given type of land surface. Throughout the nighttime, we further assume that value of k remains constant with time; thus, solution to Eq. (4) is where c is an integration constant, which represents the initial temperature difference between the air and the surface.
With the observational data from the 137 sites, here we find the value of k by the least square method, which minimizes the temperature difference between the modeled and observed values. In general, we expect k to also vary spatially but its value is assumed to be only dependent on the categorized surface type, namely "lake," "greenspace," and "urban" according to the criterion explained in Section 3.1. Figure 12 compares T 3.5 − T i of three surface types between the prediction of analytical model A, observation, and simulation, showing reasonably good agreements. The values of k for urban, greenspace, and lake surface types are −0.000493 , −0.000100 , and −0.000224s −1 , respectively.

Model B: Development of urban thermal boundary layer
Evidence from observational and simulation results shown in previous sections suggests that the air temperature and humidity vary with distance from the lake shore. Thus, we adopt a model originally developed for the simultaneous heat and water vapor transfer over natural land surfaces (Yeh and Brutsaert 1971) for urban bluespace here. Figure 13 is a schematic showing the key physical variables solved by this analytical model. The lake and urban areas have different surface temperature and humidity T ls , q ls and T us , q us , respectively, in which the surface humidity of lake q us is saturated. The mean wind u is in the direction from the lake to the urban area. An urban thermal boundary layer develops in response to the mean wind, in which the air temperature and humidity T and q will be a function of both x and y. T and q are given by the following equations (Li and Bou-Zeid 2013): where u and l are the ratios of actual surface specific humidity q us and q ls divided by the saturated surface specific humidity for the urban and lake surfaces (because lake surface is saturated, l = 1 ), respectively. T * ls = q ls ∕( , where c p is the specific heat capacity of air at constant pressure ( J∕(kg × K) ), L v is the latent heat of vaporization of water ( J∕kg ), and is the derivative of the saturated specific humidity with respect to temperature at T = T ls . The quantity g u 10 is a positive function that depends on, and decreases with u 10 ; Q u − Q l is the difference in available energy ( Q = R n − G = H + LE ) between urban and lake surfaces, where R n is net all wave radiation; G , H , and LE are ground, sensible, and latent heat flux, respectively; and f 1 (x, y) and f 2 (x, y) are functions that are only dependent on locations x and y. Derivation of the model can be found in Yeh and Brutsaert (1971) and the appendix of Li and Bou-Zeid (2013). For the nocturnal cooling considered here, Q u and Q l are zero, since the net long wave radiation is assumed to be balanced by the ground heat flux. Thus, the difference between temperature and humidity is controlled by the first term in Eqs. (6) and (7), which leads to Since model B does not consider the temporal effect, it is applied to the nighttime air temperature and humidity (6) T(x, y) − T l (y) = 1 − u l T * ls f 1 (x, y) + g u 10 Q u − Q l f 2 (x, y) averaged during 18:00LT to 23:00LT in locations close to the transition areas between different surface types because the advective effect dominates in these locations due to the large surface temperature difference. In fact, the choice of u , the ratio of actual surface specific humidity to the saturated value, can considerably impact the modeled results. The value of u may be a function of position and time but it is assumed to be a constant value (Li and Bou-Zeid 2013), which is a reasonable assumption for the lack of information on the surface humidity. Based on comparison with the observation (measurement data along sites shown in Fig. 1b), a value of u = 0.65 is adopted here. Average air temperature difference T 3.5 − T l and humidity difference q 3.5 − q l are shown in Fig. 14, along with results obtained from varying u by ± 0.1. With known T l and q ls , the model is able to predict the magnitude and variation in x for T 3.5 and q 3.5 . Some discrepancies between the model and observation exist, which are likely due to the assumed constant value of u . It could also be related to that the assumed steady-state condition in model B is not fully met during the chosen calm and clear conditions in this study. In general, model B expected to work better under steady condition (i.e., dT 3.5 ∕dt ≈ 0, dq 3.5 ∕dt ≈ 0 ), where the spatial variation of T 3.5 and q 3.5 is driven by the advective transport near lake shore region. Therefore, under the chosen meteorological condition, the local surface cooling dominates the temporal change of air temperature. Thus, compared to model B, model A better captures the underlying local negative feedback between surface and air temperatures. For humidity, weak and intermittent mean wind due to flow circulations (c.f. Figure 8) results in a lack of steady advective transport. For example, the mean streamwise velocity at 10-m height from numerical simulations is about 0.3 m/s (results not shown here). Overall, model B as a simple analytical approach is still useful when considering temporally averaged T 3.5 − T l and q 3.5 − q l and model B is able to capture the order of magnitude and spatiotemporal trends of T 3.5 − T l and q 3.5 − q l with reasonable accuracy. Compared to much more sophisticated models, e.g., the computational fluid Fig. 14 Comparison of (a) T 3.5 − T l and (b) q 3.5 − q l in analytical model B with observation averaged from 18:00LT to 23:00LT. The shaded areas are bounded by setting u between 0.55 and 0.75 dynamics model in Section 3, such analytical model can be easily implemented to estimate the influence of urban lakes to the nearby urban landscape since the lake surface temperature and u are the only two required parameters. It is therefore useful as a first-order approximation for the influence of urban lakes as a function of distance from the shore.

Conclusion
In this study, the effects of urban bluespace on the nighttime air temperature and humidity are investigated by means of a combination of observational, numerical, and analytical method in Madison, WI, under typical calm weather condition in summer. First, the observed near-surface (at 3.5 m) air temperature T 3.5 and humidity q 3.5 are analyzed. Next, semi-idealized computational fluid dynamics modeling using Fluent is conducted and the numerical results are analyzed and compared with the observation. The apparent temperature AT 3.5 , a quantity determined by the synergistic effects of wind, air temperature, and humidity, is examined. Lastly, two simple analytical models are proposed to quantify the temporal variation of air temperature over different surface types and the effect of urban lake on temperature and humidity over the adjacent urban landscape. Despite the idealized setup in the numerical simulation and the simplifying assumptions in analytical models, some insight of urban nocturnal cooling phenomenon mediated by bluespace can be learned. The key findings are summarized below.
1. Time-varying surface temperature over different surface types (i.e., a square-root time dependence) leads to distinct temporal change in meteorological quantities. Both observations and numerical simulation results demonstrate that the temporal change of the cooling rate of T 3.5 averaged over lake and greenspace surface types decreases monotonically with time, while that over urban land surface exhibits accelerated cooling in the first 2 h but decreases afterward. The apparent temperature follows a similar trend as T 3.5 . As surface temperatures further drop below the air temperature, both intensity and direction of the micro-scale air circulations change. 2. For urban areas close to open water surfaces in the studied sites, the nighttime air temperature increases, and humidity decreases with distance to water under typical calm weather condition in summer as demonstrated by both simulation and observational results along the selected sites in the measurement network. Our results further confirm that the apparent temperature in urban area decreases with distance to water, highlighting that the overall thermal comfort in urban areas close to open water surface can be compromised by higher humidity closer to the lake shore (Hu and Li 2020) (i.e., as high as 0.8 °C due to increased near-shore humidity). 3. The proposed analytical model A considers a negative feedback mechanism between air cooling rate and the differences between air and surface temperatures under the chosen urban climatological condition. With model parameter being constrained by observational data, the model results agree well with the observed average air temperature evolution for urban, greenspace, and lake surfaces. In addition, the spatial extent of lakes' influence on temperature and humidity over urban surfaces for the selected study sites can be modeled by an advection-diffusion model constrained by surface energy balance (model B, Section 4.2). Overall for both models, despite uncertainties in model parameters, they reproduce key features and give quantitative evaluation with a reasonable accuracy without requiring sophisticated computational models.
The study is limited in the following aspects. First, in the simulation, only the first-order results can be simulated by means of Reynolds-averaged Navier-Stokes modeling approach. Hence, future studies using more advanced computational fluid dynamics modeling techniques are required to understand the driving forcing of the reversed circulations. Second, the surface temperature cooling is modulated by the balance between the longwave radiation and surface heat conduction. While this is a good approximation under clear calm condition, it can be invalid under windy conditions. Furthermore, nocturnal evapotranspiration over urban surfaces and greenspace is assumed to be negligible in the numerical simulation. Although the latent heat flux in urban settings at night is small according to pervious research, the simulated 3.5-m air humidity may be underestimated. Future studies may investigate how urban surface temperature varies under other weather conditions and how advection under windy conditions can be explicitly considered in the proposed analytical model for prediction of nocturnal urban cooling mediated by urban greenspace and bluespace.