The utility of distributed hydrologic models in combination with high-resolution Weather Surveillance Radar-1988 Doppler (WSR-88D) rainfall estimates for flash flood forecasting in urban drainage basins is examined through model simulations of 10 flood events in the 14.3 km2 Dead Run watershed of Baltimore County, Maryland. The hydrologic model consists of a simple infiltration model and a geomorphological instantaneous unit hydrograph–based representation of hillslope and channel response. Analyses are based on high-resolution radar rainfall estimates from the Sterling, Virginia, WSR-88D and observations from a nested network of 6 stream gauges in the Dead Run watershed and a network of 17 rain gauge stations in Dead Run. For the three largest flood peaks in Dead Run, including the record flood on 7 July 2004, hydrologic model forecasts do not capture the pronounced attenuation of flood peaks. Hydraulic controls imposed by valley bottom constrictions associated with bridges and bridge abutments are a dominant element of the extreme flood response of small urban watersheds. Model analyses suggest that a major limitation on the accuracy of flash flood forecasting in urban watersheds is imposed by storm water management infrastructure. Model analyses also suggest that there is potential for improving model forecasts through the utilization of information on initial soil moisture storage. Errors in the rainfall field, especially those linked to bias correction, are the largest source of uncertainty in quantitative flash flood forecasting. Bias correction of radar rainfall estimates is an important element of flash flood forecasting systems.
The utility of distributed hydrologic models in combination with high-resolution Weather Surveillance Radar-1988 Doppler (WSR-88D) rainfall estimates is examined in this paper for flash flood forecasting in small urban watersheds. The hydrologic model (Morrison and Smith 2001; Giannoni et al. 2003; Smith et al. 2005a) consists of a simple infiltration model (see Ogden and Saghafian 1997) and a geomorphological instantaneous unit hydrograph (GIUH)-based representation of hillslope and channel response (see Rodriguez-Iturbe and Rinaldo 1997). The structure of the hydrologic model is dictated by computational constraints on real-time flash flood forecasting. The forecasting model has been implemented for the 14.3 km2 Dead Run watershed in Baltimore County, Maryland (see Smith et al. 2005b; Nelson et al. 2006). During the summers of 2004 and 2005, the modeling system was run in real time with a data feed of digital hybrid reflectivity (DHR; see Fulton et al. 1998) fields from the Sterling, Virginia (KLWX), WSR-88D.
Extensive development of distributed hydrologic models for flood forecasting applications has occurred during the last decade (Burnash 1995; Koren et al. 2004; Vieux 2001; Vieux et al. 2004; Julien et al. 1995; Smith et al. 2004). Similarly, there have been extensive efforts aimed at developing calibration procedures for hydrologic models (Duan et al. 1992; Gupta et al. 1998; Vrugt et al. 2005). In this paper, we assess the applicability of a distributed hydrologic model for flash flood forecasting through model analyses in an urban watershed with a nested network of stream gauging stations.
Dead Run is a 14.3 km2 tributary of the Gwynns Falls watershed in Baltimore County, Maryland (Fig. 1). Gwynns Falls is the principal study watershed of the Baltimore Ecosystem Study (BES; see Groffman et al. 2003; Smith et al. 2005b) and has a dense network of stream gauges that have been deployed in connection with BES hydrologic monitoring activities. A nested network with six stream gauges has been deployed in the Dead Run watershed at basin scales ranging from 1.2 to 14.3 km2. The stream gauging network and a dense network of rain gauges in Dead Run provide exceptional observational resources for the assessment of flash flood forecasting systems in urban watersheds. In addition, the Dead Run watershed has experienced a number of major flood events in 2004 and 2005, including the record flood in Dead Run (from a stream gauging archive that extends back to the 1950s), which occurred on 7 July 2004.
This paper is organized as follows. In section 2, a description of the Dead Run watershed and the hydrological and meteorological datasets that have been developed in Dead Run is presented. The structure and governing equations of the Network Model are discussed in section 3. Results of model analyses are presented in section 4 and in section 5 conclusions from model analyses are summarized.
2. Data and basin description
The Dead Run watershed was settled in the eighteenth century and has undergone intense urbanization since the late 1950s (Nelson et al. 2006). The basin contains a diverse mix of urban and suburban development and has an impervious fraction of approximately 36%. As a tributary of the Gwynns Falls watershed, which is the main study watershed of the BES, Dead Run has been extensively studied (see Smith et al. 2005b; Nelson et al. 2006). Hydrologic model analyses utilize observations collected during intensive summer field campaigns during the 2003–05 time period. In particular, analyses center on 10 significant flood events (Table 1) during the 3-yr period.
A nested network of five stream gauges (denoted DR1–DR5 in Fig. 1) has been installed in the Dead Run watershed to augment the U.S. Geological Survey (USGS) stream gauge at Franklintown (Fig. 1). Drainage areas for the six stream gauging stations range from 1.2 to 14.3 km2 (1.2 for DR1, 1.8 for DR2, 1.9 for DR5, 4.9 for DR3, 6.1 for DR4, and 14.3 km2 for Franklintown). Time resolution of the stage observations is 1 min for the DR1–DR5 stations and 5 min for the Franklintown USGS gauge. Rating curves for the five new stream gauges (DR1–DR5) have been developed from direct discharge measurements and hydraulic model studies (Smith et al. 2005b; Nelson et al. 2006) and are used to create 1-min discharge time series. Discharge datasets from these six stations play a central role in examining the performance of flash flood forecasting systems. Basin boundaries and stream channels were delineated from the USGS 1″ (30 m) Digital Elevation Model.
A network of 17 rain gauge stations has been deployed in the Dead Run watershed (Fig. 1). Each station has two Tru-check rain gauges that are separated by approximately 1 m and read manually after each storm event (Smith et al. 2005b, 2007). The storm total rainfall accumulation for a station is computed as the average of the accumulations from the two gauges, provided that observations from both gauges are available (occasionally, one or both gauges would be tipped or missing). The rain gauge network provides storm total accumulations that are used for bias correction of WSR-88D rainfall estimates.
Radar reflectivity observations are available from KLWX at 5–6-min time resolution and 1° azimuthal by 1-km range resolution. The polar reflectivity data are processed to 5-min and 1-km horizontal resolutions for rainfall analyses (Baeck and Smith 1998; Fulton et al. 1998; Smith et al. 2007). Reflectivity data have been utilized in two formats. For real-time implementation, the DHR product (see Fulton et al. 1998) has been utilized. An ftp protocol has been implemented for the automatic transfer of the DHR product from the National Weather Service to Princeton at the completion of each volume scan. For postanalysis applications, the archive level II volume scan reflectivity data have been obtained from the National Climatic Data Center. For rainfall estimation over Dead Run, which is approximately 70 km from the KLWX radar, the 0.5° elevation angle is used for both reflectivity products.
Rainfall rate fields, R(t, x) (mm h−1), are estimated by applying a Z–R relationship [R = aZb; Z is the radar reflectivity factor (mm6 m−3)] to the reflectivity observations. Rainfall analyses in this paper utilize the “convective” Z–R relationship with a = 0.0174 and b = 0.71. A “local” bias correction algorithm is applied to the radar rainfall estimates; the bias-corrected rainfall estimates take the form
m is the number of rain gauge locations, and xj is the spatial location of the jth rain gauge; Gij denotes the storm total accumulation (mm) for the ith storm from the jth gauge; and Ti is the duration of the ith storm. The bias correction factor is thus the ratio of the mean storm total rainfall from rain gauges to the mean storm total rainfall from radar bins containing the rain gauges (Baeck and Smith 1998). There are 17 rain gauges available for Dead Run and the bias correction for the 10 events in Table 1 ranged from 0.70 to 2.77, with 8 of the 10 having bias values greater than 1, denoting an underestimation by the radar [see Smith et al. (2007) for additional discussion of distributional properties of the multiplicative bias and strategies for the operational estimation of bias]. Rain gauge observations provide the mean storm total rainfall accumulation over the watershed and radar observations provide the temporal and spatial variabilities of the rainfall rates.
Storm total rainfall and runoff for the Dead Run basin above Franklintown (Table 1) were obtained by time integrating the rainfall rate time series and the discharge observations. Infiltration was computed as the residual between the rainfall and the runoff. The three largest events based on storm total rainfall are the 7 July 2004, 27–28 July 2004, and 28 June 2005 storms. The response time of Dead Run at Franklintown, characterized by the lag time, is typically on the order of 60 min (see also Smith et al. 2005b). The lag time is defined as the difference between the time of peak discharge and the time of peak of the basin-averaged rainfall rate. For the three largest events, the lag times are significantly longer than 60 min (80–100 min).
3. The Network Model
The Network Model is a distributed hydrologic model (Morrison and Smith 2001; Zhang et al. 2001; Zhang and Smith 2003; Giannoni et al. 2003; Turner-Gillespie et al. 2003; Smith et al. 2005a, Hicks et al. 2005; Slutzman and Smith 2006) that has been adapted for flash flood forecasting in Dead Run. The Network Model partitions the drainage basin into hillslope and channel components and represents discharge at the outlet of a drainage basin as
where Q(t) denotes the discharge (m3 s−1) at time t (s), A is the domain of the drainage basin, x is a point within A, d0(x) is the distance (m) from x to the channel network, υ0 is the overland flow velocity (m s−1), d1(x) is the distance (m) along the channel from x to the basin outlet, υ1 is the channel flow velocity (m s−1), and M(t, x) is the runoff rate (m s−1) at time t and location x. The total flow distance from x to the basin outlet is d0(x) + d1(x), the sum of the overland flow distance and the channel flow distance.
The runoff rate M(t, x) at time t and location x is computed from the rainfall rate R(t, x) using the Green–Ampt infiltration model with moisture redistribution (Ogden and Saghafian 1997). Implementation of the model requires specification of the rainfall field, soil hydraulic parameters [saturated hydraulic conductivity and tension head at the wetting front; see Ogden and Saghafian (1997)], and the overland flow distance and channel flow distance functions (which are obtained from the extracted drainage network). The saturated hydraulic conductivity Ks is the principal infiltration model parameter that is used for site-specific calibration. The initial moisture θ is a state parameter that must be specified for model implementation.
The Network Model computes discharge as the convolution of a GIUH-based scheme with the runoff field. The GIUH (Rodriguez-Iturbe and Valdes 1979; Gupta et al. 1980; Rodriguez-Iturbe and Rinaldo 1997) specifies the unit response of the drainage basin to an instantaneous pulse of runoff of unit depth. The Network Model formulation of the basin response leads to the following representation for the GIUH:
where Mδt(t, x) is equal to δt−1 for t in (0, δt] and 0 otherwise. The expression can be simplified to
The Network Model computes the travel times from each model grid cell to each location on the drainage network for which model hydrographs are desired. Computing travel times requires only the distance functions d0(x) and d1(x) and the velocity parameters υ0 and υ1. A lookup table of travel times is computed in a preprocessing step and used in the model simulations. Runoff is generated as the difference between the rainfall rate and infiltration rate, as computed from the Green–Ampt algorithm. Knowing the travel time of each grid cell to a forecast point means that the discharge at any time t can be computed as the sum of runoff contributions from individual grid cells in the basin, using the appropriate lag times computed from the lookup table. There is no need for routing flow from grid cell to grid cell. This feature of the Network Model is the key to achieving run times that are suitable for flash flood applications in urban drainage basins. Another attractive feature of the GIUH representation is that grid cell size does not have a large impact on routing and can be selected to best represent spatial heterogeneities in infiltration properties (especially if high-resolution datasets are available for specifying impervious cover).
Model analyses in the following section use a reference set of model parameters that were obtained by modeling two flood events (Fig. 2) produced by short (15 min) bursts of extreme rain rates (exceeding 100 mm h−1). The 12 June 2003 event produced a 1.31 m3 s−1 km−2 flood peak and a total runoff of 8 mm from a storm total accumulation of 17 mm (see Table 1) while the 5 November 2003 event produced a flood peak of 1.68 m3 s−1 km−2 and a total runoff of 15 mm from a storm total rainfall accumulation of 23 mm. Model parameters used for simulating the 12 June 2003 and 5 November 2003 events are υ0 = 0.088 m s−1, υ1 = 1.7 m s−1, and KS = 10 mm h−1.
Simulations were performed for the 10 events (see Table 1) to assess the potential accuracy of forecasting flash floods using WSR-88D rainfall estimates and a simple distributed hydrologic model. Three variables were chosen to compare the observed and the model discharge. These are the peak magnitude, time to peak, and the runoff volume “around the peak.” The volume around the peak is taken as the total volume during the 60-min time window centered on the peak discharge. Runoff volume is converted to depth (mm) through division by the drainage area of the basin.
Simulations were performed for the 10 events using model parameters υ0 = 0.088 m s−1, υ1 = 1.7 m s−1, and KS = 10 mm h−1. Use of the model parameters derived for the two events presented in the previous section reflects the typical setting of operational flash flood forecasting in which limited information is available for model calibration. Calibration of hydrologic models used for quantitative flash flood forecasting in urban environments requires efficient use of available stream gauging observations, driven by a solid understanding of the main sources of the regional variability in the flood response. We examine these issues in more detail below.
Two scenarios were utilized for specifying the initial moisture content. In one set of model runs, the initial moisture content was treated as a calibration parameter and each of the 10 events utilized a different initial moisture content that was calibrated for the event (results for this set of model runs are summarized in Table 2). In a second set of model runs, a constant initial moisture content with θ = 0.25 was used for each of the 10 simulations. The value of 0.25 was the average initial moisture for the 10 event simulations with variable initial moisture content (Table 2). These two sets of model runs are used below to examine the model performance, including assessments of the potential utility of information on the initial moisture content.
Model flood peaks and runoff volumes (around the peak) match the observed peaks and volumes reasonably well (Figs. 3 and 4), even under the constraint that a constant initial moisture content of 0.25 is employed. The coefficient of determination (R2), which is a measure of the correlation between the observed and simulated hydrographs, is 0.81 for the flood peak and 0.90 for the flood volumes.
Timing errors (Fig. 5) for model simulations range from −11 min (model peaks too slow) to +48 min (model peaks too fast). Errors in model forecasts of peak time are strongly dependent on discharge, with timing errors taking larger, positive values with increasing peak discharge. Timing errors depend solely on the GIUH parameters υ0 and υ1. They are largely insensitive to other model parameters and bias correction.
Sensitivity to initial soil moisture is investigated by allowing the initial moisture content to vary from event to event (results for these analyses are summarized in Table 2). The coefficient of determination increases from 0.81 to 0.96 for peak magnitude if the soil moisture is calibrated for each event. For runoff volume the coefficient of determination increases from 0.90 to 0.99. These results suggest the potential for improving the peak forecasts through information on the initial soil moisture conditions. In Smith et al. (2005b), it is shown that storm event runoff ratio in Dead Run varies systematically over the warm season (April–September), suggesting that seasonally varying initial moisture parameters could capture some of the improvement in the peak forecasts reflected in the variable initial moisture condition simulations. Utilizing continuous simulation models provides a long-term solution to the problem of representing the dependence of flash flood response on initial moisture conditions (see, e.g., Ogden and Saghafian 1997).
Rainfall estimation bias is a dominant source of error in model forecasts of flood peak magnitude. The coefficient of determination decreases from 0.81 to 0.68 for the case of no bias correction. Bias values for 8 of the 10 events are larger than 1 and Smith et al. (2007) show that bias varies systematically with event magnitude. Developing operational procedures that accommodate these systematic errors in radar rainfall estimates for flash flood events is an important element of the operational enhancement of flash flood forecasting procedures (see Smith et al. 2007 for additional discussion).
Limitations of GIUH-based models are apparent from model simulations for the 10 events. Hydraulic effects that are not accounted for in the Network Model are a source of errors in model analyses. The early rise of the simulated hydrographs compared with the observed hydrographs (see Fig. 2) is one example. The GIUH assumes that flow velocities are time invariant. In stream channels and storm sewers, however, flow velocities are not constant over time, but vary with flow magnitude, as reflected in Manning’s equation. The constant velocity assumption, which allows for computational efficiency, reflects mean conditions during the central portion of the hydrograph. Actual flow velocities are lower than model velocities during the early period of a flood, so the model hydrograph rises earlier than does the observed hydrograph.
During extreme flood events, significant flow volumes are conveyed by floodplains. For the major flood events in Dead Run, hydraulic processes associated with coupled channel–floodplain flow result in decreased flow velocities and significantly longer response times. Of particular importance in Dead Run are the effects of valley bottom constrictions associated with bridges and bridge abutments. In the DR3–DR4 reach, bridge crossings produce extensive areas of backwater for large floods (based on field observations and detailed high water maps developed from the 7 July 2004 and 28 June 2005 flood events). The storage effects of the DR3–DR4 reach explain the anomalously large flow attenuation at Franklintown for the three largest events: 7 July 2004, 27 July 2004, and 28 June 2005 (see Fig. 6). For these three events, model simulations peak 18–48 min before the observed discharge peaks. To capture these aspects of the flood response in flash flood forecasting systems requires both hydraulic models and detailed information on constrictions in valley bottoms.
Model simulations for flood events with peak discharge less than 2.0 m3 s−1 km−2 (at Franklintown) capture the basin response at the 14.3 km2 scale quite well, as illustrated by results for the 18 August 2004 flood (Fig. 7), which produced a peak unit discharge of 1.4 m3 s−1 km−2. Results for the nested set of gauges with drainage areas ranging from 1.3 to 6.1 km2, however, are mixed (see, e.g., Figs. 8 and 9). Model errors are tied both to heterogeneities of the basin response within the Dead Run basin and the difficulties of estimating the rain rate accurately as the spatial scale decreases below 10 km2 and the time scale of the basin response decreases below 60 min.
Model simulations for the 18 August 2004 event at Franklintown are among the best of the 10 events in terms of peak time, peak volume, and time to peak (Table 2), but there are large errors in peaks, volumes, and timing at each of the internal gauges. In DR3 (Fig. 8), two peaks in rainfall rate between 0215 and 0315 UTC produced two peaks in observed discharge in DR1, with the first peak the largest in the observed discharge time series. The model simulation for DR1 underestimates the first peak and overestimates the second in DR1, resulting in severe underestimation in the model simulation for DR3. The situation is reversed in the DR4 subbasin, with overestimation and a model peak that is later than the observed peak (Fig. 9). Model hydrographs for the headwater catchments DR1, DR2, and DR5 are overdispersed relative to the observed hydrographs (as illustrated in Figs. 8 and 9). The same holds for the downstream stations at DR3 and DR4, but to a lesser degree. Accurate simulation results at Franklintown (Fig. 7) are due in part to compensating errors from the DR3 and DR4 subbasins.
The interplay of storm water management infrastructure and transportation infrastructure can have a marked impact on urban flood response. In the DR4 subbasin, the response to rainfall is characterized by two distinct hydrograph peaks, even when the flood-producing rainfall is concentrated during a single period, as illustrated for the 13 August 2004 storm (Fig. 10; see also Figs. 9 and 11). This aspect of the basin response is tied to the structure of the storm drain network, and in particular the truncation of drainage into two separate units: one upstream of the Baltimore beltway (represented by the DR5 observations) and the second downstream of the beltway. The first peak in DR4 is from the local storm drain network, downstream of the beltway; the second peak is from the DR5 portion of the watershed.
There are systematic timing differences between the DR3 and DR4 subbasins of Dead Run. The most extreme case is the 13 August 2004 flood (figure not shown) for which the model hydrograph at DR3 peaks 14 min before the observed peak and the DR4 model peak is 25 min after the observed peak. The DR3 basin includes significant areas developed after the advent of storm water management regulations. As a consequence, there are numerous storm water detention structures in the DR3 basin (Smith et al. 2005b; Nelson et al. 2006), which lead to a slower response than in DR4, which has few storm water detention structures (Fig. 1). In addition, old residential development in DR4 is characterized by dense storm drain networks that rapidly concentrate runoff (Smith et al. 2005a). Capturing the impacts of storm water detention basins on flood response poses an especially difficult challenge in developing flash flood forecasting systems for urban drainage basins.
Simulation of the 28 June 2005 flood (Fig. 11; see Fig. 6c for the Franklintown simulation) produced flood peaks that are larger than the observed peak for DR4 and smaller that the observed peak at DR3. The underestimation in DR1 is even larger than at the downstream DR3 station. The 28 June 2005 storm produced rainfall accumulations over the DR1 basin that approached 100 mm in less than 30 min, based on four rain gauge observations in and adjacent to DR1 (Fig. 12). WSR-88D rainfall estimates were generally quite good for the 28 June 2005 storm, but they did not capture the peak rainfall accumulations in DR1 (Fig. 12). Consequently, the rainfall rate estimates used to produce the DR1 model simulations (Fig. 11) reflect a severe underestimation of rainfall accumulation over DR1. If the four gauges in DR1 alone are used to estimate the multiplicative bias, it increases to 1.73 from a value of 1.20 if all 17 gauges are utilized. These results highlight the dominant role of rainfall estimation error and bias correction for quantitative flash flood forecasting.
The preceding analyses highlight the difficulties inherent in quantitative flash flood forecasting in urban drainage basins that are smaller than 10 km2 in area. Problems associated with rainfall estimation error, heterogeneity of land surface properties, and the role of storm water detention structures all result in sharp increases in simulation errors as the basin scale decreases below 10 km2.
Calibration of hydrologic models for flash flood forecasting in urban regions is an important and difficult topic. There is a significant role for stream gauging observations using existing stream gauging resources, although the density of stream gauging in urban regions is often poor. There is also a role for short-term stream gauging networks like the Dead Run network. In this case, there is substantial benefit to developing nested stream gauging systems, in which it is possible to characterize the scale-dependent flood response. The transfer of information from gauged to ungauged watersheds is a key element of model implementation for urban regions. Developing quantitative measures of urban development has provided one avenue for parameter transfer from gauged to ungauged watersheds (Beighley and Moglen 2003). Simple metrics, like the fraction of impervious cover, can be useful, but results from this study also point to the utility of quantifying the extent of the storm water management infrastructure. It is also shown that the urban drainage network, comprising both the surface channel and storm drain network, is a key element of urban flood response. The increasing availability of digital datasets characterizing storm drain networks provides an important approach to improving parameterizations of hydrologic models used for flash flood forecasting.
5. Summary and conclusions
The utility of distributed hydrologic models, in combination with high-resolution WSR-88D rainfall estimates, for operational flash flood forecasting was examined through model simulations for 10 storm events in the Dead Run watershed. Analyses utilize observations from a nested set of six stream gauges in the 14.3 km2 watershed. The following are the principal observations of this study:
Distributed hydrologic models and high-resolution WSR-88D rainfall rate fields can provide important elements of site-specific flash flood forecasting systems in small urban watersheds. From the perspective of flash flood forecasting, the hydrologic response of Dead Run at the 14.3 km2 scale can be accurately modeled using the Network Model and bias-corrected WSR-88D rainfall fields, subject to the constraints discussed below. The GIUH-based framework of the Network Model provides a computationally efficient modeling system for flash flood forecasting.
Hydraulic controls imposed by valley bottom constrictions associated with bridges and bridge abutments are a dominant element of the extreme flood response of small urban watersheds. For the three largest flood peaks in Dead Run, including the record flood on 7 July 2004, hydrologic model simulations do not capture the pronounced attenuation of flood peaks. As a consequence, hydrologic model forecasts of flood peak magnitude and the timing of floods will exhibit errors that are strongly dependent on the flood magnitude.
There is potential for improving model forecasts through the utilization of information on initial soil moisture storage. A simple first step would be to incorporate a seasonally varying initial moisture content. Continuous soil moisture accounting models also provide an avenue for capturing the dependence of the flash flood response in urban drainage basins on initial moisture conditions.
A limitation on the accuracy of flash flood forecasting in urban watersheds is imposed by storm water management infrastructure. Heterogeneity of flood response in the Dead Run watershed is closely linked with storm water management infrastructure. Capturing the effects of storm water infrastructure on flash flood response requires information on the physical structures and modeling systems to deal with these features.
Decreasing the basin scale of flash flood forecasting below 10 km2 results in a sharp decrease in predictability. Difficulties in flash flood forecasting at spatial scales smaller than 10 km2 are tied both to heterogeneities of the hydrologic response and the accuracy of radar rainfall estimates.
Errors in the rainfall field produce the largest sources of uncertainty in quantitative flash flood forecasting. Bias correction of radar rainfall estimates is an important element of flash flood forecasting systems.
The authors acknowledge Gary Fisher, Peter Nelson, Elliot Holland, Matt Ballantine, Tammy Newcomer, Bernice Rosenzweig, Andy Brett, Alex Lester, Andy Lapetina, and Jane Diehl for assistance in data collection and analysis. The research was supported by the NOAA/National Weather Service Office of Hydrologic Development and the National Science Foundation (NSF Grants EAR-0208269, EAR-0409501, and ITR-0427325).
Corresponding author address: James A. Smith, Dept. of Civil and Environmental Engineering, Princeton University, Princeton, NJ 08544. Email: email@example.com