using real-time numerical weather predictions

A project established at the National Institute of Water and Atmospheric Research (NIWA) in New Zealand is aimed at developing a prototype of a real-time land- slide forecasting system. The objective is to predict tem- poral changes in landslide probability for shallow, rainfall- triggered landslides, based on quantitative weather forecasts from numerical weather prediction models. Global weather forecasts from the United Kingdom Met Office (MO) Nu- merical Weather Prediction model (NWP) are coupled with a regional data assimilating NWP model (New Zealand Lim- ited Area Model, NZLAM) to forecast atmospheric variables such as precipitation and temperature up to 48 h ahead for all of New Zealand. The weather forecasts are fed into a hy- drologic model to predict development of soil moisture and groundwater levels. The forecasted catchment-scale patterns in soil moisture and soil saturation are then downscaled us- ing topographic indices to predict soil moisture status at the local scale, and an infinite slope stability model is applied to determine the triggering soil water threshold at a local scale. The model uses uncertainty of soil parameters to pro- duce probabilistic forecasts of spatio-temporal landslide oc- currence 48 h ahead. The system was evaluated for a damag- ing landslide event in New Zealand. Comparison with land- slide densities estimated from satellite imagery resulted in hit rates of 70-90%. tralian crustal plates, is characterized by a dynamic active landscape - mountainous and hilly regions with steep slopes cover large areas of the country; uplifted and dissected sedi- mentary and volcanic rocks with low strength properties are


Introduction
Shallow landslides triggered by intense rainstorms occur in most mountainous landscapes and can cause significant damage to human lives and properties (Glade, 1998). New Zealand, located at the interface between the Pacific and Aus-Correspondence to: J. Schmidt (j.schmidt@niwa.co.nz) tralian crustal plates, is characterized by a dynamic active landscape -mountainous and hilly regions with steep slopes cover large areas of the country; uplifted and dissected sedimentary and volcanic rocks with low strength properties are widespread. New Zealand's location in the South-western Pacific means that sub-tropical cyclones hit the land mass of New Zealand on a regular basis creating heavy rainfall over large areas (McConchie, 1992). Settlement history has seen major shifts in land use and land coverage, beginning with Maori settlement around one thousand years ago and continuing with European settlement. Large land areas have been cleared of indigenous forest and converted into crop land and pasture, exposing weak parent material to the actions of weather and climate (McConchie, 1992). These factors contribute to the pronounced landslide susceptibility of New Zealand landscapes (Glade, 1998;Fig. 1) causing significant damages to economy and society (McConchie, 1992).
One of the possible ways to reduce the damaging impact of those hazards to life, infrastructure, and livestock is to increase preparedness through landslide warning systems (Wilson, 2005). Approaches in landslide warning have been based either on detailed site-specific investigations and movement monitoring programmes (Malet et al., 2005) or on statistical threshold models (Guzetti et al., 2007;Wilson, 2005). A few recent studies showed the potentials of remote sensing data (Tralli et al., 2005) and of hydrological and slope stability models to be applied within a landslide forecasting framework (Malet et al., 2005;Qui et al., 2007).
Here we aim to develop regional landslide forecasting procedures based on process-based models. To simulate and to predict rainfall-triggered landslides, three fundamental processes have to be modelled (Crozier, 1999): 1. The spatial pattern and temporal intensity variation of rainfall effects where soil can be saturated -hence accurate prediction of rainfall is necessary.
Published by Copernicus Publications on behalf of the European Geosciences Union. Extensive shallow slope failures in New Zealand's North Island soft rock country triggered from intense rainfall events (image from Hancox and Wright, 2005). Here shown landslide damage to pasture areas during the February 2004 storm event in the southwest of the North Island (Hancox and Wright, 2005;Dymond et al., 2006). Inset shows the extend of landslide damage as derived from satellite imagery (Dymond et al., 2006) and the Whangaheu catchment boundary used as a test site for this study.
2. Soil-hydrologic processes control which parts of the landscape are most affected by soil moisture increases and hence, most prone to landsliding. Therefore, understanding and simulation of soil hydrology is required. Uncertainties in models of soil hydrology are also related to the accuracy of soil and land cover properties.
3. Soil wetness changes geotechnical soil propertiesstresses and strengths are altered and -once stresses get larger than strength -the soil becomes unstable. Appropriate soil mechanical formulations in slope stability models are used to model these processes.
This paper reports initial results from a project at the National Institute of Water and Atmospheric Research (NIWA) in New Zealand aiming to develop a prototype of a realtime landslide forecasting system. The objective is to predict temporal changes in landslide probability for shallow, rainfall-triggered landslides, based on atmospheric forecasts as they are delivered by a numerical weather prediction (NWP) model. The atmospheric forecasts are fed through a model chain of hydrologic models and slope stability models to provide probabilistic forecasts of spatio-temporal landslide occurrence for vulnerable regions in New Zealand. This paper presents the technical implementation of that coupled forecasting system and some initial results for a regional rainfall event that triggered a large amount of landslides in the lower North Island of New Zealand (Fig. 1).

