***UPDATE 30 AUGUST 2012***
Just a note to mention that a revised version of this manuscript has been accepted for publication in the Springer journal Natural Hazards. The Canadian Geotechnical Journal (rightfully) rejected the paper on the grounds that it has limited direct relationship to their publicaion objectives, so I submitted it to Natural Hazards instead, in early July. I made significant changes based on some very good review comment, and it was quickly accepted upon revision.
I think this will likely hit the street within a few months, depending on how quick they are with the finer details of editing and publishing. I haven’t been through the process with this specific journal before.
I will try to place a copy of the accepted draft mauscript up shortly.
Here’s another draft journal paper, just now uploaded to the Canadian Geotechnical Journal for consideration. It consolidates previous posts about population and road density in eastern Canada and Saint Lucia.
Road density as a proxy for population density in regional scale risk modelling
BGC Engineering Inc., Victoria, BC
This paper examines the use of road network data as a proxy for interpreting population density, which is of use in regional scale qualitative risk assessment for natural hazards. Comparison of available road network and population data at various scales in Ontario and Quebec yields a best fit relationship of DP = 28 DR^2. Analysis of available high resolution topographic data for the Caribbean island nation of Saint Lucia suggests similar power law trends, with expected population density in Saint Lucia roughly half that of Canada for the same road density. Together, these findings suggest that DP ~ 10 to 30 DR^2 may represent a useful range broadly applicable for a wide variety of geographic, climactic and socio-economic settings. The Canadian relationship has been used to generate a population density model for the lowlands of eastern Ontario and southern Quebec, and this model has been compared with the spatial distribution of seismic hazard to develop a qualitative seismic risk map. The seismic risk map, presented primarily for illustrative purposes, shows elevated seismic risk in urban centres in the study area, and along a predominantly rural area east of Quebec City on both shores of the Saint Lawrence River.
Hazard, risk, population density, road density, GIS
Certain natural processes with potential to affect human health or property are widely distributed in parts of Canada. In examining hazard and risk associated with such processes, a first level assessment may be made at the regional, provincial or national scale, as a means of establishing meaningful priorities for subsequent, more detailed study. The spatial distribution of certain hazards, for example strong motions associated with earthquakes, is well understood. Other hazardous natural processes, such as flooding or landslides, may be spatially distributed more chaotically, but can in some cases be modeled to develop estimates of probable magnitude and recurrence over some meaningful area. A realistic comparison of the spatial distribution of such hazards with spatial distribution of human habitation and built infrastructure would allow the development of a risk analysis on a broad scale. Population data are not readily available in a form suitable for a meaningful comparison in risk analysis at the regional scale. This paper examines the use of road density as a proxy for population density, in an attempt to develop population density estimates with sufficient resolution for regional scale risk modeling. Strong ground motion associated with earthquake shaking in the lowlands of eastern Ontario and southern Quebec is used as a working example to illustrate the potential use of this proxy population model in regional scale qualitative risk analysis.
2 BACKGROUND: HAZARD AND RISK AT THE REGIONAL SCALE
Canada is large, with very diverse geology, geography, and physiography, and is thus affected by numerous natural hazards, which are each more or less common in certain parts of the country. This includes earthquakes, landslides, flooding, and weather hazards, as examples. Each different natural hazard is distributed differently in space and time across Canada, thus different areas are affected by a different balance of natural hazards. These hazards are a significant concern only where there is potential for them to cause damage. Rational efforts to manage the effects of these hazards should prioritize the allocation of resources, to be used perhaps for implementation of further study, construction of physical mitigative measures, or preparation of emergency response plans. Prioritization should be based on the relative value of expected loss, or risk, associated with these hazards.
Risk is a combination of hazard and consequence, where hazard is a measure of the probability of a potentially hazardous event occurring, and consequence is a measure of the probability of such an event generating harmful effects (i.e. damage, or loss) to one or more elements at risk. Consequence can be taken as the combination of likelihood of presence of a given element at risk, given occurrence of the hazard, and vulnerability of that element at risk to the given hazard. Therefore, one can develop an appreciation of risk, or expected loss, with prior knowledge of the hazard, and of the location and vulnerability of elements at risk, in both space and time.
Risk, like the natural hazards previously mentioned, is not distributed uniformly across the country, since it is a function of the presence of both hazards and elements at risk in space and time. Human settlement in Canada is concentrated relatively close to its southern border with the United States, with a relatively small number of urban centres and extensive rural settlement further north from the border. Thus human occupants and the associated built infrastructure tend to be concentrated close to the southern border and in a small number of other urban centres, rather than being distributed uniformly across the country. Since natural hazards are also not distributed evenly across the country, an estimation of the risk, or expected loss, associated with a given hazard, requires spatial comparison of the relative likelihood of both the hazard and the presence of human development. A population density map is therefore a necessary component of natural hazard risk assessment, along with a map representing spatial distribution of the hazard. The resolution of a risk map would depend on the resolution of the population and hazard models, therefore a high resolution population density map is desirable for use in regional scale qualitative risk analysis.
The present work was motivated by ongoing work in landslide susceptibility mapping in eastern Canada (Quinn, et al. 2010), and by subsequent explorations of hazard and risk associated with large landslides in sensitive clay in eastern Canada (Quinn et al. 2011). The latter paper discusses hazard and risk associated with these landslides in a conceptual way, without attempting to develop hazard or risk maps. Subsequent efforts have been made toward modeling risk, with renewed emphasis on examining risk to human habitation following a fatal landslide in Saint-Jude, Quebec in May, 2010. Initial efforts by the author to model landslide risk at the regional scale in Quebec were frustrated by a lack of meaningful population and dwelling data. This paper presents a methodology for developing rapid, meaningful population density models for Canada at the regional scale, for use in qualitative risk assessment for natural hazards.
3 POPULATION DENSITY IN A CANADIAN SETTING
Population data are recently (i.e. as of February, 2012) freely available to the general public from Statistics Canada (2007) at a variety of levels of aggregation, primarily based on political boundaries (i.e. national, provincial and municipal). Figure 1 shows municipal boundaries in most of Ontario and Quebec. The Figure also shows the approximate outline of the former Champlain Sea, which occupied a depression south of the retreating continental ice following the last glacial episode roughly 10,000 years ago. That boundary is shown to emphasize the area of interest for large landslides in eastern Canada, which is the specific geological hazard that motivated this investigation.
Figure 2 shows population density, based on the municipal boundaries, within the area of the former Champlain Sea. That Figure shows major concentrations of settlement in the Ottawa-Hull, Montreal and Quebec City areas, with a number of other smaller concentrations of population in eastern Ontario and south-central Quebec. The initial impression one obtains from this map is that the data yield a reasonable representation of population distribution, with major clusters corresponding to known urban centres. At first glance, the map therefore appears to be a suitable candidate for use in extending hazard models to risk models. However, closer examination reveals some significant issues with data resolution at larger scale, limiting the meaningfulness of resulting risk models.
Figure 3 shows the same representation of population density, enlarged to focus on the Ottawa-Hull area along the provincial border between Ontario and Quebec. This map suggests that Hull has the highest population density in this area (1630 / km^2, extracted from the population density raster), followed by Gatineau (652 / km^2) and the greater Ottawa area, including Orleans, Kanata, and other outlying areas (267 / km^2). These numbers are accurate for the indicated municipal boundaries, however they are not particularly meaningful for direct comparison, since they are obtained for municipalities with dramatically different sizes. One would expect the central core of Ottawa, on the Ontario side, to have a similar population density as the central urban core in the Hull-Gatineau area on the Quebec side, yet this is not obvious from the population densities illustrated on this basis. This anticipated similarity is indicated in Figure 4, which includes the road network overlain on top of the population density map. On the Ontario side, the greater Ottawa area clearly includes a densely developed core and dense major suburbs to the east (Orleans) and west (Kanata). However, there are also large, sparsely developed areas, including a major greenbelt with limited development between the central core and major suburbs, and significant expanses of outlying rural areas interspersed with smaller rural communities. The population density within this particular municipality is therefore not represented at a resolution suitable for meaningful risk analysis when modeled on the basis of municipal boundaries.
Population data are available for selected areas at a finer scale of subdivision in Canada. The municipal boundaries shown in Figure 1 correspond with Census Subdivisions (CSDs) used by Statistics Canada. Selected urban areas are further subdivided for census purposes into Census Tracts (CTs). Figure 5 shows the distribution of CTs in southern Ontario. CTs are present in communities with populations of 50,000 or more, and subdivide those municipalities into smaller areas with relatively stable populations between 2,500 and 8,000. Figure 6 shows CTs on the Ontario side of the Ottawa-Hull area. It is immediately evident that the CTs in this area would yield higher resolution population density information, where they exist, than that obtained from municipal boundaries alone. However, it can be seen that there are pockets of local development outside the greater Ottawa area not specifically captured by CTs, and that the density of development appears to vary within many of the CTs, in some cases significantly.
The evident disparity between local development, as inferred from the road network density, and average population density, led to the idea of using road network density as a proxy for population density on a regional, provincial or even national scale. Some work has been done in this area by Owens et al. (2010) for modeling rural population distribution to support health-related research. The present work attempts to develop approximate relationships between road and population density across a broader range of density classes, including urban areas, with consideration of different geographic settings in eastern Canada, and including comparison with related data from the Caribbean island nation of Saint Lucia.
4 ROAD DENSITY AND ITS RELATIONSHIP WITH POPULATION DENSITY
Population and dwelling counts from the 2006 Canadian census and road network data were obtained from Statistics Canada (2007). Population and dwelling data were collected at the CSD (i.e. municipality) level in Ontario, and at the CT (i.e. 2,500 to 8,000 person block) level in Quebec. The data were manipulated in a geographic information system (GIS), which was used to determine the total length of roads and total population within each census unit. These data were used with calculated surface area of each CSD and CT to obtain road (DR) and population (DP) densities, which are plotted against each other in Figure 7 and Figure 8.
Figure 7 shows all of the population versus road data for both provinces. While there is considerable scatter in the data, there is a clear trend for increasing population density with increasing road density, with them being related by a power law, with a best fit of DP = 28 DR^2, and a probable lower bound, capturing ~ 90 % of the data of DP ~ 13 DR^2. Figure 8 shows population versus road density for lower density values, where it can be seen that the general trends retain their meaning, although the lower bound trend may be a slightly better “best fit” for the Ontario CSD data, which are inferred from municipal boundaries with generally large areas. Effects of scale are examined further in Figure 9, which shows population versus road densities drawn from census divisions of different size, including the smallest (< 1 km>2), medium (> 1 km>2) and largest (> 2 km>2) CTs in Quebec, and the smallest (< 1000 km>2) and largest (> 1000 km>2) CSDs in Ontario. The best fit power law relationship maintains a good fit over several orders of magnitude, being weaker below ~ 50 persons / km^2, where there is considerably more scatter in the data. The data also appear to be trending upward above the best fit line at more than 10,000 persons / km^2. However, the range of good fit is judged to cover most of the range of practical importance. Population densities below 50 / km^2 constitute very low density, and it is not important in regional scale modeling if the true density is perhaps 5 or 100 instead of 50; similarly, densities above 10,000 / km^2 are all very high and potential errors are expected to be of limited significance in risk modeling at the regional scale.
High quality topographic data became available for the island of Saint Lucia during this work, and similar analyses were conducted for comparative purposes. The Government of Saint Lucia (personal communication, 2012) provided a copy of recent (2010) topographic data for the island, including the complete road network, and footprint of all buildings. Figure 10 shows a map of building density for Saint Lucia. The total footprint of all buildings (8.35 million m^2) can be compared with the total population (174,000) to obtain an average occupancy of 0.0208 persons / m^2 of building footprint. The building density map is converted to estimated population density by direct multiplication. Figure 11 shows the distribution of road density for Saint Lucia, determined in GIS from the available road network vector (polyline) data. One can see the building density and road density maps are qualitatively similar. They were compared statistically by generating 100,000 random points across the island, and comparing the inferred population and road densities at each point. The resulting data are plotted in Figure 12, which also shows the best fit and probable lower bound trends from the Canadian data. It can be seen that the same general power law trend is evident, although the Canadian best fit appears to be closer to an upper bound for the Saint Lucia data, and the Canadian lower bound is relatively close to a best fit for Saint Lucia.
The Canadian and Saint Lucian data are compared further in Figure 13, which shows best fit and lower bound trends for both data sets. While the trends are qualitatively similar, it is evident that the relationships between road and population density are different for Saint Lucia and eastern Canada, with population density in Canada being roughly twice that in Saint Lucia for the same road density. These differences are not surprising given the distinctly different socio-economic, climactic and geographic settings in the two countries, which would suggest very different development patterns. These two data sets examined together may suggest that a reasonable range of best fit trends, spanning a wide variety of geographic settings, could be captured by the following relationship: DP ~ 10 to 30 DR^2.
The best fit trend for eastern Canada has been used in GIS with available road network vector (polyline) data to develop a map of inferred population density within the area of the former Champlain Sea, and this is shown in Figure 14. This map may be compared to that based on municipal boundaries in Figure 2, and again shows concentrations of development around the major urban areas of Ottawa, Montreal and Quebec. Figure 15 shows the inferred population density model enlarged for the Ottawa-Hull area. This map shows the expected very high densities in the urban cores, and lower densities in the rural areas between the urban cores and outlying suburbs. Figure 16 shows the population density model with roads overlain for comparison, and it is evident the inferred population density mirrors road density, which is intuitively correlated with density of development and habitation.
5 A WORKING EXAMPLE OF USE OF THE INFERRED POPULATION DENSITY MODEL: QUALITATIVE SEISMIC RISK IN EASTERN CANADA
This section presents a brief illustrative example of use of the inferred population density model in qualitative risk assessment in eastern Canada, using the example of seismic hazard, which is widespread and relatively well understood.
Seismic hazard in Canada is defined in terms of uniform hazard spectra for different annual probabilities of occurrence. Ground motions corresponding to different return periods, or to different annual probabilities of occurrence, are provided by Halchuk and Adams (2008). These include peak ground acceleration (PGA) values, along with spectral acceleration (Sa) values for different periods, or frequencies, of ground motion. Figure 17 shows PGA values corresponding to an annual probability of 0.0021 (1 in 475 year event), the annual probability used for design of buildings in Canada prior to 2005 (National Research Council, 2005). This return period event has been selected for the qualitative risk analysis as it represents expected ground motions for a relatively low probability of occurrence. The analysis could be conducted for the current building code design motions (i.e. 1 in 2475 year event), but this would yield no substantive difference in the resulting relative risk. This map also shows recorded historical earthquakes (Natural Resources Canada, 2012) for the area, to give an indication of the distribution of past seismic activity.
It is possible to make a qualitative estimate of probable damage intensity from expected PGA by converting to expected Modified Mercalli Intensity (MMI). Various relationships are available in the literature, including the following, due to Wald et al. (1999):
 MMI = 3.66 log10(PGA) – 1.66
The MMI scale is qualitative, constituting an arbitrary ranking based on observed effects of earthquakes. Nevertheless, it is often used in practice to describe earthquake intensity, and its descriptors have tangible meaning, easily understood by the layperson, with relatively distinct breakpoints. The MMI levels of interest in eastern Canada for the 475 year seismic hazard, corresponding to PGA values up to nearly 0.6 g, range from IV to IX on the 12 level MMI scale, which ranges from I to XII. Table 1 provides abbreviated descriptors of typical damage associated with these intensity levels, along with mean PGA values inferred by Wald et al. (1999), and break points used in this paper for the subsequent analysis. The PGA limits selected for analysis vary slightly from the midpoints between Wald et al. values, and have been chosen to provide a uniform (i.e. linear) PGA scale that closely matches their best fit trend, consistent with linear trends interpreted by other authors, including Gutenberg and Richter (1942).
Estimated earthquake intensity associated with the 1 in 475 year earthquake is illustrated in Figure 18, which converts the design PGA value to an assumed MMI value using the break points listed in Table 1. It may be noted that MMI is a crude indicator of expected earthquake hazard; different structures with different dominant periods will be affected by different components of the earthquake motion, hence the use of ground motions with spectral acceleration values for several different periods. However, in examining the general potential for damage to a wide variety of structures in developed areas, PGA (and by extension the inferred MMI) is a better single indicator than any of the individual spectral accelerations. It should also be noted that the design PGA values are provided for firm ground, and that there will be local amplification or de-amplification of ground motions depending on local subsurface conditions. Use of local seismic hazard values from detailed seismic microzonation studies would be a useful advancement on the present high level analysis, particularly in areas of high relative risk. Finally, this map does not address the potential for important secondary seismic effects, such as liquefaction. Consideration of these important details is judged to be appropriate at a subsequent, more detailed level of study, following a first, high level prioritization based on qualitative risk.
The seismic intensity map shown in Figure 18 may be taken as a first approximation of seismic hazard, in that it provides a qualitative indication of the probable level of earthquake damage across the study area associated with strong earthquake motions, given the occurrence of a 475 year event and the presence of development. Risk is a combination of hazard and consequence, where consequence is a measure of the likelihood of presence and value of elements at risk, combined with consideration for their vulnerability to damage, given occurrence of the hazard. It will be assumed, for ease of analysis, that both value and vulnerability are constant across the study area for various classes of elements at risk. Therefore a qualitative risk analysis can be conducted by comparing the occurrence of the hazard with the presence of elements at risk. Subsequent, more detailed analysis of identified higher risk areas would need to take more careful consideration of value and vulnerability of elements at risk.
In this analysis, the key elements at risk are human habitants of the study area, which have been shown to be distributed spatially as illustrated in Figure 14, and physical infrastructure, including dwellings and transportation infrastructure. The primary physical infrastructure associated with human settlement is spatially distributed according to the same relationships as shown for the population density model, although with different numerical relationships. Hence, the relative findings for risk to human habitation will apply similarly for risk to physical infrastructure.
A simple algorithm has been used to determine qualitative risk values based on expected earthquake intensity and local population density, and the interpreted values are summarized in Table 2. Qualitative risk is calculated according to:
 R = D_P x 10^(MMI – 3)