Fig. 2.
Model domain of the NZLAM numerical weather prediction model for New Zealand and example of forecast output: low cloud (1000-800 hPa) and mean sea-level pressure.

Landslide forecasting technology
The aim of the study was to develop and to evaluate a landslide forecasting system based on physical principles, using output from an NWP model. The system consists of several components. A regional NWP model for New Zealand (NZLAM) predicts status of the atmosphere over New Zealand up to 48 h ahead. A model of catchment hydrology simulates soil moisture and groundwater levels over the forecast period of the weather forecasts. As the fundamental spatial resolution of the catchment model is quite large (i.e. small sub-catchments), and soil moisture varies significantly with topographic position (Beven and Kirkby, 1979), we used topographic indices to disaggregate sub-catchment average soil moisture conditions to local soil moisture conditions (Beven et al., 1995). A slope stability model calculates effects of soil moisture changes on stresses and strengths within a hillslope. Considering the uncertainties in the model parameters and forcings, uncertainty estimates of stresses and strength are used to estimate failure probabilities. The individual system components are described in the following sections.

Weather forecast components
Global weather forecasts are produced by the UK Met Office Unified Model TM (Davies et al., 2005) providing the lateral boundary conditions for the regional, high resolution New Zealand Limited Area Model (NZLAM) to forecast atmospheric variables such as precipitation and temperature out to 48 h ahead for New Zealand (Fig. 2). Both NZLAM and the global model utilise the same dynamical cores. NZLAM runs in a 6 h "warm-cycling mode", assimilating local observations, lateral boundary conditions, and initial conditions from previous runs within the model domain. NZLAM is a grid point model (324×324×38), where the distance between the grid points is 12 km. NZLAM's model structure is an implementation of the Unified Model, including three dimensional variational assimilation (Lorenc et al., 2000) of observations from land, ship, and upper air stations, drifting buoys, aircraft, and satellites. The background fields, 6 hour forecasts, are derived from the previous forecast cycle. This ensures that the model forecast is as close to reality as possible. NZLAM runs twice daily at NIWA Wellington and produces hourly model output.