D_P and MMI were defined earlier, and,
R = relative risk, a unitless parameter ranging from 1 to 10^9 (the 3 in equation  is simply a scaling factor; without it the calculated risk would instead range from 10^3 to 10^12), and which is further subdivided as follows:
Low: R < 10^3 (89.6 % of study area)
Moderate: R < 10^5 (8.9 % of study area)
High: R > 10^5 (1.5 % of study area)
A qualitative seismic risk map for the area of the former Champlain Sea in eastern Ontario and southern Quebec based on these categories is presented in Figure 19. It should be noted that these risk categories do not have a specific numeric meaning in terms of value of expected loss, but they are expected to present a strong first estimate of the spatial distribution of relative seismic risk; elevated risk is present in the map where population, seismic hazard, or both are elevated, and risk is otherwise low. The map therefore serves to focus attention spatially to areas with elevated seismic risk, to support prioritization of further study.
The qualitative risk map yields some expected results along with some unexpected findings. Elevated (moderate) risk is concentrated in all the major urban areas, including Ottawa, Montreal and Quebec, along with numerous smaller communities. This is expected, since these lowlands, contained largely in regional scale depressions of the Ottawa-Bonnechere Graben and Saint Lawrence Rift, are affected by elevated seismic hazard in comparison with most of Canada other than the west coast and high arctic. However, the highest risk is indicated along both sides of the Saint Lawrence River east of Quebec City, affecting relatively lower density, predominantly rural, areas. In these areas, risk is elevated due to significantly higher PGA values, in the range of 0.4 to 0.6 g for the 475 year event. There may be some need to refine the risk categorization algorithm to lower the relative risk in these areas, to give a better balance of seismic risk across the study area given the relatively low density of development; however, the result rightly attracts attention to this area where design motions are among the highest in Canada. The qualitative risk map in Figure 19 appears to suggest meaningful priorities for seismic risk mitigation, highlighting concentrations of development and areas with very high expected ground motions, while de-emphasizing areas with neither high seismic hazard nor elevated population density.
The preceding analysis has shown that a simple power law relationship exists between population density and road density, and that this relationship is valid over several orders of magnitude, spanning the density range of interest. For eastern Canada, we find DP ~ 28 DR^2, and comparison of data from the island nation of Saint Lucia suggests that a general relationship of DP ~ 10 to 30 DR^2 might be broadly applicable for a wide range of geographic, climactic and socio-economic settings.
The inferred relationship between population and road density has been used to generate a high resolution population density model at the regional scale, using road density in GIS as a proxy. Such models are expected to be of use in regional scale risk assessment for natural hazards, with risk being determined qualitatively, in a relative sense, to support the prioritization of further efforts in higher risk areas. This paper has presented an illustrative example of qualitative risk assessment using seismic hazard, which has well understood spatial and temporal distribution across a study area in eastern Canada. The regional scale population density map presented in this paper could be extended over the whole country for use with other natural hazards or to address seismic hazard and risk elsewhere in Canada.
BGC Engineering Inc. provided geomatics assistance in support of the analysis, which was capably undertaken by Lucy Lee, GIS analyst/developer.
Gutenberg, B., and Richter, C.F. 1942. Earthquake magnitude, intensity, energy, and acceleration. Bulletin of the Seismological Society of America. 32: 163-191.
Halchuk, S., and Adams, J. 2008. Fourth generation seismic hazard maps of Canada: maps and grid values to be used with the 2005 National Building Code of Canada. Geol. Survey of Can. Open File 5813, 32 pgs.
National Research Council of Canada. 2005. National Building Code of Canada.
Natural Resource Canada. 2012. National Earthquake Database. Available: http://www.earthquakescanada.nrcan.gc.ca/stndon/NEDB-BNDS/bull-eng.php. [Accessed March 2012].
Owens, P.M., Titus-Ernstoff, L., Gibson, L., Beach, M.L., Beauregard, S., and Dalton, M.A. 2010. Smart density: a more accurate method of measuring rural residential density for health-related research. International Journal of Health Geographics. 9: 8.
Quinn, P.E., Hutchinson, D.J. Diederichs, M.S., and Rowe, R.K. 2010. Regional scale landslide susceptibility mapping using the weights of evidence approach: an example applied to linear infrastructure. Canadian Geotechnical Journal. 47: 905-927.
Quinn, P.E., Hutchinson, D.J. Diederichs, M.S., and Rowe, R.K. 2010. Characteristics of large landslides in sensitive clay in relation to susceptibility, hazard and risk. Canadian Geotechnical Journal. 48: 1212-1232.
Statistics Canada. 2007. Population and dwelling counts: A portrait of the Canadian population. Available at: http://www12.statcan.ca/census-recensement/2006/rt-td/pd-pl-eng.cfm [Accessed February, 2012]
United States Geological Service. 2012. The modified Mercalli intensity scale. Available: http://earthquake.usgs.gov/learn/topics/mercalli.php. [Accessed March 2012].
Wald, D.J., Quitoriano, V., Heaton, T.H., and Kanamori, H. 1999. Relationships between peak ground acceleration, peak ground velocity, and modified Mercalli intensity in California. Earthquake Sectra. 15(3): 557-564.