Hydrological model
We use the spatially-distributed, physically-based model TopNet (Bandaragoda et al., 2004) to simulate basin hydrology based on weather forecast forcings. Each river basin is represented by a number of sub-catchments (each sub-basin is approximately 5 km×5 km in area), and a stream network that connects each sub-catchment to the basin outlet. TopNet has two fundamental components: a water balance model that simulates the dominant hydrological processes in each sub-catchment (snow accumulation and melt, canopy interception and throughfall, the sub-catchment water balance), and a routing model that simulates the flow of water through the river network (including reservoirs and lakes). The subbasin water balance model has its origins in the classical Top-model formulation (Beven and Kirkby 1979), which is a simplified model of topography-controlled saturated areas and baseflow from the saturated zone. TopNet extends the original Topmodel formulation to model snow accumulation and melt, and canopy interception and throughfall (Bandaragoda et al., 2004). Parameterization of the sub-basin water balance model requires data on topography (elevations, distributions of wetness index and distance to stream), soils (soil depth, porosity, conductivity), and landcover (canopy capacity, overland flow velocity) (see Bandaragoda et al., 2004 for a detailed description). These data are derived by the catchment data delivery system as described below.
For accurate modelling and forecasting of soil moisture it is critical to initialize hydrologic states (including soil moisture) as close as possible to reality. Our model is initialized once (at the "beginning of time") with arbitrary values of sub catchment model states (snow water equivalent, canopy storage, soil moisture, depth to the water table, and river level). The model is then run with a 48-h forecast from the NWP, and all model states are saved at hour 12 of the forecast. The next forecast is then initialized using the saved model states from hour 12 of the previous forecast. After several weeks to months of cycling (i.e., several wetting and drying cycles), the states in TopNet have "warmed up" to correspond to the current real-world conditions. Alternatively, we can run Top-Net with climate observations up to the start of the forecast period to estimate basin states, and then run TopNet with the NZLAM output for the 48-h forecast period. Simulations from TopNet are made using frequency distributions that describe subcatchment variability for the soil zone. Field observations show that areas of soil saturation tend to occur in convergent hollow areas. It has also been reported that landslides most commonly originate in areas of topographic convergence due to higher soil moisture (Montgomery and Dietrich, 1994). Therefore the following technique was adopted to downscale simulated soil saturation onto a finer spatial resolution using topographic variability. Beven et al. (1995) showed that local relative soil saturation m can be estimated as Where f is a parameter describing transmissivity change with soil depth and z is local soil depth. λ is the spatial average of the topographic parameter ln (a/ tan θ) for the area (the catchment) over which the average soil saturation m was estimated. This equation expresses the deviation between catchment average soil saturation m and local soil saturation m at any point in terms of the deviation of the local topographic index from its mean. We calculated the topographic parameter ln (a/ tan θ) and fine-scale soil moisture using a 30m digital elevation model.

Slope stability modelling
The geotechnical calculations applied in this study are based on the infinite slope form of the Mohr-Coulomb failure law in which the downslope component of the soil just at failurethe shear stress τ , is equal to the strength of resistance caused by cohesion c (which can consist of soil cohesion and root strength c=c s +c r ), and by frictional resistance due to the effective normal stress σ on the failure plane.
The frictional resistance term is determined by the angle of internal friction φ, the normal stress on the shear plane σ , and u -the pore pressure opposing the normal load. Downslope stress τ =w sin θ and normal stress σ =w cos θ are determined by the soil weight w and shear geometry, i.e. slope angle θ. Soil weight at the shear surface is determined by soil depth z and the depth to the saturated zone (z−h).
γ b , γ s and γ w are unit weights of moist soil, saturated soil and water, respectively. Pore pressure u=γ w h cos θ is determined by slope angle and the weight of water above the shear surface ie. the height of the water table in the soil column h. This leads to the full equation for the infinite slope model at failure.
Root cohesion varies widely in space (land cover and land use pattern) and time (growth periods) and is hard to quantitatively estimate by a sampling programme. In this study we neglected the effects of root cohesion and assumed the cohesive term of the slope is caused by soil particle cohesion only. The estimate of stability is therefore conservative -i.e. a slope considered as unstable may be stable in reality if significant amounts of root cohesion add to shear strength. However, the vulnerable pasture landscapes in New Zealand are often characterized by a lack of vegetation. Moreover, roots need to be anchored in the bedrock to have significant effect on slope stability. To find the triggering soil saturation h/z which leads to failure Eq. (1) can be rearranged to  Assuming that the complete soil column fails, m (Eq. 1) can be set to h/z, and the soil saturations derived from the hydrological model (Eq. 1) and from the slope stability model (Eq. 5) determine if the soil is unstable.

Model parameterisation
A critical part of the landslide forecasting system is to be able to parameterize the hydrology and slope stability model components for any catchment in New Zealand. We used nationally consistent databases for that purpose. Spatial information about the catchments of New Zealand's rivers has been processed in an earlier project to produce the New Zealand River Environment Classification (REC) (Snelder and Biggs, 2002). The REC includes a network of approx. 600 000 river reaches and sub-basins for New Zealand, including their reach and surface catchment geometries and topologies. Topographic reach and catchment properties have been derived by overlaying on a 30 m digital elevation model (DEM). Land cover, soil, and geological properties for all catchments in New Zealand have been derived by overlays on the New Zealand land cover database (LCDB) and the New Zealand Land Resource Inventory (LRI) -which are nationally consistent databases of land information (Newsome et al., 2000;Wilde et al., 2000). This database constitutes a large array of environmental variables for each sub-catchment in New Zealand (Table 1). Ancillary files for the individual model components (hydrology and slope stability) and for specific model domains can be extracted to generate all parameters necessary for running a forecast model for any region in New Zealand.

Uncertainty estimation
The uncertainties in landslide modelling -and especially high resolution landslide forecasting -is very high. The fol- Fig. 5. Discrimination diagram displaying the frequency distribution of predicted failure probability (Fig. 4) for areas with observed landslides and no observations of landslides (from Dymond et al., 2006) indicate that observed landslide locations correlate with higher predicted landslide probabilities (log-scale inset emphasises the differences at lower frequencies).
lowing three sources of uncertainties can be identified in the used approach. (i) Uncertainty in the weather forecasts: the rain (and other modelled atmospheric variables) will never be exactly distributed in space and time as predicted due to the inaccuracies in the weather models, parameters, and initial conditions (Simmons and Hollingworth, 2002). (ii) Uncertainty in the soil hydrology models: even if we could model the atmosphere accurately, the modelled pattern of soil saturation will never be exactly distributed as predicted due to uncertainties in the hydrology models, parameters, and initial conditions (Schaake et al., 2007). (iii) Uncertainty in the slope stability model: even if we could forecast the soil hydrology accurately, hill slopes will never exactly fail where predicted due to inaccuracies in the geotechnical models and soil parameters. Uncertainties in weather and hydrology predictions can be quantified through ensemble modelling by perturbing initial conditions and model parameters and using different model structures as ensemble members (Schaake et al., 2007). This would lead to ensembles of soil saturation scenarios for the forecast range. Uncertainties in the slope stability model can be quantified by running ensembles with different model structures and statistical realizations of geotechnical parameters.
In this study we only considered uncertainties in the parameters of the soil mechanical model. Some of the used soil parameters (e.g. soil depth, field capacity) could be extracted directly from the databases others had to be created as surrogates from lookup tables. Especially cohesion c and angle of friction φ had to be estimated from soil texture via pedotransfer functions (PTFs) (Blondeau, 1973) and contain a high degree of uncertainty. The uncertainty of those variables can be described using probability density functions (PDFs). As a simple first order approximation we described uncertainty of those variables as uniform probability distributions This leads to uncertainty bounds m min , m max in the estimated saturation degree h/z=m as: Similarly, the uncertainties in the shear stress and shear strength can be derived and the failure probability can then be calculated by the degree the modelled shear stress exceeds shear strength for each point in the model domain. It is planned to extend this uncertainty model to other types of PDFs and to explicitly quantify uncertainties in weather forecasts and simulations of soil hydrology (as described above).

Delivering forecast information
The described model components are part of a larger model suite for environmental forecasting, developed at NIWA (EcoConnect, 2008). The landslide forecasting system has been applied for two catchments in New Zealand. Relevant information is extracted from the output files to produce graphical model output twice daily (Fig. 3), including maps of landslide probability for the modelled catchment for each timestep (i.e., hourly). These maps can be animated to produce a dynamic visualization of the spatio-temporal pattern of landslide probability for the whole basin over the forecast period. Alternatively, the maximum failure probability for each point for the whole forecast period can be displayed to get a quick overview of the affected areas. As discussed above, the uncertainties in the modelled system weathersoil-hydrology -slope stability are very high -therefore displaying the results at high spatial resolutions might not be appropriate. A third forecast product can be derived from the temporal development of potential landslide-affected area for different probability classes -the product allows the forecast user to define risk and alert levels to understand the potential magnitude of the forecasted landslide hazard (Fig. 3).  Fig. 7. Relationship between predicted landslide density and observed landslide density per km 2 (Fig. 6), legend indicates number of cases. Inserted are the contingency tables (numbers in thousands) derived for different thresholds of landslide density.

Discussion and conclusions
We presented a probabilistic modelling system for forecasting shallow, rainfall-triggered landslides as they frequently occur in pastoral areas. The model has yet to be extensively calibrated and validated -here we provide an initial test using data from the extreme event in February 2004 (Schmidt et al., 2005), which caused more than 60 000 landslides in the lower North Island of New Zealand (Dymond et al., 2006). The system was run in hindcast mode for that event for the Whangaheu catchment. Figure 4 illustrates high resolution probabilistic forecasts of landslides, along with the actual occurrence from post-event satellite-based landslide mapping (Dymond et al., 2006). Qualitatively, the predicted spatial landslide distribution is comparable with the spatial distribution from satellite imagery (Fig. 4). Comparison of the high-resolution model output against the locations of landslides derived from satellite imagery in a discrimination diagram (Fig. 5;Clark and Slater, 2006) showed that areas of higher predicted landslide probability coincide with observed landslide area. However, the high-resolution model output also predicts high failure probabilities where no landslides occurred (high false alarm ratio) -as can be seen in Fig. 4, not all the predicted unstable areas fail. Predicting the exact location of landslides is extremely challenging -a more tractable problem is predicting the spatial density of landslides. The binary landslide observations can easily be aggregated to larger scales to derive landslide densities (affected area per km 2 ). Similarly, the probabilistic forecasts of landslides can be spatially aggregated to provide deterministic (single-value) forecasts of landslide density. The linkage between the high resolution probabilistic forecasts of landslide occurrence and the lower-resolution deterministic forecasts of landslide density is as follows: A single scenario of binary landslide occurrence can be derived from the probabilistic forecast by assigning an individual grid cell the value "1" if a random number from a uniform distribution, U [0, 1], is less than the forecast probability, and otherwise assigning the grid cell the value "0". The landslide density from this scenario can be computed in the same way as the observed density. Since the probabilistic forecasts define the fraction of time a landslide is expected to occur under the forecast conditions, the average spatial density from many randomized simulations will match the spatial average of the forecast probability. The accuracy of the deterministic density forecasts then depends on the statistical reliability of the probabilistic forecasts, that is, the extent to which the fraction of observed landslides matches the forecast probability. Figure 6 illustrates the spatial pattern of landslide densities derived using a 1 km 2 kernel. As expected, the spatial aggregation improves the correspondence between forecasts and observations. Figure 7 shows the joint density of predicted and observed landslide density in a scatterplot. While the scatter is quite large (42% explained variance), there is a general relationship between predictions and observations. Next we categorized the landslide densities and derived contingency measures. We achieved high hit rates (observed landslide densities match predicted probabilities) of about 70-90%. False alarm ratios (no landslides observed with predicted landslide probability) are about 30%, indicating the over prediction observed in Figs. 4 and 6. It has to be kept in mind, however, that these results are achieved without any model calibration.
The following improvements and extensions are planned to enhance the forecasting system. Given the uncertainty in the individual model components it is important to extend the current probabilistic approach to include more uncertainties in the model chain. One of the most important uncertainties is the spatial accuracy of the rainfall forecasts. All the involved models can be run in ensemble mode by perturbing the initial conditions and the model parameters. This will give multiple forecasts and enable us to quantify more uncertainties in the forecasts. Data assimilation is needed to ensure high-quality forecasts of atmospheric and hydrologic fluxes and states (Schaake et al., 2007). Currently we are assimilating atmospheric satellite observations and river discharges from recorders; we are planning to include the assimilation of soil moisture observations. We also plan to extend the forecast range by seamlessly merging short-range weather forecasts to medium-term (week, month) climate forecasts using statistical and ensemble techniques. The soil mechanical model and its parameters are based on existing (weak) data sources for New Zealand. The assumptions of the infinite slope model or "Planar Slope Model" are based on a model of a uniform, translational failure of uniform depth and slope angle. The critical failure condition is achieved only by increasing pore pressure of the soil water (and consequently decreasing shear strength of the soil) at the shear surface. This is a simplification of real-world conditions, where failure occurs in curved slip surfaces and multiple mechanisms including fluid failure. However, new techniques emerge to deal with those issues (Qui et al., 2007).
An important question is related to the communication of model outputs -including high degrees of uncertainties -to end-users using appropriate forecast products (Fig. 3). The validation results (above) showed that aggregated landslide probability on a coarser resolution was more reliable than the immediate fine-scale forecast results. The inherent uncertainties in weather simulation, hydrological modelling, and geotechnical models mean that the landslide forecast results contain high degrees of uncertainties, in particular if verified on a local scale, hence, the forecast results need to be upscaled to regional levels to be useful for applied purposes.